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

Takeshi Fukasawa The University of Tokyo<EMAIL_ADDRESS>
I thank Naoshi Doi, Hiroshi Ohashi, Suguru Otani, Katsumi Shimotsu, and
participants at JEMIOW Fall 2023 and APIOC 2023 for their helpful comments.

###### Abstract

This study investigates computationally efficient inner-loop algorithms for
estimating static / dynamic BLP models. It provides the following ideas to
reduce the number of inner-loop iterations: (1). Add a term concerning the
outside option share in the BLP contraction mapping; (2). Analytically
represent mean product utilities as a function of value functions and solve
for the value functions (for dynamic BLP); (3-1). Combine the spectral /
SQUAREM algorithms; (3-2). Choice of the step sizes. These methods are
independent and easy to implement. This study shows good performance of these
ideas by numerical experiments.

##  1 Introduction

Demand estimation is the basis of many economic analyses. The estimation
method based on the random coefficient model proposed by Berry et al., (1995)
(henceforth BLP) has been widely applied in the literature. Though Berry et
al., (1995) 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., Gowrisankaran and Rysman,
(2012), Shcherbakov, (2016)). The models are known as dynamic BLP. In these
models, products’ mean utilities $\delta$ are numerically solved in the inner
loop, and econometricians search for the parameter value minimizing the GMM
objective computed in the outer loop.

Nevertheless, one obstacle to the method is the computational time. The
computational burden can be high because of the need to repeat fixed point
iterations in the inner loop, especially when we consider estimations of
complicated demand models or estimations with large datasets. This study
investigates computationally efficient inner-loop algorithms for estimating
static and dynamic BLP models. Computationally efficient inner-loop algorithms
for estimating static BLP models have been intensively investigated in the
literature, including Kalouptsidi, (2012), Reynaerts et al., (2012), Conlon
and Gortmaker, (2020), Pál and Sándor, (2023). Nevertheless, I show there is
further room for improvement even in static BLP models, let alone dynamic BLP
models.

Regarding static BLP models, I newly introduce a new fixed-point iteration
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)$,
by slightly modifying the traditional BLP contraction mapping
$\delta^{(n+1)}=\delta^{(n)}+\left(\log\left(S_{j}^{(data)}\right)-\log\left(s_{j}(\delta^{(n)})\right)\right)$.111$S_{jt}^{(data)}$
denotes the 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 only the term
$\log\left(S_{0}^{(data)}\right)-\log\left(s_{0}(\delta^{(n)})\right)$, and it
is straightforward to code in any programming languages. When there is no
consumer heterogeneity, applying the new mapping immediately converges to the
solution after one iteration of the mapping for any initial values of
$\delta$. Though there is no guarantee that the mapping is a contraction, it
converged and reduced the number of iterations in the Monte Carlo simulation
settings experimented in the previous studies (Dubé et al., (2012), Lee and
Seo, (2015)) and Nevo, (2001) and Berry et al. (1995, 1999)’s dataset. Because
combining the spectral or SQUAREM algorithm, which has been found to
accelerate convergence in the previous studies (Reynaerts et al., (2012),
Conlon and Gortmaker, (2020), Pál and Sándor, (2023)), further reduces the
number of iterations, my recommendation is to code the new mapping first, and
then consider combining the spectral or SQUAREM algorithms if necessary.

Regarding dynamic BLP models, I propose to analytically represent mean product
utility $\delta$ as a function of value function $\check{V}$, solve for $V$ by
applying a mapping of $V$, and recover $\delta$ using the analytical formula
using $V$. Note that the previous studies (e.g., Gowrisankaran and Rysman,
(2012)) have 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$, and it leads to reducing the
number of iterations in the inner loop. When we further combine the spectral /
SQUAREM algorithm, we can further reduce the number of iterations.

Note that the spectral and SQUAREM algorithms are not only useful for the
inner-loop algorithms for static/dynamic BLP estimation, but also for
efficiently solving any type of fixed point problems222Aguirregabiria and
Marcoux, (2021) discussed that combining the spectral algorithm not only
reduces the number of fixed point iterations but also makes it possible to
find multiple equilibria in the context of nested pseudo-likelihood estimation
of dynamic discrete games., including dynamic models with value functions and
supply-side pricing equilibrium. Built-in packages of the spectral and SQUAREM
algorithms are available in some programming languages. Even when unavailable,
they are easy to code because of their simple formulas. Hence, being
accustomed to using such built-in packages or coding a function to implement
the spectral or SQUAREM algorithms tailoring practitioners’ needs would be
useful. Note that we should be careful about the choice of the step size in
the updating equations when applying the spectral or SQUAREM algorithm, though
not discussed so much in the economics literature. We discuss that the step
size
$\alpha_{S3}\equiv\frac{\left\|s^{(n)}\right\|_{2}}{\left\|y^{(n)}\right\|_{2}}$
is better than other step size choices.333See Section 5 for descriptions of
the notations. Conlon and Gortmaker, (2020) pointed out that the spectral
algorithm performed worse than the SQUAREM algorithm. However, I show
numerical results where the spectral algorithm works mostly as well as the
SQUAREM algorithm, by appropriately choosing the step size. Also, we can
accelerate the convergence by introducing variable-type-specific step sizes in
some models. We discuss it in Section 5.444Note that the SQAUREM algorithm is
a slight extension the spectral algorithm. I found that spectral and SQUAREM
algorithms performed mostly similarly, as shown in Section 6. We discuss which
algorithm (spectral or SQUAREM) practitioners should use in Remark 11 of
Section 5.2.

This study contributes to the literature by mainly providing the following
ideas or insights to reduce the number of inner-loop iterations: (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)$
for static / dynamic BLP models; (2). Analytically represent mean product
utilities $\delta$ as a function of value functions $V$ and solve for $V$ for
dynamic BLP models; (3-1). Combine the spectral or SQUAREM algorithm; (3-2).
Choice of step sizes in the spectral or SQUAREM algorithm. These methods are
independent, and relatively easy to implement. Practitioners can implement
some or all of them when they face difficulties in the computational
burden.555Fukasawa, (2024) applies 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 paper is organized as follows. Section 2 discusses the
relations to the previous studies. In Section 3, we discuss inner-loop
algorithms for estimating static BLP models. In Section 4, we discuss inner-
loop algorithms for estimating dynamic BLP models. Readers interested only in
the static BLP models can skip the section. Section 5 discusses the use of the
spectral and SQUAREM algorithms, which can be used to accelerate the
convergence of general fixed point iterations. Note that the discussion in the
section applies to any type of fixed-point iteration. Section 6 shows the
results of numerical experiments on static / dynamic BLP models. Section 7
concludes.

Appendix A shows additional results and discussions. In Appendix B, 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 has been considered in Kalouptsidi et al., (2020)
and others. Appendix C shows all the proof of the mathematical statements.

##  2 Literature

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

Conlon and Gortmaker, (2020) developed the PyBLP package written in Python for
implementing static BLP-type problems. Though the package is very convenient
to use and state-of-the-art, incorporating the knowledge found in the previous
literature, we cannot use it when we consider a BLP-type problem that cannot
be dealt with in PyBLP.666For instance, we cannot directly use PyBLP to
estimate static limited consideration set models (e.g., Goeree, (2008)) or
dynamic BLP models. In this case, we should code by ourselves, and the tricks
discussed in this study are worth considering.

Kalouptsidi, (2012) proposed an inner-loop estimation algorithm of static BLP
models, which attains faster convergence in the setting with a small number of
discrete consumer types. Doi, (2022) proposed an estimation method of static
BLP models with discrete consumer types without computationally costly fixed
point iterations, given the data availability of total sales for each consumer
type. The essential idea of these studies is that mean utilities or unobserved
product characteristics can be analytically represented as a function of
consumer-type specific value functions, given parameters and market share
data.

The idea is beneficial, especially considering dynamic BLP models. Note that
the direct application of Kalouptsidi, (2012)’s algorithm to static BLP model
does 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 3.2. We discuss these issues in detail in Appendix A.5.

Even though we restrict our attention to improving the inner-loop algorithms
for estimating static / dynamic BLP models, alternative estimation procedures
have been proposed so far.777Besides the studies mentioned above, Bonnet et
al., (2022) 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 Berry and Pakes, (2007). Salanié and
Wolak, (2022) proposed an estimator of static BLP models which 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.  Dubé et al., (2012) proposed using the MPEC
(Mathematical Programming with Equilibrium Constraint) approach for static and
dynamic BLP models. Lee and Seo, (2015) proposed the ABLP (approximate BLP)
method for static BLP models, whose idea is to iterate the process of applying
a function motivated by Newton’s method and estimating parameters. Regarding
dynamic BLP models, Sun and Ishihara, (2019) proposed a Bayesian-type
algorithm to estimate dynamic BLP models.

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
run the algorithm on a computer.888This is the point discussed in Jamshidian
and Jennrich, (1997), in the context of the acceleration of the EM algorithm.
The advantage of the methods proposed in the current paper is that we can
lower the computation costs with small additional thinking costs, because we
do not have to largely change the estimation procedure based on the standard
nested fixed point approach. Unlike the MPEC method, we do not have to install
new softwares or packages to solve constrained optimization problems.999In
MPEC, we sometimes need to specify the sparsity structure to reduce the memory
requirement. In contrast, our current method does not require such a
procedure. Unlike the ABLP method, we do not have to compute analytical
derivatives of functions101010As discussed in Miranda and Fackler, (2004) in
the context of the Newton 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 paper do not require coding derivatives, and are attractive. , and the
proposed algorithms are not restricted to static BLP models. Unlike the
Baysian-type method, we do not have to be familiarized with the using Baysian-
type techniques, such as MCMC (Markov Chain Monte Carlo) methods.

Besides, this study relates to and contributes to the economic literature
discussing the spectral-type algorithms. Reynaerts et al., (2012), Conlon and
Gortmaker, (2020), Pál and Sándor, (2023) found that combining the spectral or
SQUAREM algorithms accelerates the inner-loop convergence of static BLP
estimation. Aguirregabiria and Marcoux, (2021) proposed using the spectral
algorithm in the context of nested pseudo-likelihood estimation of dynamic
discrete games. Though not paid much attention to in the literature, the
choice of the step size largely affects the algorithm’s performance,
especially when combining with the fixed point iteration, and we discuss it in
detail in the current paper.

##  3 Static BLP model

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

###  3.1 Model

In this section, we consider a static BLP model with random coefficients. This
study proposes new algorithms whose convergence properties are better than
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 specific utility of product $j$. $\epsilon$
denotes idiosyncratic utility shocks. Let $\mathcal{J}$ be the set of
products. Then, under the assumption that $\epsilon$ follows Gumbel
distribution and the assumption that observed market share data is equal to
the market share predicted by the model, the following hold:

$\displaystyle S_{j}^{(data)}$ $\displaystyle=$
$\displaystyle\int\frac{\exp\left(\delta_{j}+\mu_{ij}\right)}{1+\sum_{k\in\mathcal{J}}\exp\left(\delta_{k}+\mu_{ik}\right)}dP(i)\approx\sum_{i\in\mathcal{I}}w_{i}\frac{\exp\left(\delta_{j}+\mu_{ij}\right)}{1+\sum_{k\in\mathcal{J}}\exp\left(\delta_{k}+\mu_{ik}\right)},$
(1) $\displaystyle S_{0}^{(data)}$ $\displaystyle=$
$\displaystyle\int\frac{\exp\left(\delta_{j}+\mu_{ij}\right)}{1+\sum_{k\in\mathcal{J}}\exp\left(\delta_{k}+\mu_{ik}\right)}dP(i)\approx\sum_{i\in\mathcal{I}}w_{i}\frac{1}{1+\sum_{k\in\mathcal{J}}\exp\left(\delta_{k}+\mu_{ik}\right)}.$
(2)

Here, we use appropriate discretizations of consumer types to numerically
solve the model, when consumer types are continuously distributed.121212See
Conlon and Gortmaker, (2020)for a review of the efficient discretization
method. $w_{i}$ denotes the fraction of type $i$ consumers, and
$\sum_{i\in\mathcal{I}}w_{i}=1$ should hold. Let $\mathcal{I}$ be the set of
discretized consumer types. If we use $I$ Monte Carlo draws to approximate the
integral, $w_{i}=\frac{1}{I}$ holds.

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

  1. 1.

Inner loop: Given $\theta_{nonlinear}$,

    1. (a)

Compute $\mu$ using $\theta_{nonlinear}$

    2. (b)

Solve for $\delta$ satisfying equations (1) and (2) given $\mu$

    3. (c)

Compute linear parameters $\theta_{linear}^{D}$ that minimize the GMM
objective given $\theta_{nonlinear}^{D}$. Let
$\theta_{linear}^{*}(\theta_{nonlinear})$ be the solution of the minimization
problem, and let $m(\theta_{nonlinear})\equiv
G^{\prime}\left(\theta_{linear}^{*}(\theta_{nonlinear}),\theta_{nonlinear}\right)WG\left(\theta_{linear}^{*}(\theta_{nonlinear}),\theta_{nonlinear}\right)$.

  2. 2.

Outer loop: Search for the value of $\theta_{nonlinear}$ minimizing the GMM
objective $m(\theta_{nonlinear})$.

Algorithm 1 Full estimation algorithm of the static BLP model

Below, we focus on Step 1(b) of Algorithm 1, namely, how to efficiently solve
for $\delta$ satisfying the constraints (1) and (2) given $\mu$.131313Conlon
and Gortmaker, (2020) 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_{linear}$ by a linear GMM
without solving the fixed point problem, as discussed in Berry, (1994).

###  3.2 Mappings of mean product utility $\delta$

To solve for $\delta$ given $\mu$ and observed market shares $S^{(data)}$,
Berry et al., (1995) proposed to iterate the update
$\delta_{j}^{(n+1)}=\delta_{j}^{(n)}+\log\left(S_{j}^{(data)}\right)-\log\left(s_{j}(\delta)\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, define a mapping
$\Phi^{\delta,\gamma}:B_{\delta}\rightarrow B_{\delta}$ for
$\gamma\in\mathbb{R}$ such that141414Though 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.:

$\displaystyle\Phi_{j}^{\delta,\gamma}\left(\delta\right)$
$\displaystyle\equiv$
$\displaystyle\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],$

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}$ corresponds to the traditional BLP contraction
mapping.

Proposition 1 guarantees 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\geq 0}\left(\delta\right)$.

######  Proposition 1.

Solution of $\delta=\Phi^{\delta,\gamma\geq 0}\left(\delta\right)$ satisfies
$S_{j}^{(data)}=s_{j}(\delta)\ \forall j\in\mathcal{J}$.

It is obvious that 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}$. Nevertheless, its converse is not trivial, without formal
proof. Proposition 1 formally shows that the converse actually holds under
$\gamma\geq 0$.

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

Set initial values of $\delta^{(0)}$. Iterate the following
$(n=0,1,2,\cdots)$:

  1. 1.

Compute $\delta_{j}^{(n+1)}=\Phi_{j}^{\delta,\gamma}(\delta^{(n)})$

  2. 2.

Exit the iteration if
$\left\|\delta^{(n+1)}-\delta^{(n)}\right\|<\epsilon_{\delta}$

Algorithm 2 Inner loop algorithm of static BLP using $\Phi^{\delta,\gamma}$

Though, 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 A.1, the size of consumer heterogeneity largely affects
the convergence speed of $\Phi^{\delta,\gamma=1}$. As the size of consumer
heterogeneity gets smaller, the convergence speed gets faster. Interestingly,
when no consumer heterogeneity exists, 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
Berry, (1994). 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 choice of $\delta^{(0)}$.

Note that, 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]-\gamma\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 Appendix A.2. Nevertheless,
the iteration converged when combining the spectral /SQUAREM algorithms, which
we discuss in detail in Section 5. In addition, as shown in Section 6, the
iteration always converged and sometimes led to drastic improvement of
convergence speed in the standard settings of Monte Carlo simulations
experimented in the previous studies (Dubé et al., (2012), Lee and Seo,
(2015)) and the estimations using Nevo, (2001) and Berry et al. (1995, 1999)’s
datasets. Because the mapping $\Phi_{j}^{\delta,\gamma=1}$ is straightforward
to implement once coding the standard BLP contraction mapping, applying the
mapping $\Phi_{j}^{\delta,\gamma=1}$ is worth considering in the static BLP
estimation.

