Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
101,635 characters · 26 sections · 83 citation commands
Model-Adaptive Approach to Dynamic Discrete Choice Models with Large State Spaces
\setcounter{page}{2}
\numberwithin{equation}{section}
\onehalfspacing
Dynamic Discrete Choice (DDC) models with large state spaces have become increasingly popular due to their ability to capture decision-making processes in complex high-dimensional settings. However, estimation and counterfactual experiments in these settings pose significant computational challenges. This paper proposes a model-adaptive approach, based on the conjugate gradient (CG) method, to solve the linear system of fixed point equations of the policy valuation operator. Our goal is to provide a fast and easily implementable method to solve the equations within a pre-specified tolerance. As the policy valuation step is fundamental to Conditional Choice Probability (CCP) estimators, full-solution estimators and counterfactual experiments with policy iteration or Newton--Kantorovich iterations, our approach offers a useful numerical tool across various empirical applications, expanding the set of complex high-dimensional settings in which DDC models can be used.
The policy valuation operator, as described in aguirregabiria2002swapping, involves solving for the value function implied by an arbitrary policy function, which may not necessarily be optimal. This value function represents the expected discounted utility if an individual behaves according to that policy function. CCP estimators (e.g., hotz1993conditional, aguirregabiria2002swapping,aguirregabiria2007sequential, pesendorfer2008asymptotic, and arcidiacono2011conditional) use the policy valuation operator to solve for the value function given a consistent estimator of CCPs. If policy iteration or Newton--Kantorovich iterations are used for full solution estimators or counterfactual experiments, each iteration employs the policy valuation operator. However, in models with large state spaces, policy valuation remains computationally demanding and requires efficient numerical methods. The accuracy of these methods is crucial for obtaining reliable estimates. dube2012improving shows that loose tolerance thresholds can lead to bias in parameter estimates. Therefore, there is a clear need for fast and accurate numerical methods for DDC models with large state spaces.
This paper proposes a model-adaptive (MA) approach to solve the linear system of fixed point equations of the policy valuation operator. The primary goal is to achieve a pre-specified tolerance while significantly reducing computational costs, enabling the use of policy valuation tools in a wide range of empirical applications. We call our approach model-adaptive as it designs the sieve space based on the model primitives (such as the transition density and utility function), and the algorithm itself selects the sieve dimension. At each step, our approach augments the sieve space with the residual from the previous iteration and projects the value function onto the augmented sieve space. It achieves a faster decay rate of the approximation error than conventional methods, such as Successive Approximation (SA), which iterates the contraction mapping, and Temporal Difference (TD), which projects the value function onto a pre-specified sieve space. Formally, we show that the approximation error decays at a superlinear rate in the sieve dimension (number of iterations), while SA achieves a linear rate and TD achieves only a sublinear rate. Furthermore, the sieve space and its dimension are automatically constructed by the algorithm, eliminating the need for researchers to design the sieve space and choose its dimension. Consequently, our method is easy to implement and converges faster than conventional methods, offering substantial computational savings.
The main computational cost depends on the number of iterations and matrix-vector multiplication operations. Our approach attains superlinear convergence on the approximation error, which substantially reduces the number of iterations required to reach a desired level of tolerance. Formally, our approach achieves an approximation error upper bound of $O((\frac{C_{1}}{\sqrt{k}})^{k})$, where $k$ is the number of iterations and $C_{1}$ is a constant. Furthermore, the bound can be improved to $O((\frac{C_{2}}{k})^{k})$ if the transition density has continuous partial derivatives. Notably, only the constants in the upper bound depend on the discount factor and properties of the transition density (such as its boundedness and smoothness). Unlike SA, whose convergence rate is governed by the contraction modulus $\beta$ and slows significantly as $\beta \to 1$, our method's superlinear convergence does not rely on the contraction property. Therefore, our method is well-suited for models with large state spaces and large discount factors, such as dynamic consumer demand models (e.g., hendel2006measuring and wang2015impact). Moreover, it works particularly well for models with smooth transition densities, such as autoregressive (AR) processes, which are commonly used in empirical applications (e.g., sweeting2013dynamic, huang2014dynamic, kalouptsidi2014time, grieco2022input, and gerarden2023demanding).
We provide implementations for both discrete and continuous state spaces. For discrete state spaces, the implementation requires only matrix operations as the system of equations takes the form of a finite-dimensional linear system. For continuous state spaces, we employ numerical integration to approximate the integral of the policy valuation step. Thus, there is a trade-off between simulation error and computational cost: increasing the number of grid points reduces simulation error at the cost of solving a larger linear system within a given tolerance. The computational cost of matrix-vector multiplication operations increases with the number of grid points. Therefore, we analyze the impact of the number of grid points on the number of iterations required to achieve a desired tolerance. We show that the number of iterations required for convergence remains approximately the same for all sufficiently large numbers of grid points. As a result, we can expect the number of iterations to be independent of the number of grid points for numerical integration. Therefore, our method allows researchers to reduce simulation error by increasing the number of grid points up to the computational limits of matrix-vector multiplication. Fast matrix-vector multiplication algorithms can further accelerate the computation (e.g., rokhlin1985rapid, greengard1987fast, and hackbusch1989fast). In addition, matrix-vector multiplication is amenable to GPU acceleration, which can further reduce the computational cost.
We illustrate the performance of our approach using three numerical experiments. We first simulate the bus engine replacement problem to visualize the convergence behavior of our method. The plot of our approximation solution shows that the method uses a few iterations to find a good sieve space. After that, it converges rapidly to the true solution.
Second, we analyze a model for dynamic consumer demand for storable goods similar to hendel2006measuring. We compare policy iteration and Newton--Kantorovich outer iterations combined with different inner solvers (MA, SA, and an exact solver), along with VFI benchmarks. We demonstrate that PI+MA is 80% faster than conventional methods such as SA and an exact equation solver. It is also more than 90 times faster than one-step value function iteration methods. These substantial computational savings open the door to the use of Bayesian MCMC estimators (chernozhukov2003mcmc), which are well-suited as the likelihood function is not differentiable in utility parameters.
Third, we examine a dynamic firm entry and exit problem in aguirregabiria2023solution. We solve the dynamic programming problem by policy iteration using MA to solve the linear system of fixed point equations. We vary the discount factor and number of grid points for numerical integration to evaluate the performance. The computational times confirm that our method improves the computational efficiency of policy iteration. The results show that the numbers of iterations required for convergence are approximately the same regardless of the numbers of grid points. Moreover, the number of iterations only slightly increases as the discount factor approaches one. Finally, we compare our method with TD and SA. The simulation results show that our method outperforms SA in terms of computational time and TD in terms of approximation error.
We apply our method to a dynamic consumer demand model for laundry detergent using Kantar's Worldpanel Take Home data. For each household size, we separately estimate the dynamic parameters using the Bayesian MCMC estimator. At each MCMC step, we solve the dynamic programming problem by policy iteration with MA. The results confirm the computational efficiency of our method in practice. We also simulate the long-run elasticities, which reveal the heterogeneous substitution patterns across different household sizes.
Sieve approximation methods have been used to solve dynamic programming problems (e.g., norets2012estimation, arcidiacono2013approximating, and wang2015impact). They have been widely used to solve the linear equations in the policy valuation step. Applications include hendel2006measuring, sweeting2013dynamic, and bodere2023dynamic. Recent work approximates the solution to the linear equation by temporal difference (see adusumilli2019temporal). However, those methods require researchers both to design the space of basis functions (e.g., polynomials, splines, or neural networks) and choose the sieve dimension (e.g., the degree of polynomials, the number of knots, and the number of hidden layers). The best choice for each application is almost always unclear. In contrast, MA constructs a model-adaptive sieve space using the algorithm itself to design that space. The sieve dimension (i.e., the number of iterations) is also determined by the algorithm. The method is guaranteed to achieve the pre-specified tolerance. Finally, we show that the approximation error decays at a superlinear rate in the sieve dimension. Other methods like TD achieve only a sublinear rate.
Successive Approximation (see kress_linear_2014), also known as fixed point iteration (judd1998numerical), is an iterative method to solve the fixed point equations of the policy valuation operator. The computational cost depends on the number of iterations and the number of matrix-vector multiplication operations. The convergence of SA relies on the $\beta$-contraction property of the transition operator $\mathcal{T}$, where $\beta$ is the discount factor. Therefore, the number of iterations of SA increases significantly as the discount factor approaches one, making the methods computationally demanding. In contrast, the superlinear convergence of our method does not rely on the contraction property, making it particularly well-suited for models with large discount factors such as consumer demand models (e.g., hendel2006measuring and wang2015impact). Moreover, our method can outperform SA even for small discount factors as SA achieves linear convergence while our method achieves superlinear convergence.
For continuous state spaces, we employ numerical integration to approximate the integral of the policy valuation operator, which implicitly discretizes the state space. Discretization is widely used in economics (e.g., sweeting2013dynamic, kalouptsidi2014time, huang2015structural, and bodere2023dynamic). rust1997using,rust1997comparison study the simulation error from numerical integration and assume the discretized equation can be solved exactly. However, there is a trade-off between simulation error and computational cost of solving the discretized equation. Increasing the number of grid points reduces the simulation error while increasing the size of the linear system to be solved; potentially making it computationally infeasible to solve the system exactly. Sieve approximation is used to solve the discretized equation (e.g., hendel2006measuring, sweeting2013dynamic, and bodere2023dynamic). Instead, we propose to solve the discretized equation within a given tolerance while minimizing computational cost. The computational cost of MA depends on the number of iterations and the number of matrix-vector multiplication operations. We can expect that the number of iterations is small and independent of the number of grid points. Therefore, MA enables researchers to reduce simulation error up to the computational limits of matrix-vector multiplication.
Outline: The remainder of the paper is organized as follows. (ref) reviews DDC models and the policy valuation operator. (ref) presents the model-adaptive approach, its implementation and computational cost. (ref) describes the theoretical properties. (ref) reports results from three numerical experiments. (ref) applies our method to a consumer demand for storable goods. (ref) concludes. The proofs are in Appendix (ref). Appendix (ref) contains details of algorithms of simulations and empirical application.
Notation: Let $\mathbb{X}$ be the support of $x$, and $\mathbb{A} := \{0,1,\cdots, A-1\}$. $\mathbb{X}$ can be discrete or continuous. For a probability distribution $\mu$ on $\mathbb{X}$, which in the continuous case is assumed to be absolutely continuous with respect to Lebesgue measure, let $L^{2}(\mathbb{X}, \mu)$ denote the space of square-integrable functions on $\mathbb{X}$. Let $\langle \cdot,\cdot \rangle_{\mu}$ and $\|\cdot\|_{\mu}$ denote the inner product and norm induced by $\mu$. Let $\nu_{Leb}$ denote the Lebesgue measure, and let $\langle \cdot,\cdot \rangle$, $\| \cdot \|$, and $L^{2}(\mathbb{X})$ denote the inner product, norm, and $L^{2}$-space, respectively. We suppress $\nu_{Leb}$ for notational simplicity. Let $\|\cdot\|_{\infty}$ denote the sup-norm of a vector. For a linear operator, $\mathcal{A}: L^{2}(\mathbb{X}) \longmapsto L^{2}(\mathbb{X})$, denote by $\|\mathcal{A}\|_{op} := \sup_{h \in L^{2}(\mathbb{X}), \|h\|=1} \|\mathcal{A} h\|$ the operator norm. Let $\mathcal{I}$ be the identity operator.
We study infinite horizon stationary dynamic discrete choice models as in rust1994structural. In each discrete period $t = 0, 1, 2, \ldots$, an individual chooses $a_{t} \in \mathbb{A}$ to maximize her discounted expected utility:
where $\beta \in (0,1)$ is the discount factor, $u(x_{t},a_{t})$ is the period utility, $x_{t} \in \mathbb{X}$ is the observable state (to researchers) that follows a first-order Markov process with a transition density $f(x_{t+1} | x_{t}, a_{t})$, $\varepsilon_{t} \in \mathbb{R}^A$ is a vector of unobservable i.i.d. type I extreme value shocks with Lebesgue density $g(\varepsilon_{t})$, and $\varepsilon_t(a_t)$ is the element of $\varepsilon_t$ corresponding to $a_t$.
Under regularity conditions (see rust1994structural), the utility maximization problem has a solution and the optimal value function $V_{opt}(x,\varepsilon)$ is the unique solution to the Bellman equation:
where $(x',\varepsilon')$ denotes the next period's state and utility shock. Under the i.i.d. assumption on utility shocks, integrating the utility shocks out, the integrated Bellman equation has the following form:
where $V_{opt}(x)$ is the integrated value function and $v^{*}(x,a)$ is the conditional value function defined as:
The conditional choice probability (CCP) is the probability that action $a$ is optimal conditional on observable state $x$ defined by:
where $\mathbbm{1}\{\cdot\}$ is the indicator function. Under the distributional assumption on utility shocks, CCPs take the following form:
The following map is derived in arcidiacono2011conditional as a corollary of Hotz--Miller Inversion Lemma by hotz1993conditional, $\forall \ (x,a) \in \mathbb{X} \times \mathbb{A}$:
where $\kappa$ is the Euler constant. The equation (ref) establishes a crucial link between the integrated value function $V_{opt}(x)$, conditional value function $v^{*}(x,a)$, and conditional choice probabilities $p^{*}(a|x)$. It provides a powerful tool for the policy valuation step, which is essential for both policy iteration and CCP estimators.
The policy valuation operator (see aguirregabiria2002swapping) maps an arbitrary policy function to the value function using (ref). The value function represents the expected discounted utility of an individual if she behaves today and in the future according to that policy, which is not necessarily optimal. aguirregabiria2002swapping shows that the value function is obtained by solving the following equation for $V$ given a policy function $p$:
In this paper, we focus on solving (ref) for a given policy function. Therefore, we introduce the following notations and suppress the dependence on $p$:
Using this notation, we can rewrite (ref) as a linear system of fixed point equations:
Both CCP estimators and the policy iteration method (see howard1960dynamic) require solving equation (ref) for $V$. For CCP estimators, $p$ is replaced with its consistent estimator $\hat{p}$. For policy iteration, at iteration $i$, we solve for $V_{i}$ associated with $p_{i}$ from the previous iteration. Subsequently, the policy improvement step updates the policy function as follows:
where $v_{i}(x,a) = u(x,a) + \beta \mathbb{E}_{x'|x,a} V_{i}(x')$. This process iterates until convergence of the policy function is achieved.
An alternative full-solution method is the Newton--Kantorovich (NK) iteration applied to the integrated Bellman equation, as considered in rust1987optimal. Define the Bellman operator $\Gamma: \mathbb{R}^{|\mathbb{X}|} \to \mathbb{R}^{|\mathbb{X}|}$ by:
The fixed point of $\Gamma$ is the integrated value function $V_{opt}$. The Jacobian of $\Gamma$ at $V_{k}$ is the matrix $\mathcal{T}_{k}$ with entries:
where $p_{k}(a|x) := \frac{\exp(v_{k}(x,a))}{\sum_{a'}\exp(v_{k}(x,a'))}$ is the CCP induced by $V_{k}$. The NK step is:
Note that $\mathcal{T}_{k}$ has the same structure as $\mathcal{T}$ in (ref). In addition, the NK step coincides with the policy iteration step under the i.i.d. extreme value assumption. Therefore, (ref) is a linear system of the form $(\mathcal{I} - \mathcal{T})V = u$, where the operator $\mathcal{T}_{k}$ and right-hand side $\Gamma(V_{k}) - V_{k}$ change at each NK iteration. The MA approach developed in the next section applies directly to solving (ref). While NK iterations enjoy quadratic convergence to the fixed point of $\Gamma$ (see rust1988maximum), each NK step requires solving the linear system (ref), which becomes computationally demanding for large state spaces. Our method provides an efficient solver for this inner linear system.
This section presents the model-adaptive (MA) approach, its implementation and computational cost. Our method employs the Conjugate Gradient (CG) method, an iterative approach for solving large linear systems. The CG method is commonly attributed to hestenes1952methods. For comprehensive textbooks, see kelley1995iterative and han2009theoretical.
The CG was originally designed for solving linear systems with self-adjoint operators. However, the operator $\mathcal{T}$ is not necessarily self-adjoint. If $\mathcal{T}$ were self-adjoint with respect to the inner product space $\langle \cdot, \cdot \rangle_{\mu}$, then the Markov chain would be time-reversible, which is a strong assumption in many practical settings. Therefore, instead of solving (ref) for $V$ directly, we propose to solve the following equation for $y$:
and set $V = (\mathcal{I} - \mathcal{T}^{*})y$ to solve (ref), where $\mathcal{T}^{*}$ is the adjoint operator of $\mathcal{T}$.\footnote{The adjoint operator is similar to matrix transpose in the finite-dimensional case. For formal definition, see for example han2009theoretical Chapter 2.6. For the specific adjoint operator used in our model-adaptive approach see (ref) below.}
For discrete state spaces, (ref) boils down to a finite-dimensional linear system. For continuous state spaces, we will use deterministic numerical integration to approximate the integral in (ref).
The adjoint operator $\mathcal{T}^{*}$ is defined with respect to an inner product space. The choice of the inner product does affect the convergence rate of the approximation error. Nevertheless, approximation solutions on different spaces all converge superlinearly under regularity conditions. To achieve the fastest decay of the approximation error when using CG, we propose to solve (ref) on $L^{2}(\mathbb{X})$ and define $\mathcal{T}^{*}$ by the inner product $\langle \cdot, \cdot \rangle$. (ref) formally discusses the convergence rate. (ref) shows the existence and uniqueness of the solution to (ref) on $L^{2}(\mathbb{X})$, denoted by $y^{*}$. Moreover, (ref) proves $V^{*} = (\mathcal{I} - \mathcal{T}^{*}) y^{*}$ where $V^{*}$ is the solution to (ref) on $L^{2}(\mathbb{X}, \mu)$.
The key idea of the conjugate gradient method is to iteratively build a solution by searching along directions that are conjugate (i.e., orthogonal with respect to the operator $(\mathcal{I} - \mathcal{T})(\mathcal{I} - \mathcal{T}^{*})$). At each iteration, the algorithm: (i) chooses a step size $\alpha_{k-1}$ that minimizes the residual along the current search direction $s_{k-1}$; (ii) computes the new residual $r_k$; and (iii) constructs the next search direction $s_k$ by combining the new residual with the previous direction, where the coefficient $\beta_{k-1}$ ensures conjugacy. The residuals $r_0, r_1, \ldots$ are mutually orthogonal by construction, and they form the basis of the model-adaptive sieve space. The model-adaptive approach is as follows:
The core idea behind our approach is, at each iteration, to augment the sieve space with the residual from the previous iteration. By construction, the updates are orthogonal to previous updates. And, it is easy to show that $y_{k} \in \text{span}\{r_{0}, r_{1}, \cdots, r_{k-1}\}$ where $\{r_{i}\}_{i \leq k-1}$ is the sequence of residuals produced by previous iterations.\footnote{The space is also called the Krylov subspace. The Krylov subspace of order $k$ generated by a matrix $\mathcal{A}$ and vector $b$ is $\mathcal{K}_{k}(\mathcal{A}, b) := \text{span}\{b, \mathcal{A} b, \mathcal{A}^{2} b, \ldots, \mathcal{A}^{k-1} b\}$ (See han2009theoretical Page 251).} In other words, $\left\{r_{i} \right\}_{i\leq k-1}$ is the model-adaptive sieve space after iteration $k$. Thus, after $k$ iterations the sieve dimension equals $k$, and the number of iterations directly determines the quality of the approximation. In (ref), we show that $y_{k}$ minimizes $\|(\mathcal{I} - \mathcal{T}^{*})y - V^{*}\|$ over the model-adaptive sieve space, and (ref) shows the superlinear convergence of the approximation error. The next section discusses the implementation and computational cost of our method.
Discrete State Spaces: For discrete state spaces, (ref) reduces to a finite-dimensional linear system as:
where $\mathcal{I}$ is the identity matrix and $\mathcal{T}$ is the discounted transition matrix. Therefore, (ref) only involves matrix-vector multiplications. The algorithm is as follows\footnote{We refer to judd1998numerical for other iterative methods such as Gauss--Jacobi and Gauss--Seidel algorithms.}:
Continuous State Spaces: For continuous state spaces, our method has to use numerical integration. We propose to use a deterministic numerical integration rule such as a Quasi-Monte Carlo rule to approximate the integral in (ref). Let $\mathbb{M} := \{x_{1},\cdots,x_{M} \}$ be the set of deterministic grid points used to approximate the integral. Note that implementing (ref) on $\mathbb{M}$ (with the transition density normalized) implicitly solves the following equation:
where $\hat{\mathcal{T}}_{M}$ is the matrix whose $(i, j)$-th element is $\beta \tilde{f}(x_{j}|x_{i})$, $\tilde{f}(x_{i}|x) := \frac{f(x_{i}|x)}{\sum_{j}f(x_{j}|x)}$ is the normalized transition density assuming the denominator is non-zero, $u_{M}$ is an $M$-dimensional vector with $u_{i} = u(x_{i})$ and $\mathcal{I}_{M}$ is an $M \times M$ identity matrix. Note that after solving (ref) at iteration $\hat{k}_{M}$ for all $x \in \mathbb{M}$, we can interpolate $\widehat{V}_{\widehat{k}_M}^{ma}(x)$ for $x \in \mathbb{X} \backslash \mathbb{M} $ using:
where $\widehat{V}^{ma}_{\widehat{k}_{M}}(x_{i}) := [(\mathcal{I}_{M} - \hat{\mathcal{T}}_{M}^{T}) y_{\widehat{k}_{M}}]_{i}$ and $y_{\widehat{k}_{M}}$ is the approximate solution to (ref).
The continuous state-space algorithm is as follows:
For discrete state spaces, the total computational cost of our method is $O(\hat{k}|\mathbb{X}|^{2})$ where $\hat{k}$ is the number of iterations required for convergence and $O(|\mathbb{X}|^{2})$ is the cost of matrix-vector multiplication. As written in (ref), each iteration requires three matrix-vector multiplications: two for $(\mathcal{I} - \mathcal{T})(\mathcal{I} - \mathcal{T}^{T}) y_{k}$ and one for $(\mathcal{I} - \mathcal{T}^{T})s_{k}$. In practice, the cost can be reduced to two matrix-vector multiplications per iteration by using the recurrence $r_{k} = r_{k-1} - \alpha_{k-1} (\mathcal{I} - \mathcal{T})(\mathcal{I} - \mathcal{T}^{T}) s_{k-1}$ to update the residual, since $(\mathcal{I} - \mathcal{T}^{T})s_{k-1}$ is already computed for $\alpha_{k-1}$ and one additional application of $(\mathcal{I} - \mathcal{T})$ yields $(\mathcal{I} - \mathcal{T})(\mathcal{I} - \mathcal{T}^{T}) s_{k-1}$. Due to the superlinear convergence, the number of iterations is expected to be small.
For continuous state spaces, the total computational cost is $O(\hat{k}_{M}M^{2})$ where $\hat{k}_{M}$ can vary with $M$. As discussed before, there is a trade-off between simulation error and computational time. A large $M$ leads to a small simulation error but a higher computational cost of solving the equation within the same tolerance. For different $M$, the cost of matrix-vector multiplication is determined by $O(M^{2})$. The algorithm still converges superlinearly as shown in (ref). Therefore, $\hat{k}_{M}$ is still expected to be small. Moreover, we will show that $\hat{k}_{M}$ is approximately the same for all sufficiently large $M$, which suggests that $\hat{k}_{M}$ is independent of $M$. Therefore, increasing $M$ primarily affects the computational cost through matrix-vector multiplication rather than $\hat{k}_{M}$. This property offers a significant computational advantage as it primarily relies on matrix-vector multiplication that is amenable to GPU acceleration and fast matrix-vector multiplication methods mentioned in the introduction.
This section discusses the theoretical properties of the model-adaptive approach. We will show the superlinear convergence of MA. We compare MA with TD and SA. For continuous state spaces, we will consider the simulation error from the numerical integration and prove the number of iterations is approximately the same for all sufficiently large numbers of grid points. We impose the following regularity conditions:
Under (ref), $\mathcal{T}$ maps $L^{2}(\mathbb{X},\mu)$ to itself. Moreover, $\mathcal{T}$ is a $\beta$-contraction with respect to $\|\cdot\|_{\mu}$ (see for example bertsekas2015dynamic). Consequently, (ref) has a unique solution on $L^{2}(\mathbb{X}, \mu)$. (ref) ensures that the solution is uniformly bounded by $\frac{C_{u}}{1-\beta}$. To achieve the fastest convergence rate of the approximation error when using CG, we propose to solve (ref) on $L^{2}(\mathbb{X})$. Assumptions (ref) and (ref) are used to prove the existence and uniqueness of the solution to (ref) on $L^{2}(\mathbb{X})$. Moreover, it $\mu$-almost surely equals the solution on $L^{2}(\mathbb{X},\mu)$, which are summarized in the following theorem:
Our approach first enjoys the following nice property:
(ref) shows that $V^{ma}_{k} := (\mathcal{I} - \mathcal{T}^{*})y_{k}$ minimizes the approximation error over $y \in \text{span}\{r_{0}, r_{1}, \cdots, r_{k-1}\}$. The orthogonality of the basis functions implies that the approximation error $\|V^{ma}_{k} - V^{*}\|$ decreases monotonically. As the residuals are informative about the solution, the projection onto the adaptive-sieve space can lead to a faster convergence rate of the approximation error than conventional methods.
Before establishing the convergence rate of our method, we define the concept of superlinear convergence. The concept of $R$-convergence, which is in analogy to the Cauchy root test for the convergence of series, quantifies the convergence rate of a sequence of approximation solutions:
(ref) establishes the superlinear convergence of the residual and the approximation error. It also shows the monotonic convergence of the approximation error. The decay rate of $c_{k}$ is at most $\frac{1}{k}$ and at least $\frac{1}{\sqrt{k}}$. Those two bounds hold for all inner product spaces under regularity conditions, while the lower bound is achieved by solving the equation on $L^{2}(\mathbb{X})$ with the continuity assumption on the partial derivatives of the transition density. As $c_{k} \to 0$, it is straightforward to show that algorithm (ref) will achieve the tolerance after a finite number of iterations. See the proof of (ref) for explicit expressions of the constants.
This section compares the convergence rate of the approximation error of our method with TD and SA. Denote by $V^{sa}_{k}$ and $V^{td}_{k}$ the SA and TD approximation solutions.\footnote{For TD, see tsitsiklis1996analysis and dann2014policy. For SA, see kress_linear_2014.} Let $\mathbb{S}_{k} := \text{span} \{\phi_{1}, \cdots, \phi_{k}\}$ be a sieve space where $(\phi_{1}, \cdots, \phi_{k})$ are basis functions, and let $\Pi_{\mathbb{S}_{k}}(V) := \operatornamewithlimits{argmin}_{h \in \mathbb{S}_{k}} \| V - h \|_{\mu}$ be the projection operator onto $\mathbb{S}_{k}$. We impose the following assumption on the projection bias of TD:
(ref) imposes upper and lower bounds on the projection bias. The upper bound is standard in the literature. Similar lower bounds on function approximation by neural nets can be found in yarotsky2017error Lemma 3. The following theorem establishes the convergence rate of SA and TD:
(ref) shows that the SA method converges linearly and the TD method converges sublinearly. It implies that to achieve a pre-specified tolerance, the number of iterations required for SA and the sieve dimension for TD can be much larger than for our method as it converges superlinearly.
For continuous state spaces, numerical integration introduces two sources of error. The main results of this section are: (i) the total error decomposes into a simulation error (from discretization) and an approximation error (from the MA iterations), and (ii) the number of MA iterations required for convergence is approximately independent of the number of grid points $M$. This is the as mesh independence principle. Together, these results imply that increasing $M$ improves accuracy without increasing the number of iterations.
Under smoothness and regularity conditions on the transition density and utility function (detailed in Appendix (ref)), these conditions are natural in practice, as the transition density is often very smooth; for example, autoregressive processes are commonly used to model the transition of state variables (e.g., erdem2003brand, hendel2006measuring, aguirregabiria2007sequential, aw2011r, and gowrisankaran2012dynamics). The approximation error converges superlinearly, and the number of iterations required for a given tolerance is approximately independent of $M$ --- a property known as the mesh independence principle (see atkinson1997numerical). The formal statement is as follows:
This section presents three numerical experiments. First, we simulate a bus engine replacement model to visualize the convergence of MA. Second, we analyze a model of consumer demand for storable goods similar to hendel2006measuring. We compare policy iteration and Newton--Kantorovich outer iterations combined with different inner solvers (MA, SA, and an exact solver), along with VFI benchmarks. The results show that MA opens the door to the use of Bayesian MCMC estimators for such models. Finally, we examine a single-firm entry and exit problem described in aguirregabiria2023solution. We show that MA can improve the computational efficiency of policy iteration. We also compare the performance of MA against: SA and TD.
This section simulates a bus engine replacement problem to visualize the convergence behavior of MA. We adapt the setting in arcidiacono2011conditional.
At each period $t \leq \infty$, an agent chooses to maintain $ a_{t} = 1$ or replace $a_{t} = 0$ the engine. The replacement cost is $RC$. The maintenance cost is linear in mileage with accumulated mileage up to 25, i.e., $ u(x_{t},1) = \theta_{1} \min\{x_{t},25\}$, where $x_{t}$ is the mileage of the engine. Moreover, mileage accumulates in increments of 0.125. The period utility of the agent is $u(x_{t}, a_{t}, \varepsilon_{t}) = (1 - a_{t}) (RC + \varepsilon_{t}(0)) + a_{t} (\theta_{1} \min\{x_{t},25\} + \varepsilon_{t}(1))$ where $(\varepsilon_{t}(0), \varepsilon_{t}(1))$ are i.i.d extreme value type I distributed shocks. The transition probability of $x_{t}$ is specified as:
where we set $\theta_{1} = -0.15, \theta_{2} = 1, RC = -2$, and $\beta = 0.9$. To visualize, we solve the DP problem and use the true CCPs to construct the linear systems.
Our method takes 15 iterations to solve the equation. (ref) visualizes our approximation solution for the first 9 iterations and the basis functions. Each panel shows the impact on the approximate solution of adding one additional sieve basis function to the previous value.
(ref) plots the $L_{2}$-norm of the residuals and the sup-norm of the approximation error. The approximation solution is close to the true solution after iteration 6, which implies that our algorithm constructs a good sieve space using 6 iterations. After that, the approximation converges to the true solution rapidly. This finding aligns with (ref) as the approximation error dramatically decreases at iteration 6-7. Other norms also decay dramatically after 6 iterations. Moreover, the $L_{2}$-norm of the approximation errors decreases monotonically, consistent with (ref).
In this section, we analyze a model of consumer demand for storable goods similar to hendel2006measuring. We aim to show that our method opens the door to the use of Bayesian MCMC estimation methods for such models.
At each period, given prices $\bm{p}$ and inventories $I$, a household of size $n$ decides which brand $k$ to purchase, how much to purchase $j$ and how much to consume $C$.\footnote{Further details of the model are presented in (ref). For each household size, parameter values are set equal to the estimated values reported in (ref).} Let $\omega_{j}\left(\bm{p} \right)$ be the indirect utility from brand choice and $x := \{\bm{\omega},I\}$. Given the assumptions in hendel2006measuring, brand choice is purely a static problem and the consumer's value function $V(x)$ for the dynamic problem is the solution to:
where $C(x,j) := \operatornamewithlimits{argmax}_{c_{min} \leq c \leq c_{max}} \Bigl[ U(c,I,j;\theta) + \omega_{j} \left(\bm{p} \right) + \beta \mathbb{E}[V(x') | x, c,j] \Bigr]$.
We compare several solution methods combining two outer iteration frameworks (PI and NK) with different inner solvers (MA, SA, and an exact linear equation solver\footnote{We use scipy.sparse.linalg.spsolve.}), along with VFI and one-step VFI as benchmarks. The pseudocode for PI+MA and its VFI and one-step VFI variants is given in Appendix (ref).\footnote{For PI and NK algorithms, the stopping rules are $\|p^{i+1} - p^{i}\|_{\infty} \leq 10^{-4}$, $\|r_{k}\|_{\infty} \leq 10^{-8}$, and $\|C^{i+1} - C^{i}\|_{\infty} = 0 $. For VFI and one-step VFI, the stopping rules are $\|V^{i+1}(x) - V^{i}(x)\|_{\infty} \leq 10^{-8}$, and $\|C^{i+1} - C^{i}\|_{\infty} = 0 $.} For a given consumption function, we update the value function using either PI or NK outer iterations with the chosen inner solver, then update the consumption function and iterate until convergence. The one-step updating was used in osborne2018approximating to implement the Bayesian estimator proposed by imai2009bayesian—henceforth, IJC. At each MCMC step, rather than iterate to convergence, the algorithm updates both the value function and the consumption function only by one-step.
(ref) compares the computational time and number of iterations required by each algorithm across simulations for different household sizes. NK+MA is the fastest method, requiring approximately 0.5 minutes across all household sizes, followed by PI+MA at around 0.7 minutes. NK+SA requires about 1.0 minutes and PI+SA about 2.4 minutes. Exact solution takes about 4.6 minutes. VFI converges in about 1.5 minutes. These results demonstrate that combining MA with NK outer iterations substantially improves computational efficiency. One-step VFI requires about 64 minutes, more than 100 times slower than NK+MA.
Our method opens the door to the use of Bayesian MCMC estimators. To simulate 10,000 MCMC steps, (ref) shows that PI+MA requires between 1.6 and 3.3 days depending on household size.\footnote{At each MCMC step, the value function obtained from the previous MCMC step is used as the initial guess for the current dynamic programming problem, thereby reducing the computational time.}
An alternative estimator to the Bayesian MCMC estimator is the IJC approach. However, the performance of the one-step VFI suggests that the IJC approach may not be the most suitable estimator for this problem. In (ref), the one-step VFI approach requires more than 1600 iterations to converge even when the true parameters are known. Consequently, due to the extremely large number of iterations required, we expect the total computational time using the IJC approach would be substantially higher than our proposed approach.
This section examines a single-firm entry and exit problem described in aguirregabiria2023solution. We compare the performance of the model-adaptive approach (MA) against successive approximation (SA) and temporal difference (TD) as inner solvers for the linear system (ref). These inner solvers are embedded within policy iteration (PI) and Newton--Kantorovich (NK) outer iterations. We also include value function iteration (VFI) as a benchmark and examine how computational time scales with the state space size.
At each period $t \leq + \infty$, a firm decides whether to exit ($a_{t}=0$) or enter ($a_{t}=1$) the market. For an active firm, the profit is $\pi(1,x_{t}) + \varepsilon_{t}(1)$. $\pi(1,x_{t})$ equals the variable profit $VP_{t}$ minus fixed cost $FC_{t}$, and minus entry cost $EC_{t}$. For an inactive firm, $\pi(0,x_{t})$ is normalized to be 0, and the profit is $\varepsilon_{t}(0)$. The variable profit is $VP_{t} = (\theta_{0}^{VP} + \theta_{1}^{VP} z_{1t} + \theta_{2}^{VP} z_{2t}) \exp(w_{t}) $ where $w_{t}$ is the productivity shock, $z_{1t}$ and $z_{2t}$ are exogenous state variables that affect price-cost margin. The fixed cost is $FC_{t} = \theta_{0}^{FC} + \theta_{1}^{FC} z_{3t}$, and the entry cost is $EC_{t} = (1-a_{t-1}) (\theta_{0}^{EC} + \theta_{1}^{EC} z_{4t})$ where $(1-a_{t-1})$ indicates that the entry cost is paid if the firm is inactive at the previous period ($a_{t-1} = 0$). Continuous state variables follow AR(1) process: $z_{jt} = 0.6 z_{jt-1} + \epsilon_{jt}$, $w_{t} = 0.2 + 0.6 w_{t-1} + \epsilon_{t}$, where $\epsilon_{jt}$, $\epsilon_{t}$ follows i.i.d standard normal. The true parameters $\theta^{*}$ are chosen to be $\theta_{0}^{VP} = 0.5, \theta_{1}^{VP} = 1, \theta_{2}^{VP} = -1, \theta_{0}^{FC} = 1.5, \theta_{1}^{FC} = 1, \theta_{0}^{EC} = 1, \theta_{1}^{EC} = 1$. We use $M$-point Tauchen's method (tauchen1986finite) to discretize each of 5 continuous state variables and obtain the transition matrix where $M = \{6,7,8,9,10\}$. Moreover, we set the discount factor $\beta = \{0.95, 0.975, 0.980, 0.985, 0.990, 0.995, 0.999\}$.
The transition matrix $F$ admits a Kronecker product factorization $F = F_{a} \otimes F_{1} \otimes F_{2} \otimes \cdots \otimes F_{d}$, where $F_{a}$ corresponds to the discrete action state and each $F_{j}$ is the $M \times M$ Tauchen transition matrix for the $j$-th continuous variable. This structure allows computing the matrix-vector product $Fv$ without ever forming or storing $F$ explicitly. The cost of one matrix-vector product is $O(d M |\mathcal{X}_{M}|)$, compared with $O(|\mathcal{X}_{M}|^{2})$ for a dense matrix-vector product. Since each MA iteration requires two such products (one with $(\mathcal{I} - \beta F)$ and one with its adjoint), the total cost of $k$ MA iterations is $O(k \, d \, M \, |\mathcal{X}_{M}|)$. The memory requirement is $O(d M^{2} + |\mathcal{X}_{M}|)$—storing $d$ matrices of size $M \times M$ plus the iterate—rather than $O(|\mathcal{X}_{M}|^{2})$ for the dense matrix.
We compare the computational performance of the above methods for solving the dynamic programming problem with $M=6$ grid points per dimension ($|\mathcal{X}_{M}| = 15{,}552$) and $\beta = 0.95$. All methods exploit the Kronecker product structure of the transition matrix for efficient matrix-vector multiplication.\footnote{The stopping rules are $\|p_{k+1} - p_{k}\|_{\infty} \leq 10^{-8}$ for PI, $\|r_{k}\|_{\infty} \leq 10^{-8}$ for the inner solvers, and $\|\Gamma(V) - V\|_{\infty} \leq 10^{-8}$ for NK outer iterations. Initial guesses: $p_{0} = \frac{1}{2}, y_{0} = 0$.}
(ref) presents the results. Panel A reports the solve time for the linear system (ref) given the true CCP, which is obtained by solving the DP problem to convergence using VFI. MA converges in 120 iterations (0.16 seconds), outperforming SA which requires 389 iterations (0.21 seconds). TD produces a large residual of $5.31 \times 10^{2}$, indicating that the polynomial sieve space is inadequate for approximating the value function in this model.
Panel B compares full-solution methods that solve the DP problem directly. VFI converges in 479 iterations (1.06 seconds). Among PI-based methods, PI+MA (0.72 seconds) outperforms PI+SA (0.91 seconds). NK+MA achieves the fastest total time (0.44 seconds). The speed advantage of NK over PI reflects two factors: NK converges in 4 outer iterations versus 5 for PI, and its quadratic outer convergence yields fewer total inner iterations (307 versus 514).
In (ref), we visualize the convergence of MA for $M=6$ given the true CCP. We plot the sup-norm of residuals, $\|r_{k}\|_{\infty}$, $L_2$-norm of the residuals, $\|r_{k}\|_{2}$, the sup-norm of the approximation error, $\|V_{k} - V^{*}\|_{\infty}$, and $L_2$ norm of the approximation error, $\|V_{k} - V^{*}\|_{2}$. All figures are plotted as functions of the iteration count $k$ for different $\beta$. (ref) exhibit a similar pattern: they increase initially, reach a peak, and then decrease rapidly. (ref) shows the $L_{2}$-norm of approximation error decreases monotonically consistent with (ref). Notably, while $\beta = 0.999$ initially results in larger residuals and approximation errors compared to smaller discount factors, it still achieves convergence within a comparable number of iterations. This demonstrates the ability of MA to handle large discount factors, which is challenging for SA.
We examine how the computational cost of PI+MA scales with the state space size.\footnote{An alternative approach to solve DP is the Euler-Equation method (aguirregabiria2023solution). However, it works for models where the only endogenous state variable is the previous action (see their Definition 1). This feature is satisfied in the simulation, though it is restrictive in general.} (ref) presents the average number of MA iterations and average computational time per policy iteration step for $M = 6, 7, 8, 9, 10$. PI converges after 5 iterations for all $M$ and $\beta$.
These results provide empirical support for (ref). As $M$ increases, the average number of iterations remains relatively stable for all $\beta$, consistent with the mesh independence property. For instance, at $\beta = 0.95$, the average number of iterations is 103 for $M=6$ and 96 for $M=10$. The computational time scales approximately linearly in $|\mathcal{X}_{M}|$, consistent with the $O(d \, M \, |\mathcal{X}_{M}|)$ cost of Kronecker matrix-vector multiplication described in (ref). Additionally, there is a slight increase in the average number of iterations (from around 96 to 120) as $\beta$ increases from 0.95 to 0.999 at $M=10$. As it only takes up to 1.5 seconds per policy iteration step to solve the linear equation with $|\mathcal{X}_{M}| = 200{,}000$, MA improves the computational efficiency of policy iteration where the main computational cost is solving the linear system of equations.
This section applies our method to a dynamic consumer demand model for laundry detergent using Kantar's Worldpanel Take Home data. We first describe the data, model, and estimation procedure. Then, we discuss the implication of the results.
The analysis of the Great Britain laundry detergent industry is based on the data from 1st January 2017 until 31st December 2019. It captures detailed information on a representative sample of British households' purchases of fast-moving products, including food, drink, and laundry detergents. The data has been used in previous studies such as dubois2014prices,dubois2020well. Households use barcode scanners to record all their grocery purchases. For each purchase, the data includes key information such as price, quantity, product characteristics, and purchase date.
We consider the market for laundry detergent. A laundry detergent product is defined by its quantity, brand, and chemical properties (bio/non-bio). Laundry detergents are available in various formats such as liquid, powder, and gel, each with different dosage metrics. To standardize quantity across formats, we define the quantity purchased by the number of washes. (ref) presents the top 10 brands and bio/non-bio combinations, which account for 76.82% of the total observed purchases. We restrict our analysis to these top 10 brand and bio/non-bio combinations, and assume they are available to all consumers.
In Figure (ref) we show the histogram of quantities purchased\footnote{We restricted our analysis to products with number of washes between 10 and 100.}, which suggests three natural clusters corresponding to small, medium, and large sizes. Therefore, we use $3$-means clustering to aggregate the quantities for each brand and bio/non-bio combination into small (23), medium (40), and large (64) sizes. (ref) summarizes the number of observations for each cluster. We then calculate the weekly average transaction price for each product cluster, and use it as the price index in our model.
In (ref) we describe purchase statistics by household size. We restrict our analysis to households who purchase laundry detergent at least 3 times and at most 15 times each year for household sizes 1-4, and 3 to 20 times each year for household size 5. The average size of laundry detergent purchased increases with household size, ranging from 32.2 washes per purchase for single-person households to 38.1 washes for households with four or five. The average price per wash, conditional on purchase, shows a slight decreasing trend as household size increases, from £0.152 for single-person households to £0.142 for the largest households. The average number of weeks between purchases decreases as household size increases, with single-person households making purchases every 9.7 weeks on average, while households with five members purchase every 7.4 weeks. The average quantity of laundry detergent purchased per year increases with household size, ranging from 165 washes per year for single-person households to 256 washes for the largest households. (ref) indicates heterogeneity in purchasing behavior across different household sizes. Therefore, we separately estimate the model for each household size.
In the model, laundry detergent is storable, and households derive utility from both consumption and purchase. As in hendel2006measuring, we assume that laundry detergents are perfect substitutes in consumption, which implies that the unobserved state variable, inventory, is one-dimensional. At week $t \leq + \infty$, consumer $i$ chooses discrete consumption $c_{t}$ and purchase decision $a_{jklt}$, where $j \in \{0,23,40,64\}$ is the quantity measured by the number of washes, $k$ refers to the brand, and $l$ is bio/non-bio. Let $\bm{a}_{t}$ be a vector of purchase decisions with $\sum_{j,k,l} a_{jklt} = 1$, where $a_{jklt} = 1$ indicates the purchase of quantity $j$ of brand $k$ and bio/non-bio $l$, and $a_{jklt} = 0$ otherwise. The period utility\footnote{Our model differs from hendel2006measuring as we do not include a taste shock for the consumption level $c_{t}$. osborne2018approximating also found that it is difficult to identify its distribution.} is given by:
where $I_{t}$ is the inventory at the beginning of $t$ and $\bm{p}_{t}$ is a vector of prices. $\varepsilon_{jklt}$ is the choice-specific utility shock following an i.i.d. Type I Extreme Value distribution. $(\theta_{1},\theta_{2})$ capture the marginal utility of consumption. $\theta_{3}$ captures the carrying cost that depends on the inventory at next period: $I_{t+1} := I_{t} + j_{t} - c_{t}$. The fixed cost of making a purchase is $\theta_{4}$. We impose an upper bound on the inventory $ I_{t} \in \{0 ,1, \cdots, I_{max}\}$ where\footnote{We impose the upper bound of 80 washes because it approximately equals one-third of the average annual washes (256) for the largest household size. This bound helps maintain computational feasibility while still capturing consumers' stockpiling.} $I_{max} = 80$. The consumption is bounded from below and above by $c_{min,t} := \max\{0, I_{t} + j_{t} - I_{max}\}$ and $c_{max,t} := I_{t} + j_{t}$, respectively. $\theta_{5}$ captures the price sensitivity. $\xi_{k}$ and $\delta_{l}$ are brand and bio/non-bio fixed effects per wash.
A household maximizes expected life-time utility:
where $\beta=0.99$. The state space is defined by $s_{t} := (I_{t},\bm{p}_{t})$.
This section presents the estimation procedure. We mainly follow the three-step procedure in hendel2006measuring. We assume $\bm{p}_{t}$ follows an exogenous first-order Markov process. Therefore, inventory $I_{t}$ is the only endogenous state variable. This allows us to decompose the decisions into dynamic decisions $(c_{t},j_{t})$ and static decisions $(k_{t},l_{t})$, as the evolution of $I_{t}$ is determined only by $(c_{t},j_{t})$. Moreover, as shown in hendel2006measuring, the choice probability of purchasing brand $k$ and bio/non-bio $l$ conditional on purchasing quantity $j$ has the conditional logit form:
allowing for simple estimation of $(\theta_{5},\xi_{k},\delta_{l})$.
Then, the inclusive value of purchasing quantity $j$ is defined by:
where $\omega_{j,t}$ is the indirect utility of purchasing quantity $j$ at period $t$. We impose the inclusive value sufficiency (IVS) assumption as in hendel2006measuring:
Given (ref), the lagged inclusive value $\bm{\omega}_{t-1}$ is a sufficient statistic to forecast $\bm{\omega}_{t}$. As a result, the state space can be reduced from $s_{t} := (I_{t},\bm{p}_{t})$ to $x_{t} := \{I_{t}, \bm{\omega}_{t}\}$. We forecast $\bm{\omega}_{t}$ using a VAR(1) and then discretize it into $5^{3}$ bins. The simplified dynamic programming problem based on the IVS assumption is:
where $U(c,I,j;\theta) = \theta_{1} c + \theta_{2} c^{2} + \theta_{3} I'^{2} + \theta_{4}\mathbbm{1}(j > 0)$, $I' := I + j - c$ is next-period inventory, $\varepsilon$ is a 4-dimensional vector of i.i.d Type I Extreme Value distribution, and $V(x) = \mathbb{E}_{\varepsilon} V(x,\varepsilon)$.
We decompose the joint optimization problem into sequential optimization problems: first choose purchase quantity $j$, then choose consumption. For the second problem, given state $x$, purchase quantity $j$, and $V(x)$, define the consumption function as:
For the first problem, the choice probability of purchasing quantity $j$ is:
where $v(x,j) := U(C(x,j),I,j;\theta) + \omega_{j} + \beta \mathbb{E}[V(x') | x, C(x,j),j]$. Combining these two steps, we have:
For a given consumption function, we solve (ref) by policy iteration with our model-adaptive approach. Then, we update the consumption function using (ref), and iterate until convergence (see (ref)).
As the argmax operator in (ref) for discrete outcomes is not differentiable, we estimate the dynamic parameters by the Bayesian MCMC estimator. The estimation details are described in Appendix (ref).
We separately estimate the model for different household sizes. (ref) presents the results of the conditional logit model. The coefficient $\theta_{5}$ measures the price sensitivity, which is negative and statistically significant at 5% level across all household sizes. Single-person households exhibit the highest price sensitivity at -0.501, followed by households with three members at -0.416, households with two members at -0.345, and households with four members at -0.307. Five-person households have the lowest price sensitivity at -0.223. The estimates also suggest the heterogeneity in price sensitivity across different household sizes.
In (ref) we report the estimates of dynamic parameters and computational times.\footnote{At each MCMC step, the value function obtained from the previous MCMC step is used as the initial guess for the current dynamic programming problem.} All estimates have the expected signs and are statistically significant at the 5% level. The utility from consumption, determined jointly by the linear term ($\theta_1$) and quadratic term ($\theta_2$), shows varying patterns across household types. The linear coefficient $\theta_1$ ranges from its lowest value of 1.878 for four-person households to its peak of 3.583 for three-person households. Five-person households have the lowest magnitude for the quadratic term $\theta_2$ (-8.423). Single-person households have the lowest fixed cost of making a purchase ($\theta_4 = -4.195$), while five-person households have a relatively high fixed cost ($\theta_4 = -5.349$). The average time ranges from 0.23 to 0.48 minutes per Metropolis--Hastings step. The total time of the MCMC estimator ranges from 1.6 to 3.3 days, which confirms the computational efficiency of PI+MA in practice.
To investigate the model fit, we compare the simulated and observed purchase behavior. For each household size, we simulate consumption and purchase behavior for the same number of households as observed in the data over 156 weeks (3 years).\footnote{The first 30% periods of both the simulated and observed data are discarded as it is used to simulate the initial inventory for the estimation.} The simulated market shares in (ref) reasonably match the observed data. Figure (ref) plots the hazard rates and their confidence intervals. We also estimate the static model without the carrying cost.\footnote{For the static model, households consume the entire pack when a purchase is made. Therefore, $U(\omega_{t},j_{t},\varepsilon_{t}):= \theta_{1}^{static} j_{t} + \theta_{2}^{static} j_{t}^{2} + \omega_{jt} + \theta_{4}^{static} \mathbbm{1}(j_{t} > 0) + \varepsilon_{jt}$. The hazard rate is then the probability of making a purchase.} All household sizes exhibit a similar pattern, with the hazard rates initially increasing until stabilizing at the highest level. For single-person and two-person households, the hazard rates increase until around 8 weeks, while the hazard rates increase until around 5 weeks for larger households. The hazard rates all peak around 15%. The dynamic model outperforms the static model in capturing the hazard rates as the static model, by construction, cannot capture the increasing hazard rate observed in the data.
Long-run elasticities measure the effects of permanent price changes on market shares. (ref) simulates the long-run elasticities for different household sizes. Own-price elasticities are larger in absolute value for larger pack sizes. They are all greater than 1 in absolute value, except the small pack for household size 5 (-0.742), indicating that the demand is generally elastic. For single-person households, the own-price elasticity for large packs is -4.067, meaning that if the prices of all large packs increase by 1%, the market shares of large packs will decrease by 4.067%. As household size increases, there is a trend towards lower price elasticities, suggesting that larger households are less sensitive to price changes. This trend can partially be explained by the price coefficients ($\theta_5$) as larger households are less price sensitive, except for household size 3, as shown in (ref).
The demand for each size is elastic in general, indicating consumers are likely to substitute out of the pack size that has increased in price and to other pack sizes or not purchase at all. For single-person households, the cross-price elasticity between medium and large packs is 0.601, higher than the cross-price elasticity between medium and small packs (0.256). This pattern continues across household sizes - for 2-person households, the medium-large cross-price elasticity is 0.378, for 3-person households it increases to 0.628, showing stronger substitution effects for larger sizes. The cross-elasticities between medium and small packs remain relatively stable across household sizes (ranging from 0.256 to 0.268 for households of 1-3 persons), while the substitution effects involving large packs tend to be stronger. These findings suggest that consumers are more likely to substitute between larger pack sizes when prices change, particularly in households with 2-3 members, though this effect diminishes somewhat for the largest households. Overall, the magnitudes of our long-run elasticities are consistent with prior estimates in the literature: hendel2006measuring report long-run own-price elasticities between $-2.5$ and $-4.3$ for 128 oz.\ liquid detergent, which aligns well with our estimates for larger pack sizes.
In this paper, we propose a model-adaptive approach to solve the linear system of fixed point equations of the policy valuation operator. Our method adaptively constructs the sieve space and chooses its dimension. It converges superlinearly, while conventional methods do not. We demonstrate through simulations that our model-adaptive sieves dramatically improve the computational efficiency of policy iteration and Newton--Kantorovich iterations, and open the door to the use of Bayesian MCMC estimators for DDC models.
We apply our method to analyze consumer demand for laundry detergent using Kantar's Worldpanel Take Home data. The empirical application confirms that our approach improves computational efficiency of policy iteration in practice, which opens the door to the estimation of DDC models by Bayesian MCMC estimators. To investigate the model fit, we simulate market shares and hazard rates for different household sizes. The model achieves a reasonable model fit. We also simulate the long-run elasticities. The results show the heterogeneous substitution patterns across different household sizes.