Fast and simple inner-loop algorithms of static / dynamic BLP estimations

Fast and simple inner-loop algorithms of static / dynamic BLP estimations

Introduction

Demand estimation is the basis for many economic analyses. The estimation method based on the random coefficient logit model proposed by berry1995automobile (henceforth BLP) has been widely used in the literature. Although berry1995automobile initially considered a static demand model, the model and estimation procedure have been extended to dynamic demand models, such as durable goods and goods with switching costs (e.g., gowrisankaran2012dynamics; shcherbakov2016measuring). These models are known as dynamic BLP. In these models, the products' mean utilities $\delta$ are numerically solved in the inner loop, and econometricians find the parameter values that minimize the GMM objective in the outer loop.

Nevertheless, significant challenge of this method is the computational time, as fixed-point iterations in the inner loop can be burdensome, especially for complex demand model estimations or large datasets. This study explores computationally efficient inner-loop algorithms for estimating both static and dynamic BLP models. While computationally efficient inner-loop algorithms for static BLP models have been extensively studied in the literature, including kalouptsidi2012market, reynaerts2012enhencing, conlon2020best, and pal2023comparing, there remains considerable potential for improvement in static BLP models, and even more so in dynamic BLP models.

Regarding static BLP models, I introduce a new fixed-point iteration mapping $\delta_{j}^{(n+1)}=\Phi_{j}^{\delta,\gamma=1}(\delta_{j}^{(n)})\equiv \delta_{j}^{(n)}+\left(\log\left(S_{j}^{(data)}\right)-\log\left(s_{j}(\delta^{(n)})\right)\right)-\left(\log\left(S_{0}^{(data)}\right)-\log\left(s_{0}(\delta^{(n)})\right)\right)$. This slightly modifies the traditional BLP contraction mapping $\delta_{j}^{(n+1)}=\delta_{j}^{(n)}+\left(\log\left(S_{j}^{(data)}\right)-\log\left(s_{j}(\delta^{(n)})\right)\right)$.\footnote{$S_{j}^{(data)}$ denotes product $j$'s observed market share, and $s_{j}(\delta)$ denotes the market share predicted by the structural model.} The difference with the traditional BLP contraction mapping is the term $\log\left(S_{0}^{(data)}\right)-\log\left(s_{0}(\delta^{(n)})\right)$, which is straightforward to implement in any programming language. Interestingly, when consumer heterogeneity is absent, the iteration converges after one iteration, regardless of initial values $\delta^{(0)}$. In such cases, $\delta_{j}=\log(S_{j}^{(data)})-\log(S_{0}^{(data)})$ holds, as shown by berry1994estimating. For the mapping $\Phi^{\delta,\gamma=1}$, $\Phi_{j}^{\delta,\gamma=1}\left(\delta\right)=\log\left(S_{j}^{(data)}\right)-\log\left(S_{0}^{(data)}\right)$ holds, and the output of $\Phi^{\delta,\gamma=1}$ is equal to the true $\delta$ for any input $\delta$. When the consumer heterogeneity exists, the iteration may not converge after one iteration. However, it still inherits the good convergence property. This study finds that the new mapping converges and reduces the number of iterations in the Monte Carlo simulation settings experimented in the previous studies (dube2012improving; lee2015computationally) and datasets from nevo2001measuring and Berry et al. berry1995automobile,berry1999voluntary. The speed-up is especially prominent when the consumer heterogeneity is relatively small, as theoretically shown in Appendix (ref).\footnote{Though there is no guarantee that the mapping is a contraction, we can guarantee the global convergence of the iterations using the new mapping by adding a few lines in the programming code, by using the fact that the traditional BLP contraction mapping is a contraction. For details, see Appendix (ref). If practitioners are conservative concerning the convergence, the procedure is worth considering.}

This study also compares various fixed-point iterations acceleration methods (Anderson acceleration, Spectral, SQUAREM), and finds that Anderson acceleration, which has been understudied in the literature\footnote{duch2023evaluating is the exception. As discussed in the online appendix of their paper, they applied the stabilized version of the Anderson acceleration method proposed by zhang2020globally to estimate a random coefficient nested logit model with a large nest parameter.}\footnote{conlon2020best pointed out in footnote 81 of their paper that “Anderson acceleration” using “anderson” function in Python Scipy package was too slow and unreliable to be worth considering. The current study verify it by using PyBLP. However, the “Anderson acceleration” method in Python Scipy package is aimed at solving general nonlinear equations, and the algorithm is not the same as the “Anderson acceleration” for fixed point iterations discussed in the current paper.}, performs better than the others, and further reduces both the number of iterations and computation time. The strategy remains effective even when the size of consumer heterogeneity is large. As the new mapping can be easily coded by adding a few lines to the traditional BLP contraction mapping, the current study suggests coding the new mapping first, and then considering combining the acceleration of fixed point iterations, especially Anderson acceleration, if needed. Anderson acceleration method is intuitive and relatively easy to code.

Regarding dynamic BLP models, I propose analytically representing the mean product utility $\delta$ as a function of value function $V$, solveing for $V$ by applying a mapping of $V$, and recovering $\delta$ using the analytical formula and $V$ values. Previous studies (e.g., gowrisankaran2012dynamics) applied an inner loop algorithm essentially consisted of solving for two types of variables $\delta$ and $V$. In the current algorithm, we only have to solve for only one type of variable $V$, thus reducing inner loop iterations. Combining fixed-point iteration acceleration methods, especially Anderson acceleration further reduces the number of iterations.

Fixed-point iteration acceleration methods are useful not only for the inner-loop algorithms of static/dynamic BLP estimations, but also for efficiently solving various fixed-point problems, including dynamic models with value functions and supply-side pricing equilibrium. Built-in packages of the acceleration methods, including Anderson acceleration, Spectral, and SQUAREM, are available in some programming languages. Even when unavailable, they are easy to code because of their simple formulas. Hence, being familializing oneself with accustomed to using such built-in packages or developing custom coding to implement the acceleration methods to meet practitioners' needs would be useful. Note that when we apply the spectral or SQUAREM algorithm, we should be careful about the choice of step sizes $\alpha$, though not paid much attention to in the previous economics literature. When we do not use appropriate step sizes, the algorithms may cause divergence. This is briefly discussed in Section (ref).

This study contributes to the literature by presenting key insights aimed at reducing inner-loop iterations in static/dynamic BLP estimations: (1). New mapping $\delta_{j}^{(n+1)}=\delta_{j}^{(n)}+\left(\log\left(S_{j}^{(data)}\right)-\log\left(s_{j}(\delta^{(n)})\right)\right)-\left(\log\left(S_{0}^{(data)}\right)-\log\left(s_{0}(\delta^{(n)})\right)\right)$; (2). Analytically represent mean product utilities $\delta$ as a function of value functions $V$ and solve for $V$ (for dynamic BLP models); (3). Combine the acceleration method of the fixed-point iterations, especially the Anderson acceleration. These methods are independent, and relatively easy to implement. Practitioners can selectively implement some when they face computational challenges.\footnote{fukasawa2024lightbulb applies some of these ideas to the dynamic demand estimation of the light bulb market by specifying a durable goods model with forward-looking consumers and replacement demand. The computation time is reduced by more than 10 times with these strategies. Without them, the expected estimation time would be more than 10 hours. }

The rest of this article is organized as follows. Section (ref) examines the relationships with the previous studies. Section (ref) discusses inner-loop algorithms for static BLP. Section (ref) explores inner-loop algorithms for dynamic BLP models. Readers focusing on static BLP models can skip this section. Section (ref) discusses the acceleration methods of fixed point iterations, including Anderson acceleration, Spectral, and SQUAREM. Note that the discussion in this section applies to any fixed-point iteration. Section (ref) presents the results of numerical experiments on static/dynamic BLP models. Finally, Section (ref) presents the conclusions.

Appendix (ref) shows additional results and discussions. Appendix (ref) contains all the proofs of the mathematical statements. The Supplemental Appendix shows further results and discussions. In the Supplemental Appendix, we discuss that the idea of the proposed inner-loop algorithm for estimating dynamic BLP models can be used in estimating dynamic discrete choice models with unobserved payoffs, which have been considered in kalouptsidi2020linear and others.

Literature

First, this study relates to the literature on computationally efficient estimation of static/dynamic BLP models.

conlon2020best developed the PyBLP package written in Python for implementing static BLP models' estimations and simulation. Although the package is highly convenient to use and state-of-the-art, incorporating the knowledge from existing literature, it is not suitable for BLP-type problems that PyBLP cannot address.\footnote{For instance, we cannot directly use PyBLP to estimate static limited consideration set models (e.g., goeree2008limited) or dynamic BLP models.} In such cases, we should write our own code, and the techniques discussed in the current study should be considered.

kalouptsidi2012market proposed inner-loop estimation algorithms of static BLP models that achieve faster convergence in the setting with a small number of discrete consumer types. doi2022simple proposed an estimation method for static BLP models with discrete consumer types that eliminates the need for computationally costly fixed-point iterations, given the total sales data availability for each consumer type. The main concept of these studies is that mean utilities or unobserved product characteristics can be analytically represented as a function of value functions specific to consumer types, given parameters and market share data.

The idea is beneficial, especially in the context of dynamic BLP models. While the direct application of kalouptsidi2012market's algorithm to static BLP models do not work well when the number of consumer types is large. Nevertheless, I show that it works well by simplifying and slightly modifying the original algorithm. The new algorithm corresponds to the mapping of value function $V$ discussed in Section (ref). We discuss these issues in detail in the Supplemental Appendix.

Although we focus on improving the inner-loop algorithms for estimating static/dynamic BLP models through Nested fixed-point (NFXP) approach, several alternative estimation procedures have been proposed so far.\footnote{Besides the studies mentioned above, bonnet2022yogurts proposed an inner-loop estimation algorithm of static demand models, using the idea of two-sided matching. The algorithm is applicable to static demand models other than static BLP models, including the pure characteristics model considered in berry2007pure. salanie2022fast proposed an estimator of static BLP models that does not require solving fixed-point problems and is approximately correct. The estimator can be used as the initial values of parameters in the standard nested fixed point estimation. } dube2012improving proposed the MPEC (Mathematical Programming with Equilibrium Constraint) approach for static and dynamic BLP models. lee2015computationally proposed the ABLP (approximate BLP) method for static BLP models, which iterates the process inspired by Newton's method to estimate parameters. Regarding dynamic BLP, Sun2019 proposed a Bayesian-type algorithm.

Although promising, recent studies on the MPEC have indicated that its performance is not necessarily as good as that of the NFXP approach. pal2023comparing showed in static BLP models that the computational times using the MPEC are longer than those with the NFXP approach using the traditional BLP contraction mapping and the spectral/SQUAREM algorithm as the acceleration method. In addition, the frequency of reaching the global optimum of the optimization problem is lower for MPEC.\footnote{In general, there is no guarantee that the GMM objective function in the BLP estimation is convex, and we may reach a local minimum of the GMM objective function.} Sun2019 also showed that the MPEC performed worse than the NFXP. As mentioned in dube2012improving, it is known that the MPEC is computationally costly when the Jacobian matrix of the constraint of the constrained optimization problem is dense, and the approach may not be always practical.

Furthermore, in general, there are two types of costs for any algorithm. The first one is “thinking costs”, which are associated with thinking about the validity of an algorithm for practitioners' needs and coding the algorithm. The second one is “computer costs”, which are associated with the computation time required to execute the algorithm on a computer.\footnote{jamshidian1997acceleration discuss this point in the context of the acceleration of the EM algorithm.} The proposed methods in the current article reduce the computation costs with small additional thinking costs, as they require only minor changes to the estimation procedure based on the standard NFXP approach widely applied by practitioners. Unlike the MPEC method, we do not have to install new software or packages to solve constrained optimization problems.\footnote{In MPEC, sometimes we must specify a sparsity structure to reduce the memory requirements. Our current method, however, does not require this procedure.} Unlike the ABLP method\footnote{pal2023comparing showed numerical results where the ABLP performs worse than the NFXP method using the traditional BLP contraction mapping and the spectral/SQUAREM algorithm as the acceleration method.}, the proposed methods do not require computing the analytical derivatives of functions\footnote{As discussed in miranda2004applied in the context of the Newton's method for solving nonlinear equations, analysts might make coding errors in coding derivatives of functions, and such a procedure should be avoided if possible. The algorithms proposed in the current article do not require coding derivatives, and are attractive.}, and the proposed algorithms are not restricted to static BLP models. Moreover, unlike the Baysian-type methods, there is no need to introduce Baysian-type techniques, such as Markov Chain Monte Carlo (MCMC) methods.\footnote{The method proposed by Sun2019 relies on a fixed-point mapping whose convergence is not necessarily fast in the NFXP approach, as shown in the current paper. The insights presented in the current paper might be useful for improving the performance of Sun2019's method.}

Finally, this study contributes to the literature discussing acceleration methods of fixed point iterations for BLP estimations. reynaerts2012enhencing, conlon2020best, and pal2023comparing found that combining the spectral or SQUAREM algorithms accelerates the inner-loop convergence of the static BLP estimation. The current study finds that the Anderson acceleration method, which is not discussed in these studies, outperforms the spectral and SQUAREM algorithms. Note that the choice of the step sizes in the spectral and SQUAREM algorithms largely affects the algorithm's performance, though not paid much attention to in the literature, and the current study briefly discusses it.

Static BLP model

This section focuses on the standard static BLP models without a nest structure (Random Coefficient Logit (RCL) model). We discuss dynamic BLP models without nest structure in Section (ref), and static BLP models with nest structure (Random Coefficient Nested Logit (RCNL) model) in Appendix (ref). Let $\left\Vert \cdot\right\Vert _{\infty}$ be the sup norm, and let $\left\Vert \cdot\right\Vert _{2}$ be the L2 norm. $\left\Vert \cdot\right\Vert $ also denotes a norm. $\epsilon_{\delta},\epsilon_{V},\epsilon$ denote values of the inner-loop tolerance levels.\footnote{conlon2020best recommended setting the values between 1E-12 and 1E-14 when using 64-bit computers and applying the standard BLP contraction mapping.}

Model

This section considers a static BLP model with random coefficients. This study proposes new algorithms whose convergence properties outperform those of the traditional BLP contraction mapping, and also generalizable to dynamic BLP models. Let consumer $i$'s utility when buying product $j$ be $v_{ij}=\delta_{j}+\mu_{ij}+\epsilon_{ij}$, and utility when buying nothing be $v_{i0}=\epsilon_{i0}$. $\delta_{j}$ denotes product $j$'s mean utility, and $\mu_{ij}$ denotes consumer $i$-specific utility of product $j$. $\epsilon$ denotes the idiosyncratic utility shocks. Let $\mathcal{J}$ be the set of products. Then, under the assumption that $\epsilon$ follows Gumbel distribution and that the observed market share data match the market share predicted by the model, the following equations hold:

Here, $dP(i)$ denotes the density of consumer $i$. We use appropriate discretizations of the consumer types to solve the model numerically, when consumer types are continuously distributed.\footnote{See conlon2020bestfor a review of the efficient discretization methods.} Let $\mathcal{I}$ be the set of discretized consumer types. $w_{i}$ denotes the fraction of type $i$ consumers, and $\sum_{i\in\mathcal{I}}w_{i}=1$ hold. If we use $I$ Monte Carlo draws to approximate the integral, $w_{i}=\frac{1}{I}$ holds. Note that equation ((ref)) can be derived from ((ref)).

In empirical applications, the distribution of $\mu_{ij}$ is parameterized by $\theta_{n}$. For instance, $\mu_{ij}=\theta_{n}X_{j}\nu_{i}$ where $\nu_{i}\sim N(0,\Sigma)$. The mean product utility $\delta$ is represented as $\delta_{j}=X_{j}\theta_{l}+\xi_{j}$, where $\xi_{j}$ denotes the unobserved product characteristics of product $j$. Econometricians seek the values of $(\theta_{l},\theta_{n})$ minimizing the GMM objective $G(\theta_{l},\theta_{n})WG(\theta_{l},\theta_{n})$ where $G(\theta_{l},\theta_{n})\equiv Z\xi(\theta_{l},\theta_{n})=Z\left(\delta(\theta_{n})-X\theta_{l}\right)$, $W$ denotes the weight matrix, and $Z$ denotes the appropriate instrumental variables. $\delta(\theta_{n})$ denotes the values of $\delta$ satisfying equations ((ref)) and ((ref)). Algorithm (ref) shows the formal steps.

Below, we focus on Step 1(b) of Algorithm (ref), namely, how to efficiently solve for $\delta$ satisfying the constraints ((ref)) and ((ref)) given $\mu$ and observed market shares $S^{(data)}$.\footnote{conlon2020best discuss the details of overall estimation procedures of static BLP models. Note that when no consumer heterogeneity exists, we can recover the value of $\delta$ by $\delta_{j}=\log\left(S_{j}^{(data)}\right)-\log\left(S_{0}^{(data)}\right)$, and we can estimate the linear parameters $\theta_{l}$ by a linear GMM without solving the fixed-point problem, as discussed in berry1994estimating.}

Mappings of mean product utility $\delta$

To solve for $\delta$ given $\mu$ and observed market shares $S^{(data)}\equiv\left(S_{j}^{(data)}\right)_{j\in\mathcal{J}}$, berry1995automobile proposed iterating the update $\delta_{j}^{(n+1)}=\delta_{j}^{(n)}+\log\left(S_{j}^{(data)}\right)-\log\left(s_{j}\left(\delta^{(n)}\right)\right)$ with appropriate initial values $\delta^{(0)}$ until convergence. Here, $s_{j}(\delta)\equiv\sum_{i}w_{i}\frac{\exp\left(\delta_{j}+\mu_{ij}\right)}{1+\sum_{k\in\mathcal{J}}\exp\left(\delta_{k}+\mu_{ik}\right)}$ denotes the market share of product $j$ predicted by the structural model. The mapping is known as the BLP contraction mapping.

Motivated by the BLP contraction mapping, we define a mapping $\Phi^{\delta,\gamma}:B_{\delta}\rightarrow B_{\delta}$ for $\gamma\in\mathbb{R}$ such that\footnote{Although the mapping $\Phi_{j}^{\delta,\gamma}$ implicitly depends on the values of $S^{(data)}$ and $\mu$, we omit expressing them as arguments of the mapping to simplify the exposition.}:

where $s_{0}(\delta)=1-\sum_{k\in\mathcal{J}}s_{k}(\delta)$. $B_{\delta}$ denotes the space of $\delta\equiv\left\{ \delta_{j}\right\} _{j\in\mathcal{J}}$. Obviously, $\Phi^{\delta,\gamma=0}$ matches the traditional BLP contraction mapping.

Proposition (ref) ensures that we can find the solution $\delta$ satisfying $S_{j}^{(data)}=s_{j}(\delta)\ (j\in\mathcal{J})$ by alternatively solving the fixed-point constraint $\delta=\Phi^{\delta,\gamma\geq0}\left(\delta\right)$.

Clearly, the solution of $S_{j}^{(data)}=s_{j}(\delta)\ (j\in\mathcal{J})$ satisfies $\delta_{j}=\Phi_{j}^{\delta,\gamma}\left(\delta\right)=\delta_{j}+\left[\log\left(S_{j}^{(data)}\right)-\log\left(s_{j}(\delta)\right)\right]-\gamma\left[\log\left(S_{0}^{(data)}\right)-\log\left(s_{0}(\delta)\right)\right]\ j\in\mathcal{J}$. However, its converse is not trivial, even without formal proof. Proposition (ref) formally shows that the converse actually holds under $\gamma\geq0$.

Algorithm (ref) shows the algorithm to solve for $\delta$ using the mapping $\Phi^{\delta,\gamma}$.

Although, in principle, we can choose any real value of $\gamma$, we mainly consider $\gamma=1$ because of its good convergence properties. As discussed in detail in Appendix (ref), the consumer heterogeneity size largely affects the convergence speed of $\Phi^{\delta,\gamma=1}$. As consumer heterogeneity decreases, the convergence speed increases. In the absence of consumer heterogeneity, the iteration immediately converges after one iteration, regardless of the choice of initial values $\delta^{(0)}$. When no consumer heterogeneity exists, $\delta_{j}=\log(S_{j}^{(data)})-\log(S_{0}^{(data)})$ holds, as shown in berry1994estimating. Regarding the mapping $\Phi^{\delta,\gamma=1}$, $\Phi_{j}^{\delta,\gamma}\left(\delta\right)=\log\left(S_{j}^{(data)}\right)-\log\left(S_{0}^{(data)}\right)$ holds, and the output of $\Phi^{\delta,\gamma}$ is equal to the true $\delta$ for any input $\delta$.

In general, there is no guarantee that $\Phi^{\delta,\gamma=1}$ is a contraction. It implies there is no guarantee that the iteration with $\delta_{j}^{(n+1)}=\delta_{j}^{(n)}+\left[\log\left(S_{j}^{(data)}\right)-\log\left(s_{j}(\delta^{(n)})\right)\right]-\left[\log\left(S_{0}^{(data)}\right)-\log\left(s_{0}(\delta^{(n)})\right)\right]$ converges. In fact, we can construct a simple numerical example where $\Phi_{j}^{\delta,\gamma=1}$ is not a contraction in an extreme setting where consumer heterogeneity is too large, as shown in the Supplemental Appendix. However, the iteration converges when combining the spectral/SQUAREM algorithm, which we discuss in detail in Section (ref). In addition, as shown in Section (ref), the iteration always converges and sometimes leads to drastic improvement of convergence speed in the standard settings of Monte Carlo simulations experimented in the previous studies (dube2012improving; lee2015computationally) and the estimations using datasets used by nevo2001measuring and Berry et al. berry1995automobile,berry1999voluntary. As mapping $\Phi^{\delta,\gamma=1}$ is straightforward to implement once coding the standard BLP contraction mapping, applying the mapping $\Phi^{\delta,\gamma=1}$ is worth considering in the static BLP estimation.

Please note that the global convergence of iterations using $\Phi^{\delta,\gamma=1}$can be ensured by adding a few lines in the programming code, by using the fact that $\Phi^{\delta,\gamma=0}$, traditional BLP contraction mapping, is a contraction. For details, see Appendix (ref). Conservative practitioners can consider this process for convergence.

Concerning the convergence speed of the BLP contraction mapping ($\Phi^{\delta,\gamma=0}$), it is slow when the outside option share is small (dube2012improving). In contrast, regarding the new mapping $\Phi^{\delta,\gamma=1}$, the convergence speed is less sensitive to the outside option share, as formally discussed in Appendix (ref) and numerically shown in Section (ref).

Mappings of value functions $V$

Although it may seem intuitive to consider mappings of $\delta$ to solve for $\delta$, we can alternatively solve for consumers' value functions $V$ using a mapping of $V$, and then recover $\delta$ using an analytical formula. Here, consumer $i$'s value function, or the inclusive value, is defined as follows:

Below, we discuss the approach solving for $V$. It seems that applying the approach is not necessary, when only considering the static BLP models. The approach directly solving for $\delta$ is more intuitive and standard. However, when examining dynamic BLP models, introducing the approach solving for $V$ is attractive, because it closely relates with the value function iterations typically used in solving dynamic models, as discussed in Section (ref).

kalouptsidi2012market first proposed solving for consumer-type specific variables in the inner loop of static BLP estimations. However, the original algorithms have convergence problems. The following algorithms overcome the issues by choosing appropriate mappings that have dualistic relations with the mappings of $\delta$. For further discussion on the kalouptsidi2012market's algorithms, see the Supplemental Appendix.

First, because $S_{j}^{(data)}=\exp\left(\delta_{j}\right)\cdot\sum_{i}w_{i}\frac{\exp\left(\mu_{ij}\right)}{1+\sum_{k\in\mathcal{J}}\exp\left(\delta_{k}+\mu_{ik}\right)}=\exp\left(\delta_{j}\right)\cdot\sum_{i}w_{i}\frac{\exp\left(\mu_{ij}\right)}{\exp\left(V_{i}\right)}$ holds by ((ref)) and ((ref)), the following equation should hold:

Based on this formula, we define a function $\Phi^{V,\gamma}:B_{V}\rightarrow B_{V}$ for $\gamma\in\mathbb{R}$ such that:

where $\iota_{\delta\rightarrow V}:B_{\delta}\rightarrow B_{V}$ is a mapping such that:

and $\iota_{V\rightarrow\delta}^{\gamma}:B_{V}\rightarrow B_{\delta}$ is a mapping such that:

$B_{V}$ denotes the space of $V\equiv\left\{ V_{i}\right\} _{i\in\mathcal{I}}$. Note that $\Phi_{i}^{V,\gamma}\left(V\right)=\log\left(1+\left(\sum_{j\in\mathcal{J}}S_{j}^{(data)}\frac{\exp(\mu_{ij})}{\sum_{i}w_{i}\exp(\mu_{ij})\exp(-V_{i})}\right)\cdot\left(\frac{\sum_{i}w_{i}\exp(-V_{i})}{S_{0}^{(data)}}\right)^{\gamma}\right)$ holds.

The following proposition shows the validity of using the mapping $\Phi^{V,\gamma\geq0}$:

Hence, to solve for $\delta$, two algorithms can be considered. The first is to apply $\Phi^{\delta,\gamma}$ iteratively and solve for $\delta$ with appropriate initial values of $\delta$, as discussed in Section (ref). The second one is to apply $\Phi^{V,\gamma}$ iteratively with appropriate initial values of $V$, and obtain $\delta$, which is analytically computed using $V$ in the iteration. Algorithm (ref) outlines this second approach.

We can easily prove a dualistic relationship between $\Phi^{\delta,\gamma}$ and $\Phi^{V,\gamma}$. As shown in the Monte Carlo simulation (Section (ref)), their convergence speed is generally the same, provided that they take a common value $\gamma$. The following proposition is a formal statement:

Figure (ref) graphically shows the relationships between the mappings of $\delta$ and $V$. Here, $\delta$-($\gamma$) denotes the algorithm using the mapping $\Phi^{\delta,\gamma}$, and $V$-($\gamma$) denotes the algorithm using the mapping $\Phi^{V,\gamma}$.

In Appendix (ref), we discuss the convergence properties of the mapping $\Phi^{V,\gamma}$. As in the case of $\Phi^{\delta,\gamma}$, the size of the consumer heterogeneity affects the convergence speed.

Dynamic BLP model

Model

Next, we consider a dynamic BLP model, where consumers' forward-looking decisions are integrated into the static BLP models. Similar to the static BLP models, we assume that consumer types are appropriately discretized to solve the model numerically. Let $\mathcal{I}$ be the set of consumer types, and let $w_{i}$ be the fraction of the type $i$ consumers.

The discounted sum of utility of consumer $i$ at state $(x_{it},\Omega_{t})$ when buying product $j$ at time $t$ is:

The discounted sum of utility of consumer $i$ at state $(x_{it},\Omega_{t})$ when not buying any product at time $t$ is:

Here, $\beta$ denotes consumers' discount factor. $x_{it}$ denotes individual state variables, such as durable product holdings (for durable goods) or the product purchased in the previous period (for goods with switching costs). $\Omega_{t}$ denotes market-level state variables, including characteristics of products sold in the market. $\Omega_{t}$ may include the mean product utilities $\delta$. $\epsilon$ denotes i.i.d. idiosyncratic utility shocks, and $V$ denotes the (integrated) value function, defined by $V_{it}(x_{it},\Omega_{t})\equiv E_{\epsilon}\left[\max_{j\in\mathcal{A}_{t}(x_{it})}v_{ijt}(x_{it},\Omega_{t})\right]$, where $E_{\epsilon}$ denotes the expectation operator concerning $\epsilon$. $E_{t}$ denotes the expectation operator concerning the realizations of the future states $(x_{it+1},\Omega_{t+1})$.

Assuming that $\epsilon$ follows i.i.d. mean zero type-I extreme value distribution and that the observed market shares match the predicted market shares, $\delta$ and $V$ satisfy the following equations:

where

Here, $S_{jt}^{(data)}$ denotes product $j$'s market share at time $t$, and $S_{0t}^{(data)}$ denotes the fraction of consumers not buying anything at time $t$. $Pr_{it}(x_{it})$ denotes the fraction of type $i$ consumers at state $x_{it}$. For example, for durables, $Pr_{it}(x_{it}=\emptyset)$ denotes the fraction of type $i$ consumers not holding any durable product at time $t$, if we let $x_{it}=\emptyset$ be the state where consumers do not hold any durable products. $\mathcal{A}_{t}(x_{it})\subset\mathcal{J}_{t}\cup\{0\}$ denotes the consideration set of consumers at state $x_{it}$ at time $t$. For instance, in the durable goods example above, $\mathcal{A}_{t}(x_{it}\neq\emptyset)=\{0\}$ holds if consumers holding any durable product do not buy additional products. Here, $\mathcal{J}_{t}$ denotes the set of products sold at time $t$. The term $\frac{\exp\left(\delta_{jt}+\mu_{ijt}(x_{it},\Omega_{t})+\beta E_{t}\left[V_{it+1}(x_{it+1},\Omega_{t+1})|x_{it},\Omega_{t},a_{it}=j\right]\right)}{\exp\left(V_{it}(x_{it},\Omega_{t})\right)}\equiv s_{ijt}^{(ccp)}(x_{it},\Omega_{t})$ denotes the conditional choice probability (CCP) of type $i$ consumers at state $(x_{it},\Omega_{t})$ choosing product $j$ at time $t$. $\frac{\exp\left(\mu_{i0t}(x_{it},\Omega_{t})+\beta E_{t}\left[V_{it+1}(x_{it+1},\Omega_{t+1})|x_{it},\Omega_{t},a_{it}=0\right]\right)}{\exp\left(V_{it}(x_{it},\Omega_{t})\right)}\equiv s_{i0t}^{(ccp)}(x_{it},\Omega_{t})$ denotes the CCP of type $i$ consumers at state $(x_{it},\Omega_{t})$ choosing the outside option at time $t$.

Algorithm

The overall estimation algorithm is similar to the static BLP: As for dynamic BLP, we should solve for $\delta$ satisfying equations ((ref))-((ref)) given $\mu$.\footnote{Unlike the static BLP models, there is no guarantee that $\delta$ satisfying ((ref))-((ref)) is unique. If multiple $\delta$ satisfying these equations exist, we should adopt $\delta$ with the lowest GMM objective.} As in the static BLP, we focus on efficient algorithms for solving for $\delta$ satisfying ((ref))-((ref)) given $\mu$.

Here, we examine a nonstationary environment, where consumers have perfect foresight on the state transition of $\Omega_{t}$, and $\Omega_{t}$ remains constant after the terminal period $T$.\footnote{Similar specification was considered in gowrisankaran2012dynamics and conlon2012dynamic. Analogous specification was used in igami2017estimating in the context of dynamic discrete games. Besides, regarding consumer expectations except for the terminal period, we can relax the assumption of perfect foresight to rational expectations where consumer expectations are on average correct. For details, see the discussion in fukasawa2024lightbulb. The assumption of rational expectation has been applied in kalouptsidi2020linear and empirical studies cited therein. } Although this specification differs from the one under inclusive value sufficiency, widely applied in earlier studies (e.g., gowrisankaran2012dynamics), the algorithm under the current specification is relatively simple and is the basis of algorithms under alternative specifications. Besides, we consider a setting where the values of $\left\{ Pr_{it}(x_{it})\right\} _{x_{it}\in\chi,t=1,\cdots,T}$ are known, to clarify the algorithm's essence. In Section (ref), we address cases where the values of $\left\{ Pr_{it}(x_{it})\right\} _{x_{it}\in\chi,t=1,\cdots,T}$ are unknown.

Equations ((ref)), ((ref)), and ((ref)) imply $V$ is the solution of $V=\Phi^{V,\gamma}(V)\ \ \forall\gamma\in\mathbb{R}$, where

Here, we define a function $\Psi_{V\delta\rightarrow V}^{\gamma}:B_{V}\times B_{\delta}\rightarrow B_{V}$ such that:

We also define a function $\iota_{V\rightarrow\delta}:B_{V}\rightarrow B_{\delta}$ such that:

{

}

$B_{V}$ is a space of $V\equiv\left\{ V_{it}(x_{it},\Omega_{t})\right\} _{x_{it}\in\mathcal{\chi},i\in\mathcal{I},t=1,\cdots,T}$, and $B_{\delta}$ is a space of $\delta\equiv\left\{ \delta_{jt}\right\} _{j\in\mathcal{J}_{t},t=1,\cdots,T}$. Note that $V_{iT+1}(x_{iT+1},\Omega_{T+1}^{(data)})=V_{iT}(x_{iT},\Omega_{T}^{(data)})$ holds by the assumption that consumers have perfect foresight on the state transition of $\Omega_{t}$, and $\Omega_{t}$ remains constant after the terminal period $T$.

Therefore, when the values of $Pr_{it}(x_{it})$ are known, we can solve for $\delta$ by iteratively applying $\Phi^{V,\gamma}$ and obtain $\delta$ in the process. Algorithm (ref) shows the steps:

Algorithms used in the previous studies

The upper part of Algorithm (ref) was applied in earlier studies. The algorithm is simple: Researchers first solve for the value function $V$ given the values of mean product utilities $\delta$, then calculate the difference between the observed and predicted market shares based on the structural model. If the norm of the difference is larger than a tolerance level, the BLP contraction mapping $\delta_{jt}^{(n+1)}=\delta_{jt}^{(n)}+\phi\left[\log(S_{jt}^{(data)})-\log(s_{jt}(\delta))\right]$ is applied\footnote{In dynamic BLP models, there is no guarantee that the iteration $\delta_{jt}^{(n+1)}=\delta_{jt}^{(n)}+\left[\log(S_{jt}^{(data)})-\log(s_{jt}(\delta^{(n)}))\right]$ converges. Hence, introducing a dampening parameter $\phi\in(0,1]$ is sometimes necessary to stabilize the convergence. Proposition 2 of Sun2019 shows that the mapping $\delta_{jt}^{(n+1)}=\delta_{jt}^{(n)}+\phi\left[\log(S_{jt}^{(data)})-\log(s_{jt}(\delta^{(n)}))\right]$ is a contraction if $\phi$ is “sufficiently small” and additional conditions hold.} and repeat the process again. However, the nested loops increase the computational burden.

We can think of an alternative algorithm shown in the lower part of Algorithm (ref), which avoids the nested loops by jointly updating the values of $\delta$ and $V$ based on the fixed-point constraint.\footnote{The idea of jointly updating two types of variables can also be found in Pakes_McGuire1994 algorithm to solve dynamic games. In the algorithm, essentially, all the firms' value functions and strategic variables are jointly updated.} Note that we should solve for two different types of variables $\delta$ and $V$ in the algorithm. In contrast, in the proposed algorithm (Algorithm (ref)), $\delta$ is analytically represented as a function of $V$, and we must only solve for $V$, which speeds up convergence. The superior performance of the proposed algorithm is shown in Section (ref).

Acceleration methods of fixed point iterations

This section briefly describes three acceleration methods of fixed point iterations: Anderson acceleration, spectral algorithm, and SQUAREM algorithm.

Anderson acceleration

Anderson acceleration, originally proposed by anderson1965iterative, is a method to accelerate the convergence of fixed-point iterations.

Suppose that we want to solve $x=\Phi(x)$ for $x$. The Anderson acceleration method solves the problem by the steps shown in Algorithm (ref).

Intuitively, $x^{(n+1)}=\sum_{l=0}^{m_{n}}\theta_{l}^{(n)}x^{(n-l)}$ is chosen to assign a larger weight on $x^{(n-l)}$ with smaller absolute values of $f_{n-l}\equiv\Phi\left(x^{(n-l)}\right)-x^{(n-l)}$. As noted in fang2009two, this method closely relates to the quasi-Newton method. $m$ or $m_{n}$ is sometimes called memory size, because it represents how long we retain the information of past iterations $\left(f_{n-m_{n}},\cdots,f_{n}\right)$. In the numerical experiments shown in Section (ref), $m=5$ works well.

In Step 3(c) of the algorithm, we must solve a linear constrained optimization problem. The problem can be easily solved by alternatively solving the linear least squares problem\footnote{Strictly speaking, we cannot rule out the possibility that the values of $\left(f_{n-m_{n}},\cdots,f_{n}\right)$ are close to collinear, and the matrix $\left(\mathcal{F}_{n}^{T}\mathcal{F}_{n}\right)$ is close to singular. In that case, we choose one of $\theta$ minimizing $\left\Vert \sum_{l=0}^{m_{n}}\theta_{l}^{(n)}f_{n-l}\right\Vert _{2}^{2}$. In the numerical experiments in Section (ref), using the mldivide function suppressing any warnings on singular matrices in MATLAB and the lstsq function in numpy package in Python works well.} outlined in the note of Algorithm (ref).

Spectral algorithm

The spectral algorithm is designed to solve nonlinear equations and nonlinear continuous optimization problems. To solve a nonlinear equation $F(x)=0$, $x^{(n)}$ is iteratively updated as follows, given the initial values of $x^{(0)}$:\footnote{Newton's method, which uses the updating equation $x^{(n+1)}=x^{(n)}-\left(\nabla F(x^{(n)})\right)^{-1}F\left(x^{(n)}\right)\ (n=0,1,2,\cdots)$, attains fast convergence around the solution. Although Newton's method can be applied to solve the equation, computing $\nabla F(x^{(n)})$ requires coding analytical first derivatives of $F$, which is not an easy task, especially when the function $F$ is complicated. Moreover, especially when $n_{x}$, the dimension of $x$, is large, computing the inverse of $n_{x}\times n_{x}$ matrix $\nabla F(x^{(n)})$ is computationally costly. Hence, the use of the spectral algorithm is attractive from the perspective of simplicity and computational cost.}

\[ x^{(n+1)}=x^{(n)}+\alpha^{(n)}F\left(x^{(n)}\right)\ (n=0,1,2,\cdots) \]

Here, $\alpha^{(n)}$ is typically based on the values of $s^{(n)}\equiv x^{(n)}-x^{(n-1)}$ and $y^{(n)}\equiv F(x^{(n)})-F(x^{(n-1)})$, which can vary across $n=0,1,2,\cdots$.

The idea of the spectral algorithm can be utilized in the context of fixed point iterations. Suppose we want to solve a fixed point constraint $x=\Phi(x)$. Then, by letting $F(x)=\Phi(x)-x$, we can solve $x=\Phi(x)$ by iteratively updating the values of $x$ as follows:\footnote{Spectral algorithm can be thought of as a generalization of extrapolation method (cf. judd1998numerical) with varying values of $\alpha$.}

As discussed in detail in varadhan2008simple,\footnote{Step size $\alpha_{S3}\equiv\frac{\left\Vert s^{(n)}\right\Vert _{2}}{\left\Vert y^{(n)}\right\Vert _{2}}$ can be derived from a simple optimization problem $\min_{\alpha^{(n)}}\frac{\left\Vert x^{(n+1)}-x^{(n)}\right\Vert _{2}^{2}}{\left|\alpha^{(n)}\right|}=\frac{\left\Vert s^{(n)}+\alpha^{(n)}y^{(n)}\right\Vert _{2}^{2}}{\left|\alpha^{(n)}\right|}$, as discussed in varadhan2008simple.} the step size $\alpha^{(n)}=\alpha_{S3}\equiv\frac{\left\Vert s^{(n)}\right\Vert _{2}}{\left\Vert y^{(n)}\right\Vert _{2}}>0$ works well. The iteration choosing negative values of $\alpha^{(n)}$ leads to farther point from the solution when the mapping $\Phi$ is a contraction, as discussed in the Supplemental Appendix. Although conlon2020best and pal2023comparing mentioned $\alpha_{S1}\equiv-\frac{s^{(n)\prime}y^{(n)}}{y^{(n)\prime}y^{(n)}}$, they may take negative values, which can destabilize the convergence. In contrast, $\alpha_{S3}\equiv\frac{\left\Vert s^{(n)}\right\Vert _{2}}{\left\Vert y^{(n)}\right\Vert _{2}}$ always takes positive values, and it is preferable to apply $\alpha_{S3}$ as the step size. Algorithm (ref) shows the detailed steps in the spectral algorithm under the step size $\alpha_{S3}$.

Variable-type-specific step size

Although the standard spectral algorithm uses a scalar $\alpha^{(n)}$, I newly introduce the idea of variable-type-specific step sizes to the spectral algorithm.

To introduce the idea, suppose we would like to solve a nonlinear equation $F(x)=\left(

\right)=0\in\mathbb{R}^{n_{1}+n_{2}}$, where $x=\left(

\right)\in\mathbb{R}^{n_{1}+n_{2}}$. When using the standard step size $\alpha_{S3}$, we set $\alpha^{(n)}=\frac{\left\Vert s^{(n)}\right\Vert _{2}}{\left\Vert y^{(n)}\right\Vert _{2}}$, where $s^{(n)}\equiv\left(

\right)-\left(

\right)$, $y^{(n)}\equiv\left(

\right)-\left(

\right)$, and update $x$ by $x^{(n+1)}=x^{(n)}+\alpha^{(n)}F\left(x^{(n)}\right)$.

Here, we can alternatively update the values of $x=\left(

\right)$ by $x^{(n+1)}=\left(

\right)=\left(

\right)+\left(

\right),$ where $\alpha_{m}^{(n)}=\frac{\left\Vert s_{m}^{(n)}\right\Vert _{2}}{\left\Vert y_{m}^{(n)}\right\Vert _{2}}$, $s_{m}^{(n)}\equiv x_{m}^{(n)}-x_{m}^{(n-1)}$, $y_{m}^{(n)}\equiv F_{m}(x^{(n)})-F_{m}(x^{(n-1)})\ (m=1,2)$. The idea can be generalized to the case where $N\in\mathbb{N}$ types of variables exist.

As discussed earlier, $\alpha^{(n)}$ is chosen to accelerate the convergence. In principle, any choices of $\alpha^{(n)}$ are allowed in the spectral algorithm provided it is effective, and its numerical performance is more critical.\footnote{We can justify the strategy with a simple thought experiment. Suppose $x_{1}^{(n)}$ is fixed at the true value $x_{1}^{*}$, and $x_{2}^{(n)}$ is not. Then, we should solve the equation $F_{2}(x_{2};x_{1}^{*})=0$ as an equation of $x_{2}$, and setting $\alpha_{2}^{(n)}=\frac{\left\Vert s_{2}^{(n)}\right\Vert _{2}}{\left\Vert y_{2}^{(n)}\right\Vert _{2}}$ is desirable. If we assume $\alpha_{1}^{(n)}=\alpha_{2}^{(n)}=\alpha^{(n)}$, we should set $\alpha^{(n)}=\frac{\left\Vert s^{(n)}\right\Vert _{2}}{\left\Vert y^{(n)}\right\Vert _{2}}=\frac{\sqrt{\left\Vert s_{1}^{(n)}\right\Vert _{2}^{2}+\left\Vert s_{2}^{(n)}\right\Vert _{2}^{2}}}{\sqrt{\left\Vert y_{1}^{(n)}\right\Vert _{2}^{2}+\left\Vert y_{2}^{(n)}\right\Vert _{2}^{2}}}$, which may not be equal to $\frac{\left\Vert s_{2}^{(n)}\right\Vert _{2}}{\left\Vert y_{2}^{(n)}\right\Vert _{2}}$.}

SQUAREM algorithm

varadhan2008simple proposed the SQUAREM (Squared Polynomial Extrapolation Methods) algorithm, which solves the fixed-point problem $x=\Phi(x)$, and the updating equation is as follows:

where $s^{(n)}\equiv\Phi\left(x^{(n)}\right)-x^{(n)}$, $y^{(n)}\equiv F\left(\Phi\left(x^{(n)}\right)\right)-F\left(x^{(n)}\right)=\Phi^{2}\left(x^{(n)}\right)-2\Phi\left(x^{(n)}\right)+x^{(n)}$, and $F(x)=\Phi(x)-x$. Step size $\alpha^{(n)}$ is chosen based on the values of $s^{(n)}$ and $y^{(n)}$.

The SQUAREM algorithm is closely related to the spectral algorithm. Updating equation ((ref)) can be reformulated as:

Here, we define a spectral update function, defined by $\Psi_{spectral}\left(x;\Phi,\alpha\right)\equiv(1-\alpha)x+\alpha\Phi(x)$. The right-hand side can be regarded as the output of the spectral update $\Psi_{spectral}\left(x^{(n)};\Phi,\alpha^{(n)}\right)$ and $\Psi_{spectral}\left(\Phi\left(x^{(n)}\right);\Phi,\alpha^{(n)}\right)$.

Algorithm (ref) shows the steps when we the SQUAREM algorithm with step size $\alpha_{S3}\equiv\frac{\left\Vert s^{(n)}\right\Vert _{2}}{\left\Vert y^{(n)}\right\Vert _{2}}$ is applied:

Choice of acceleration methods

As numerically shown in Section (ref), the Anderson acceleration outperforms both the spectral and SQUAREM algorithms in static and dynamic BLP applications, and the current study recommends using the Anderson acceleration as the acceleration method. Because even the spectral and SQUAREM algorithms work fairly well, they can be good alternatives in case the Anderson acceleration does not work well.

Regarding convergence, Anderson acceleration is mainly designed to accelerate the convergence of a contraction mapping, and the iteration may not converge when the mapping does not have contraction properties.\footnote{zhang2020globally developed a globalization strategy for nonexpansive mappings, which is a slight extension of contraction mappings, in Anderson acceleration. However, the strategy may not work when there is no guarantee that the iteration is nonexpansive.} In contrast, regarding the spectral algorithm, there are many studies on globalization strategies to stabilize convergence (e.g., la2006spectral; huang2017new), and the spectral algorithm can be a good alternative for practitioners prioritizing convergence. Note that, in static BLP models, the global convergence can be ensured regardless of acceleration methods used, simplify by adding a few lines in the programming code, utilizing the fact that the BLP contraction mapping is a contraction. For details, see Appendix (ref).

Numerical Experiments

This section shows the results of numerical experiments for static and dynamic BLP models. All the experiments were run on a laptop computer with the CPU AMD Ryzen 5 6600H 3.30 GHz, 16.0 GB of RAM, Windows 11 64-bit, and MATLAB 2022b. In this study, we assume an algorithm solving a fixed-point constraint $x=\Phi(x)$ does not converge, when the value of $\left\Vert \Phi(x)-x\right\Vert _{\infty}$ takes an infinite or NaN value, or the number of function evaluations reaches a pre-determined maximum.

Static BLP model

Here, we define $\delta$-(0) and $\delta$-(1) as the algorithms using mapping $\Phi^{\delta,\gamma=0}$ and $\Phi^{\delta,\gamma=1}$, respectively. Similarly, we denote $V$-(0) and $V$-(1) as the algorithms using mapping $\Phi^{V,\gamma=0}$ and $\Phi^{V,\gamma=1}$, respectively.

Monte Carlo simulation

To compare the performance of the algorithms, I conducted Monte Carlo simulation. The settings are the same as those of lee2015computationally and dube2012improving, and similar to those of pal2023comparing. Consumer $i$'s utility for choosing product $j$ is $U_{ij}=X_{j}\theta_{i}+\xi_{j}+\epsilon_{ij}\ (j=1,\cdots,J)$, and utility for not buying is $U_{i0}=\epsilon_{i0}$. Let $\delta_{j}\equiv X_{j}\cdot E[\theta_{i}]+\xi_{j}$. For details on the data-generating process, see the Supplemental Appendix. Number of simulation draws is set to 1000. The tolerance level of the convergence is set to 1E-13. The initial value $\delta_{j}^{(0)}$ is set to $\log(S_{j}^{(data)})-\log(S_{0}^{(data)})$, and the initial value of $V_{i}^{(0)}$ is set to 0.

We evaluate the performance of the BLP inner-loop algorithms in solving for $\delta$ given parameter values $\theta_{n}$ (standard deviation of random coefficients $\theta_{i}$). First, market share data are generated based on the true nonlinear parameter values $\theta_{n}^{*}$ and true $\delta^{*}$. As analysts lack knowledge of the true nonlinear parameter values during estimation process, we imitate the setting by drawing parameter values $\theta_{n}\sim U[0,2\theta_{n}^{*}]$, and measuring the performance of each algorithm under the parameters. This process is repeated 50 times.

Here, we mainly focus on the number of function evaluations, rather than CPU time. Although the CPU time indicates the exact duration required for iteration termination, generally it largely depends on the algorithm's syntax and the characteristics of each programming language. In contrast, the number of function evaluations remains unaffected by these factors as we focus on the values.\footnote{In the case of the standard fixed-point iterations, Anderson acceleration, and the spectral algorithm, the number of function evaluations equals the number of iterations. Regarding the SQUAREM algorithm, the number of function evaluations is roughly half of the number of iterations, and we focus on the number of function evaluations as the indicator. }

Tables (ref) and (ref) show the results for $J=25$ and $J=250$.\footnote{The current study also conducted numerical experiments in settings with $J=25$ and larger constant term of the utility than the setting in Table (ref). The results are also consistent with the discussions in this section. They are available upon request.} The mappings $\delta$-(1) and $V$-(1) require significantly fewer iterations than $\delta$-(0) and $V$-(0), with their superior performance more evident for $J=250$. As discussed in Section (ref) and detailed in Appendix (ref), the convergence speed of the BLP contraction mapping ($\delta$-(0)) can be slow when the outside option share is small, whereas $\delta$-(1) remains unaffected. As shown in the tables' notes, the mean outside option share is 0.850 for $J=25$, and 0.307 for $J=250$. The observations are in line with the discussion. Moreover, the numerical results highlight the duality between the mappings of $\delta$ and $V$: the numbers of iterations are mostly the same for $\delta$-($\gamma$) and $V$-($\gamma$) ($\gamma=0,1$).

Combining the acceleration methods further reduces the number of function evaluations. Anderson acceleration is the fastest, and outperforms both the spectral and SQUAREM algorithms, which show similar performance. Note that the spectral and SQUAREM algorithms perform similarly. Here, we use the step size $\alpha_{S3}$. When we use the other step sizes ($\alpha_{S1}$,$\alpha_{S2}$), the iterations sometimes do not converge. Detailed results are provided in the Supplemental Appendix.

Note that the value of $\left\Vert \log(S_{j}^{(data)})-\log(s_{j})\right\Vert _{\infty}$ might not be close to 0, even when the iteration converges. For instance, the convergence of the iteration based on the mapping $\delta$-(1) is determined by $\left\Vert \Phi^{\delta,\gamma=1}(\delta)-\delta\right\Vert =\left\Vert \left(\log(S_{j}^{(data)})-\log(s_{j})\right)-\left(\log(S_{0}^{(data)})-\log(s_{0})\right)\right\Vert <\epsilon$, not $\left\Vert \log(S_{j}^{(data)})-\log(s_{j})\right\Vert <\epsilon$. Hence, to ensure that the iteration yields a solution with sufficiently small $\left\Vert \log(S_{j}^{(data)})-\log(s_{j})\right\Vert $, the tables display the mean values of $DIST\equiv\left\Vert \log(S_{j}^{(data)})-\log(s_{j})\right\Vert _{\infty}$ computed using the iteration results. They also show the percentage of settings where the value of $DIST$ is smaller than 1E-12. The results indicate that all algorithms except for the slow $\delta$-(0), $V$-(0) algorithms satisfy the precision constraint.

Replications of nevo2001measuring and Berry et al. berry1995automobile,berry1999voluntary

To further examine the performance of inner-loop algorithms, we use the datasets from nevo2001measuring and Berry et al. berry1995automobile,berry1999voluntary to estimate demand parameters under different inner-loop algorithms.

In the numerical experiments, the PyBLP package (version 0.13.0) in Python (conlon2020best) was used. Because the package itself does not support the new mapping $\delta$-(1), the original PyBLP source code was modified to introduce the new mapping $\delta$-(1). In addition, the steps of the Anderson acceleration were newly coded. Regarding the spectral algorithm, we cannot choose step sizes other than $\alpha_{S2}$, because the dfsane function in the SciPy package, which the PyBLP relies on, exclusively supports $\alpha_{S2}$. Hence, the dfsane function in the SciPy package was modified to allow different step sizes. Moreover, because the current dfsane function does not support such iterations without globalization strategies, the code was adjusted to implement such iterations. Notably, it rarely affects the results, as shown in the Supplemental Appendix.

\paragraph*{nevo2001measuring's dataset}

Settings largely reflect those specified in Figure 5 of conlon2020best, including initial parameter values and tolerance levels, except for the choice of step sizes, the existence of globalization strategies, and the choice of mappings.

Table (ref) shows that the algorithm $\delta$-(1) is more than twice as fast as the algorithm $\delta$-(0). Incorporating Anderson acceleration further reduces total function evaluations by approximately 75%. In particular, the values of total objective evaluations and the objective remain consistent across all inner-loop algorithms, indicating that the choice of the inner-loop algorithm does not affect the outer-loop estimation in this setting.

\paragraph*{Berry et al.berry1995automobile,berry1999voluntary's dataset}

The same settings as the one specified in Figure 6 of conlon2020best were applied. \footnote{Besides globalization strategies and the choice of mappings, the key difference lies in the optimization tolerance level, which is set to 1E-4 here, instead of the default setting 1E-8. The results under the optimization tolerance 1E-8 are shown in the Supplemental Appendix. As discussed, the choice of the inner-loop algorithms has a slight impact on the inner-loop numerical errors, and may largely affect the convergence of the optimization, especially with a tight optimization tolerance level. We discuss that the convergence speed of outer-loop iterations can be sensitive to the inner-loop numerical errors, even when estimated parameters are not sensitive to inner-loop numerical errors. For details, see the Supplemental Appendix.}

Table (ref) shows that the algorithm $\delta$-(1) reduces the number of function evaluations by more than 10%. Combining the spectral/Anderson algorithm further reduces the number of function evaluations by more than 90%.

Convergence speed of $\delta$-(1)

In the estimation using the Berry et al.berry1995automobile,berry1999voluntary's dataset, applying $\delta$-(1) rather than $\delta$-(0) had minimal impact on convergence speed, differing from the Monte Carlo experiments and the estimation using the nevo2001measuring's dataset. However, the results align with the convergence properties of $\delta$-(1). As briefly discussed in Section (ref) and discussed in detail in Appendix (ref), the convergence speed of $\delta$-(1) is significantly influenced by the size of the consumer heterogeneity. As shown in Proposition (ref) in Appendix (ref), $\left\Vert \Phi^{\delta,\gamma=1}(\delta_{1})-\Phi^{\delta,\gamma=1}(\delta_{2})\right\Vert _{\infty}\leq\widetilde{c_{\gamma=1}}\left(\left\Vert \max\{\delta_{1},\delta_{2}\}-\delta_{1}\right\Vert _{\infty}+\left\Vert \max\{\delta_{1},\delta_{2}\}-\delta_{2}\right\Vert _{\infty}\right)\leq2\widetilde{c_{\gamma=1}}\cdot\left\Vert \delta_{1}-\delta_{2}\right\Vert _{\infty}$ holds, where $\widetilde{c_{\gamma=1}}\equiv\sup_{i\in\mathcal{I},J^{*}\subset\mathcal{J},\delta}s_{iJ^{*}}(\delta)-\inf_{i\in\mathcal{I},J^{*}\subset\mathcal{J},\delta}s_{iJ^{*}}(\delta)\leq1$. Because the value of $\widetilde{c_{\gamma=1}}$ itself is not easy to compute when the number of products is large, we alternatively consider $\max_{i\in\mathcal{I}}s_{i0t}(\delta)-\min_{i\in\mathcal{I}}s_{i0t}(\delta)$ at the converged $\delta$, as a rough approximation of the term. By the proposition, $\delta$-(1) is expected to converge slowly when consumer heterogeneity, measured by $\max_{i\in\mathcal{I}}s_{i0t}(\delta)-\min_{i\in\mathcal{I}}s_{i0t}(\delta)\leq1$, is close to 1.

Table (ref) shows the values of $\max_{i\in\mathcal{I}}s_{i0t}(\delta)-\min_{i\in\mathcal{I}}s_{i0t}(\delta)$ from the numerical experiments presented earlier. Using the Berry et al.berry1995automobile,berry1999voluntary's dataset, the value of $\max_{i}s_{i0t}(\delta)-\min_{i}s_{i0t}(\delta)$ is nearly 1. In contrast, in the Monte Carlo experiment with $J=250$ and the estimation using the nevo2001measuring's dataset, the values are far from 1.

From the numerical experiments, we can obtain the following findings:

Considering these observations, it is preferable to first consider applying the new algorithm $\delta$-(1) in the standard fixed-point iteration, and then consider combining the Anderson acceleration.

Dynamic BLP model

Next, the numerical results of inner-loop algorithms for estimating the dynamic BLP models are shown. Here, we consider the model of perfectly durable goods, following the setting considered in Section 5.2 of Sun2019. Consumers are forward-looking, and do not buy a product once they purchase any product.

When consumer $i$ purchases product $j$ at time $t$, they receive utility $U_{ijt}=X_{jt}^{\prime}\theta_{i}+\xi_{jt}+\epsilon_{ijt}$. If no product is purchased given that they do not own any durable product, the utility is $U_{i0t}=\beta E_{t}[V_{it+1}(\Omega_{t+1})|\Omega_{t}]+\epsilon_{i0t}$. Here, $X_{jt},\xi_{jt}$ denotes the observed / unobserved product characteristics of product $j$ at time $t$, and let $\delta_{jt}\equiv X_{jt}^{\prime}E[\theta_{i}]+\xi_{jt}$. $V_{it}(\Omega_{t})$ denotes the (integrated) value function of consumer $i$ not owning any product at time $t$. Besides, let $\mu_{ijt}\equiv X_{jt}^{\prime}\left(\theta_{i}-E[\theta_{i}]\right)$.

Assuming that $\epsilon$ follows i.i.d. type-I extreme value distribution, the CCP of product $j$ for consumer $i$ not holding any product at time $t$ is $s_{ijt}^{(ccp)}=\frac{\exp\left(\delta_{jt}+\mu_{ijt}\right)}{\exp(V_{it}(\Omega_{t}))}$. Moreover, the CCP of not buying any product for consumer $i$ not holding any product at time $t$ is $s_{i0t}^{(ccp)}=\frac{\exp\left(\beta E_{t}[V_{it+1}(\Omega_{t+1})|\Omega_{t}]\right)}{\exp(V_{it}(\Omega_{t}))}.$

Because consumers exit the market after a purchase, $Pr0_{it}$, the fraction of type $i$ consumers not owning any durable product, satisfies $Pr0_{it+1}=Pr0_{it}\cdot s_{i0t}^{(ccp)}$. Note that $Pr0_{it}$ corresponds to $Pr_{it}(x_{it})$ in the general dynamic BLP model discussed in Section (ref).

As in the Monte Carlo experiments of static BLP models, I draw the values of $\theta_{n}$ (standard deviation of random coefficients) by $\theta_{n}\sim U[0,2\theta_{n}^{*}]$, where $\theta_{n}^{*}$ denotes the true parameters. The number of simulation draws to approximate continuous consumer distribution is set to 50, and the process is repeated 50 times. Further details on the data-generating process are provided in the Supplemental Appendix. Regarding future expectations, we consider two specifications: Perfect foresight, and Inclusive value sufficiency. We discuss them below.

Model under perfect foresight

First, consider the setting where consumers have perfect foresight regarding the transition of $\Omega_{t}$, and $\Omega_{t}$ remaining constant after the terminal period $T$, as described in Section (ref). Algorithm (ref) shows the proposed algorithm, which we denote $V$-($\gamma$).

The tolerance level of the convergence is set to 1E-12. Observations from the current data-generating process indicate that the simulated market stabilizes after $t=50$. Consequently, the value of $T$ is set to 50.

To compare the performance of the proposed and the traditional algorithms, both algorithms are run. In the traditional algorithm, the values of $\delta$ and $V$ are jointly updated until convergence. Unlike Algorithm (ref) , Step 1(a) is omitted and instead add a step $\delta_{jt}^{(n+1)}=\delta_{jt}^{(n)}+\log\left(S_{jt}^{(data)}\right)-\log\left(s_{jt}(\delta^{(n)},V^{(n)})\right)-\gamma\left[\log\left(S_{0t}^{(data)}\right)-\log\left(s_{0t}\left(\delta^{(n)},V^{(n)}\right)\right)\right]\ (\gamma=0,1)$ between Steps 2 and 3. In addition, the termination of the iteration is determined by $\left\Vert V^{(n+1)}-V^{(n)}\right\Vert <\epsilon_{V},\left\Vert \delta^{(n+1)}-\delta^{(n)}\right\Vert <\epsilon_{\delta}$. We denote the algorithm as $V\delta$-($\gamma$) (joint). Note that the nested version of the algorithm (upper part of Algorithm (ref)) exhibits poor performance, as shown in the Supplemental Appendix. Thus, the results are excluded in this section.

As in the static BLP models, the results of the algorithms combining the acceleration methods are also presented. For the spectral and SQUAREM methods, “time-dependent step sizes” are introduced, based on the discussion of the variable-type-specific step sizes in Section (ref).\footnote{Supplemental Appendix of the current paper shows the performance of the time-dependent step sizes compared to the case of time-independent step sizes. The results imply that introducing time-dependent step sizes leads to faster convergence by several times.} This specification is motivated by the nonstationarity of the dynamic model. In the case of Algorithm (ref), $V$ is updated by $V_{it}^{(n+1)}\leftarrow\alpha_{t}^{(n)}V_{it}^{(n+1)}+(1-\alpha_{t}^{(n)})V_{it}^{(n)}$, where

$\alpha_{t}^{(n)}\equiv\frac{\left\Vert s_{t}^{(n)}\right\Vert _{2}}{\left\Vert y_{t}^{(n)}\right\Vert _{2}}$, $s_{t}^{(n)}\equiv V_{t}^{(n)}-V_{t}^{(n-1)},y_{t}^{(n)}\equiv F_{t}\left(V^{(n)}\right)-F_{t}\left(V^{(n-1)}\right),F_{t}(V)\equiv\text{\ensuremath{\Phi_{t}(V)-V_{t}}}$. As $\alpha_{t}^{(n)}$ occasionally take too large values (e.g., 100), which might lead to unstable convergence, the maximum value of $\alpha_{t}$ is set to 10.\footnote{Coding time-specific step sizes is not difficult in any programming language for any $T$. In the setting above, $V\equiv\left(V_{t=1}\cdots,V_{t=T}\right)$ is a $|\mathcal{I}|\times T$ dimensional variable. Then, $s^{(n)}\equiv V^{(n)}-V^{(n-1)}$ and $y^{(n)}\equiv F\left(V^{(n)}\right)-F\left(V^{(n-1)}\right)$ are also $|\mathcal{I}|\times T$ dimensional variables. Then, when applying the $\alpha_{S3}$-type step size, we can update the values of $V$ by $V^{(n+1)}\leftarrow\alpha^{(n)}V^{(n+1)}+(1-\alpha^{(n)})V^{(n)}$, where $\alpha^{(n)}\equiv\left(

\right)=\left(

\right)./\left(

\right)$. Here, ``$.$'' denotes the element-wise operation of matrices. $\left(

\right)$ can be computed by $\left(

\right)=$sqrt.$\left(colsums\left(s^{(n)}.\textasciicircum2\right)\right)$. Here, “colsums” denotes the sum of the values in each column of the matrix, and “sqrt” denotes the squared root.}

Table (ref) shows that combining the acceleration methods significantly reduces the number of function evaluations and accelerates convergence. Among these acceleration methods, the Anderson acceleration method outperforms the spectral and SQUAREM. Furthermore, the new algorithms $V$-(0) and $V$-(1) require significantly fewer function evaluations than the traditional algorithms $V\delta$-(0) (joint) and $V\delta$-(1) (joint). In addition, the traditional algorithms $V\delta$-(0) (joint) and $V\delta$-(1) (joint) are not so stable, when the spectral/SQUAREM algorithms are combined.

Note that we cannot clearly see the prominently superior performance of $V$-(1) relative to $V$-(0) in the table. In the current setting, the values of the outside option CCPs are relatively large as shown in the table note, and the advantage of introducing the outside option shares in the updating equations is not necessarily very large. In contrast, the Supplemental Appendix additionally shows numerical results in a setting where the values of the outside option CCPs are smaller. It shows results where $V$-(1) with Anderson acceleration is on average 5$\sim$10% faster than $V$-(0) with Anderson acceleration. Intuitively, the mapping $V$-(0) shares similarity with both the BLP contraction mapping and the value function iteration mapping in dynamic discrete choice models. The convergence speed of the former can be slow as the discount factor gets closer to 1. Hence, the introduction of the outside option shares in the updating equation may not drastically reduce the number of iterations when the discount factor is large, though it mitigates the convergence speed problem associated with the BLP contraction mapping. However, it can reduce the number of iterations, and it is helpful because dynamic BLP estimations are typically computationally burdensome.

Model under Inclusive value sufficiency

Next, we examine a setting where consumers form expectations based on their inclusive values. Generally, $\Omega_{t}$,\-\- which includes product characteristics, is a high-dimensional variable, and fully specifying their stochastic state transitions becomes impractical, especially when dealing with tens of products. Hence, many previous studies (e.g., hendel2006measuring; gowrisankaran2012dynamics) have adopted the idea of inclusive value sufficiency. Specifically, we alternatively use inclusive value $\omega_{it}\equiv\log\left(\sum_{k\in\mathcal{J}_{t}}\exp(\delta_{kt}+\mu_{ikt})\right)$ as the state of the dynamic problem. Under this setting, consumers form expectations based on the transitions of the state $\omega_{it}$.

In this study, the value of $T$ is set to 25, as in Sun2019. Besides, unlike the perfect foresight scenario, the model under inclusive value sufficiency operates within a mostly stationary framework. A common scalar $\alpha^{(n)}$ is used in the spectral/SQUAREM algorithms.

Because the algorithm becomes more complex due to the need to take grid points, the complete steps of the algorithm are detailed in the Supplemental Appendix. Nevertheless, the idea of analytically representing $\delta$ as a function of $V$ remains unchanged.

Table (ref) shows the results. As in the case under perfect foresight, combining the acceleration methods speeds up convergence, and the Anderson acceleration outperforms the others. In addition, the new algorithms $V$-(0) and $V$-(1) are faster than $V\delta$-(0) and $V\delta$-(1), especially when combining the spectral or SQUAREM algorithms.

Conclusion

This study examined computationally fast inner-loop algorithms for estimating static and dynamic BLP models. To minimize the number of inner-loop iterations, the following ideas and insights are proposed: (1). New mapping $\delta^{(n+1)}=\delta^{(n)}+\left(\log\left(S_{j}^{(data)}\right)-\log\left(s_{j}(\delta^{(n)})\right)\right)-\left(\log\left(S_{0}^{(data)}\right)-\log\left(s_{0}(\delta^{(n)})\right)\right)$; (2). Analytically represent the mean product utilities $\delta$ as a function of value functions $V$ and solve for $V$ (for dynamic BLP); (3). Combine an acceleration method of fixed point iterations, especially Anderson acceleration. They are independent and easy to implement. These proposed methods would ease empirical analyses of markets with large datasets or complex demand models under static/dynamic BLP framework, facing problems of computational burden.

Although BLP demand models are considered, whether the ideas and insights can also be applied to other demand models, such as pure characteristics models in berry2007pure, is an interesting topic for further research.

Finally, although this study focuses on improving the convergence speed of the inner-loop algorithms, the mappings with fast convergence could potentially applied to other estimation procedures, such as the MPEC method (e.g., dube2012improving) and sequential estimation algorithms (e.g., lee2015computationally). Future research could examine these possibilities.