Note that the convergence speed of the BLP contraction mapping
($\Phi^{\delta,\gamma=0}$) is slow when the outside option share is small, as
pointed out in Dubé et al., (2012). In contrast, regarding the new mapping
$\Phi^{\delta,\gamma=1}$, the convergence speed is not so sensitive to the
outside option share, as formally discussed in Appendix A.1 and numerically
shown in Section 6.

######  Remark 1.

Applying the mapping $\Phi^{\delta,\gamma=1}$ requires additionally
subtracting the term
$\left[\log\left(S_{0}^{(data)}\right)-\log\left(s_{0}(\delta)\right)\right]$,
compared to the case with traditional BLP contraction mapping
$\Phi^{\delta,\gamma=0}$. It implies the computational cost of applying a
mapping per iteration increases when we apply $\Phi^{\delta,\gamma=1}$ rather
than $\Phi^{\delta,\gamma=0}$. Nevertheless, the additional computational cost
is generally much smaller than the one for applying $\Phi^{\delta,\gamma=0}$.
Regarding $\Phi^{\delta,\gamma=0}$, the computationally costly part is the
computation of
$s_{j}(\delta)=\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)}$,
which requires the following operations: (1). Summation of a $\delta_{j}$ and
$\mu_{ij}$ for all $i\in\mathcal{I}$ and $j\in\mathcal{J}$; (2). Computation
of exponential of $\delta_{j}+\mu_{ij}$; (3). Summation of a
$|\mathcal{J}|$-dimensional vector
$\left\\{\exp\left(\delta_{j}+\mu_{ij}\right)\right\\}_{j\in\mathcal{J}}$ for
all $i\in\mathcal{I}$; (4). Division of $\exp\left(\delta_{j}+\mu_{ij}\right)$
by $1+\sum_{k\in\mathcal{J}}\exp\left(\delta_{k}+\mu_{ik}\right)$ for all
$i\in\mathcal{I}$ and $j\in\mathcal{J}$; (5). Summation of
$\left\\{w_{i}\frac{\exp\left(\delta_{j}+\mu_{ij}\right)}{1+\sum_{k\in\mathcal{J}}\exp\left(\delta_{k}+\mu_{ik}\right)}\right\\}_{i\in\mathcal{I}}$
for all $j\in\mathcal{J}$. In contrast, computing a scalar
$s_{0}(\delta)=1-\sum_{k\in\mathcal{J}}s_{k}(\delta)$ requires the summation
of $\sum_{k\in\mathcal{J}}s_{k}(\delta)$ using already computed
$s_{k}(\delta)$, and it is much less costly compared to the computation of
$s_{j}(\delta)$. Consequently, applying the mapping $\Phi^{\delta,\gamma=1}$
is not so costly to apply compare to the case of apply the BLP contraction
mapping $\Phi^{\delta,\gamma=0}$.

###  3.3 Mappings of value functions $V$

Though it would be natural to consider mappings of $\delta$ to solve for
$\delta$, we can alternatively solve for consumers’ value function $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:

$\displaystyle V_{i}$ $\displaystyle\equiv$
$\displaystyle\log\left(1+\sum_{k\in\mathcal{J}}\exp\left(\delta_{k}+\mu_{ik}\right)\right)$
(3)

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.
Nevertheless, when considering dynamic BLP models, introducing the approach
solving for $V$ is attractive, because it has a close relationship with value
function iterations used to solve dynamic models, as discussed in Section 4.

The idea of solving for consumer-type specific variables in the inner loop of
static BLP estimations was first proposed in Kalouptsidi, (2012).
Nevertheless, the original algorithms proposed by Kalouptsidi, (2012) have
convergence problems. The algorithms we discuss below overcome the issues by
choosing appropriate mappings that have dualistic connections with the
mappings of $\delta$. For further discussion on Kalouptsidi, (2012)’s
algorithm, see Appendix A.5.

First, by equations (1) and (3), the following equation should hold:

$\displaystyle\delta_{j}$ $\displaystyle=$
$\displaystyle\log\left(S_{j}^{(data)}\right)-\log\left(\sum_{i}w_{i}\exp\left(\mu_{ij}-V_{i}\right)\right)$

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

$\displaystyle\Phi^{V,\gamma}(V)$ $\displaystyle\equiv$
$\displaystyle\iota_{\delta\rightarrow
V}\left(\iota_{V\rightarrow\delta}^{\gamma}\left(V\right)\right),$

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

$\displaystyle\iota_{\delta\rightarrow V,i}\left(\delta\right)$
$\displaystyle\equiv$
$\displaystyle\log\left(1+\sum_{j\in\mathcal{J}}\exp\left(\delta_{j}+\mu_{ij}\right)\right),$

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

$\displaystyle\iota_{V\rightarrow\delta,j}^{\gamma}\left(V\right)$
$\displaystyle\equiv$
$\displaystyle\log\left(S_{j}^{(data)}\right)-\log\left(\sum_{i}w_{i}\exp\left(\mu_{ij}-V_{i}\right)\right)-\gamma\log\left(\frac{S_{0}^{(data)}}{\sum_{i}w_{i}\exp(-V_{i})}\right)$

$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\geq 0}$:

######  Proposition 2.

$\delta$ such that $V=\Phi^{V,\gamma\geq 0}(V),\
\delta=\iota_{V\rightarrow\delta}^{\gamma}(V)$ satisfies
$S_{j}^{(data)}=s_{j}(\delta)\ \forall j\in\mathcal{J}$.

Hence, to solve for $\delta$, we can think of two algorithms. The first is to
apply $\Phi^{\delta,\gamma}$ iteratively and solve for $\delta$ with
appropriate initial values of $\delta$. 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 3 shows the algorithm based on the latter approach.

Set initial values of $V^{(0)}$. Iterate the following $(n=0,1,2,\cdots)$:

  1. 1.

Compute $\delta^{(n)}=\iota_{V\rightarrow\delta}^{\gamma}\left(V^{(n)}\right)$

  2. 2.

Update $V$ by $V^{(n+1)}=\iota_{\delta\rightarrow V}\left(\delta^{(n)}\right)$

  3. 3.

Exit the iteration if $\left\|V^{(n+1)}-V^{(n)}\right\|<\epsilon_{V}$

Algorithm 3 Inner loop algorithm of static BLP using $\Phi^{V,\gamma}$

######  Remark 2.

As in the case of $\Phi^{\delta,\gamma=1}$ and $\Phi^{\delta,\gamma=0}$(cf.
Remark 1), computational cost of applying $\Phi^{V,\gamma=1}$ is not largely
different from $\Phi^{V,\gamma=0}$, because computing a scalar
$s_{0}(V)=\sum_{i}w_{i}\exp(-V_{i})$ is much less costly than computing
$\sum_{j\in\mathcal{J}}\exp\left(\delta_{j}+\mu_{ij}\right)$.

######  Remark 3.

If we define a mapping $\Psi_{V\delta\rightarrow
V}^{\gamma}\left(V,\delta\right)=\log\left(1+\left(\sum_{j\in\mathcal{J}}\exp(\delta_{j}+\mu_{ij})\right)\cdot\left(\frac{\sum_{i}w_{i}\exp(-V_{i})}{S_{0}^{(data)}}\right)^{\gamma}\right)$,
$\Phi^{V,\gamma}(V)=\Psi_{V\delta\rightarrow
V}^{\gamma}\left(V,\iota_{V\rightarrow\delta}\left(V\right)\right)$ holds.
Hence, we can alternatively compute
$\delta^{(n)}=\iota_{V\rightarrow\delta}\left(V^{(n)}\right)$, and update $V$
by $V^{(n+1)}=\Psi_{V\delta\rightarrow
V}^{\gamma}\left(V^{(n)},\delta^{(n)}\right)$. This type of algorithm is used
in the dynamic BLP model discussed in the next section.

We can easily show that there is a dualistic relation between
$\Phi^{\delta,\gamma}$ and $\Phi^{V,\gamma}$. As shown in the Monte Carlo
simulation (Section 6), their convergence speed is mostly the same, as long as
they take a common value $\gamma$. The following proposition shows the formal
statement:

######  Proposition 3.

(Duality of mappings of $\delta$ and $V$) The following equations hold for all
$\gamma\in\mathbb{R}$:

$\displaystyle\Phi^{V,\gamma}$ $\displaystyle=$
$\displaystyle\iota_{\delta\rightarrow
V}\circ\iota_{V\rightarrow\delta}^{\gamma},$
$\displaystyle\Phi^{\delta,\gamma}$ $\displaystyle=$
$\displaystyle\iota_{V\rightarrow\delta}^{\gamma}\circ\iota_{\delta\rightarrow
V}.$

Figure 1 graphically shows the relations 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}$.

Figure 1: Duality between the mappings of $\delta$ and $V$

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

##  4 Dynamic BLP model

###  4.1 Model

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

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

$\displaystyle v_{ijt}(x_{it},\Omega_{t})$ $\displaystyle=$
$\displaystyle\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]+\epsilon_{ijt}(a_{it}).$

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

$\displaystyle v_{i0t}(x_{it},\Omega_{t})$ $\displaystyle=$
$\displaystyle\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]+\epsilon_{i0t}(a_{it}).$

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, such as the characteristics
of products sold in the market. $\Omega_{t}$ may include 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
states $(x_{it+1},\Omega_{t+1})$.

Under the assumption that $\epsilon(a)$ follows i.i.d. mean zero type-I
extreme value distribution and observed market shares are equal to the
predicted market shares, $\delta$ and $V$ satisfy the following equations:

$\displaystyle S_{jt}^{(data)}$ $\displaystyle=$ $\displaystyle
s_{jt}(V,\delta),$ (4) $\displaystyle S_{0t}^{(data)}$ $\displaystyle=$
$\displaystyle s_{0t}(V,\delta),$ (5) $\displaystyle
V_{it}(x_{it},\Omega_{t})$ $\displaystyle=$
$\displaystyle\log\left(\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)+\right.$
$\displaystyle\ \ \ \
\left.\sum_{j\in\mathcal{A}_{t}(x_{it})-\\{0\\}}\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)\right),$

where

$\displaystyle s_{jt}(V,\delta)$ $\displaystyle\equiv$
$\displaystyle\sum_{i\in\mathcal{I}}w_{i}\sum_{x_{it}\in\chi}Pr_{it}(x_{it})\cdot\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})\right)},$
$\displaystyle s_{0t}(V,\delta)$ $\displaystyle\equiv$
$\displaystyle\sum_{i\in\mathcal{I}}w_{i}\sum_{x_{it}\in\chi}Pr_{it}(x_{it})\cdot\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)}.$

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 instance, in the case of 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 product.
$\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 products do not buy any 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})\right)}\equiv
s_{ijt}^{(ccp)}(x_{it},\Omega_{t})$ denotes the conditional choice probability
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 conditional choice probability
of type $i$ consumers at state $(x_{it},\Omega_{t})$ choosing the outside
option at time $t$.

###  4.2 Algorithm

The full estimation algorithm is similar to the static BLP: In the case of
dynamic BLP, we should solve for $\delta$ satisfying equations (4)-(4.1) given
$\mu$.151515Unlike static BLP models, there is no guarantee that $\delta$
satisfying (4)-(4.1) 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 (4)-(4.1) given $\mu$.

Here, we consider a setting of 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$. 161616The
specification was first considered in Conlon, (2012). Analogous specification
was used in Igami, (2017) 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
Fukasawa, (2024). The assumption of rational expectation has been applied in
Kalouptsidi et al., (2020) and empirical studies cited therein.  Though the
specification differs from the one under inclusive value sufficiency, which
has been applied in many previous studies (e.g., Gowrisankaran and Rysman,
(2012)), the algorithm under the current specification is relatively simple
and is the basis of algorithms under alternative specifications. Besides, we
consider the setting where the values of $Pr_{it}(x_{it})$ are known, to make
clear the essence of the algorithm. In Section 6.2, we allow for the case
where the values of $Pr_{it}(x_{it})$ are unknown.

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

$\displaystyle\Phi^{V,\gamma}(V)$ $\displaystyle\equiv$
$\displaystyle\Psi_{V\delta\rightarrow
V}^{\gamma}\left(V,\iota_{V\rightarrow\delta}\left(V\right)\right)$ (7)

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

$\displaystyle\Psi_{V\delta\rightarrow\delta,it}^{\gamma}(V,\delta)$
$\displaystyle\equiv$
$\displaystyle\log\left(\exp\left(\mu_{i0t}(x_{it},\Omega_{t}^{(data)})+\beta
E_{t}\left[V_{it+1}(x_{it+1},\Omega_{t+1}^{(data)})|x_{it},a_{it}=0\right]\right)+\right.$
$\displaystyle\left.\sum_{j\in\mathcal{A}_{t}(x_{it})-\\{0\\}}\exp\left(\delta_{jt}+\mu_{ijt}(x_{it},\Omega_{t}^{(data)})+\beta
E_{t}\left[V_{it+1}(x_{it+1},\Omega_{t+1}^{(data)})|x_{it},a_{it}=j\right]\right)\cdot\left(\frac{s_{0t}(V,\delta)}{S_{0t}^{(data)}}\right)^{\gamma}\right)$

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

$\displaystyle\iota_{V\rightarrow\delta,jt}(V)$ $\displaystyle\equiv$
$\displaystyle\log\left(S_{jt}^{(data)}\right)-\log\left(\sum_{i\in\mathcal{I}}w_{i}\sum_{x_{it}\in\chi}Pr_{it}(x_{it})\cdot\frac{\exp\left(\mu_{ijt}(x_{it},\Omega_{t}^{(data)})+\beta
E_{t}\left[V_{it+1}(x_{it+1},\Omega_{t+1}^{(data)})|x_{it},a_{it}=j\right]\right)}{\exp\left(V_{it}(x_{it},\Omega_{t}^{(data)})\right)}\right)$

