Discrete Choice and Competitive Reactions: End-to-End Simulation with the R Package cash

Discrete Choice and Competitive Reactions: End-to-End Simulation with the R Package cash

Introduction

As much as product design is a creative task, most of the time it is also an optimization problem in terms of market share, revenue or profit. To maximize the latter objective, which encompasses the former two as special cases, it is necessary to determine the respective products' unit contribution margins as well as their expected (relative or absolute) demand. Since the pioneering works on conjoint and discrete choice analysis by 1971_Green_Rao, later Nobel laureate 1974_McFadden and 1983_Louviere_Woodworth, conjoint choice experiments have matured into a default for the elicitation of consumer preferences 2000_Louviere_et_al,2007_Gustafsson_et_al,2009_Train,2014_Rao, 2015_Hensher_et_al,2021_Baier_Brusch,2024_Hess_Daly, which allow us to forecast a product's share of choice (i.e., the expected relative demand) within a certain competitive scenario. More specifically, in such experiments, respondents are repeatedly asked to choose their preferred option from varying sets of alternatives resembling real purchasing situations. The resulting discrete choices are then decomposed into the respondents' individual utilities for the presented feature levels of the product by training correspondingly defined random utility models on the data. With these individual feature level utilities, it is possible to predict the respondents' choice behavior in a competitive scenario of alternatives and therefore the shares of choice required for product (line) optimization.

Now, after (successfully) launching a new or modified product (line), maybe due to the help of the above mentioned product (line) optimization, competitors will most likely react to not lose market share, revenue or profit. Anticipating such reactions has the potential of sustainably strengthening a firm's competitive position. While there exists a substantial body of research on product (line) optimization 2008_Belloni_et_al, 2024_Baier_Voekler, this does not hold true for the simulation of competitive reactions bringing together discrete choice models and non-cooperative game theory 2025_Dressler_et_al, especially the 1951_Nash and 1934_vonStackelberg game-theoretic strategies. As far as we are aware, there is neither a free and open-source nor a commercial software at hand in this regard. Although the existing studies are inevitably built on code, no one seems to have developed and published an official package. Furthermore, only a fraction of them concentrated on the firms' strategic product policy, i.e., the more complex long-run competition on both product price and design (i.e., non-price features), none of which integrated fully Bayesian choice models into their simulation, let alone under different types of choice behavior (probabilistic versus deterministic) 2025_Dressler_et_al.