$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},\Omega_{t}^{(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$.

Hence, in the case where 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 4 shows the steps:

Set initial values of $V^{(0)}$. Iterate the following $(n=0,1,2,\cdots)$:

  1. 1.

Compute $\delta^{(n)}=\iota_{V\rightarrow\delta}\left(V^{(n)}\right)$

  2. 2.

Update $V$ by $V^{(n+1)}=\Psi_{V\delta\rightarrow
V}^{\gamma}\left(V^{(n)},\delta^{(n)}\right)$

  3. 3.

Exit the iteration if $\left\|V^{(n+1)}-V^{(n)}\right\|<\epsilon_{V}$

Algorithm 4 Proposed dynamic BLP inner-loop algorithm

Note. After the convergence, we should verify that
$\left\|\log(S_{jt}^{(data)})-\log(s_{jt}(\delta,V))\right\|$ is sufficiently
small. See also Remark 4.

######  Remark 4.

Though $\delta$ satisfying equations (4)-(4.1) satisfies
$V=\Phi^{V,\gamma}(V)$, I have not shown that the solution of
$V=\Phi^{V,\gamma}(V)$ actually satisfies equations (4)-(4.1), unlike the
static BLP. If multiple solutions exist in $V=\Phi^{V,\gamma}(V)$, we cannot
rule out the possibility that the solution does not satisfy
$S_{j}^{(data)}=s_{j}(\delta,V)$. Hence, we should verify that
$\left\|\log(S_{jt}^{(data)})-\log(s_{jt}(\delta,V))\right\|$ is sufficiently
small after the convergence of the algorithm, using the solution $(\delta,V)$.
Note that Bellman equation (4.1) holds for any $\gamma\in\mathbb{R}$ under
$S_{jt}^{(data)}=s_{jt}(\delta,V)\ \forall j,t$.

#### Algorithms used in the previous studies

In the previous studies, essentially upper part of Algorithm 5 has been
applied. The idea of the algorithm is simple: Researchers first solve for the
value function $V$ given the values of product mean product utilities
$\delta$. Then, they compute the difference between the observed market shares
and predicted market shares based on the structural model. If the norm of the
difference is larger than the tolerance level, apply the BLP contraction
mapping
$\delta_{jt}^{(n+1)}=\delta_{jt}^{(n)}+\phi\left[\log(S_{j}^{(data)})-\log(s_{j}(\delta))\right]$,171717In
dynamic BLP models, there is no guarantee that the iteration
$\delta_{jt}^{(n+1)}=\delta_{jt}^{(n)}+\left[\log(S_{j}^{(data)})-\log(s_{j}(\delta^{(n)}))\right]$
converges. Hence, introducing a dampening parameter $\phi\in(0,1]$ is
sometimes necessary to stabilize the convergence. Proposition 2 of Sun and
Ishihara, (2019) shows that the mapping
$\delta_{jt}^{(n+1)}=\delta_{jt}^{(n)}+\phi\left[\log(S_{j}^{(data)})-\log(s_{j}(\delta^{(n)}))\right]$
is a contraction if $\phi$ is "sufficiently small" and additional conditions
hold. and repeat the process again. Nevertheless, the computational burden can
be high because of the existence of nested loops.

We can think of an alternative algorithm, which is shown in the lower part of
Algorithm 5 . The algorithm avoids the nested loops by jointly updating the
values of $\delta$ and $V$ based on the fixed point constraint.181818The idea
of jointly updating two types of variables can also be found in Pakes and
McGuire, (1994) 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 4 ), $\delta$ is analytically represented as a function of $V$, and
we only have to solve for $V$, which leads to faster convergence. I show
superior performance of the proposed algorithm in Section 6.2.

  * •

The case with nested loops:

Set initial values of $\delta^{(0)}$. Iterate the following
$(n=0,1,2,\cdots)$:

    1. 1.

Set initial values of $V^{(0,n)}$. Iterate the following $(m=0,1,2,\cdots)$:

      1. (a)

Update $V$ by $V^{(m+1,n)}=\Psi_{V\delta\rightarrow
V}^{\gamma=0}\left(V^{(m,n)},\delta^{(n)}\right)$

      2. (b)

If $\left\|V^{(m+1,n)}-V^{(m,n)}\right\|<\epsilon$, let $V^{(n)*}$ be
$V^{(m+1,n)}$ and exit

    2. 2.

Update $\delta$ by $\delta_{jt}^{(n+1)}=\Phi_{\delta
V\rightarrow\delta}^{\gamma=0,\phi}(\text{$\delta^{(n)},V^{(n)*}$})\equiv\delta_{jt}^{(n)}+\phi\left[\log\left(S_{jt}^{(data)}\right)-\log\left(s_{jt}(\delta^{(n)},V^{(n)*})\right)\right]$,
where $\phi\in(0,1]$ is a dampening parameter.

    3. 3.

Exit the iteration if
$\left\|\delta^{(n+1)}-\delta^{(n)}\right\|<\epsilon_{\delta}$

  * •

The case with one loop:

Set initial values of $\delta^{(0)},V^{(0)}$. Iterate the following
$(n=0,1,2,\cdots)$:

    1. 1.

Update $\delta$ by
$\delta^{(n+1)}=\delta_{jt}^{(n)}+\phi\left[\log\left(S_{jt}^{(data)}\right)-\log\left(s_{jt}(\delta^{(n)},V^{(n)})\right)\right]$

    2. 2.

Update $V$ by
$V^{(n+1)}=\Psi_{V\delta\rightarrow\delta}^{\gamma=0,\phi}\left(V^{(n)},\delta^{(n)}\right)$

    3. 3.

Exit the iteration if
$\left\|\delta^{(n+1)}-\delta^{(n)}\right\|<\epsilon_{\delta}$ and
$\left\|V^{(n+1)}-V^{(n)}\right\|<\epsilon_{V}$

Algorithm 5 Traditional dynamic BLP inner-loop Algorithm

######  Remark 5.

Schiraldi, (2011), Gowrisankaran and Rysman, (2012) considered durable goods
models with replacement demand, where $\mu_{i0t}(x_{it},\Omega_{t})$ depends
on $\left\\{\delta_{kt}\right\\}_{k\in\mathcal{J}_{t}}$, depending on the
product holding represented by $x_{it}$. Even in such settings, equation (4)
still holds, and we can derive an equation on $V$ without using $\delta$.

######  Remark 6.

Typically, the values of
$\left\\{Pr_{it}(x_{it})\right\\}_{i\in\mathcal{I},t=1,\cdots,T,x_{it}\in\chi}$
are unknown. Hence, we need to impose assumptions on the form of
$\left\\{Pr_{it}(x_{it})\right\\}_{i\in\mathcal{I},t=1,\cdots,T,x_{it}\in\chi}$
at the initial period and also solve for the variables. Here, suppose that the
values of
$\left\\{Pr_{it=1}(x_{it=1})\right\\}_{i\in\mathcal{I},x_{it}\in\chi}$ are
known.191919We can alternatively assume that $t=1$ is at the stationary state,
as in Fukasawa, (2024). Kasahara and Shimotsu, (2009) discuss the issue in the
context of standard dynamic discrete choice model estimation with unobserved
heterogeneity. Generally,
$\left\\{Pr_{it+1}(x_{it+1})\right\\}_{i\in\mathcal{I},t=1,\cdots,T,x_{it+1}\in\chi}$
satisfies the following:

$\displaystyle Pr_{it+1}(x_{it+1})$ $\displaystyle=$
$\displaystyle\sum_{x_{it}}Pr{}_{it}(x_{it})\cdot\sum_{j\in\mathcal{A}_{t}(x_{it})}F(x_{it+1}|x_{it},a_{it}=j)\cdot
s_{ijt}^{(ccp)}(x_{it},\Omega_{t}),$ (8)

where $F(x_{it+1}|x_{it},a_{it}=j)$ denotes the state transition probability
of state $x_{it}$. Then, as discussed in more specific model settings in
Section 6.2, $Pr_{it}(x_{it})$ can be represented as a function of V. Hence,
even when the values of $Pr_{it}(x_{it})$ are unknown, we can slightly alter
Algorithm 4 to accommodate unknown $Pr_{it}(x_{it})$.202020Intuitively,
$s_{ijt}^{(ccp)}(x_{it},\Omega_{t})$ is a function of $V$ and $\delta$.
Because delta can be analytically represented as a function of $V$,
$Pr_{it}(x_{it})\ (t=1\cdots,T)$ can be sequentially represented as a function
of $V$.

##  5 Spectral algorithm and SQUAREM algorithm

Here, we discuss the spectral algorithm (Section 5.1) and the SQUAREM
algorithm (Section 5.2), which can be used to accelerate fixed-point
iterations. Note that the discussion here applies to any fixed-point
iteration, including BLP contraction mapping, value function iterations and
others. Though not discussed much in the economics literature, the choice of
step size plays a vital role in the algorithm’s convergence, and we discuss it
in detail in Section 5.1. Besides, this study proposes to apply variable-type-
specific step sizes in the spectral / SQUAREM algorithms. The idea is
effective, especially in nonstationary dynamic models. We discuss it in
Section 5.3.

We discuss that the step size $\alpha_{S3}$ is better than other step size
choices. In the context of BLP inner-loop algorithms, Conlon and Gortmaker,
(2020) compared the spectral algorithm with step size $\alpha_{S2}$ and the
SQUAREM algorithm with step size $\alpha_{S3}$. (See also Remark 11 in Section
5.2.) They reported that the SQUAREM algorithm performs better than the
spectral algorithm. Nevertheless, when we apply the step size $\alpha_{S3}$ in
the spectral algorithm, the spectral and SQUAREM algorithms performed mostly
similarly, as shown in Section 5.1. Hence, it seems that what is more
important is the choice of the step size, rather than the type of the
algorithm (spectral or SQUAREM).

The discussions in Sections 5.1 and 5.2 largely rely on Varadhan and Roland,
(2008), discussing the application of the spectral and SQUAREM algorithm for
the EM algorithm, which requires fixed point iterations. In this paper, we
discuss these algorithms by comparing the concepts which would be familiar to
economists, such as the extrapolation / relaxation methods for fixed point
iterations. We discuss them in detail, because a deeper understanding of the
algorithms is essential for choosing appropriate step sizes.212121To make the
comparison easier, I modify the notations in Varadhan and Roland, (2008),
including the sign of the step size. Because some of the previous studies and
built-in packages use the sign of the step size differently, readers should be
careful about it. In the economics literature, Aguirregabiria and Marcoux,
(2021) also describes the spectral algorithm in the context of nested pseudo-
likelihood estimation of dynamic discrete choice models. This study
additionally discuss the role of the step sizes.

###  5.1 Spectral algorithm

In this subsection, we first discuss the idea of the spectral algorithm with
fixed-point iteration, and then discuss the relations with the
extrapolation/relaxation methods, which would be more familiar to economists.
Then, we discuss desirable the sign of the step size and the choice of the
step size.

Spectral algorithm is an algorithm designed to solve nonlinear equations. When
we want to solve a nonlinear equation $F(x)=0$, $x^{(n)}$ is iteratively
updated in the following way given initial values of $x^{(0)}$ and the choice
of a scalar $\alpha^{(n)}$:

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

Here, generally $\alpha^{(n)}$ is taken based on the values of $s^{(n)}\equiv
x^{(n)}-x^{(n-1)}$ and $y^{(n)}\equiv F(x^{(n)})-F(x^{(n-1)})$, and it can
vary across $n=0,1,2,\cdots$. It is known that the 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, and the spectral algorithm chooses the scalar
$\alpha^{(n)}$ is chosen so that $\alpha^{(n)}$ “approximates” the matrix
$-\left(\nabla F(x^{(n)})\right)^{-1}$ well. Even though it would be possible
to apply the Newton’s method to solve nonlinear equations, 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. In addition,
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, using the spectral algorithm is attractive from the viewpoint of
simplicity and computational cost. We will discuss the details of the choice
of $\alpha^{(n)}$ later.

The idea of the spectral algorithm can be utilized in the context of fixed
point iteration. Suppose we want to solve a fixed point constraint
$x=\Phi(x)$. Then, by letting $F(x)=\Phi(x)-x$, we should iteratively update
the values of $x$ to solve $x=\Phi(x)$ in the following way:

$\displaystyle x^{(n+1)}$ $\displaystyle=$ $\displaystyle
x^{(n)}+\alpha^{(n)}\left(\Phi\left(x^{(n)}\right)-x^{(n)}\right)\
(n=0,1,2,\cdots)$ (9)

#### Relations with the Extrapolation method / Relaxation method

As discussed on page 166 of Judd, (1998), one way to accelerate or stabilize
fixed point iteration is the extrapolation method. The method is also known as
the relaxation method.222222In the industrial organization literature, the
relaxation method was used in Kasahara and Shimotsu, (2012) to stabilize the
convergence of nested pseudo-likelihood estimation of dynamic discrete games.
Besides, the idea of introducing a dampening parameter $\phi\in(0,1]$ in the
BLP contraction mapping of the traditional dynamic BLP inner-loop algorithm
(cf. Algorithm 5) is similar to the relaxation method. If a variable $x$
satisfies a fixed point constraint $x=\Phi(x)$, we can iteratively update the
value of $x$ by $x^{(n+1)}=\alpha\Phi(x^{(n)})+(1-\alpha)x^{(n)}$ given a pre-
specified scalar $\alpha$. If the iteration under $\alpha=1$ converges slowly,
we should try a larger value $\alpha>1$, hoping to accelerate the convergence.
If the iteration under $\alpha=1$ is unstable, we should try a smaller value
$0<\alpha<1$ to stabilize the convergence. One problem of the extrapolation
method combined with the fixed point iterations is that the appropriate value
of $\alpha$ is unclear. We typically have to manually try several values of
$\alpha$ to efficiently solve the problem.

Regarding the spectral algorithm with the updating equation (9), the updating
equation can be reformulated as follows:

$\displaystyle x^{(n+1)}$ $\displaystyle=$
$\displaystyle\alpha^{(n)}\Phi(x^{(n)})+(1-\alpha^{(n)})x^{(n)}\
(n=0,1,2,\cdots)$

The only difference with the extrapolation method is that a scalar $\alpha$
can vary, and the value is determined flexibly based on the predetermined
formula and the convergence process. The spectral algorithm avoids manually
finding an appropriate value of $\alpha$, as long as the appropriate
predetermined formula is chosen.

#### Desirable sign of the step size $\alpha^{(n)}$

As discussed above, $\alpha^{(n)}$ should be chosen to accelerate convergence
unless it gets unstable. Generally, choosing $\alpha^{(n)}<0$ is not a good
choice, because it might lead to the opposite direction of the original fixed
point iteration $x^{(n+1)}=\Phi(x^{(n)})$. Figure 2 illustrates it, and the
following proposition shows the formal statement:232323The proposition is not
shown in Varadhan and Roland, (2008), and I show it to clarify the point.

######  Proposition 4.

Suppose that $\Phi:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is a contraction
mapping of modulus $K\in[0,1)$ on $L^{2}$ metric space. Let $x^{*}$ be the
unique solution of $x=\Phi(x)$. Then, the following holds for $\alpha<0$:

$\displaystyle\left\|\left(\alpha\Phi(x)+(1-\alpha)x\right)-x^{*}\right\|_{2}$
$\displaystyle>$ $\displaystyle\left\|x-x^{*}\right\|_{2}$

Proposition 4 shows that $\alpha^{(n)}\Phi(x)+(1-\alpha^{(n)})x$ is further
from the solution $x^{*}$ than $x$, and it might lead to divergence of the
iterations, when $\Phi$ is a contraction mapping and $\alpha^{(n)}$ takes a
negative value242424Note that we cannot deny the possibility that the
iteration converges faster when setting $\alpha^{(n)}<0$ for some $n$.
Nevertheless, in general, it would be natural to expect that setting
$\alpha^{(n)}<0$ for some $n$ would make the convergence slower and more
unstable.. Even when there is no guarantee that $\Phi$ is a contraction,
choosing a negative $\alpha$ might not be adequate when $\Phi$ is “close to” a
contraction.

Figure 2: Idea of the extrapolation method and the spectral algorithm

Notes.

$x^{*}$ denotes the solution of the fixed point problem $x=\Phi(x)$.

Proposition 4 implies
$x^{(n+1)}=\alpha^{(n)}\Phi(x^{(n)})+(1-\alpha^{(n)})x^{(n)}$ is further than
$x^{(n)}$ from $x^{*}$, when $\alpha^{(n)}<0$ holds. In contrast, if a
positive value $\alpha^{(n)}$ is appropriately chosen, $x^{(n+1)}$ can be
closer to $x^{*}$ than the case of $x^{(n)}$.

#### Choice of the step size

##### Step size $\alpha_{S1}$

In the spectral algorithm, $x^{(n+1)}=x^{(n)}+\alpha^{(n)}F(x^{(n)})$ and
$x^{(n)}=x^{(n-1)}+\alpha^{(n-1)}F(x^{(n-1)})$ hold. Because we would like to
compute values of $x^{(n+1)}$ and $x^{(n)}$ close to the solution $x^{*}$,
$x^{(n+1)}$ and $x^{(n)}$ should be close as much as possible. Hence, we
consider the following minimization problem:

$\displaystyle\min_{\alpha^{(n)}}\left\|x^{(n+1)}-x^{(n)}\right\|_{2}^{2}$
$\displaystyle=$
$\displaystyle\left\|s^{(n)}+\alpha^{(n)}y^{(n)}\right\|_{2}^{2},$

where $s^{(n)}=x^{(n)}-x^{(n-1)}$ and $y^{(n)}=F(x^{(n)})-F(x^{(n-1)})$. The
solution of the minimization problem is
$\grave{\alpha_{S1}^{(n)}=-\frac{s^{(n)\prime}y^{(n)}}{y^{(n)\prime}y^{(n)}}}$.
This is the first choice of $\alpha^{(n)}$ considered in Barzilai and Borwein,
(1988). Varadhan and Roland, (2008) pointed out that $\alpha_{S1}$ might take
a value sufficiently close to zero when $s^{(n)\prime}y^{(n)}\approx 0$. In
this case, $x^{(n+1)}\approx x^{(n)}$ holds, and the iteration might suspend
at the point far from $x^{*}$, and the specification is problematic.

##### Step size $\alpha_{S2}$

We can alternatively consider the following problem:

$\displaystyle\min_{\alpha^{(n)}}\frac{\left\|x^{(n+1)}-x^{(n)}\right\|_{2}^{2}}{\left(\alpha^{(n)}\right)^{2}}$
$\displaystyle=$
$\displaystyle\frac{\left\|s^{(n)}+\alpha^{(n)}y^{(n)}\right\|_{2}^{2}}{\left(\alpha^{(n)}\right)^{2}}.$

The solution is
$\grave{\alpha_{S2}^{(n)}=-\frac{s^{(n)\prime}s^{(n)}}{s^{(n)\prime}y^{(n)}}}$.
This is the second choice of $\alpha^{(n)}$ considered in Barzilai and
Borwein, (1988), and used in La Cruz et al., (2006). Varadhan and Roland,
(2008) pointed out that $\alpha_{S2}$ might take a too large value when
$s^{(n)\prime}y^{(n)}\approx 0$. In this case, the iteration might diverge,
and the specification is also problematic.

##### Step size $\alpha_{S3}$

We can alternatively consider the following problem:

$\displaystyle\min_{\alpha^{(n)}}\frac{\left\|x^{(n+1)}-x^{(n)}\right\|_{2}^{2}}{\left|\alpha^{(n)}\right|}$
$\displaystyle=$
$\displaystyle\frac{\left\|s^{(n)}+\alpha^{(n)}y^{(n)}\right\|_{2}^{2}}{\left|\alpha^{(n)}\right|}.$

If we restrict the solution to be positive, the solution is
$\alpha_{S3}^{(n)}=\frac{\left\|s^{(n)}\right\|_{2}}{\left\|y^{(n)}\right\|_{2}}$.
Note that $\alpha_{S3}=\sqrt{\alpha_{S1}\alpha_{S2}}$ holds under
$s^{(n)\prime}y^{(n)}<0$. This is the specification proposed in Varadhan and
Roland, (2008), and the problems associated with the case of
$s^{(n)\prime}y^{(n)}\approx 0$ would not occur when we use
$\alpha_{S3}$252525In the numerical literature, Dai et al., (2015) and Huang
and Wan, (2017) also discussed the effectiveness of the step size
$\alpha_{S3}$..

These step sizes are regarded as the approximation of $-\left(\nabla
F(x^{(n)})\right)^{-1}$ that appears in the Newton update. To see this
intuitively, consider a setting where $x$ is a scalar. Then, $y$ and $s$ are
also scalars, and $\alpha_{S1}=\alpha_{S2}=-\frac{s^{(n)}}{y^{(n)}}$ holds.
Regarding
$\alpha_{S3}\equiv\frac{\left\|s^{(n)}\right\|_{2}}{\left\|y^{(n)}\right\|_{2}}=\frac{|s^{(n)}|}{|y^{(n)}|}$,
we can show that $-\frac{s^{(n)}}{y^{(n)}}$ is positive for $F(x)=\Phi(x)-x$
when $|\Phi^{\prime}(x)|<1$ holds, which is the sufficient condition for
$\Phi$ to be a contraction mapping,262626We can show
$-\frac{s^{(n)}}{y^{(n)}}=-\frac{F\left(x^{(n+1)}\right)-F\left(x^{(n)}\right)}{x^{(n+1)}-x^{(n)}}=1-\frac{\Phi\left(x^{(n+1)}\right)-\Phi\left(x^{(n)}\right)}{x^{(n+1)}-x^{(n)}}>0$
under the condition. To show the claim, It is sufficient to show that
$\frac{\Phi\left(x\right)-\Phi\left(y\right)}{x-y}<1\ \forall
x,y\in\mathbb{R}$ if $\left|\Phi^{\prime}\left(x\right)\right|<1\ \forall
x\in\mathbb{R}$. Let $x<y$ without loss of generality. Then,
$\left|\frac{\Phi\left(x\right)-\Phi\left(y\right)}{x-y}\right|=\frac{\left|\int_{u=0}^{u=x-y}\Phi^{\prime}\left(x+u\right)du\right|}{x-y}=\frac{\int_{u=0}^{u=x-y}\left|\Phi^{\prime}\left(x+u\right)\right|du}{x-y}<\frac{\int_{u=0}^{u=x-y}1du}{x-y}=1$
holds. and
$\alpha_{S3}=\alpha_{S1}=\alpha_{S2}=\frac{s^{(n)}}{y^{(n)}}=-\text{$\left(\frac{F\left(x^{(n+1)}\right)-F\left(x^{(n)}\right)}{x^{(n+1)}-x^{(n)}}\right)$}^{-1}>0$
holds in the setting.
$\left(\frac{F\left(x^{(n+1)}\right)-F\left(x^{(n)}\right)}{x^{(n+1)}-x^{(n)}}\right)^{-1}$
can be regarded as the approximation of the derivative $\left(\nabla
F(x^{(n)})\right)^{-1}$ in the Newton’s method, and the spectral update
“approximates” the Newton update.272727In the case $x$ is a scalar, the method
is equivalent to the secant method. See Section 3.4 of Miranda and Fackler,
(2004). Besides, Section 5 of Varadhan and Roland, (2008) discusses why the
spectral-type algorithms work well by giving a simple example.

In summary, when we apply the spectral algorithm with step size
$\alpha_{S3}\equiv\frac{\left\|s^{(n)}\right\|_{2}}{\left\|y^{(n)}\right\|_{2}}$,
the detailed algorithm is as follows:

  1. 1.

Set initial values of $x^{(0)}$. Choose $\alpha^{(0)}$ and tolerance level
$\epsilon$.

  2. 2.

Iterate the following $(n=0,1,2,\cdots)$:

    1. (a)

Compute $\Phi\left(x^{(n)}\right)$ and
$F\left(x^{(n)}\right)=\Phi\left(x^{(n)}\right)-x^{(n)}$.

    2. (b)

Compute $s^{(n)}\equiv x^{(n)}-x^{(n-1)},y^{(n)}\equiv
F(x^{(n)})-F(x^{(n-1)}),\text{and
}\alpha^{(n)}=\frac{\left\|s^{(n)}\right\|_{2}}{\left\|y^{(n)}\right\|_{2}}\text{\
if\ }n\geq 1.$

    3. (c)

Compute $x^{(n+1)}=x^{(n)}+\alpha^{(n)}F(x^{(n)})$

    4. (d)

If $\left\|\Phi(x^{(n)})-x^{(n)}\right\|<\epsilon$, exit the iteration.
Otherwise, go back to Step 2(a).

Algorithm 6 Spectral algorithm with fixed point mapping

######  Remark 7.

In general, there is no guarantee that the iteration converges. Hence,
globalization strategies, which guarantee convergence regardless of the choice
of any initial points $x^{(0)}$, are proposed in the literature on numerical
analysis. The idea is as follows. When using the step size
$\alpha^{(n)}=\alpha_{S3}^{(n)}=\frac{\left\|s^{(n)}\right\|_{2}}{\left\|y^{(n)}\right\|_{2}}$,
we first try
$x^{(n+1)}=\alpha_{S3}^{(n)}\Phi(x^{(n)})+(1-\alpha_{S3}^{(n)})x^{(n)}$, and
if $\left\|F(x^{(n+1)})\right\|^{2}$ is “too large”, we next try
$\alpha^{(n)}=\pm\lambda\alpha_{S3}^{(n)}$ ($\lambda\in(0,1)$) to compute
$x^{(n+1)}=\alpha_{S3}^{(n)}\Phi(x^{(n)})+(1-\alpha_{S3}^{(n)})x^{(n)}$ and
$\left\|F(x^{(n+1)})\right\|^{2}$. We repeat the process until we can find a
plausible value $\alpha^{(n)}$ where $\left\|F(x^{(n+1)})\right\|^{2}$ is
“relatively small”.282828The process is known as backtracing in the numerical
literature. For instance, La Cruz et al., (2006) proposed a nonmonotone
globalization strategy, which applies to the spectral algorithm without
computing the derivatives of $F$. The algorithm is known as DF-SANE
(Derivative-Free Spectral Algorithm for Nonlinear Equation). Note that the
globalization strategy is not restricted to a specific choice of
$\alpha^{(n)}$. In the algorithm above, the step should be introduced between
Step 2(c) and 2(d).

Nevertheless, in our BLP estimation, the algorithm without globalization
strategies always converged, at least as far as I experimented. The same
experience was also reported in Pál and Sándor, (2023), and it seems not
essential to introduce globalization strategies, at least in the BLP
estimation.

######  Remark 8.

We cannot deny the possibility that $\alpha^{(n)}\rightarrow 0\
(n\rightarrow\infty)$. In that case,
$\left\|x^{(n+1)}-x^{(n)}\right\|=\alpha^{(n)}\left\|\Phi(x^{(n)})-x^{(n)}\right\|\rightarrow
0$ might occur even when $\left\|\Phi(x^{(n)})-x^{(n)}\right\|$ is not close
to zero. Hence, checking $\left\|\Phi(x^{(n)})-x^{(n)}\right\|<\epsilon$,
rather than $\left\|x^{(n+1)}-x^{(n)}\right\|<\epsilon$, is essential.

###  5.2 SQUAREM algorithm

Varadhan and Roland, (2008) proposed the SQUAREM (Squared Polynomial
Extrapolation Methods) algorithm. The algorithm is designed to solve a fixed
point problem $x=\Phi(x)$, and the updating equation is as follows:

$\displaystyle x^{(n+1)}$ $\displaystyle=$ $\displaystyle
x^{(n)}+2\alpha^{(n)}s^{(n)}+\left(\alpha^{(n)}\right)^{2}y^{(n)}$ (10)

where $s^{(n)}\equiv\Phi\left(\delta^{(n)}\right)-\delta^{(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)}$.

Note that the SQUAREM algorithm is closely related to the spectral algorithm.
Updating equation (10) can be reformulated as:

$\displaystyle x^{(n+1)}$ $\displaystyle=$
$\displaystyle\left(1-\alpha^{(n)}\right)\left[(1-\alpha^{(n)})x^{(n)}+\alpha^{(n)}\Phi\left(x^{(n)}\right)\right]+\alpha^{(n)}\left[(1-\alpha^{(n)})\Phi\left(x^{(n)}\right)+\alpha^{(n)}\Phi^{2}\left(x^{(n)}\right)\right]$
$\displaystyle=$
$\displaystyle\left(1-\alpha^{(n)}\right)\Psi_{spectral}\left(x^{(n)};\Phi,\alpha^{(n)}\right)+\alpha^{(n)}\Psi_{spectral}\left(\Phi\left(x^{(n)}\right);\Phi,\alpha^{(n)}\right).$

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 outputs of spectral update
$\Psi_{spectral}\left(x^{(n)};\Phi,\alpha^{(n)}\right)$ and
$\Psi_{spectral}\left(\Phi\left(x^{(n)}\right);\Phi,\alpha^{(n)}\right)$.

Based on the similar ideas, we can choose the step size
$\alpha_{S1}\equiv-\frac{s^{(n)\prime}y^{(n)}}{y^{(n)\prime}y^{(n)}}$,
$\alpha_{S2}\equiv-\frac{s^{(n)\prime}s^{(n)}}{s^{(n)\prime}y^{(n)}}$, and
$\alpha_{S3}\equiv\frac{\left\|s^{(n)}\right\|_{2}}{\left\|y^{(n)}\right\|_{2}}$.
Regarding the SQUAREM algorithm, it is unclear whether choosing a negative
step size $\alpha$ leads to a further point, unlike Proposition 4 in the case
of the spectral algorithm. Nevertheless, because the SQUAREM algorithm is
closely related to the spectral algorithm as discussed above, choosing a
positive step size $\alpha$ would be desirable.

In summary, when we apply the SQUAREM algorithm with step size
$\alpha_{S3}\equiv\frac{\left\|s^{(n)}\right\|_{2}}{\left\|y^{(n)}\right\|_{2}}$,
the detailed algorithm is as follows:

  1. 1.

Set initial values of $x^{(0)}$ and tolerance level $\epsilon$.

  2. 2.

Iterate the following $(n=0,1,2,\cdots)$:

    1. (a)

Compute $\Phi(x^{(n)})$ and
$\Phi^{2}(x^{(n)})\equiv\Phi\left(\Phi(x^{(n)})\right)$

    2. (b)

Compute
$s^{(n)}\equiv\Phi\left(x^{(n)}\right)-x^{(n)},y^{(n)}\equiv\Phi^{2}\left(x^{(n)}\right)-2\Phi\left(x^{(n)}\right)+x^{(n)}$,
and
$\alpha^{(n)}=\frac{\left\|s^{(n)}\right\|_{2}}{\left\|y^{(n)}\right\|_{2}}$.

    3. (c)

Compute
$x^{(n+1)}=x^{(n)}+2\alpha^{(n)}r^{(n)}+\left(\alpha^{(n)}\right)^{2}v^{(n)}$.

    4. (d)

If $\left\|\Phi(x^{(n)})-x^{(n)}\right\|<\epsilon$, exit the iteration.
Otherwise, go back to step 2(a).

Algorithm 7 SQUAREM algorithm

######  Remark 9.

As in the spectral algorithm, we cannot deny the possibility that
$\alpha^{(n)}\rightarrow 0\ (n\rightarrow\infty)$. In that case,
$\left\|x^{(n+1)}-x^{(n)}\right\|=\alpha^{(n)}\left\|\Phi(x^{(n)})-x^{(n)}\right\|\rightarrow
0$ might occur even when $\left\|\Phi(x^{(n)})-x^{(n)}\right\|$ is not close
to zero. Hence, checking $\left\|\Phi(x^{(n)})-x^{(n)}\right\|<\epsilon$,
rather than $\left\|x^{(n+1)}-x^{(n)}\right\|<\epsilon$, is essential.

######  Remark 10.

Because the previous studies investigating the spectral and SQUAREM algorithms
applied different step sizes in their numerical experiments, some of which are
confusing, I summarize the settings here.

Reynaerts et al., (2012) chose the step size $\alpha_{S1}$ for the spectral
algorithm, and $\alpha_{S3}$ for the SQUAREM algorithm.292929They mentioned
that they chose the default step size of squarem function provided by the
SQUAREM package in R. Because the default step size is $\alpha_{S3}$, it seems
they chose the step size
$\alpha_{S3}=\frac{\left\|s^{(n)}\right\|_{2}}{\left\|y^{(n)}\right\|_{2}}$,
rather than
$sgn\left(s^{(n)\prime}y^{(n)}\right)\frac{\left\|s^{(n)}\right\|_{2}}{\left\|y^{(n)}\right\|_{2}}$
in equation (15) of the paper.  Conlon and Gortmaker, (2020) mentioned the
step size $\alpha_{S1}$ in equation (22) for the SQUAREM algorithm, but the
default step size in PyBLP they developed is $\alpha_{S3}$. Regarding the
spectral algorithm, they implicitly chose the step size $\alpha_{S2}$, because
the dfsane function provided in SciPy package PyBLP package relies on only
allows the step size $\alpha_{S2}$. Pál and Sándor, (2023) chose the step size
$\alpha_{S1}$ for the spectral and the SQUAREM algorithms.

Reynaerts et al., (2012) and Conlon and Gortmaker, (2020) discussed that the
SQUAREM algorithm performs better than the spectral algorithm. Nevertheless,
they did not use the step size $\alpha_{S3}$ in the spectral algorithm, might
lead to the worse performance of the spectral algorithm. As shown in Section
6, the spectral and the SQUAREM algorithms perform similarly in the Monte
Carlo estimation experimented in the previous studies, and estimations using
Nevo, (2001) and Berry et al. (1995, 1999)’s datasets.

######  Remark 11.

As numerically shown in Section 6, the spectral and SQUAREM algorithms perform
similarly in the BLP inner-loop algorithms, when we apply the step size
$\alpha_{S3}$. Hence, both algorithms are worth considering because of their
easy coding and application. In contrast, Varadhan and Roland, (2008) showed a
numerical example where SQUAREM works much better than the spectral, even when
choosing the step size $\alpha_{S3}$. Hence, which algorithm practitioners
should use depends on their preference. Regarding the spectral algorithm,
there is a large number of studies on the globalization strategies to attain
global convergence (e.g., La Cruz et al., (2006), Huang and Wan, (2017)), and
we can utilize these strategies when we make much of global convergence. Also,
as shown in Table 13 of Appendix A.3, built-in packages of the spectral
algorithm are available in many programming languages, though the step size
$\alpha_{S3}$ is not available. In contrast, regarding the SQUAREM algorithm,
globalization strategies do not exist for general fixed point mappings. Also,
built-in packages of the SQUAREM algorithm with step size $\alpha_{S3}$ are
available only in PyBLP (Python) and R language. If practitioners make much of
global convergence, though not introduced in the numerical experiments in
Section 6, they should use the spectral algorithm.

###  5.3 Variable-type-specific step size

Even though the standard spectral /SQUAREM algorithms use a scalar
$\alpha^{(n)}$, I newly introduce the idea of variable-type-specific step
sizes to the spectral / SQUAREM algorithms,303030The idea of variable-specific
step size can also be found in the diagonal Broyden Jacobian approximation
method (Waziri et al., (2010)), though it is motivated by the quasi-Newton
method rather than the spectral algorithm.  which is sometimes useful to
accelerate the convergence, especially nonstationary dynamic models.313131See
also the results in Appendix A.8. Note that readers only interested in
standard static BLP models can skip the section so far. Though I discuss it in
the context of the spectral algorithm, the idea can be used in the SQUAREM
algorithm.

To introduce the idea, suppose we would like to solve a nonlinear equation
$F(x)=\left(\begin{array}[]{c}F_{1}(x)\\\
F_{2}(x)\end{array}\right)=0\in\mathbb{R}^{n_{1}+n_{2}}$, where
$x=\left(\begin{array}[]{c}x_{1}\\\
x_{2}\end{array}\right)\in\mathbb{R}^{n_{1}+n_{2}}$. When using the standard
step size $\alpha_{S3}$, we set
$\alpha^{(n)}=\frac{\left\|s^{(n)}\right\|_{2}}{\left\|y^{(n)}\right\|_{2}}$,
where $s^{(n)}\equiv\left(\begin{array}[]{c}x_{1}^{(n)}\\\
x_{2}^{(n)}\end{array}\right)-\left(\begin{array}[]{c}x_{1}^{(n-1)}\\\
x_{2}^{(n-1)}\end{array}\right)$,
$y^{(n)}\equiv\left(\begin{array}[]{c}F_{1}(x^{(n)})\\\
F_{2}(x^{(n)})\end{array}\right)-\left(\begin{array}[]{c}F_{1}(x^{(n-1)})\\\
F_{2}(x^{(n-1)})\end{array}\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(\begin{array}[]{c}x_{1}\\\ x_{2}\end{array}\right)$ by
$x^{(n+1)}=\left(\begin{array}[]{c}x_{1}^{(n+1)}\\\
x_{2}^{(n+1)}\end{array}\right)=\left(\begin{array}[]{c}x_{1}^{(n)}\\\
x_{2}^{(n)}\end{array}\right)+\left(\begin{array}[]{c}\alpha_{1}^{(n)}F_{1}(x^{(n)})\\\
\alpha_{2}^{(n)}F_{2}(x^{(n)})\end{array}\right),$ where
$\alpha_{m}^{(n)}=\frac{\left\|s_{m}^{(n)}\right\|_{2}}{\left\|y_{m}^{(n)}\right\|_{2}}$,
$s_{m}^{(n)}\equiv x_{m}^{(n)}-x_{m}^{(n-1)}$, $y_{m}^{(n)}=\equiv
F(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 in Section 5.1, $\alpha^{(n)}$ is chosen to accelerate the
convergence. In principle, any choices of $\alpha^{(n)}$ are allowed in the
spectral algorithm as long as it works, and its numerical performance is more
critical. In Appendix A.6, we discuss the analytical justifications of
variable-type-specific step sizes.

##  6 Numerical Experiments

In this section, I show the results of numerical experiments for static /
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\|\Phi(x)-x\right\|_{\infty}$ takes an infinite or NaN value, or the
number of function evaluations reaches a pre-determined maximum.

###  6.1 Static BLP model

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

####  6.1.1 Monte Carlo simulation

To compare the performance of the algorithms, I conducted the Monte Carlo
simulation. The setting is the same as Lee and Seo, (2015) and Dubé et al.,
(2012), and similar to Pál and Sándor, (2023). Consumer $i$’s utility when
choosing product $j$ is $U_{ij}=X_{j}\theta_{i}+\xi_{j}+\epsilon_{ij}\
(j=1,\cdots,J)$, and utility when not buying anything is
$U_{i0}=\epsilon_{i0}$. Let $\delta_{j}\equiv X_{j}\cdot
E[\theta_{i}]+\xi_{j}$. For details of the data-generating process, see
Appendix A.10.1. 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.

Because we focus on the performance of BLP inner-loop algorithms, we measure
the performance of each algorithm to solve for $\delta$ given parameter values
$\theta_{nonlinear}$ (standard deviation of random coefficients $\theta_{i}$).
I first construct a market share data based on the true nonlinear parameter
values $\theta_{nonlinear}^{*}$ and true $\delta^{*}$. In the estimation
process, analysts do not know the true nonlinear parameter values. To imitate
the setting, I draw the parameter values $\theta_{nonlinear}\sim
U[0,2\theta_{nonlinear}^{*}]$, and measure the performance of each algorithm
under the parameters. I repeat the process 50 times.

Here, we mainly focus on the number of function evaluations, rather than the
CPU time. Though the CPU time represents the exact time required until the
termination of the iterations, generally it largely depends on the syntax of
each algorithm and characteristics of each programming language. In contrast,
the numbers of function evaluations are not sensitive to the syntaxes and
types of programming languages, and we mainly focus on the values.323232In the
case of the standard fixed point iterations 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 1 and 2 shows the results in the settings where $J=25$ and $J=250$. The
results shows that the mappings $\delta$-(1) and $V$-(1) require much fewer
iterations than $\delta$-(0) and $V$-(0). The superior performance of
$\delta$-(1) and $V$-(1) is prominent under $J=250$. As shortly discussed in
Section 3.2 and formally discussed in Appendix A.1, the convergence speed of
the BLP contraction mapping ($\delta$-(0)) can be slow when the outside option
share is small, but $\delta$-(1) is not. As shown in the tables’ notes, the
mean outside option share is 0.850 under $J=25$, and 0.307 under $J=250$. The
observations are in line with the discussion. Besides, we can numerically see
the duality between the mappings of $\delta$ and $V$: the numbers of
iterations are mostly the same for $\delta$-($\gamma$) and $V$-($\gamma$).

When we further combine the spectral / SQUAREM algorithms, we can further
reduce the number of function evaluations. 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. For the results, see Appendix A.3.

Note that the value of
$\left\|\log(S_{j}^{(data)})-\log(s_{j})\right\|_{\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\|\Phi^{\delta,\gamma=1}(\delta)-\delta\right\|=\left\|\left(\log(S_{j}^{(data)})-\log(s_{j})\right)-\left(\log(S_{0}^{(data)})-\log(s_{0})\right)\right\|<\epsilon$,
not $\left\|\log(S_{j}^{(data)})-\log(s_{j})\right\|<\epsilon$. Hence, to
confirm that the iteration yielded the solution with sufficiently small
$\left\|\log(S_{j}^{(data)})-\log(s_{j})\right\|$, in the tables I show the
mean values of
$DIST\equiv\left\|\log(S_{j}^{(data)})-\log(s_{j})\right\|_{\infty}$ computed
using the solution of the iterations. I also show the percentage of settings
where the value of $DIST$ is smaller than 1E-12. The results show that all the
algorithms except for the slow $\delta$-(0), $V$-(0) algorithms satisfy the
precision constraint.

Table 1: Results of Monte Carlo simulation (Static BLP model; Continuous
consumer types; $J=25$)

| $J$ | Func. Eval. | Mean | Conv. | Mean | $DIST<\epsilon_{tol}$
---|---|---|---|---|---|---
| Mean | Min. | 25th | Median. | 75th | Max. | CPU time (s) | (%) | $\log_{10}\left(DIST\right)$ | (%)
$\delta$-(0) (BLP) | 25 | 41.28 | 5 | 8 | 12 | 31 | 582 | 0.0087 | 100 | -15.6 | 100
$\delta$-(0) (BLP) + Spectral | 25 | 11.54 | 5 | 7 | 8 | 13 | 39 | 0.00298 | 100 | -16.4 | 100
$\delta$-(0) (BLP) + SQUAREM | 25 | 12.66 | 5 | 8 | 9 | 14 | 52 | 0.003 | 100 | -16.5 | 100
$\delta$-(1) | 25 | 13.7 | 5 | 7 | 9 | 19 | 37 | 0.00324 | 100 | -15.9 | 100
$\delta$-(1) + Spectral | 25 | 8.7 | 5 | 7 | 7 | 10 | 18 | 0.00242 | 100 | -16.6 | 100
$\delta$-(1) + SQUAREM | 25 | 8.8 | 5 | 6 | 8 | 10 | 18 | 0.00218 | 100 | -16.7 | 100
$V$-(0) | 25 | 42.4 | 4 | 8 | 11 | 31 | 593 | 0.01162 | 100 | -15.8 | 100
$V$-(0) + Spectral | 25 | 10.48 | 4 | 6 | 8 | 12 | 33 | 0.00354 | 100 | -16.6 | 100
$V$-(0) + SQUAREM | 25 | 11.38 | 4 | 7 | 8.5 | 13 | 48 | 0.00344 | 100 | -16.4 | 100
$V$-(1) | 25 | 13.72 | 4 | 7 | 8.5 | 18 | 40 | 0.00416 | 100 | -16 | 100
$V$-(1) + Spectral | 25 | 9.02 | 4 | 6 | 7 | 12 | 21 | 0.00298 | 100 | -16.5 | 100
$V$-(1) + SQUAREM | 25 | 9.62 | 4 | 6 | 7 | 13 | 24 | 0.003 | 100 | -16.5 | 100

Notes. $DIST\equiv\left\|\log(S_{j}^{(data)})-\log(s_{j})\right\|_{\infty}$,
$\epsilon_{tol}=$1E-12.

Number of simulation draws is set to 1000.

Maximum number of function evaluations is set to 1000.

Mean Outside option share is 0.850.

Table 2: Results of Monte Carlo simulation (Static BLP model; Continuous
consumer types; $J=250$)

| $J$ | Func. Evals. | Mean | Conv. | Mean | $DIST<\epsilon_{tol}$
---|---|---|---|---|---|---
| Mean | Min. | 25th | Median. | 75th | Max. | CPU time (s) | (%) | $\log_{10}\left(DIST\right)$ | (%)
$\delta$-(0) (BLP) | 250 | 205.1 | 26 | 60 | 109.5 | 192 | 1000 | 0.77564 | 94 | -13.3 | 96
$\delta$-(0) (BLP) + Spectral | 250 | 36.2 | 15 | 23 | 29.5 | 38 | 150 | 0.14018 | 100 | -14.6 | 100
$\delta$-(0) (BLP) + SQUAREM | 250 | 40.3 | 15 | 24 | 30 | 42 | 155 | 0.15276 | 100 | -14.4 | 100
$\delta$-(1) | 250 | 25.04 | 11 | 19 | 25 | 28 | 58 | 0.09596 | 100 | -14.6 | 100
$\delta$-(1) + Spectral | 250 | 14.44 | 8 | 12 | 14.5 | 16 | 23 | 0.05616 | 100 | -15.2 | 100
$\delta$-(1) + SQUAREM | 250 | 14.78 | 9 | 12 | 15 | 16 | 24 | 0.05644 | 100 | -15.4 | 100
$V$-(0) | 250 | 214.32 | 25 | 62 | 117 | 214 | 1000 | 0.95844 | 94 | -13.4 | 96
$V$-(0) + Spectral | 250 | 32.1 | 13 | 22 | 27 | 36 | 89 | 0.14696 | 100 | -14.8 | 100
$V$-(0) + SQUAREM | 250 | 31.38 | 15 | 21 | 28 | 37 | 78 | 0.14162 | 100 | -14.7 | 100
$V$-(1) | 250 | 26.32 | 11 | 19 | 26 | 29 | 62 | 0.11768 | 100 | -14.8 | 100
$V$-(1) + Spectral | 250 | 17.64 | 10 | 14 | 16 | 21 | 31 | 0.07654 | 100 | -15.1 | 100
$V$-(1) + SQUAREM | 250 | 17.24 | 11 | 14 | 16 | 20 | 30 | 0.07592 | 100 | -15.2 | 100

Notes. $DIST\equiv\left\|\log(S_{j}^{(data)})-\log(s_{j})\right\|_{\infty}$,
$\epsilon_{tol}=$1E-12.

Number of simulation draws is set to 1000.

Maximum number of function evaluations is set to 1000.

Mean Outside option share is 0.307.

####  6.1.2 Replications of Nevo, (2001) and Berry et al. (1995, 1999)

To further experiment the performance of the inner-loop algorithms, we use
Nevo, (2001) and Berry et al.(1995, 1999)’s dataset and estimate demand
parameters under different inner-loop algorithms.

In the numerical experiments presented here, I utilized the PyBLP package
(version 0.13.0) in Python (Conlon and Gortmaker, (2020)). Because the package
itself does not allow the new mapping $\delta$-(1), I modify the original
PyBLP source code to introduce the new mapping $\delta$-(1). In addition,
regarding the spectral algorithm, we cannot choose step sizes other than
$\alpha_{S2}$, because the dfsane function in the SciPy package the PyBLP
relies on only allows $\alpha_{S2}$. Hence, I also modify the dfsane function
in the SciPy package to allow different step sizes. Besides, because the
current dfsane function does not allow the iterations without globalization
strategies, I modify the code to implement iterations without globalization
strategies. Note that it rarely changes the results, as shown in Appendix A.3.

#####  Nevo, (2001)’s dataset

I apply mostly the same settings as the one specified in Figure 5 of Conlon
and Gortmaker, (2020), 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 3 shows the results. The results show that applying the algorithm
$\delta$-(1) is more than two times faster than the algorithm $\delta$-(0). If
we further combine the spectral or SQUAREM algorithms, we can further reduce
the number of total function evaluations by more than 2 times. Note that the
values of the total objective evaluations and the objective are common,
regardless of the choice of the inner-loop algorithm. It implies the choice of
the inner-loop algorithm did not affect the outer-loop estimation.

Table 3: Estimation results using Nevo, (2001)’s dataset

| Mean feval | Total obj eval | Total feval | Objective
---|---|---|---|---
$\delta$-(1) + Spectral | 19.209 | 57 | 102921 | 4.562
$\delta$-(1) + SQUAREM | 19.911 | 57 | 106683 | 4.562
$\delta$-(1) | 43.288 | 57 | 231937 | 4.562
$\delta$-(0) + Spectral | 27.796 | 57 | 148931 | 4.562
$\delta$-(0) + SQUAREM | 26.866 | 57 | 143949 | 4.562
$\delta$-(0) | 94.714 | 57 | 507475 | 4.562

Notes.

“Total obj eval” denotes the total number of objective evaluations in the GMM
estimation.

“Total feval” denotes the total number of function evaluations in the GMM
estimation.

“Objective” denotes the GMM objective value.

“Mean feval” denotes the mean number of function evaluations, defined by Total
feval / (Number of markets $\times$ Total objective evalutions).

##### Berry et al.(1995, 1999)’s dataset

I apply mostly the same settings as the one specified in Figure 6 of Conlon
and Gortmaker, (2020). Other than the existence of globalization strategies,
and the choice of mappings, the difference is the optimization tolerance
level. I set it to 1E-4 here, rather than the default setting 1E-8. The
results under the optimization tolerance 1E-8 is shown in Appendix A.4. As
discussed in the section, the choice of inner-loop algorithms slightly affects
inner-loop numerical errors, and it might largely affect the convergence of
the optimization, especially when setting 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 Section
A.4.

Table 4 shows the results. The results show that applying the algorithm
$\delta$-(1) reduces the number of function evaluations by more than 10%. If
we further combine the spectral / SQUAREM algorithm, we can further reduce the
number of function evaluations by more than 70%.

Table 4: Estimation results using Berry et al.(1995, 1999)’s dataset

| Mean feval | Total obj eval | Total feval | Objective
---|---|---|---|---
$\delta$-(1) + Spectral | 44.153 | 70 | 61814 | 497.336
$\delta$-(1) + SQUAREM | 48.229 | 70 | 67520 | 497.336
$\delta$-(1) | 216.938 | 70 | 303713 | 497.336
$\delta$-(0) + Spectral | 47.518 | 70 | 66525 | 497.336
$\delta$-(0) + SQUAREM | 46.7 | 70 | 65380 | 497.336
$\delta$-(0) | 246.358 | 70 | 344901 | 497.336

Notes.

“Total obj eval” denotes the total number of objective evaluations in the GMM
estimation.

“Total feval” denotes the total number of function evaluations in the GMM
estimation.

“Objective” denotes the GMM objective value.

“Mean feval” denotes the mean number of function evaluations, defined by Total
feval / (Number of markets $\times$ Total objective evalutions).

Outer loop tolerance is set to 1E-4.

Overall, we can obtain the following findings from the results of the
numerical experiments:

  * •

Algorithm $\delta$-(1) works better than $\delta$-(0), especially when the
outside option share is small

  * •

Combining the spectral or SQUAREM algorithms accelerates convergence

  * •

The performance of the spectral and the SQUAREM algorithms under the step
sizes $\alpha_{S3}$ are not largely different

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

###  6.2 Dynamic BLP model

Next, I show the numerical results of inner-loop algorithms for estimating
dynamic BLP models. Here, we consider the model of perfectly durable goods,
following the setting used in Section 5.2 of Sun and Ishihara, (2019).
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 the utility
$U_{ijt}=x_{jt}^{\prime}\theta_{i}+\xi_{jt}+\epsilon_{ijt}$. When not buying a
product given that they do not own any durable product, they receive a utility
$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}$ denotes the integrated value
function of consumer $i$ at time $t$. Besides, let $\mu_{ijt}\equiv
x_{jt}^{\prime}\left(\theta_{i}-E[\theta_{i}]\right)$.

Under the assumption that $\epsilon$ follows i.i.d. type-I extreme value
distribution, conditional choice probability of choosing 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})}$,
and conditional choice probability 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})}.$

Because consumers exit the market once a purchase is made, $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 4.

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

####  6.2.1 Model under perfect foresight

First, we consider the setting where consumers form perfect foresight on the
transition of $\Omega_{t}$, and $\Omega_{t}$ is constant after the terminal
period $T$, as considered in Section 4. Algorithm 8 shows the proposed
algorithm, which we denote $V$-($\gamma$).

The tolerance level of the convergence is set to 1E-12. Because I found that
the simulated market became mostly stable after $t=50$ under the current data-
generating process, I set the value of $T$ to 50.

Set initial values of $V^{(0)}$. Iterate the following $(n=0,1,2,\cdots)$:

  1. 1.

For $t=1:T$,

    1. (a)

Using $Pr0_{it}$ and $V_{it}^{(n)}$, compute
$\delta_{jt}^{(n)}=\iota_{V\rightarrow\delta,jt}^{\gamma}\left(V^{(n)}\right)=\log\left(S_{jt}^{(data)}\right)-\log\left(\sum_{i}w_{i}Pr0_{it}\cdot\frac{\exp\left(\mu_{ijt}\right)}{\exp\left(V_{it}\right)}\right)$

    2. (b)

Compute
$s_{ijt}^{(ccp)}=\frac{\exp\left(\delta_{jt}^{(n)}+\mu_{ijt}\right)}{\exp(V_{it}^{(n)})}$
for $i\in\mathcal{I},j\in\mathcal{J}$

    3. (c)