For this reason, we introduce the \proglang{R} 2026_RCoreTeam package \pkg{cash} (conjoint $+$ Nash), which is available from the Comprehensive R Archive Network (CRAN) at \url{https://CRAN.R-project.org/package=cash} and can be installed and loaded via

\pkg{cash} contributes functions for simulating Nash equilibria based on discrete choice models and additionally for the steps to reach the latter, thereby offering a tool for end-to-end simulation as well as simulation on the basis of external data, whether it be real or artificial. To combine these elements, it partly draws from established discrete choice \proglang{R} packages.

Doing justice to the structure of \pkg{cash}, i.e., the viable end-to-end simulation and simultaneous option to leverage the functions individually (or in parts), the rest of this paper consists of ten sections, seven of them following the order of \pkg{cash}'s eight functions if employed as an ensemble (the last two functions share a section), each presenting the respective function's theoretical background and practical use as independently as possible without compromising their cohesiveness. As \pkg{cash} has recently been put to the test in a large-scale Monte Carlo study by us demonstrating its potential, the theoretical background of the functions is taken from the methodological framework of the corresponding paper 2025_Dressler_et_al. In particular, the functions are concerned with preference generation ((ref)), choice design ((ref)), error and response simulation ((ref)), hierarchical Bayesian estimation of mixed logit models ((ref)), design assessment ((ref)), convergence and model assessment ((ref)), and pre-computations and Nash competition ((ref)). Before diving into the functions, (ref) gives some preliminary information. (ref) discusses the adjustments required to consider competition on market share or revenue instead of profit. (ref) addresses \pkg{cash}'s limitations and closes with ideas for future development.

Preliminaries

In this preliminary section, we would like to briefly explain two overarching aspects of the package that were initially intended for internal use only but are now still present because they might be helpful to the user.

1. The functions (as well as their output files) are labelled from \code{B_} to \code{I_} ($+$ name) for reasons of clarity regarding their logical order. \code{A_} is reserved for the (potential) file calling the functions.

2. Instead of returning its output to be stored in a variable, each function outputs multiple .rds files (change each function's argument \code{temp} to \code{FALSE} to not save into the session's temporary but the current working directory) to be re-imported by subsequent functions. These output files will be described in the respective sections, but there exists a single one that all functions share, called \code{A-I_finaloutput.rds}. It is a \code{list} of \code{lists} which is propagated through the functions, ultimately containing eight \code{list} objects, one for each of the eight functions: \code{$fun1}, ..., \code{$fun8}. The elements of each of those eight \code{list} objects are \setlist{nolistsep}

In other words, \code{A-I_finaloutput.rds} merges the functions' in- and outputs. As such a file takes up a lot of additional storage depending on the functions' inputs (and therefore outputs), many elements are deliberately set to \code{NULL} because they are less important, too large and/or already included. This way, \code{A-I_finaloutput.rds} still outlines everything that exists in the \pkg{cash} environment to the current point in time, but the file size stays negligibly small. If the information covered by \code{A-I_finaloutput.rds} suits the user's needs, some or all of the individual .rds files may be deleted to save storage and neaten the directory. This can be quite handy if the functions are run hundreds or thousands of times as an ensemble 2025_Dressler_et_al. However, the first of the eight functions must always be executed once to initiate \code{A-I_finaloutput.rds} and avoid subsequent importing errors. Note that if a function is called twice in the same directory, its .rds files including the respective \code{$fun} in \code{A-I_finaloutput.rds} are overwritten.

Preference generation

Theoretical background

As in many fields of science, there are numerous unsolved questions in the domain of discrete choice analysis and its extensions which are impossible to answer without simulated data. In 2025_Dressler_et_al, for example, we tackled the issue to which degree state-of-the-art discrete choice models are capable of uncovering the true Nash equilibria arising under the true consumer preferences. Obviously, true consumer preferences are not known unless they are simulated. To still obtain valid outcomes in such an endeavor, we think it is of great importance to generate synthetic utilities approximating estimates seen in empirical research and commercial applications. Derived from the latter, 2010a_Wirth,2010b_Wirth proposed an approach for sampling means and variances of a multivariate normal preference structure that has recently been employed by 2019_Hein_et_al,2020_Hein_et_al,2022_Hein_et_al, 2024_Goeken_et_al,2025_Dressler_et_al after re-examination of resulting densities based on extensive practical experience.

For the theoretical explanation, let $l\in L=\{1,...,\mathcal{l}\}$ be the feature index, $m\in M=\{1,...,\mathcal{m}\}$ be the level index, $o\in O=\{1,...,\mathcal{o}\}$ be the parameter index such that $\mathcal{o}=\mathcal{l}(\mathcal{m}-1)$, and $\lfloor{\cdot}\rceil$ be the rounding to the nearest integer function. The definition of $\mathcal{o}$ implies modelling of main effects only and having discrete features (including price) with a constant number of levels (not a prerequisite). Since firms often only use certain price points, it seems reasonable to assume discreteness for the price feature as well. Following 2010a_Wirth,2010b_Wirth, we draw $\lfloor{0.1\mathcal{o}}\rceil$, $\mathcal{o}-2\lfloor{0.1\mathcal{o}}\rceil$, $\lfloor{0.1\mathcal{o}}\rceil$ random samples from the continuous uniform densities $U(-5,-2)$, $U(-2,2)$, $U(2,5)$, respectively, and concatenate them to a vector $\bar{\bm{\upbeta}}$ of preference means. Simultaneously, $\mathcal{o}$ realizations of a random variable $X$ are taken as preference variances, where

The Gamma density $\Gamma$, with shape $\kappa$ and scale $\theta$, was fit to empirical data by 2010a_Wirth,2010b_Wirth and truncated using realizations of the uniformly distributed continuous random variables $Z_{(\cdot)}$. Two resulting tuples of distributional parameters allow for a less and a more heterogeneous structure of preference variances to emerge ($0.7,1.5,0.08,0.4,9,11$ and $0.7,4.5,0.2,2,13,18$ for $\kappa,\theta,\zeta_1,\xi_1,\zeta_2,\xi_2$). With the realizations of $X$ stacked into a vector $\bm{\upsigma}^2$, each pair of components from $\bar{\bm{\upbeta}}$ and $\bm{\upsigma}^2$ is subsequently used to draw a vector $\bm{\upbeta}_o$ of $\mathcal{i}$ random samples from a normal density, where $i\in I=\{1,...,\mathcal{i}\}$ denotes the respondent index. This leads to the preferences being dispersed differently within each feature level, which is arguably more realistic than setting a constant variance.

We then choose to assign the vectors $\bm{\upbeta}_o \forall o$ to the feature levels by combining them to an $\mathcal{i}\times\mathcal{o}$ preference matrix $\bm{\mathrm{B}}$ according to the order given by $\mathcal{o}$ random samples without replacement from a discrete uniform density $D(1,\mathcal{o})$. If we wish to define a feature $l$ with monotonously changing part-worths corresponding to a specific sorting of its levels (e.g., price), we actually do the preceding steps for $\mathcal{m}$ instead of $\mathcal{m}-1$ levels of $l$, rearrange the $\mathcal{m}$ entries of $l$ in $\bm{\mathrm{B}}$ in order for each $i$, shift them so that the first column of $l$ becomes a null vector, and lastly remove the latter from $\bm{\mathrm{B}}$. This is to avoid any violation of such a monotonicity condition by the reference category when the coding scheme of the design matrix is expected to be dummy in the response simulation and the re-estimation.

(ref) gives a typical example of homogeneous and heterogeneous synthetic utilities generated with the above method. Precisely, it displays histograms and kernel density estimators (using Gaussian kernels and Silverman's 1986_Silverman bandwidth) for $\mathcal{l}(\mathcal{m}-1)=2(5-1)=8$ parameters based on $\mathcal{i}=500$ individuals, showing plausible differences in heterogeneity and preference order without violating a monotonicity constraint exemplarily imposed on $l=1$ (violet densities). It can also be seen that the mixtures of normals resulting from this interference are not a matter of concern.

Practical use

Function and input \newline The corresponding function from \pkg{cash} is \code{B_prefgen()}:

\code{resps} defines the number of respondents and \code{lvls} the number of features and levels for which we would like to generate synthetic utilities. While \code{resps} takes an integer, \code{lvls} has to be given a vector of integers. The length of the vector is interpreted as the number of features and the values as the number of levels per feature. The function automatically drops a level per feature, as the reference level does not have to be simulated, so it is necessary to specify the actual number of levels per feature. Except for \code{seed} (always addresses \code{set.seed()} in \pkg{cash}), the other arguments can be either \code{TRUE} or \code{FALSE}. If \code{het = TRUE}, the utilities of each feature level are more heterogeneous across respondents. If \code{randomizecols = TRUE}, the generated column vectors containing the utilities for all respondents for each feature level are randomized (recall that the means for the feature levels are drawn in order from the continuous uniform densities $U(-5,-2)$, $U(-2,2)$, $U(2,5)$). If \code{sortfirstattdesc = TRUE}, the utilities of the first feature are sorted descendingly (taking the reference level (first level) utility of zero into account). We primarily implemented the latter to be able to impose a monotonicity constraint on the price feature. Since the first feature is expected to be price within \pkg{cash}, ordered from lowest to highest price, with \code{sortfirstattdesc = TRUE}, the highest utility corresponds to the lowest price and so forth. If \code{saveplot = TRUE}, a very basic \code{ggplot} 2016_Wickham of the final densities is saved.

Output and example \newline Besides \code{A-I_finaloutput.rds} and the optional \code{ggplot}, \code{B_prefgen()} outputs the two .rds files \setlist{nolistsep}

In the following example, the former is a $50\times 8$ (50 respondents $\times$ 2 features with 5 levels each minus the reference level for each feature, first feature is price with a monotonicity constraint) heterogeneous preference matrix and the latter is the vector \code{c(5, 5)} given to the argument \code{lvls}.

Choice design

Theoretical background

With the ongoing objective of preserving generality, in our view, the generation of choice designs for the estimation of main effects based on synthetic utilities needs to strive for optimality from both the statistical and practical perspective as usual, although artificial respondents certainly are quite forgiving with respect to some of the otherwise highly important decisions to be made at this stage. Therefore, the selection of the numbers of levels, alternatives and choice sets should serve to facilitate meeting respondents' capabilities (hypothetically) as well as three well-known principles of efficient choice designs, namely orthogonality, level balance and minimal overlap 1996_Huber_Zwerina,2010_Zwerina_et_al. While staying in the commonly recommended ranges considering cognitive, behavioral and cost aspects 1996_Johnson_Orme,1997_Pinnell_Englert,2001_Hensher_et_al, 2005_Caussade_et_al,2006_Hoogerbrugge_vanderWagt,2007_Haaijer_Wedel, 2012_Kurz_Binner,2013_Louviere_et_al,2019_Hein_et_al, in 2025_Dressler_et_al, for example, we defined the number of levels to be constant across features, the number of alternatives to be equal to the number of levels and the number of choice sets to be a multiple of the latter (which, of course, is often not feasible and sometimes not even desirable in real-world settings). Symmetric designs (i.e., a constant number of levels across features) may also enable (hypothetical) control of a possible number of levels effect 1981_Currim_et_al,1982_Wittink_et_al,1990_Wittink_et_al, 1990_Green_Srinivasan.

After these pre-adjustments (which are the responsibility of the user), an established option is to apply the modified Fedorov algorithm 1972_Fedorov,1980_Cook_Nachtsheim to generate choice designs by maximizing D-efficiency, the most frequently employed statistical optimality criterion therefor 2019_Street_Viney. As it cannot be ruled out that the candidate set itself prohibits reaching optimal D-efficiency in case an orthogonal array is used, the candidate set of alternatives from which the algorithm swaps may be the full factorial if computational feasibility is given. Additionally, initiating from multiple random start designs is recommended 2010_Kuhfeld. Further, note that the Fisher information matrix depends on the parameters in logit models 2019_Street_Viney such that a prior assumption has to be made for maximizing D-efficiency. If the prior is set to zero, utility balance as the fourth principle of efficient choice designs is only trivially satisfied 2010_Zwerina_et_al. Specifying plausible non-zero (Bayesian) priors and blocking larger designs 2009_Rose_Bliemer,2018_Walker_et_al may be considered to improve efficiency. Lastly, and motivated solely by the practical perspective, the avoidance of set duplicates seems sensible.

Designs for predictive testing (hold-outs) may be generated using the same procedure, with the uniqueness constraint ideally being extended to the union of training and test sets. Given the artificial nature of the respondents, there is no need for randomizing the sets.

Practical use

Function and input \newline The corresponding function from \pkg{cash} is \code{C_choicedes()}:

First, we need to define the number of features and levels via the argument \code{lvls}, equivalent to \code{B_prefgen()}, but this time automatically through the corresponding .rds file coming from \code{B_prefgen()} (or an own one, of course) containing the vector. Second, it is necessary to specify the desired number of choice sets and alternatives per choice set via \code{setstraining}, \code{altstraining} (for estimation), \code{setstest}, \code{altstest} (for predictive testing). Third, \code{prior} always has to be a Bayesian prior initially, i.e., an $\mathcal{i}\times\mathcal{o}$ preference matrix $\bm{\mathrm{B}}$ as outputted by \code{B_prefgen()}. If \code{zeroprior = TRUE}, the entries of this matrix are set to zero, and if \code{bayesian = FALSE}, its column means are taken as the (non-Bayesian) prior.

For the actual generation of the choice designs (\code{createdesign = TRUE}), \code{C_choicedes()} borrows functions from the \proglang{R} package \pkg{idefix} 2020_Traets_et_al in reduced form. As in \pkg{idefix}, we can choose between \code{algorithm = "modifiedfedorov"} and (faster but possibly less efficient) \code{"coordinateexchange"} (\code{trace = TRUE} enables tracing for both algorithms) 2020_Traets_et_al. The number of random start designs to test can be given to \code{startdesignstotest}, and by setting \code{orthoarrayascandmodfed = TRUE}, an orthogonal array generated with the \proglang{R} package \pkg{DoE.base} 2018_Groemping is employed as the candidate set instead of the full factorial in case of the modified Fedorov algorithm.

The user is explicitly invited to make use of \pkg{idefix}'s full potential (e.g., blocking larger designs) or other software 2010_Kuhfeld,2010_Zwerina_et_al by importing externally generated designs. To do so, they must come as .csv files with each row being an alternative (i.e., in long-format), row names as first column and column names as first row (the value separator is controllable by \code{csvseparator} and the names themselves are irrelevant because they are redone internally). The arguments for the paths to the .csv files are \code{externaltrainingdesign} and \code{externaltestdesign}. If the external designs are already dummy-coded, \code{codingofexternal} must be set to \code{"dummy"}, and if no coding scheme has been applied (i.e., each column corresponds to a feature and the entries are integers representing the levels), \code{codingofexternal} must be set to \code{"none"} to convert them to dummy-coding (first level becomes reference level).

As \code{C_choicedes()} needs the information for the above and other internal processes (e.g., the creation of respondent and choice IDs), the correct specification of \code{lvls}, \code{setstraining}, \code{altstraining}, \code{setstest}, \code{altstest} and \code{prior} is also mandatory if \code{createdesign = FALSE} and external designs are imported. Only for the sake of completeness, efficiency measures of the external designs may be stored in the arguments \code{efficiencyof[...]}.

If not realized by the external designs (since it is also not included in \code{C_choicedes()} through \pkg{idefix}), the exclusion of set duplicates within and between designs may be checked with \code{F_designeval()} (see (ref)). There is a reason for \code{F_} not being named \code{D_}, which will be elucidated in the respective section, but note that \code{F_} is an exception in regard to the logical order from the practical perspective and can be directly called after \code{C_} as one would expect. \newline

Output and example \newline \code{C_choicedes()} outputs the two .rds files \setlist{nolistsep}

each being a \code{list} with the design matrix stored in \code{$design} in long-format and the D(B)-error in \code{$error} (among other information).

In the first of the two examples below, designs are generated with \code{C_choicedes()} for 2 features with 5 levels each (see \code{B_levels.rds} from \code{B_prefgen()} in (ref)), 15 choice sets and 5 alternatives per choice set.

For the second example, let us import the following external .csv designs, ext1.csv and ext2.csv 2010_Kuhfeld,2010_Zwerina_et_al\footnote{\copyright\ 2024 SAS Institute Inc. SAS and all other SAS Institute Inc. product or service names are registered trademarks or trademarks of SAS Institute Inc., Cary, NC, USA.}, where ext1 (on the left, for training) has the same specifications as the first example and ext2 (on the right, for testing) only consists of 5 choice sets:

Here, the row and column names are meaningful, but that does not matter.

Error and response simulation

Theoretical background

Next, following random utility theory, noise must be added to the individual deterministic portion of utility of the artificial respondents across alternatives and choice sets prior to the simulation of choices. For the generation of these random error terms, previous simulation studies fixed a scalar with which to multiply the standard variance of the chosen error density 2002a_Andrews_et_al,2003_Andrews_Currim,2010a_Wirth, 2010b_Wirth,2020_Hein_et_al or set the variance of the latter to a percentage of the combined variance of the deterministic and the stochastic portion of utility 1975_Srinivasan,1981_Wittink_Cattin, 1989_Wedel_Steenkamp,1996_Vriens_et_al,2002b_Andrews_et_al. The second approach is far more flexible as it automatically adjusts for changes in the variability of the deterministic utility (e.g., when the number of features serves as an experimental factor). Nonetheless, explicitly defining a ratio of measures of dispersion still seems rather intangible because firstly, the inherent emphasis on larger values has to be taken into account when using variance or standard deviation. Secondly, expectation, skewness and kurtosis of both distributions (deterministic utility and error) are at least computationally neglected, and thirdly, the stochasticity in the subsequent assignment to the deterministic utilities is completely unrestricted.

In our view, it is preferable to directly control the proportion of error to deterministic utility, but to do so, the error has to be independent from the expectation of the deterministic utilities. Fortunately, in the case of zero-centering, the deterministic utilities (as well as the errors) themselves are deviations in their entirety, which makes the absence of translational invariance irrelevant. Otherwise, proportionate errors would actually be too large, while not zero-centering the errors might even make it impossible to sample proportionate and bidirectional ones. Therefore, zero-centering is a prerequisite for both.

Based on these thoughts, we propose controlling post hoc for the median of proportions of error to deterministic utility, which we call the median relative Gumbel error (MRGE). The distributional choice is arbitrary and just driven by the nature of the logit model used in \pkg{cash}. Let $\bm{\mathrm{x}}_{jk}$ be a $1\times \mathcal{o}$ vector of dummy variables defining the alternative $j\in J=\{1,...,\mathcal{j}\}$ in choice set $k\in K=\{1,...,\mathcal{k}\}$ such that $\bm{\upbeta}_i\bm{\mathrm{x}}_{jk}^T$ is the deterministic utility of respondent $i$ for $j$ in $k$ (linear-additive part-worth model), and let $\varepsilon_{ijk}$ be a realization of the random error $E\sim \mathrm{Gumbel}(\lambda,s)$ with location $\lambda$ and scale $s$ 2002_Stephenson. As usual, an individual's total utility is assumed to be $\bm{\upbeta}_i\bm{\mathrm{x}}_{jk}^T+\varepsilon_{ijk}$. Now, if $\bm{\mathrm{v}}_{\mathrm{abs}}$ and $\bm{\mathrm{\varepsilon}}_{\mathrm{abs}}$ denote the vectors containing $|\bm{\upbeta}_i\bm{\mathrm{x}}_{jk}^T$| and $|\varepsilon_{ijk}|$ $\forall i,j,k$, respectively, and $\oslash$ represents the Hadamard division, then the MRGE can be expressed as

The MRGE attempts to enhance interpretability by employing a more accessible fraction (e.g., let the error be 50% of the deterministic utility on average) and to increase robustness through its post hoc inspection. The median is used to prevent the errors from being scaled down without the suppression of outliers.

Applying the tuning procedure below ((ref)), the standard deviation $\sigma$ of the error density is iteratively adjusted until a required MRGE target value holds. Besides $\bm{\mathrm{v}}_{\mathrm{abs}}$, its mean $\bar{v}_{\mathrm{abs}}$ and the MRGE target value, it takes a learning rate $d$, a tolerance $t$ and a maximum number of iterations $r_\mathrm{{max}}$ as inputs, which do not have to be precisely calibrated to guarantee functionality.

After having generated the errors, the choices are simulated based on the resulting total utilities and the first choice rule,

Practical use

Function and input \newline The corresponding function from \pkg{cash} is \code{D_responsesim()}:

As the goal is to let our artificial respondents answer our choice designs, \code{D_responsesim()}, by default, imports the synthetic utilities from \code{B_} as well as the choice designs from \code{C_} through the accordingly named first three arguments. Although this data does not necessarily have to originate from \code{B_} and \code{C_}, it must always be in .rds format and exhibit the identical structure. Fortunately, the option of end-to-end simulation allows the user to take a look at such .rds defaults by simply generating them with the preceding functions.

While $\bm{\mathrm{v}}_{\mathrm{abs}}$ and $\bar{v}_{\mathrm{abs}}$ as inputs for the error simulation (see (ref)) are automatically determined via the loaded preference matrix, MRGEtarget, $d$, $t$ and $r_\mathrm{{max}}$ must be defined explicitly. We preset well-functioning values for $d$ (\code{learningrate}), $t$ (\code{stoppingcriterion}) and $r_\mathrm{{max}}$ (\code{maxiterations}) employed in 2025_Dressler_et_al such that these arguments may be left unchanged, but MRGEtarget (\code{MRGEtarget}), i.e., the desired magnitude of the error, has to be set to a percentage of the deterministic utility (so to speak).

If \code{transform = TRUE}, \code{D_responsesim()} stores the training design and responses in an additional way such that they can be processed by the second and last external discrete choice \proglang{R} package we draw from in \pkg{cash}, namely \pkg{bayesm} 2025_Rossi, which is used for model estimation in \code{E_} (see (ref)). With \code{simulate = FALSE}, it is possible to skip the error and response simulation to prepare external (experimental or differently simulated) data for model estimation. To do so, the imported design .rds files have to be identical to the ones generated by \code{C_} (\code{$design}, \code{$resps}, \code{$betas}, \code{$lvls}, \code{$sets...}, \code{$alts...}, \code{$respondentid...}, \code{$choiceid...}, rest can be \code{NA}), but \code{$design} must comprise the designs of all respondents with the responses in the last column of this matrix. \newline

Output and example \newline \code{D_responsesim()} therefore outputs up to three .rds files, which are \setlist{nolistsep}

The first two are the \code{$design} elements from the \code{C_} outputs, replicated for each respondent and extended by two columns, the simulated total utilities and the first choices. The last file is similar to the first one but tailored to the \pkg{bayesm} requirements.

To give an example, the following call of \code{D_responsesim()} generates answers to the choice designs from \code{C_} through the use of the deterministic utilities from \code{B_} and an additional MRGE of 30%.

To be clear with regard to \code{simulate}, we would get the same \code{D_bayesmdatatraining.rds} as with the previous settings if we now imported an .rds with the above matrix (without \code{sim u}) in \code{$design} and set \code{simulate = FALSE}.

Hierarchical Bayesian estimation of mixed logit models

Theoretical background

As stated in (ref), we let the stochastic portion of utility be i.i.d. Gumbel such that the individual's choice probability $\eta_{ijk}$ of choosing alternative $j$ in choice set $k$ given $\bm{\upbeta}_i$ is defined by the conditional logit model 1974_McFadden

For the coefficients $\bm{\upbeta}_i$, a multivariate normal prior $\phi$ with $\mathcal{o}\times1$ mean vector $\bar{\bm{\upbeta}}$ and $\mathcal{o}\times\mathcal{o}$ covariance matrix $\bm{\Sigma}$ is specified which describes the population's unimodal preference heterogeneity and leads to the (unconditional) mixed logit probability 1999_Brownstone_Train

Based on this well-established single component approach for the distribution of $\bm{\upbeta}_i$, the likelihood of $i$'s choice sequence $\bm{\mathrm{y}}_i$ then becomes

For the Bayesian estimation of the hyperparameters $\bar{\bm{\upbeta}}$ and $\bm{\Sigma}$ 1999_Allenby_Rossi, we extend the hierarchical model by the conjugate hyperprior distributions $\phi(\bar{\bm{\upbeta}};\bm{\upmu},\bm{\Omega})$ (normal) and $\omega(\bm{\Sigma};\mathcal{v},\bm{\Psi})$ (inverse Wishart) with $\mathcal{o}\times1$ mean vector $\bm{\upmu}$, $\mathcal{o}\times\mathcal{o}$ covariance matrix $\bm{\Omega}$, degrees of freedom $\mathcal{v}$ $(\geq\mathcal{o})$ and $\mathcal{o}\times\mathcal{o}$ scale matrix $\bm{\Psi}$. Using Bayes' theorem and omitting the marginalization over $\bm{\upbeta}_i$, the joint posterior of this 3-stage framework can be written as

To obtain draws from $P$, we use a hybrid Gibbs sampler with a random walk Metropolis-Hastings step, which allows sampling of $\bm{\upbeta}_i\forall i$ alongside the hyperparameters while avoiding simulation of the intractable $\mathcal{o}$-dimensional integral from (ref) and (ref) 2005_Rossi_et_al,2009_Train. As can be seen in 2025_Dressler_et_al, we refrain from imposing constraints on the parameter estimation and always keep the entire Markov chain for evaluation purposes. However, we control for monotonicity post hoc if necessary (see (ref)). In the rest of this paper, we will denote the hyperparameters as $\Theta$ and refer to posterior draws of $\bm{\upbeta}_i$ as just draws.

Practical use

Function and input \newline The corresponding function from \pkg{cash} is \code{E_hbmxlest()}:

It employs the sampler implemented by the \code{rhierMnlRwMixture()} function from the \proglang{R} package \pkg{bayesm} 2025_Rossi.

If not specified differently, \code{E_hbmxlest()} automatically imports the design and response .rds files from \code{C_} and \code{D_} through its first two arguments. The length of the Markov chain has to be defined by an integer given to \code{mcmciterations}. Pre-set to every 5,000 iterations, the progress update frequency of the sampler can be controlled via \code{nthiterationtoprint}. As a single component approach for the distribution of $\bm{\upbeta}_i$ is assumed, the argument \code{moncomponents} may be ignored for now. Although we prefer to save the entire chain and avoid constraints in the estimation, pre-thinning the chain and sign-restricting the draws of the $\mathcal{m}-1$ price levels (first feature) can be done with \code{nthdrawtokeep} (e.g., \code{nthdrawtokeep = 2} leads to deletion of every other draw) and \code{signresfirstatt}, respectively. \code{compresslist} is just for reducing the size of the output.

Following the findings of 2019_Hein_et_al, we leave the prior specifications at default, which is why there are no respective arguments. Here, too, the user is invited to make use of the package's or other software's full potential by importing her own data into upcoming functions. Compared to \code{C_} and \code{D_}, no external data preparation is included in \code{E_}. As with all the functions in general, the user may look at and mimic the structure of the outputs of \code{E_}/inputs of the upcoming functions in such a case.

Output and example \newline \code{E_hbmxlest()} only outputs the .rds file \setlist{nolistsep}

which is a \code{list} carrying (among other information) the model tensor in \code{$betadraw}.

As an example, estimating a model given the (training) design from \code{C_} and the responses from \code{D_} can be achieved through

Based on the 50 respondents $\times$ 8 parameters defined in \code{B_} and the Markov chain length of 20000, the model tensor is of size 50 $\times$ 8 $\times$ 20000.

Looking at the 19860th and 19960th draw here (the two matrices below) for the first two respondents (rows) and first four parameters (columns, price coefficients), we can see that in the 19860th draw (top matrix) the first respondent's (first row) (unconstrained) utility for the third price level (second column) violates the monotonicity introduced to the true preferences in \code{B_}. This issue will be addressed in (ref).

If \code{$betadraw} is required to hold more than four billion elements, the \pkg{RcppArmadillo} error \code{Cube::init(): requested size is too large; suggest to enable ARMA_64BIT_WORD} \newline might show up. To resolve this, create .R/Makevars (if not already present) in the home directory, add the compiler flag \code{CXXFLAGS = -DARMA_64BIT_WORD=1} to the Makevars file and reinstall \pkg{bayesm}.

Design assessment

Theoretical background

In (ref), we already mentioned that the design assessment is an exception to our naming convention reflecting the order within the ensemble of functions. The primary reason for this is to be able to presuppose (ref)'s theoretical underpinnings here as well, thereby avoiding redundancy in explanation.

Recall the three well-known principles of efficient choice designs from (ref). We assess them by computing the correlation matrix for each design (orthogonality), the marginal frequencies of the feature levels in each design (level balance) and the % of equal feature levels for each pair of alternatives within each set within the training and the test designs (minimal overlap). As promised, we also calculate the % of equal alternatives for each pair of sets within and between the training and the test designs, allowing to check if there exist any set duplicates.

Furthermore, a hierarchical Bayesian mixed logit model (see (ref)) is estimated based on random responses, which should lead to parameter estimates close to zero for an efficient choice design. As a recovery measure, we determine the root mean squared error (RMSE) between the true (zero) and re-estimated individual part-worths on parameter level for each draw as well as across parameters for each draw. Averaging them across draws then gives us two different ARMSE results (parameter level values, total value). The same is done with the mean squared error (MSE) and the mean absolute error (MAE).

The statistical criterion that was optimized for in the design generation is not re-calculated, as it is usually provided there.

Practical use

Function and input \newline The corresponding function from \pkg{cash} is \code{F_designeval()}:

The designs are loaded through the first two arguments (if the designs have been replicated and response columns have been added to prepare external data in \code{D_} for model estimation but the design evaluation is of interest although the choice experiment must have already taken place, \code{$design} has to be changed to single designs without responses). \code{burninpercentage} is for discarding burn-in iterations of the Markov chain before calculating the recovery measures. Model estimation for the test design based on random responses might not be sensible, as it is not used for training, and preventing it by setting \code{recoverycheckfortestdesign = FALSE} may even be required to avoid Cholesky decomposition errors in the Metropolis-Hastings step if the matrix has been altered (e.g., by keeping only few sets). The other arguments are identical to \code{E_} (see (ref)). \newline

Output and example \newline Like \code{E_}, \code{F_designeval()} only outputs a single .rds file, namely \setlist{nolistsep}

which is a \code{list} containing all the information described above.

Given the defaults and available outputs of previous functions, the call may just be

In view of the foregoing definitions, the list elements are mostly self-explanatory, but let us peek at the evaluation of our training design. It exhibits near orthogonality (lower left/upper right quadrant of correlation matrix shows small linear dependencies), perfect level balance (equal marginal frequencies) and perfect minimal overlap (\code{max overlap within sets} is zero). As \code{max overlap between sets} is 60%, there are no set duplicates. It is important to mention that, in the field of discrete choice analysis, the term (minimal) overlap is usually connected to the comparison of alternatives within sets (i.e., the assessment of feature level equality between product configurations). Here, we also use overlap to refer to the equality of alternatives between sets (necessary to identify set duplicates).

Note that the complete \code{overlap within sets} matrices list every comparison of alternatives independent of the resulting % of equal feature levels (i.e., also if two alternatives show 0% overlap with regard to their feature levels as seen below).

In contrast, the \code{overlap between sets} matrices (check for duplicate sets) only show the comparisons leading to at least one match of two alternatives. If there are multiple equal alternatives between two sets, all matches except for the last one carry \code{NA} in the column \code{overlap}. For example, row 2 of the \code{overlap between sets within training} matrix below expresses that \code{set 1 alt 1} is equal to \code{set 4 alt 1} (ignore the numbers in the column header, as they only serve to indicate that each row is a comparison of two alternatives from two sets), but because row 3 provides that \code{set 1 alt 3} equal to \code{set 4 alt 4}, row 2 carries \code{NA} and row 3 0.4. 2 out of 5 alts are equal between those two sets, i.e., 40%.

Convergence and model assessment

Theoretical background

Guided by the recommendations of 2013_Gelman_et_al regarding the diagnosis of convergence, we simulate an independent second chain with identical settings except for the seed, discard burn-in iterations, split each chain in half, and assess mixing and stationarity by computation of the uni- and multivariate potential scale reduction factor ((M)PSRF) 1992_Gelman_Rubin, 1998_Brooks_Gelman using the within- and between-sequence variances 2006_Plummer_et_al. If numerous experimental conditions are to be studied, we suggest to reverse-engineer a conservative number of iterations to hold constant for burn-in and assessment of convergence based on the most demanding treatments 2025_Dressler_et_al.

When monotonicity constraints are imposed on the synthetic utilities, it may be reasonable to remove violating draws (incorrect order and signs, see (ref)) from the primary chain, and if the lowest remaining number across individuals is insufficient, to estimate the necessary primary chain length by extrapolation. For each individual, the required number of acceptable draws is then collected starting from the end of the (potentially new) primary chain. Finally, to reduce serial correlation, the draws should be thinned 2009_Train,2013_Gelman_et_al independently of the removal 2000_Orme_Baker.

In line with the Bayesian estimation, the resulting $\mathcal{i}\times\mathcal{o}\times\mathcal{n}$ tensor of draws $n\in N=\{1,...,\mathcal{n}\}$ is validated by computing credible intervals for common measures of parameter recovery and predictive accuracy. If the credibility level $1-\alpha$ points at decimal indices for the limits, the weighted mean of the two respective values of the measure is taken at both limits (i.e., at the lower limit, $(\frac{\alpha}{2}\mathcal{n})\,\mathrm{mod}\,1$ and $1-(\frac{\alpha}{2}\mathcal{n})\,\mathrm{mod}\,1$ as weights for the lower and upper value, respectively). In the following, we restrict the mathematical representations to the draws, but the measures are provided for point estimates (posterior means), too.

To assess the parameter recovery, for each draw we calculate the RMSE between the true and re-estimated individual part-worths, $\beta_{io}$ and $\hat{\beta}_{nio}$, as well as the average Pearson correlation across the $\mathcal{i}$ individual part-worth vector pairs, $\bm{\upbeta}_{i}$ and $\hat{\bm{\upbeta}}_{ni}$ 1996_Vriens_et_al, 2002a_Andrews_et_al,2002b_Andrews_et_al,2019_Hein_et_al,2020_Hein_et_al, 2024_Goeken_et_al. Following 2002a_Andrews_et_al, we assume that the Gumbel's scale $s$, known to be inextricably included in $\hat{\bm{\upbeta}}_{ni}$ as divisor of the unscaled part-worths, is estimated correctly as the square root of the factor by which the standard variance of $\frac{\pi^2}{6}$ is scaled during the error generation in the response simulation (see (ref) in (ref)) and can be cancelled out by multiplying the estimated part-worths with the true $s$ such that\footnote{To be precise, 2002a_Andrews_et_al introduced the scaling to the simulated part-worths and re-scaled the RMSE.}

In regard to the Pearson correlation, which may be the more appropriate measure of parameter recovery given its scaling invariance, we think it is quite desirable that the average of coefficients tends to underestimate the true correlation. Monte Carlo studies also show that there is a risk of inflating positive bias further by employing the Fisher transformation for correction 2015_Bishara_Hittner. In contrast to previous papers, we therefore calculate the mean conservatively as

The out-of-sample predictive accuracy is evaluated through an individual as well as an aggregate measure for each draw, namely the hit rate across hold-out sets averaged over individuals (Equation (ref)) and the shares of (first) choice RMSE over alternatives and hold-out sets (Equation (ref)) 1996_Vriens_et_al,2002a_Andrews_et_al,2002b_Andrews_et_al, 2019_Hein_et_al,2020_Hein_et_al,2022_Hein_et_al,2024_Goeken_et_al.

The rationale behind the use of the first choice rule here is rooted in the central part of the methodological framework, i.e., (ref). Given the broader scope of required considerations in comparison to the logit choice rule, it is just more suitable for the technical explanations, but analogous to the parameter sets (draws, means, simulated/true), the $\mathrm{RMSE}^{\mathrm{soc}}$ is also reported under the logit rule (Equation (ref)).

In the above, the focus is not only limited to the draws and first choice rule but generally variants of the measures favored by us. To be exhaustive, we provide \setlist{nolistsep}

Practical use

Function and input \newline The corresponding function from \pkg{cash} is \code{G_modeleval()}:

The default imports of previous outputs should be clear at this point. If an internal model based on external data or an external model is to be evaluated, set \code{experimental = TRUE}, ignore \code{simulatedindividualbetas} (and therefore recovery), let the .rds for \code{trainingdesign} and \code{testdesign} be identical to the ones generated by \code{C_} except that \code{$design} must be a matrix containing the designs of all respondents with responses in the last column (as explained in \Cref{sec:responsesim}). The \code{...simuresponses} arguments then have to be .rds comprising just those matrices (see \Cref{sec:responsesim} as well). The .rds for \code{bayesmdata} may be generated with \code{D_} for external data, and concerning \code{model} (if external), mimic the output of \code{E_} (\code{$nmix} and \code{$loglike} can be disregarded). For model evaluation of point estimates only, the matrix has to be duplicated to hand over the tensor \code{G_} expects from \code{E_} in \code{$betadraw}.

The next four arguments refer to convergence. If the value of \code{burninpercentage} is $<1$, it is interpreted as a percentage, and if it is $\geq1$, it is interpreted as an absolute number of burn-in iterations to be discarded. \code{chainlengthpsrf} serves to define the absolute number of iterations after burn-in to be used for computation of the (M)PSRF and therefore to avoid calculating the latter based on potentially divergent chain lengths across multiple experimental conditions. Although not recommended, \code{secondmodelforgelman = FALSE} allows to skip simulation of an independent second chain, leading to the primary chain being split into two (artificial) chains. If an independent second (real) chain shall be simulated, the respective seed can be given to \code{seedsecondmodel}. In the latter case, \code{G_} always provides the (M)PSRF for split and unsplit chains. If \code{saveplot = TRUE}, the last respondent's primary chain trace plots for each parameter are saved.

\code{removereversals = TRUE} and \code{removepositives = TRUE} result in the removal of draws (remaining after burn-in) that violate price monotonicity (incorrect order and signs). Note that \code{removepositives} is not considered if \code{removereversals = FALSE}. The desired number of draws after removal before thinning has to be specified via \code{drawspostremoval}. If this threshold is not met, \code{G_} stops to execute and tells the lowest remaining number across individuals. To resolve this, either \code{drawspostremoval} has to be set to this number, a new primary chain has to be simulated with higher \code{mcmciterations} (see (ref)) or \code{burninpercentage} has to be reduced. If the threshold is met, the draws get truncated to \code{drawspostremoval} and thinned according to \code{thinning} (as with \code{nthdrawtokeep} in \code{E_}, e.g., \code{thinning = 2} leads to deletion of every other draw).

\code{alphanormality} and \code{crediblelevel} take the significance level for the frequentist normality tests and the probability coverage of the Bayesian credible intervals for the parameter recovery and predictive accuracy measures, respectively.

Besides the trace plots mentioned above, \code{saveplot} is somewhat similar to \code{saveplot} from \code{B_}, but in contrast to the simulated preferences, the estimated ones are available as posterior draws and means. Therefore, \code{saveplot = TRUE} not only produces a very basic \code{ggplot} of the final posterior means with identical settings to the graphic from \code{B_} to facilitate their comparison but also visualizes the draws. For each parameter, the population's distribution is shown for a selected number of the final $\mathcal{n}$ draws in equal distance (controlled by \code{show}, e.g., every 10th draw if \code{show = 50} and $\mathcal{n} = 500$), as well as the draws' distribution for a selected number of respondents. \newline

Output and example \newline \code{G_modeleval()} outputs the three .rds files \setlist{nolistsep}

The first one is a \code{list} carrying the assessment results and the other two are the final preferences arrays to work with.

As an example, let us discard half of the Markov chain (10000 out of 20000 iterations specified in \code{E_}), use the full other half for the convergence checks, set our minimum for the remaining draws after removal of draws violating price monotonicity to 500 and the thinning to 5 (every 5th draw is kept).

This leaves us with a tensor of size 50 $\times$ 8 $\times$ 100 for the draws.

As in \code{F_designeval.rds}, the list elements of \code{G_modeleval.rds} are self-explanatory given the information above. Therefore, we will just pick a few of them to assess our model, leaving the rest for interested users to explore.

With regard to convergence, the average MPSRF is close enough to perfect convergence ($=1$) for our purposes here 2013_Gelman_et_al.

The highest number of draws per respondent violating price monotonicity is 9077 out of 10000, which means if we had specified \code{drawspostremoval} > 923, \code{G_} would have stopped execution and told us this number via error message.

Concerning recovery and predictive accuracy, let us take a look at the 95% credible intervals of the last three measures for which the mathematical representations have been shown in the theoretical background.

To prove that these results indicate good performance of the final model, it is sensible (aside from theoretical baselines of measures and experience) to compare them with results based on a random preference tensor (or any other model of your choice). Although not vital here, note that whenever \code{removereversals = FALSE} (which has to be the case because if we removed reversals, we would have to create a non-violating (non-random) tensor to avoid errors), \code{drawspostremoval} is not considered and \code{thinning} has to be increased to still be left with the same target number of final draws. In our case, \code{drawspostremoval = 500} and \code{thinning = 5} lead to $\mathcal{n}=100$. To get the same $\mathcal{n}$ with \code{removereversals = FALSE}, specify \code{thinning = 100} to thin 10000 draws after burn-in down to 100.

The credible intervals of the correlation and hit rate span their baseline expectations of 0 and 20% (5 alternatives), respectively, given random parameters (i.e., random answers). In contrast to the correlation and hit rate, for which a theoretical baseline always exists independent of the actual preferences, the baseline of the shares of choice error is highly dependent on the true distribution of choices (e.g., error will be low if the alternatives are chosen equally frequent by the respondents because that is also the result of random guessing). As with the first two measures, the shares of choice error improves drastically when using the original model.

We now have to re-run \code{G_} again with the original model to re-overwrite \code{G_}'s outputs which are currently based on the random tensor.

Pre-computations and Nash competition

Theoretical background

Before undertaking game-theoretic simulations based on the estimates from (ref), it is necessary to define the competing firms' objective. Following the seminal conjoint-based research on non-cooperative competitive reactions from the long-run perspective 1993_Choi_DeSarbo, 1997_Green_Krieger, a firm $w\in W=\{1,...,\mathcal{w}\}$ is assumed to search for a product (line) $a\in A=\{1,...,\mathcal{a}\}$ that maximizes the total contribution margin $\pi_{wak^-}$ given the partial (excluding $w$) competitive scenario $k^-\in K^-=\{1,...,\mathcal{k}^-\}$ by varying price and design (i.e., non-price feature levels) simultaneously.

To mathematically formulate and extend this optimization problem, also in view of the more recent works in this field 2025_Dressler_et_al, let $J$ now be the set of indices for products in the complete (including $w$) competitive scenario $k$ (equivalent to a choice set before), let $q\in Q=\{1,...,\mathcal{q}\}$ be the product index in the optimizing firm's line such that $Q\subset J$, and take $|Q|(\geq1)$ as exogeneously fixed. Furthermore, let $\bm{\mathrm{p}}$ be a $1\times\mathcal{m}$ price vector, and let $\bm{\mathrm{c}}$ be a $1\times\mathcal{m}(\mathcal{l}-1)$ vector that contains the cost of the non-price feature levels. Because the row vector $\bm{\mathrm{x}}_{(\cdot)}$ represents a complete product configuration (see (ref)), we are able to describe price and design separately by splitting $\bm{\mathrm{x}}_{(\cdot)}$ into the subvectors $\bm{\mathrm{x}}_{(\cdot)p}$ and $\bm{\mathrm{x}}_{(\cdot)c}$, respectively. If we expand $\bm{\mathrm{x}}_{(\cdot)}$ and $\hat{\bm{\upbeta}}_{(\cdot)}$ by the reference category of each feature beforehand, i.e., change the coding scheme of $\bm{\mathrm{x}}_{(\cdot)}$ and add zeros to $\hat{\bm{\upbeta}}_{(\cdot)}$, the optimization problem can be written as

For each product $q$, which has to differ from the other (if existing) products $q^{*}\in Q(q^*\neq q)$ in at least one feature (Equation (ref)), the demand is calculated by taking the sum of first choices $\hat{f}_{niq}$ (Equation (ref)) over individuals $i$ and averaging it across draws $n$. This simulation almost surely converges to the expected demand, which is the integral over the distribution of individual preferences $\bm{\upbeta}$ and the posterior of hyperparameters $\Theta$ with regard to the first choices (scaled by our volume $\mathcal{i}$ here). Compared with logit probabilities, the first choices are scaling invariant (i.e., not influenced by the scale factor $s^{-1}$ of the part-worths) and immune to the property of independence of irrelevant alternatives (i.e., not prone to share inflation for similar products). Their possible downside of unrealistic determinism 1979_Shocker_Srinivasan is reduced by implicitly obtaining the individual-level share of draws with maximum utility for the respective product (note that $\frac{1}{\mathcal{n}}\sum_{n,i}\hat{f}_{niq}$ is the same as $\sum_{i}(\frac{1}{\mathcal{n}}\sum_{n}\hat{f}_{niq})$). However, $\hat{f}_{niq}$ and $\hat{\bm{\upbeta}}_{ni}$ may be replaced to test the effect of various combinations between choice rules (first and logit choice) and parameter sets (draws vs. point estimates (posterior means) vs. true preferences) on the optimization and equilibria outcomes. We could as well incorporate the preference uncertainty by calculating the demand (and optimizing the objective function) for every draw of the hyperparameters,

and subsequently simulating the equilibria for every draw like 2014_Allenby_et_al to build up posterior distributions of equilibria. Though, in our case, the equilibrium quantity is multidimensional 2014_Allenby_et_al, which greatly restricts the manageability of such distributions with respect to interpretability and comparability.

If multiple first choices are present within a complete competitive scenario, indicated by the cardinality of set $S$, here, the 100% probability is equally divided between them (Equations (ref), (ref)). Despite being identical with high enough frequency, sampling produces slightly different total contribution margins for complete competitive scenarios consisting of the same set of winning products, which may be unsuitable for the structural analysis of equilibria (the commonness of recurrences and ties will become apparent in the following). It is worth mentioning that we do not (need to) implement such a tie-breaking strategy in the response simulation and the model assessment, as a tie can only arise there in the highly improbable case of two distinct products showing the identical and at the same time largest total utility (safety mechanism: if this happened in the response simulation, \pkg{bayesm} would throw the error \code{nrow(X) ne p*length(yi); exception at unit...}, telling us that there are more choices than sets). Furthermore, we do not restrict competition by forcing a lower bound on the share of choice 2012_Kuzmanovic_Martic.

The expected demand is then multiplied with the corresponding contribution margin of a single unit, which is computed by subtracting $q$'s cost $\bm{\mathrm{c}}\bm{\mathrm{x}}_{qc}^T$ as well as a scalar $\delta$ from $q$'s price $\bm{\mathrm{p}}\bm{\mathrm{x}}_{qp}^T$ 2019_Kuzmanovic_et_al. $\delta$ should not only comprise a base cost term but also the cost for features that are assumed to be unchangeable from the firms' perspective or irrelevant for the consumers' choice. Note that the accuracy of the (relative) cost structure most definitely is a critical factor affecting the validity of the simulation outcomes 1990_Choi_et_al,1994_Choi_DeSarbo. The marginal cost of production is assumed to be constant 2014_Allenby_et_al, as is the cost of repositioning.

Due to the discrete domain (i.e., set of binary integers) on which the nonlinear objective function (Equation (ref)) is defined, the optimization problem is of combinatorial nature and a solution cannot be derived analytically. The number of theoretically possible product configurations $\tau=\mathcal{m}^{\mathcal{l}}$ (continuing the example of symmetric designs) and the number of theoretically possible product line configurations $\mathcal{a}=\binom{\tau}{\mathcal{q}}$ grow exponentially with $\mathcal{l}$ and $\mathcal{q}$, respectively. As is the case with numerous combinatorial optimization problems, there also exists no exact numeric algorithm capable of solving it in polynomial time, making it NP-hard 1990_Kohli_Sukumar. In line with the effort of eliminating uncontrolled systematic influences, we nevertheless optimize using complete enumeration such that the structural properties of the equilibria can be fully captured and are not biased by artefacts from heuristics.

Apart from that, complete enumeration's time complexity does not always have to be disadvantageous, which is why we also refrain from implementing other exact methods for now. Compared with procedures guaranteeing a solution's global optimality, even the full potential of heuristics 2008_Belloni_et_al,2024_Baier_Voekler in terms of runtime superiority does not come into play when the solution space $A$ of the optimization problem is rather small. In the simulation of competitive reactions the latter can be computationally very limited if there is a large number of optimization problems to be solved consecutively. This will be elaborated upon in the subsequent.

Given the discrete domain of the objective function, our game-theoretic solution concept of interest cannot be derived analytically either. Hence, we simulate dynamic closed-loop games (i.e., multi-stage games with mutually observable past actions) of myopic best responses in a sequential manner to obtain the fixed points known as pure strategy Nash equilibria 1838_Cournot,1951_Nash, 1991_Fudenberg_Tirole.

More precisely, and again closely following 1993_Choi_DeSarbo and 1997_Green_Krieger as well as 1995_Gutsche and 2000_Steiner_Hruschka, firms take turns in maximizing their total contribution margin $\pi_{wak^-}$ depending on the others' product lines until no firm can benefit from unilaterally changing its product line. In accordance with the aforementioned articles, the competitors are all assumed to be active and to be symmetric in regard to prices, cost structure, estimated consumer preferences and number of products. Moreover, the number of competitors is expected to remain constant throughout a game 2014_Allenby_et_al.

If the sequence of each firm optimizing once,

is denoted a round $b \in B=\{1,...,\mathcal{b}\}$ and $k^-_0\in K^-$ is the index of the initial state (i.e., the partial competitive scenario at the beginning of the $\mathrm{t\hat{a}tonnement}$), here, a Nash equilibrium can be formally expressed as the singleton

$\mathcal{T}^{k^-_0}$ contains the index $k^{k^-_0}_{b}$ referencing the complete competitive scenario $k$ which is present at the end of round $b$ and (first) remained unchanged for two consecutive rounds after starting from $k^-_0$. It is crucial to set an upper limit for the rounds in order to prevent a game from running infinitely in the absence of an equilibrium, especially if no comprehensive detection mechanism is implemented that checks for all different types of cycles. The latter requires the (partial) comparison between the current $k_{b}$ and each of the preceding $k_{b'}$ with $b' \in B'=\{1,...,b-2\}$. We decided to just look for the shortest possible cycle (2-round cycle) and this only once when the upper limit of the rounds is reached (this is more efficient if the expected number of cyclic games is rather low because then this examination is not even triggered once in the majority of games),

On the other hand, at least two rounds have to be played to see if $\mathcal{T}^{k^-_0}\neq\emptyset$ because regardless of the order of movement of the firms, there is no initial state $k_0$ to meaningfully compare $k_1$ with. In other words, a computational dependence on an initial product line configuration of the firm that comes first in the reaction sequence does not exist (and we do not choose one at random), as each game starts with this firm optimizing over its $\mathcal{a}$ possibilities.

Thus, there are at most $\mathcal{k}^-$ initial states ($\mathcal{a}^\mathcal{w-1}$ theoretically possible competitive scenarios of $\mathcal{w}-1$ firms) to start a game from, which we exhaustively go through for two reasons. Firstly, the effect of each unique initial competitive scenario as well as its order of movement variants (thanks to symmetric competitors) on the equilibria can be observed, and secondly, every existing equilibrium is guaranteed to be found since they are inevitably represented in the initial states.

It is noteworthy that with complete enumeration the firm reacting first already needs to calculate $\pi_{wak^-}$ once for all $\mathcal{k}$ (i.e., $\mathcal{a}^\mathcal{w}$) theoretically possible complete competitive scenarios when optimizing over its $\mathcal{a}$ product line configurations at the start of each of the $\mathcal{k}^-$ games. If enough memory is available, there are major advantages to pre-computing a matrix

\AtBeginEnvironment{bNiceMatrix}{\everymath}

comprising the $\mathcal{k}$ complete competitive scenarios as well as the corresponding total contribution margins $\pi_{wk}\forall k$ (i.e., $\pi_{wak^-}\forall a,k^-$) from an arbitrary but constant viewpoint $w$ (e.g., $w=1$). After pre-optimizing $\mathcal{k}^-$-times over $\mathcal{a}$ product line configurations in $\bm{\mathcal{M}}$, the resulting matrix

can easily be utilized as a look-up table for the best response $\operatorname*{arg\,max}_a\pi_{wak^-}$ to a given partial competitive scenario $k^-$. Recalling the property of symmetry, it is evident that $\bm{\mathcal{M}}^{\mathrm{opt}}$ allows to circumvent the repeated calculation of identical total contribution margins and optima across firms and rounds in the games. Consequently, the runtime of the $\mathcal{k}^-$ games to be played in an experimental condition becomes neglectable. Since increasing compute is usually less feasible than increasing memory, we chose to take this memory-heavy path.\footnote{For the large-scale Monte Carlo study in 2025_Dressler_et_al, we used three rack servers with 104 physical cores and 2,560 GiB RAM in total: 1x Dell\textsuperscript{\texttrademark} PowerEdge\textsuperscript{\texttrademark} R450 with 2x Intel\textsuperscript{\textregistered} Xeon\textsuperscript{\textregistered} Silver 4316 CPUs and 16x 64 GiB DDR4 3200 MT/s RDIMMs, 2x Dell\textsuperscript{\texttrademark} PowerEdge\textsuperscript{\texttrademark} R440s with 2x Intel\textsuperscript{\textregistered} Xeon\textsuperscript{\textregistered} Silver 4216 CPUs and 12x 64 GiB DDR4 3200 MT/s RDIMMs each.}

As implicitly stated above, the $\mathcal{k}^-$ optimization problems to be solved ($\mathcal{k}^-\mathcal{w}\mathcal{b}$ without pre-computation) pose the main computational bottleneck. Due to the exponential growth with $\mathcal{w}$, the solution space $A$ of a single optimization problem can be so small that both heuristics and exact methods are comparably fast but the vast initial state space $K^-$ precludes computational feasibility. Thus, if $A$ calls for heuristics, $K^-$ certainly will too. The primary way of restoring feasibility therefore is to prune $K^-$. To the best of our knowledge, such an approach has not yet been developed, but even if it were to exist, it would only be employed here with the presence of a mathematical proof demonstrating the equivalence of results.

Apart from the latter, the runtime of the pre-computation of $\bm{\mathcal{M}}$ is substantially decreased by implementing the workhorse functions in C++ 2026_Eddelbuettel_et_al, and fortunately, the task is embarrassingly parallel (for reproducible parallel computing, we use the \proglang{R} packages \pkg{doParallel} 2022a_Microsoft_Weston, \pkg{foreach} 2022b_Microsoft_Weston and \pkg{doRNG} 2026_Gaujoux). Additionally, we pre-compute the possible line and product configurations of a firm, the $\tau\times\mathcal{i}\times\mathcal{n}$ tensor of (exponentiated) product utilities as well as the $\tau\times 1$ vector of product contribution margins to serve as look-ups 2008_Belloni_et_al. This avoids unnecessary re-computations also within the main pre-computation process of $\bm{\mathcal{M}}$ itself. With regard to memory, the number of elements in $\bm{\mathcal{M}}^{({\mathrm{opt}})}$ is minimized by mapping the pre-computed extended line and product configurations of each firm to a single integer analogous to (ref).

We were given the opportunity to test our implementation (modified to GPU computing) on a blade server of an exascale supercomputer currently under development. Note that an extrapolation of the results indicated that even on entire machines leading the TOP500 list, simulation of many interesting, yet moderate scenarios would still be far out of reach, which casts a different light on computational limitations the user might encounter on usual machines.

Practical use

Function and input \newline The corresponding functions from \pkg{cash} are \code{H_precomputeM()}:

and \code{I_preoptandnashgame()}:

As its name reveals, \code{H_precomputeM()} is responsible for the pre-computation of $\bm{\mathcal{M}}$, which makes it the computationally most demanding function in \pkg{cash} (depending on the settings, of course). If it is to be used independently (i.e., with external data), it may be helpful to know that from \code{C_trainingdesign.rds} it only extracts \code{$resps}, \code{$betas} and \code{$lvls}.

Its argument \code{estimatedindividualbetas} takes the final preference array of draws, point estimates/posterior means or simulated/true utilities.

Price and cost levels have to be imported from two .csv files through \code{pricecostmatrix} (a .csv only comprising the features that have been defined in \code{B_} as being part of the conjoint choice experiment) and \code{pricecostmatrixfull} (a .csv containing all features relevant to the contribution margin). These .csv have to have row names as first column, column names as first row, as many rows as there are features + 1, as many columns as there are (max) levels + 1, and price as first feature. Undefined cells (if the number of levels differs across features) must be \code{NA}, and a zero indicates that the level has no cost. Obviously, the order of the features and levels has to match the order of the estimates in the preference array, e.g., the price levels have to be sorted in ascending order if a monotonicity constraint was introduced to the price feature in \code{B_} and enforced in \code{G_}. Via \code{basecostfix}, the fixed part of the base cost term $\delta$ can be set, while \code{basecostvar = TRUE} leads to the levels of the features included in the full but not partial price cost matrix being randomly sampled to then add their corresponding costs to $\delta$. If $\delta$ shall only be a fixed term ($\geq0$), the same .csv can be given to both price cost matrix arguments, but \code{basecostvar} must be \code{FALSE}.

If \code{"first"} as the default of \code{choicerule} is to be applied, \code{ruletosolvefcutilitytie} allows to choose between \code{"splitting"} and \code{"sampling"} as a tie-breaking strategy for the case of multiple first choices being present within a complete competitive scenario. \code{prodsperline} and \code{competitors} both have to be given an integer specifying the number of products per line and number of competitors, respectively.

\code{addnullparameters = FALSE} suppresses the expansion of the preference array by zeros for the reference category, which may be necessary in case of an external preference array. The progress of the pre-computation of $\bm{\mathcal{M}}$ is saved to a .txt file. Its update frequency can be controlled via a numeric value (minutes) given to \code{progressupdatedelay}.

\code{I_preoptandnashgame()}, on the other hand, pre-computes $\bm{\mathcal{M}}^{{\mathrm{opt}}}$ and uses it for fast simulation of the Nash games.

To do so, it needs $\bm{\mathcal{M}}$ from \code{H_}, as well as additional information to, e.g., convert the single-integer line-identifiers back to readable product feature levels. The corresponding .rds files for these first two arguments (\code{scenariosmatrix} and -\code{info}) differ in results and name depending on the choice rule and parameter settings (so, the defaults may need to be changed). The other two, which contain the pre-computed product and line configuration matrices, do not. \code{maxrounds} defines the the upper limit $\mathcal{b}$ of the game rounds and has to be an even integer (for internal reasons). \newline

Output and example \newline \code{H_precomputeM()} outputs the five .rds files \setlist{nolistsep}

We have already explained the first four above, as they serve as inputs for \code{I_}. The last one holds the first 10 rows of $\bm{\mathcal{M}}$ for quick checks.

\code{I_preoptandnashgame()} only outputs the .rds file \setlist{nolistsep}

which includes all relevant information about the equilibrium simulations. Note that the terminal \code{A-I_finaloutput.rds} also gets this choice rule and parameter set name extension.

Both functions additionally generate a similar \code{[]_}progress.txt file.

As an example, we let \code {H_} import the final preference array of draws from \code{G_}, which is the default, and the following two .csv, pcm2atts.csv and pcm6atts.csv, with the price and cost levels for \code{pricecostmatrix} (recall that we defined 2 features with 5 levels each in \code{B_}) and \code{pricecostmatrixfull} (6 features with 5 levels each), respectively:

These price and cost levels (in €) are taken from 2025_Dressler_et_al, where we were able to gain insight into the (relative) cost structure of a well-known computer manufacturer and modeled competition between firms in this sector. In addition to price, we determined five modifiable, discriminating features driving the consumers' choice and the firms' contribution margin for notebooks, namely display size, central processing unit (CPU), solid state drive (SSD) capacity, battery life and random access memory (RAM). For the purposes of our example, further elaboration on the study's details is not necessary. Specific feature levels will be provided as needed.

Let us leave the two base cost arguments at preset, which means ignoring the fixed part of $\delta$ (no cost for, e.g., housing and mainboard, and no explicitly chosen cost levels for CPU, SSD capacity, battery life and RAM) but sampling cost levels for the CPU, SSD capacity, battery life and RAM from the full price cost matrix. Let us also use the defaults for the choice rule, i.e., first choice, as well as the tie-breaking strategy, i.e., splitting, and simulate competition between 2 firms, each offering a product line consisting of 2 products. \code{progressupdatedelay} is set to a few seconds to force H\code{_}progress.txt to be generated.

As the .rds defaults in \code{I_} correspond to the first choice and draws settings employed here in \code{H_}, there is no need to adjust any file names when calling \code{I_} in this case.

Because \code{I_equilibriainfo[][].rds} comprises all contents of \code{H_scenariosinfo[][].rds} and more, we will skip the latter and jump directly to the former and its most important components. But before we do so, it seems sensible to practically explain the structure of $\bm{\mathcal{M}}$ and the progress file by taking a look at \code{H_scenariosmatrixhead[][].rds} and H\code{_}progress.txt.

Like (ref) showed theoretically, $\bm{\mathcal{M}}$ consists of as many columns as there are competitors plus one for the contribution margin. The integers in the competitors' columns refer to rows in the pre-computed line configuration matrix which contain integers referring to rows in the pre-computed product configuration matrix which contain interpretable feature levels. The unit of the contribution margin is defined by the unit of the price cost matrix entries.

The progress file reports the following information:

If the specified delay is very short, some oscillation of the time estimates is normal due to performance fluctuation across cores in parallelization.

Now, we finally arrive at what we have been working towards throughout this paper and package, namely the output of \code{I_}. Its \code{list} begins with the provision of a few metrics describing $\bm{\mathcal{M}}$ and $\bm{\mathcal{M}}^{{\mathrm{opt}}}$. Having defined 2 features with 5 levels each, the number of theoretically possible product configurations is 25.

By choosing competition on a product line with 2 products, we get 300 theoretically possible line configurations ($\binom{25}{2}$, i.e., 25 $\times$ 25 minus 25 lines of identical products, divided by two to get rid of the flips).

Specifying 2 competitors then leads to 90,000 theoretically possible complete competitive scenarios, i.e., 90,000 contribution margins to be pre-computed (rows of $\bm{\mathcal{M}}$) and 300 optimization problems to be pre-solved (rows of $\bm{\mathcal{M}}^{{\mathrm{opt}}}$).

Furthermore, we can access the sampled variable base cost levels (matching the order of the (excess) features in the price cost matrix).

Extending the \code{[]_}progress.txt, both functions output a run time and memory usage matrix. While for \code{I_}, it reports just the run time (mins) and object size (Mb) for different parts of the function, for \code{H_}, it also calculates the percentage of free system memory taken up (currently only on Windows and Linux).

With regard to the equilibria, \code{I_} presents information about them in three primary ways. As there are $\mathcal{k}^-$ initial states ($\mathcal{a}^\mathcal{w-1}$ theoretically possible (partial) competitive scenarios of $\mathcal{w}-1$ firms) to start a game from (here, 300), identical equilibria may arise multiple times. Such duplicates are generally only counted once per run (i.e., per function call). Flips, which are equilibria consisting of the same products or lines but swapped between firms (recall the competitors’ symmetry), are not treated as duplicates. Note that flips are only possible for differentiated equilibria, i.e., equilibria in which not all competitors share the same product (line).

In the following first matrix, the equilibria are listed in form of their single-integer line-identifiers, i.e., as rows of $\bm{\mathcal{M}}$. Identical equilibria are included once except if the number of rounds to reach them differs between initial states (there exists a reduced form in \code{$`unique equilibria in lines`} without rounds and therefore without duplicates due to round differences).

It also saves the initial states (instead of the last product lines present) from which no equilibrium has been reached for further investigation. If a game in a run does not end in an equilibrium, either more rounds have to be played or a cycle is present. As written in the theoretical background, we are only certain about the 2-round cycles. To distinguish between these two cases, the first competitor (first column) is \code{NA} if a 2-round cycle has been detected and zero otherwise (larger cycle present or round limit too low).

In our example, three equilibria have been found. \code{54-51} is a flip of \code{51-54}, and \code{130-130} shows up three times because it took a different number of rounds to reach it from some of the underlying initial states. There exist no initial states from which no equilibrium was reached. Otherwise, there would have been rows with \code{NA} or \code{0} in column \code{linecomp1} and corresponding initial states (instead of the last product lines present) in column \code{linecomp2}. \code{130-130} is the pareto-superior equilibrium, as it results in much higher contribution margins for both competitors.

The second matrix contains as many rows as there are (non-duplicate) equilibria plus as many rows as necessary to view each of those equilibria from the perspective of each competitor to calculate the corresponding contribution margins. In \code{130-130}, there is no differentiation between the competitors, which leaves us with the same contribution margin for both of them (rows 1, 2). In \code{51-54}, on the other hand, the two competitors differ in their product lines, leading to two different contribution margins (rows 3, 4). Note that the column names have to be ignored, as the lines of each equilibrium get exchanged to be able to always treat the first column as the competitor from whose perspective the contribution margin is calculated. E.g., in the original second equilibrium \code{51-54} (row 3), competitor 2 chooses line 54. In row 4, we act as if competitor 1 chooses line 54 to be able to calculate the contribution margin the same way. In row 4, competitor 1 is therefore to be interpreted as competitor 2. The two resulting rows for the flip of \code{51-54} may be ignored, as they do not provide any additional insight.

Ultimately, the third way of presenting the equilibria is a list of matrices in which the equilibria are converted back to a readable feature level format. Each matrix has as many rows as there are products per line times competitors and as many columns as there are product features. In our example, each equilibrium consists of 2 products with 2 features for each of the 2 competitors.

Through comparison of the numbers with the actual feature levels 2025_Dressler_et_al, we can conclude that in \code{130-130} (first matrix) both firms offer two notebooks at 599 € (second price level), one with a 14" and one with a 15" display (second and third display level).

In \code{51-54} and its flip (second and third matrix) both competitors offer one notebook at 299 € (first price level) with a 15" display (third display level). For the second product in their lines, both offer a notebook at 599 € (second price level), but one offers it with a 14" display (second display level) and the other with a 17" display (fifth display level).

Market share and revenue as special cases of profit optimization

Due to the fact that market share and revenue are both implicitly embedded in the objective function of the profit optimization problem presented herein (Equation (ref) in Section (ref)),

there is no need to alter the framework to be able to simulate corresponding equilibria. For competition on market share, prices and costs must not enter the objective function. If the price vector $\bm{\mathrm{p}}$ and the cost vector $\bm{\mathrm{c}}$ are replaced with the all-ones vector $\mathbb{1}$ and the null vector $\mathbb{0}$, respectively, and the base cost $\delta$ is set to zero as well,

each product's demand gets multiplied with the same contribution margin of one and therefore stays unchanged. For competition on revenue, on the other hand, only costs must not enter the objective function,

which results in each product's demand being multiplied with just price. To achieve that in \pkg{cash}, the entries of the price cost matrices have to be modified accordingly ((ref)), and, for the market share objective, the first feature can be treated like any other non-price feature throughout the framework.

Limitations and future development

In this paper, we illustrated the implementation and application details of the \proglang{R} package \pkg{cash}. To the best of our knowledge, \pkg{cash} is the first package to contribute a framework for (end-to-end) simulation of competitive reactions based on discrete choice analysis.

To close off, we would like to address \pkg{cash}'s primary limitations, which are congruent with potential avenues for future development. As in our view they would be worthwhile, we currently think about extensions to 2025_Dressler_et_al \setlist{nolistsep}

Finally, an approach for pruning the vast initial state space (see (ref)) without loss of equilibrium information would be a significant milestone, as the simulation of more complex scenarios might thereby come within reach.

Acknowledgments

The authors thank Philipp Aschersleben for his helpful comments on the pre-optimization described in (ref).