Compute $s_{i0t}^{(ccp)}=1-\sum_{j\in\mathcal{J}_{t}}s_{ijt}^{(ccp)}$

    4. (d)

Update $Pr0_{it+1}=Pr0_{it}\cdot s_{i0t}^{(ccp)}$

  2. 2.

Update $V$ by $V_{it}^{(n+1)}=\Psi_{V\delta\rightarrow
V,it}^{\gamma}\left(V^{(n)},\delta^{(n)}\right)=\log\left(\exp\left(\beta
V_{it+1}^{(n)}\right)+\sum_{j\in\mathcal{J}_{t}}\exp\left(\delta_{jt}^{(n)}+\mu_{ijt}\right)\cdot\left(\frac{s_{0t}(V^{(n)})}{S_{0t}^{(data)}}\right)^{\gamma}\right)$
for $t=1,\cdots,T$. Here, let $V_{iT+1}^{(n)}=V_{iT}^{(n)}$. $s_{0t}(V^{(n)})$
is computed by $s_{0t}(V^{(n)})=\sum_{i}w_{i}\frac{\exp\left(\beta
V_{it+1}^{(n)}\right)}{\exp\left(V_{it}\right)}$.

  3. 3.

Exit the iteration if $\left\|V^{(n+1)}-V^{(n)}\right\|<\epsilon_{V}$

Algorithm 8 Inner-loop Algorithm of dynamic BLP (Perfectly durable goods;
Perfect foresight)

To compare the performance of the proposed algorithm with the traditional
algorithm, I run both algorithms. Regarding the traditional algorithm, we
jointly update the values of $\delta$ and until convergence. Unlike Algorithm
8 , we skip Step 1(a), 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. Also, the termination of the iteration is
determined by
$\left\|V^{(n+1)}-V^{(n)}\right\|<\epsilon_{V},\left\|\delta^{(n+1)}-\delta^{(n)}\right\|<\epsilon_{\delta}$.
We denote the algorithm as $V\delta$-($\gamma$) (joint). Note that the nested
version of the algorithm (upper part of Algorithm 5) performs poorly as shown
in Appendix A.9, and I do not show the results in this section.

As in static BLP models, I also show the results of the algorithms combining
the spectral / SQUAREM algorithms. Here, we introduce "time-dependent step
sizes", based on the discussion of variable-type-specific step sizes in
Section 5.3.333333Appendix A.8 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 8, $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\|s_{t}^{(n)}\right\|_{2}}{\left\|y_{t}^{(n)}\right\|_{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{$\Phi_{t}(V)-V_{t}$}$.
Because I find $\alpha_{t}^{(n)}$ sometimes take too large values (e.g., 100),
which might lead to unstable convergence, I set the maximum value of
$\alpha_{t}$ to 10. 343434Coding time-specific step sizes is not difficult in
any programming language even when $T$ is large. 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(\begin{array}[]{ccc}\alpha_{t=1}^{(n)}&\cdots&\alpha_{t=T}^{(n)}\end{array}\right)=\left(\begin{array}[]{ccc}\left\|s_{t=1}^{(n)}\right\|_{2}&\cdots&\left\|s_{t=T}^{(n)}\right\|_{2}\end{array}\right)./\left(\begin{array}[]{ccc}\left\|y_{t=1}^{(n)}\right\|_{2}&\cdots&\left\|y_{t=T}^{(n)}\right\|_{2}\end{array}\right)$.
Here, “.” denotes the element-wise operation of matrices.
$\left(\begin{array}[]{ccc}\left\|s_{t=1}^{(n)}\right\|_{2}&\cdots&\left\|s_{t=T}^{(n)}\right\|_{2}\end{array}\right)$
can be computed by
$\left(\begin{array}[]{ccc}\left\|s_{t=1}^{(n)}\right\|_{2}&\cdots&\left\|s_{t=T}^{(n)}\right\|_{2}\end{array}\right)=$sqrt.$\left(\text{colsums}\left(s^{(n)}\text{.\textasciicircum
2}\right)\right)$. Here, “colsums” denotes the sum of the values in each
column of the matrix, and “sqrt” denotes the squared root.

Table 5 shows the results. The results show that combining the spectral or
SQUAREM algorithms leads to a much smaller number of function evaluations and
faster convergence. Also, the new algorithms $V$-(0) and $V$-(1) require much
fewer function evaluations. We can also see that the traditional algorithms
$V\delta$-(0) (joint) and $V\delta$-(1) (joint) are not so stable, especially
when combining the spectral / SQUAREM algorithms.

Note that we cannot clearly see the superior performance of $V$-(1) relative
to $V$-(0) in the table. In the current setting, the values of outside option
CCPs are roughly between 0.95 and 0.99, and the advantage of introducing the
outside option shares in the updating equations is not so large.

Table 5: Results of Dynamic BLP Monte Carlo simulation (Perfectly durable
goods; Perfect foresight)

| $J$ | Func. Evals. | Mean | Conv. | Mean | $DIST<\epsilon_{tol}$
---|---|---|---|---|---|---
| Mean | Min. | 25th | Median. | 75th | Max. | CPU time (s) | (%) | $\log_{10}\left(DIST\right)$ | (%)
$V$-(0) | 25 | 2598.65 | 2513 | 2559.5 | 2611.5 | 2638.5 | 2652 | 3.535 | 100 | -14.2 | 100
$V$-(0) + Spectral | 25 | 491.95 | 413 | 456.5 | 471 | 490.5 | 732 | 0.73355 | 100 | -16.8 | 100
$V$-(0) + SQUAREM | 25 | 379.35 | 351 | 362 | 374.5 | 383 | 492 | 0.6032 | 100 | -16.3 | 100
$V$-(1) | 25 | 2539.5 | 2461 | 2498 | 2547.5 | 2572.5 | 2607 | 3.4376 | 100 | -14.2 | 100
$V$-(1) + Spectral | 25 | 486.15 | 416 | 459.5 | 473 | 491 | 700 | 0.70785 | 100 | -16.9 | 100
$V$-(1) + SQUAREM | 25 | 358.95 | 82 | 353.5 | 367 | 380 | 490 | 0.46135 | 95 | -15.5 | 95
$V\delta$-(0) (joint) | 25 | 2787.65 | 2578 | 2660.5 | 2796.5 | 2877 | 3000 | 2.5892 | 90 | -14 | 100
$V\delta$-(0) (joint) + Spectral | 25 | 1574.65 | 31 | 1286 | 1735.5 | 2049 | 3000 | 1.49735 | 80 | -12.7 | 85
$V\delta$-(0) (joint) + SQUAREM | 25 | 1471.4 | 982 | 1174 | 1439 | 1647 | 2334 | 1.44845 | 100 | -14.6 | 100
$V\delta$-(1) (joint) | 25 | 2781.45 | 2575 | 2656 | 2788.5 | 2867 | 3000 | 2.69165 | 90 | -14 | 100
$V\delta$-(1) (joint) + Spectral | 25 | 1914.95 | 1239 | 1495.5 | 1861 | 2138 | 3000 | 2.0351 | 90 | -14.8 | 100
$V\delta$-(1) (joint) + SQUAREM | 25 | 1395.6 | 12 | 1166 | 1442 | 1650 | 2354 | 1.31035 | 95 | -13.9 | 95

Notes. $DIST\equiv\left\|\log(S_{j}^{(data)})-\log(s_{j})\right\|_{\infty}$,
$\epsilon_{tol}=$1E-12.

Maximum number of function evaluations is set to 3000.

####  6.2.2 Model under Inclusive value sufficiency

Next, we consider the 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 is not so realistic, especially when there are
tens of products. Hence, many previous studies (e.g., Hendel and Nevo, (2006),
Gowrisankaran and Rysman, (2012)) have employed the idea of inclusive value
sufficiency. More 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 the setting, consumers form
expectations based on the transitions of the state $\omega_{it}$.

In this study, I set the value of $T$ to 25, as in Sun and Ishihara, (2019).
Besides, unlike the case under perfect foresight, the model under inclusive
value sufficiency is confined to a mostly stationary framework, and I use a
common scalar $\alpha^{(n)}$ in the spectral / SQUAREM algorithms.

Because the algorithm is more complicated because we need to take grid points,
I show the complete steps of the algorithm in Appendix A.10.2. Nevertheless,
the idea of analytically representing $\delta$ as a function of $V$ is the
same.

Table 6 shows the results. As in the case under perfect foresight, combining
the spectral / SQUAREM algorithms leads to much faster convergence. 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 / SQUAREM
algorithms.

Table 6: Results of Dynamic BLP Monte Carlo simulation (Perfectly durable
goods; Inclusive value sufficiency)

| $J$ | Func. Evals. | Mean | Conv. | Mean | $DIST<\epsilon_{tol}$
---|---|---|---|---|---|---
| Mean | Min. | 25th | Median. | 75th | Max. | CPU time (s) | (%) | $\log_{10}\left(DIST\right)$ | (%)
$V$-(0) | 25 | 1347.25 | 755 | 1264 | 1368 | 1480.5 | 1799 | 2.3244 | 100 | -14.9 | 100
$V$-(0) + Spectral | 25 | 229.15 | 167 | 209 | 232.5 | 246.5 | 300 | 0.40905 | 100 | -15.7 | 100
$V$-(0) + SQUAREM | 25 | 224.35 | 178 | 198 | 223.5 | 241.5 | 309 | 0.39385 | 100 | -15.5 | 100
$V$-(1) | 25 | 1340.55 | 772 | 1253 | 1362 | 1478.5 | 1799 | 2.3489 | 100 | -15 | 100
$V$-(1) + Spectral | 25 | 232.2 | 168 | 215 | 234.5 | 250.5 | 292 | 0.41725 | 100 | -15.6 | 100
$V$-(1) + SQUAREM | 25 | 219.15 | 176 | 206.5 | 218.5 | 229.5 | 310 | 0.38325 | 100 | -15.7 | 100
$V\delta$-(0) (joint) | 25 | 1374.7 | 705 | 1297.5 | 1399 | 1496.5 | 1806 | 2.11635 | 100 | -14.9 | 100
$V\delta$-(0) (joint) + Spectral | 25 | 565.15 | 314 | 507 | 583 | 635 | 706 | 0.9322 | 100 | -15.2 | 100
$V\delta$-(0) (joint) + SQUAREM | 25 | 562.4 | 370 | 533 | 569 | 617 | 720 | 0.86865 | 100 | -15.1 | 100
$V\delta$-(1) (joint) | 25 | 1371.3 | 712 | 1294.5 | 1394.5 | 1495 | 1805 | 2.165 | 100 | -14.9 | 100
$V\delta$-(1) (joint) + Spectral | 25 | 569.7 | 280 | 523.5 | 575.5 | 625 | 795 | 0.93755 | 100 | -15.2 | 100
$V\delta$-(1) (joint) + SQUAREM | 25 | 563.7 | 358 | 525 | 559 | 615 | 746 | 0.8755 | 100 | -15 | 100

Notes. $DIST\equiv\left\|\log(S_{j}^{(data)})-\log(s_{j})\right\|_{\infty}$,
$\epsilon_{tol}=$1E-12.

Maximum number of function evaluations is set to 3000.

##  7 Conclusion

This study has investigated computationally fast inner-loop algorithms for
estimating static / dynamic BLP models. I provide the following ideas and
insights to reduce the number of inner-loop iterations: (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)$
for static / dynamic BLP models; (2). Analytically represent mean product
utilities $\delta$ as a function of value functions $V$ and solve for $V$ for
dynamic BLP models; (3-1). Combine the spectral or SQUAREM algorithm; (3-2).
Choice of step sizes in the spectral or SQUAREM algorithm. They are
independent and easy to implement. The proposed methods would ease empirical
analyses of markets with large datasets or complex demand models using the
framework of static / dynamic BLP, facing problems of computational burden.

Though we consider BLP demand models, whether the ideas and insights can also
be applied to other demand models, such as pure characteristics models in
Berry and Pakes, (2007), is an interesting topic for further research.

Finally, though we focus on improving the convergence speed of the inner-loop
algorithms, we might be able to utilize the mappings with fast convergence in
other estimation procedures, such as MPEC method (e.g., Dubé et al., (2012))
and sequential estimation algorithms (e.g., Lee and Seo, (2015)). I leave it
for further research.

## References

  * Aguirregabiria and Marcoux, (2021) Aguirregabiria, V. and Marcoux, M. (2021).  Imposing equilibrium restrictions in the estimation of dynamic discrete games.  Quantitative Economics, 12(4):1223–1271.
  * Arcidiacono and Miller, (2011) Arcidiacono, P. and Miller, R. A. (2011).  Conditional choice probability estimation of dynamic discrete choice models with unobserved heterogeneity.  Econometrica, 79(6):1823–1867.
  * Barzilai and Borwein, (1988) Barzilai, J. and Borwein, J. M. (1988).  Two-point step size gradient methods.  IMA journal of numerical analysis, 8(1):141–148.
  * Berry et al., (1995) Berry, S., Levinsohn, J., and Pakes, A. (1995).  Automobile prices in market equilibrium.  Econometrica, 63(4):841–890.
  * Berry et al., (1999) Berry, S., Levinsohn, J., and Pakes, A. (1999).  Voluntary export restraints on automobiles: Evaluating a trade policy.  American Economic Review, 89(3):400–431.
  * Berry and Pakes, (2007) Berry, S. and Pakes, A. (2007).  The pure characteristics demand model.  International Economic Review, 48(4):1193–1225.
  * Berry, (1994) Berry, S. T. (1994).  Estimating discrete-choice models of product differentiation.  The RAND Journal of Economics, 25(2):242–262.
  * Bonnet et al., (2022) Bonnet, O., Galichon, A., Hsieh, Y.-W., O’hara, K., and Shum, M. (2022).  Yogurts choose consumers? estimation of random-utility models via two-sided matching.  The Review of Economic Studies, 89(6):3085–3114.
  * Conlon and Gortmaker, (2020) Conlon, C. and Gortmaker, J. (2020).  Best practices for differentiated products demand estimation with pyblp.  The RAND Journal of Economics, 51(4):1108–1161.
  * Conlon, (2012) Conlon, C. T. (2012).  A dynamic model of prices and margins in the lcd tv industry.  mimeo, Columbia University, 80:1433–1504.
  * Dai et al., (2015) Dai, Y.-H., Al-Baali, M., and Yang, X. (2015).  A positive barzilai–borwein-like stepsize and an extension for symmetric linear systems.  In Numerical Analysis and Optimization: NAO-III, Muscat, Oman, January 2014, pages 59–75. Springer.
  * Doi, (2022) Doi, N. (2022).  A simple method to estimate discrete-type random coefficients logit models.  International Journal of Industrial Organization, 81:102825.
  * Du and Varadhan, (2018) Du, Y. and Varadhan, R. (2018).  Squarem: An r package for off-the-shelf acceleration of em, mm and other em-like monotone algorithms.  arXiv preprint arXiv:1810.11163.
  * Dubé et al., (2012) Dubé, J.-P., Fox, J. T., and Su, C.-L. (2012).  Improving the numerical performance of static and dynamic aggregate discrete choice random coefficients demand estimation.  Econometrica, 80(5):2231–2267.
  * Fukasawa, (2024) Fukasawa, T. (2024).  When do firms sell high durability products? The case of Light Bulb Industry.  mimeo.
  * Goeree, (2008) Goeree, M. S. (2008).  Limited information and advertising in the us personal computer industry.  Econometrica, 76(5):1017–1074.
  * Gowrisankaran and Rysman, (2012) Gowrisankaran, G. and Rysman, M. (2012).  Dynamics of consumer demand for new durable goods.  Journal of Political Economy, 120(6):1173–1219.
  * Grigolon and Verboven, (2014) Grigolon, L. and Verboven, F. (2014).  Nested logit or random coefficients logit? A comparison of alternative discrete choice models of product differentiation.  Review of Economics and Statistics, 96(5):916–935.
  * Hendel and Nevo, (2006) Hendel, I. and Nevo, A. (2006).  Measuring the implications of sales and consumer inventory behavior.  Econometrica, 74(6):1637–1673.
  * Huang and Wan, (2017) Huang, S. and Wan, Z. (2017).  A new nonmonotone spectral residual method for nonsmooth nonlinear equations.  Journal of Computational and Applied Mathematics, 313:82–101.
  * Igami, (2017) Igami, M. (2017).  Estimating the innovator’s dilemma: Structural analysis of creative destruction in the hard disk drive industry, 1981–1998.  Journal of Political Economy, 125(3):798–847.
  * Iizuka, (2007) Iizuka, T. (2007).  Experts’ agency problems: evidence from the prescription drug market in japan.  The RAND journal of economics, 38(3):844–862.
  * Jamshidian and Jennrich, (1997) Jamshidian, M. and Jennrich, R. I. (1997).  Acceleration of the EM algorithm by using quasi-Newton methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 59(3):569–587.
  * Judd, (1998) Judd, K. L. (1998).  Numerical methods in economics.  MIT press.
  * Kalouptsidi, (2012) Kalouptsidi, M. (2012).  From market shares to consumer types: Duality in differentiated product demand estimation.  Journal of Applied Econometrics, 27(2):333–342.
  * Kalouptsidi et al., (2020) Kalouptsidi, M., Scott, P. T., and Souza-Rodrigues, E. (2020).  Linear iv regression estimators for structural dynamic discrete choice models.  Journal of Econometrics.
  * Kasahara and Shimotsu, (2009) Kasahara, H. and Shimotsu, K. (2009).  Nonparametric identification of finite mixture models of dynamic discrete choices.  Econometrica, 77(1):135–175.
  * Kasahara and Shimotsu, (2012) Kasahara, H. and Shimotsu, K. (2012).  Sequential estimation of structural models with a fixed point constraint.  Econometrica, 80(5):2303–2319.
  * La Cruz et al., (2006) La Cruz, W., Martínez, J., and Raydan, M. (2006).  Spectral residual method without gradient information for solving large-scale nonlinear systems of equations.  Mathematics of computation, 75(255):1429–1448.
  * Lee and Seo, (2015) Lee, J. and Seo, K. (2015).  A computationally fast estimator for random coefficients logit demand models using aggregate data.  The RAND Journal of Economics, 46(1):86–102.
  * Lee and Seo, (2016) Lee, J. and Seo, K. (2016).  Revisiting the nested fixed-point algorithm in blp random coefficients demand estimation.  Economics Letters, 149:67–70.
  * Miranda and Fackler, (2004) Miranda, M. J. and Fackler, P. L. (2004).  Applied computational economics and finance.  MIT press.
  * Nevo, (2001) Nevo, A. (2001).  Measuring market power in the ready-to-eat cereal industry.  Econometrica, 69(2):307–342.
  * Pakes and McGuire, (1994) Pakes, A. and McGuire, P. (1994).  Computing Markov-Perfect Nash Equilibria: Numerical Implications of a Dynamic Differentiated Product Model.  The RAND Journal of Economics, 25(4):555.
  * Pál and Sándor, (2023) Pál, L. and Sándor, Z. (2023).  Comparing procedures for estimating random coefficient logit demand models with a special focus on obtaining global optima.  International Journal of Industrial Organization, 88:102950.
  * Reynaerts et al., (2012) Reynaerts, J., Varadha, R., and Nash, J. C. (2012).  Enhencing the convergence properties of the blp (1995) contraction mapping.
  * Salanié and Wolak, (2022) Salanié, B. and Wolak, F. A. (2022).  Fast, detail-free, and approximately correct: Estimating mixed demand systems.  Technical report, Working paper.
  * Schiraldi, (2011) Schiraldi, P. (2011).  Automobile replacement: a dynamic structural approach.  The RAND journal of economics, 42(2):266–291.
  * Shcherbakov, (2016) Shcherbakov, O. (2016).  Measuring consumer switching costs in the television industry.  The RAND Journal of Economics, 47(2):366–393.
  * Sun and Ishihara, (2019) Sun, Y. and Ishihara, M. (2019).  A computationally efficient fixed point approach to dynamic structural demand estimation.  Journal of Econometrics, 208(2):563–584.
  * Varadhan and Gilbert, (2010) Varadhan, R. and Gilbert, P. (2010).  Bb: An r package for solving a large system of nonlinear equations and for optimizing a high-dimensional nonlinear objective function.  Journal of statistical software, 32:1–26.
  * Varadhan and Roland, (2008) Varadhan, R. and Roland, C. (2008).  Simple and globally convergent methods for accelerating the convergence of any em algorithm.  Scandinavian Journal of Statistics, 35(2):335–353.
  * Waziri et al., (2010) Waziri, M., Leong, W., Hassan, M., and Monsi, M. (2010).  Jacobian computation-free newton method for systems of non-linear equations.  Journal of numerical Mathematics and stochastics, 2(1):54–63.

##  Appendix A Additional results and discussions

###  A.1 Convergence properties of the mappings (Static BLP)

####  A.1.1 Convergence properties of $\Phi^{\delta,\gamma}$

The following proposition shows the convergence properties of the mapping
$\Phi^{\delta,\gamma}$.

######  Proposition 5.

(a). Under $\gamma\in[0,1]$, the following inequalities hold:

$\displaystyle\left\|\Phi^{\delta,\gamma}(\delta_{1})-\Phi^{\delta,\gamma}(\delta_{2})\right\|_{\infty}$
$\displaystyle\leq$ $\displaystyle
c_{\gamma}\left\|\delta_{1}-\delta_{2}\right\|_{\infty}$
$\displaystyle\left\|\Phi^{\delta,\gamma}(\delta_{1})-\Phi^{\delta,\gamma}(\delta_{2})\right\|_{\infty}$
$\displaystyle\leq$
$\displaystyle\widetilde{c_{\gamma}}\left(\left\|\max\\{\delta_{1},\delta_{2}\\}-\delta_{1}\right\|_{\infty}+\left\|\max\\{\delta_{1},\delta_{2}\\}-\delta_{2}\right\|_{\infty}\right)$
$\displaystyle\leq$ $\displaystyle
2\widetilde{c_{\gamma}}\cdot\left\|\delta_{1}-\delta_{2}\right\|_{\infty}$

where $c_{\gamma}\equiv\sup_{j\in\mathcal{J},\delta\in
B_{\delta}}\sum_{k\in\mathcal{J}}\left|\sum_{i}s_{ik}\left(\frac{w_{i}s_{ij}(\delta)}{s_{j}(\delta)}-\gamma\frac{w_{i}s_{i0}(\delta)}{s_{0}(\delta)}\right)\right|$
and

$\widetilde{c_{\gamma}}\equiv\begin{cases}\sup_{i\in\mathcal{I},J^{*}\subset\mathcal{J},\delta}s_{iJ^{*}}(\delta)=\sup_{i\in\mathcal{I},\delta}\left(1-s_{i0}(\delta)\right)\leq
1&\text{if}\ \gamma=0\\\
\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)\leq
1&\text{if}\ \gamma=1\\\
\max\left\\{\sup_{i\in\mathcal{I},J^{*}\subset\mathcal{J},\delta\in
B_{\delta}}s_{iJ^{*}}(\delta)-\gamma\inf_{i\in\mathcal{I},J^{*}\subset\mathcal{J},\delta\in
B_{\delta}}s_{iJ^{*}}(\delta),\right.\\\ \left.\ \ \
\gamma\sup_{i\in\mathcal{I},J^{*}\subset\mathcal{J},\delta\in
B_{\delta}}s_{iJ^{*}}(\delta)-\inf_{i\in\mathcal{I},J^{*}\subset\mathcal{J},\delta\in
B_{\delta}}s_{iJ^{*}}(\delta)\right\\}&\text{Otherwise}\end{cases}$

(b). Under $\gamma=0$,

$\left\|\Phi^{\delta,\gamma=0}(\delta_{1})-\Phi^{\delta,\gamma=0}(\delta_{2})\right\|_{\infty}\leq\sup_{i\in\mathcal{I},\delta\in
B_{\delta}}\left(1-s_{i0}(\delta)\right)\cdot\left\|\delta_{1}-\delta_{2}\right\|_{\infty}\leq\left\|\delta_{1}-\delta_{2}\right\|_{\infty}$
holds.

(c). If $|\mathcal{J}|=1,$

$\displaystyle\left\|\Phi^{\delta,\gamma}(\delta_{1})-\Phi^{\delta,\gamma}(\delta_{2})\right\|_{\infty}$
$\displaystyle\leq$
$\displaystyle\widetilde{c_{\gamma}}\cdot\left\|\delta_{1}-\delta_{2}\right\|_{\infty}\leq\left\|\delta_{1}-\delta_{2}\right\|_{\infty}$

First, we consider the case of $|\mathcal{J}|=1$ (only one product). The
proposition implies $\Phi^{\delta,\gamma=0}$ (BLP contraction mapping) has a
contraction property353535See also Remark 12. with modulus
$\sup_{i\in\mathcal{I},J^{*}\subset\mathcal{J}}s_{iJ^{*}}(\delta)=1-\inf_{i\in\mathcal{I}}s_{i0}(\delta)$.
In contrast, $\Phi^{\delta,\gamma=1}$ has a contraction property with modulus
$\sup_{i\in\mathcal{I},J^{*}\subset\mathcal{J}}s_{iJ^{*}}(\delta)-\inf_{i\in\mathcal{I},J^{*}\subset\mathcal{J}}s_{iJ^{*}}(\delta)<\sup_{i\in\mathcal{I},J^{*}\subset\mathcal{J}}s_{iJ^{*}}(\delta)$.
Hence, we can expect that the latter attains faster convergence.

When $|\mathcal{J}|>1$ holds, there is no guarantee that $c_{\gamma}$ is less
than 1. Hence, there is no guarantee that $\Phi^{\delta,\gamma=1}$ is a
contraction. Nevertheless, $c_{\gamma=1}$ takes a small value when the
consumer heterogeneity, measured by
$\frac{s_{ij}}{s_{j}}-\frac{s_{i0}}{s_{0}}$, is small. If there is no consumer
heterogeneity, $\frac{s_{ij}}{s_{j}}-\frac{s_{i0}}{s_{0}}=0$ holds, and the
modulus is equal to 0. It implies the iteration using the mapping immediately
converges to the solution. The proposition also implies that the magnitude of
consumer heterogeneity is important for the convergence speed of
$\Phi^{\delta,\gamma=1}$.

Besides, the inequality
$\left\|\Phi^{\delta,\gamma}(\delta_{1})-\Phi^{\delta,\gamma}(\delta_{2})\right\|_{\infty}\leq\widetilde{c_{\gamma}}\left(\left\|\max\\{\delta_{1},\delta_{2}\\}-\delta_{1}\right\|_{\infty}+\left\|\max\\{\delta_{1},\delta_{2}\\}-\delta_{2}\right\|_{\infty}\right)$
for $|\mathcal{J}|\in\mathbb{N}$ implies

$\left\|\Phi^{\delta,\gamma}(\delta_{1})-\Phi^{\delta,\gamma}(\delta_{2})\right\|_{\infty}\leq\widetilde{c_{\gamma}}\left\|\delta_{1}-\delta_{2}\right\|_{\infty}=\left(\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)\right)\left\|\delta_{1}-\delta_{2}\right\|_{\infty}\leq\left\|\delta_{1}-\delta_{2}\right\|_{\infty}$
under $\delta_{1}\geq\delta_{2}$. Hence, $\Phi^{V,\gamma=1}$ has a property
similar to contractions.

######  Remark 12.

If we assume $B_{\delta}=\mathbb{R}^{|\delta|}$, the value of
$\sup_{i\in\mathcal{I},\delta\in B_{\delta}}\left(1-s_{i0}(\delta)\right)$ can
be 1, which implies $\Phi^{\delta,\gamma=0}$ is not a contraction, strictly
speaking. Hence, we should restrict the domain or the range of the mapping
$\Phi^{\delta,\gamma=0}$ to guarantee that the mapping is a contraction, as
discussed in the Appendix of Berry et al., (1995).

####  A.1.2 Convergence properties of $\Phi^{\delta,\gamma}$

We define the following terms:+

$\displaystyle prob_{0|i_{2}}(V)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{1+\sum_{j\in\mathcal{J}}S_{j}^{(data)}\frac{\exp(\mu_{i_{2}j})}{\sum_{i\in\mathcal{I}}w_{i}\exp(\mu_{ij})\exp(-V_{i})}},$
$\displaystyle prob_{I^{*}|j}(V)$ $\displaystyle\equiv$
$\displaystyle\frac{\sum_{i\in
I^{*}\subset\mathcal{I}}w_{i}\exp(\mu_{ij})\exp(-V_{i})}{\sum_{i\in\mathcal{I}}w_{i}\exp(\mu_{ij})\exp(-V_{i})},$
$\displaystyle prob_{I^{*}|0}(V)$ $\displaystyle\equiv$
$\displaystyle\frac{\sum_{i\in
I^{*}\subset\mathcal{I}}w_{i}\exp(-V_{i})}{\sum_{i\in\mathcal{I}}w_{i}\exp(-V_{i})}.$

Intuitively, $prob_{0|i}$ denotes the probability that consumer $i$ chooses
alternative $j$. $prob_{I^{*}\subset\mathcal{I}|j}$ denotes the fraction of
consumers in group $I^{*}\subset\mathcal{I}$ choosing alternative $j$ among
those choosing $j$.

Then, we obtain the following results regarding the convergence properties of
$\Phi^{V,\gamma}$.

######  Proposition 6.

(a). Under $\gamma\in[0,1]$, the following inequalities hold:

$\displaystyle\left\|\Phi^{V,\gamma}(V_{1})-\Phi^{V,\gamma}(V_{2})\right\|_{\infty}$
$\displaystyle\leq$ $\displaystyle 2\sup_{i_{2}\in\mathcal{I},V\in
B_{V}}\left((1-prob_{0|i_{2}}(V))\cdot\left(\frac{\sum_{i\in\mathcal{I}}w_{i}\exp(-V_{i})}{S_{0}^{(data)}}\right)^{\gamma}\right)\cdot\left\|V_{1}-V_{2}\right\|_{\infty},$
$\displaystyle\left\|\Phi^{V,\gamma}(V_{1})-\Phi^{V,\gamma}(V_{2})\right\|_{\infty}$
$\displaystyle\leq$ $\displaystyle\sup_{i\in\mathcal{I},V\in
B_{V}}\left((1-prob_{0|i}(V))\cdot\left(\frac{\sum_{i\in\mathcal{I}}w_{i}\exp(-V_{i})}{S_{0}^{(data)}}\right)^{\gamma}\right)\cdot$
$\displaystyle\sup_{I^{*}\subset\mathcal{I},j\in\mathcal{J},V\in
B_{V}}\left|prob_{I^{*}|j}(V)-\gamma prob_{I^{*}|0}(V)\right|\cdot$
$\displaystyle\left(\left\|\max\\{V_{1},V_{2}\\}-V_{1}\right\|_{\infty}+\left\|\max\\{V_{1},V_{2}\\}-V_{2}\right\|_{\infty}\right).$

(b). Under $\gamma=0$,
$\left\|\Phi^{V,\gamma}(V_{1})-\Phi^{V,\gamma}(V_{2})\right\|_{\infty}\leq\sup_{V\in
B_{V},i_{2}\in\mathcal{I}}\left((1-prob_{0|i_{2}}(V))\right)\cdot\left\|V_{1}-V_{2}\right\|_{\infty}\leq\left\|V_{1}-V_{2}\right\|_{\infty}$
holds.

The result in (b) corresponds to the contraction property of the BLP
contraction mapping $\Phi^{\delta,\gamma=0}$. Regarding the results in (a),
the first inequality implies $\Phi^{V,\gamma}$ is not always a contraction
mapping if $2\sup_{i_{2}\in\mathcal{I},V\in
B_{V}}\left((1-prob_{0|i_{2}}(V))\cdot\left(\frac{\sum_{i\in\mathcal{I}}w_{i}\exp(-V_{i})}{S_{0}^{(data)}}\right)^{\gamma}\right)>1$.
Nevertheless, the second inequality shows it has a property similar to
contraction mappings and has good convergence properties. Here, consider the
setting where
$\frac{\sum_{i\in\mathcal{I}}w_{i}\exp(-V_{i})}{S_{0}^{(data)}}\approx 1$,
$V^{*}>V$, and $\gamma=1$. Then,
$\left\|\Phi^{V,\gamma}(V)-\Phi^{V,\gamma}(V^{*})\right\|_{\infty}\leq\sup_{i\in\mathcal{I},B\in
B_{V}}\left((1-prob_{0|i}(V))\right)\cdot\sup_{I^{*}\subset\mathcal{I},j\in\mathcal{J},V\in
B_{V}}\left|prob_{I^{*}|j}(V)-prob_{I^{*}|0}(V)\right|\cdot\left\|V^{*}-V\right\|_{\infty}$
holds. $\sup_{I^{*}\subset\mathcal{I},j\in\mathcal{J},V\in
B_{V}}\left|prob_{I^{*}|j}(V)-prob_{I^{*}|0}(V)\right|$ is the maximum
absolute value of the difference between $prob_{I^{*}|0}(V)$ and
$prob_{I^{*}|j}(V)$, and it is less than 1. Hence,
$\left\|\Phi^{V,\gamma}(V)-\Phi^{V,\gamma}(V^{*})\right\|_{\infty}\leq\exists
K\left\|V^{*}-V\right\|_{\infty}\leq\left\|V^{*}-V\right\|_{\infty}\ (K\leq
1)$ holds under the setting.

Besides, $prob_{I^{*}|j}(V)-prob_{I^{*}|0}(V)$ take small values when the
consumer heterogeneity is small, implying we can attain high convergence speed
in the model with small consumer heterogeneity. If the size of consumer
heterogeneity is zero as in logit models without random coefficients, the
mapping is equivalent to
$\Phi^{V,\gamma=1}(V)=\log\left(1+\left(\sum_{j\in\mathcal{J}}S_{j}^{(data)}\right)\cdot\left(\frac{1}{S_{0}^{(data)}}\right)\right)=-\log\left(S_{0}^{(data)}\right)$,
whose right hand side is the solution of the fixed point problem, and it
immediately converges to the solution after applying the mapping
$\Phi^{V,\gamma=1}$ once.

###  A.2 The case with large consumer heterogeneity

As discussed in Appendix A.1, the size of consumer heterogeneity is critical
for the performance of the new mappings
$\Phi^{\delta,\gamma=1},\Phi^{V,\gamma=1}$. This section shows results under
large consumer heterogeneity. Unlike the standard setting with relatively
small consumer heterogeneity experimented in Dubé et al., (2012), Lee and Seo,
(2015) and others, applying the mapping
$\Phi^{\delta,\gamma=1},\Phi^{V,\gamma=1}$ did not lead to convergence.
Nevertheless, the spectral algorithm works well even under the setting.

#### Settings

Consumer $i$’s utility when choosing product $j$ is
$U_{ij}=\delta_{j}+\mu_{ij}+\epsilon_{ij}$, and utility when not buying
anything is $U_{i0}=\epsilon_{i0}$. Suppose $\epsilon$ follows Gumbel
distribution. Let $|\mathcal{J}|=2$ (2 products), $|\mathcal{I}|=2$ (2
consumer types), $\delta_{j=1}=0$, $\delta_{j=2}=-1,$ and
$\mu_{i=1,j=1}=\mu_{i=2,j=2}=10,$ $\mu_{i=1,j=2}=\mu_{i=2,j=1}=0$,
$w_{i=1}=0.1,w_{i=2}=0.9$.

Table 7 shows the choice probabilities for each consumer type. As the table
shows, there is large consumer heterogeneity: Type 1 consumers strongly prefer
product 1, and type 2 consumers strongly prefer product 2.

Table 7: Choice probabilities by consumer types | $i=1$ | $i=2$
---|---|---
$j=1$ | 0.9999 | 0.0001
$j=2$ | 0.0000 | 0.9998
$j=0$ | 0.0000 | 0.0001

Unlike the setting in Section 6, I use the true $\mu_{ij}(i=1,2;j=1,2)$ to
solve for $\delta$, to make the setting replicable.

#### Results

Table 8 shows the results. When we apply the algorithm $\delta$-(1) (mapping
$\Phi^{\delta,\gamma=1}$), the algorithm does not converge. Regarding
$\delta$-(1), Figure 3 shows the trend of the norm
$\left\|\delta_{j}^{(n+1)}-\delta_{j}^{(n)}\right\|_{\infty}$. If the mapping
$\Phi^{\delta,\gamma=1}$ is a contraction, the value of
$\left\|\delta_{j}^{(n+1)}-\delta_{j}^{(n)}\right\|_{\infty}$ should decrease
as the number of iterations increases, and
$\left\|\delta_{j}^{(n+2)}-\delta_{j}^{(n+1)}\right\|_{\infty}-\left\|\delta_{j}^{(n+1)}-\delta_{j}^{(n)}\right\|_{\infty}<0$
should hold. Nevertheless, they do not hold in the figure. This implies that
$\delta$-(1) is not a contraction in the current setting.

Nevertheless, even when using $\delta$-(1), combining the spectral algorithm
works well, as shown in the number of function evaluations.

Table 8: Results of Monte Carlo simulation (Static BLP model; Too large
consumer heterogeneity)

| $J$ | Func. Evals. | Mean | Conv. | $\log_{10}\left(DIST\right)$ | $DIST<\epsilon_{tol}$
---|---|---|---|---|---|---
| Mean | Min. | 25th | Median. | 75th | Max. | CPU time (s) | (%) | (%)
$\delta$-(0) (BLP) | 25 | 2000 | 2000 | 2000 | 2000 | 2000 | 2000 | 0.099 | 0 | -3.8 | 0
$\delta$-(0) (BLP) + Spectral | 25 | 41 | 41 | 41 | 41 | 41 | 41 | 0.005 | 100 | -13.7 | 100
$\delta$-(0) (BLP) + SQUAREM | 25 | 119 | 119 | 119 | 119 | 119 | 119 | 0.006 | 100 | -13.1 | 100
$\delta$-(1) | 25 | 2000 | 2000 | 2000 | 2000 | 2000 | 2000 | 0.073 | 0 | -4.2 | 0
$\delta$-(1) + Spectral | 25 | 98 | 98 | 98 | 98 | 98 | 98 | 0.005 | 100 | -14.9 | 100
$\delta$-(1) + SQUAREM | 25 | 35 | 35 | 35 | 35 | 35 | 35 | 0.001 | 100 | -14.2 | 100
$V$-(0) | 25 | 2000 | 2000 | 2000 | 2000 | 2000 | 2000 | 0.103 | 0 | -3.4 | 0
$V$-(0) + Spectral | 25 | 36 | 36 | 36 | 36 | 36 | 36 | 0.004 | 100 | -14.9 | 100
$V$-(0) + SQUAREM | 25 | 1013 | 1013 | 1013 | 1013 | 1013 | 1013 | 0.052 | 100 | -13.1 | 100
$V$-(1) | 25 | 2000 | 2000 | 2000 | 2000 | 2000 | 2000 | 0.095 | 0 | -4.9 | 0
$V$-(1) + Spectral | 25 | 36 | 36 | 36 | 36 | 36 | 36 | 0.003 | 100 | -15.5 | 100
$V$-(1) + SQUAREM | 25 | 58 | 58 | 58 | 58 | 58 | 58 | 0.003 | 100 | -14 | 100

Notes.

$DIST\equiv\left\|\log(S_{j}^{(data)})-\log(s_{j})\right\|_{\infty},$$\epsilon_{tol}=$1E-12.

Maximum number of function evaluations is set to 2000.

Figure 3: Example of non-contraction mapping $\delta$-(1)

Note.

$\Delta\log_{10}\left(\left\|\delta_{j}^{(n+1)}-\delta_{j}^{(n)}\right\|_{\infty}\right)$
is defined by
$\log_{10}\left(\left\|\delta_{j}^{(n+2)}-\delta_{j}^{(n+1)}\right\|_{\infty}\right)-\log_{10}\left(\left\|\delta_{j}^{(n+1)}-\delta_{j}^{(n)}\right\|_{\infty}\right)$.

###  A.3 Choice of step sizes in the spectral / SQUAREM algorithm

Here, I show the results of static BLP Monte Carlo simulations under
alternative choices of step sizes ($\alpha_{S1},\alpha_{S2}$). In addition to
these step sizes, I also experiment the step size $\alpha_{S3}^{\prime}\equiv
sgn\left(s^{(n)\prime}y^{(n)}\right)\frac{\left\|s^{(n)}\right\|_{2}}{\left\|y^{(n)}\right\|_{2}}$,
which introduced in the BB package of R language. Here,
$sgn\left(s^{(n)\prime}y^{(n)}\right)$ denotes the sign of
$s^{(n)\prime}y^{(n)}$. Results are shown in Tables 9, 10, and 11. Unlike the
case of $\alpha_{S3}$, the algorithms sometimes do not converge, and lead to
slower convergence under these specifications. The results imply choosing the
step size $\alpha_{S3}$ is better from the viewpoint of stable and fast
convergence.

Note that using the step sizes $\alpha_{S1},\alpha_{S2}$ do not always lead to
much worse performance. Table 12 shows the performance of each algorithm under
different choice of the step size, using Nevo, (2001)’s dataset. In the case
of the estimation using Nevo, (2001)’s dataset, the total number of function
evaluation are mostly the same for all the step sizes.

In Table 13, I show the list of built-in packages of the spectral and SQUAREM
algorithms. The step size $\alpha_{S3}$ is not available in all the built-in
packages of the spectral / SQUAREM algorithms. Though using $\alpha_{S3}$ is
desirable, they are not available in all the packages. Using built-in packages
with step size $\alpha_{S1},\alpha_{S2}$ might be sufficient for practitioners
who would like to use built-in packages rather than delving into coding
issues.

Table 9: Results of Monte Carlo simulation (Static BLP model; Continuous
consumer types; $\alpha_{S1}$)

| $J$ | Func. Evals. | Mean | Conv. | Mean | $DIST<\epsilon_{tol}$
---|---|---|---|---|---|---
| Mean | Min. | 25th | Median. | 75th | Max. | CPU time (s) | (%) | $\log_{10}\left(DIST\right)$ | (%)
$\delta$-(0) (BLP) + Spectral | 250 | 135.1 | 14 | 22 | 29.5 | 40 | 1000 | 0.44626 | 90 | -13.6 | 92
$\delta$-(0) (BLP) + SQUAREM | 250 | 223.3 | 15 | 23 | 32 | 41 | 1000 | 0.70986 | 80 | -12.6 | 82
$\delta$-(1) + Spectral | 250 | 14.4 | 8 | 12 | 14 | 16 | 24 | 0.04972 | 100 | -15.2 | 100
$\delta$-(1) + SQUAREM | 250 | 14.64 | 9 | 12 | 15 | 16 | 23 | 0.04886 | 100 | -15.4 | 100
$V$-(0) + Spectral | 250 | 30.14 | 13 | 21 | 26.5 | 38 | 64 | 0.11384 | 100 | -14.8 | 100
$V$-(0) + SQUAREM | 250 | 33.5 | 15 | 22 | 28 | 36 | 158 | 0.13004 | 100 | -14.7 | 100
$V$-(1) + Spectral | 250 | 17.36 | 10 | 13 | 16 | 21 | 31 | 0.06754 | 100 | -15.1 | 100
$V$-(1) + SQUAREM | 250 | 36.46 | 11 | 14 | 16 | 18 | 1000 | 0.1545 | 98 | -14.8 | 98

Notes. $DIST\equiv\left\|\log(S_{j}^{(data)})-\log(s_{j})\right\|_{\infty}$,
$\epsilon_{tol}=$1E-12.

Number of simulation draws is set to 1000.

Maximum number of function evaluations is set to 1000.

Mean Outside option share is 0.307.