Solving Dynamic Discrete Choice Models Using Smoothing and Sieve Methods

Solving Dynamic Discrete Choice Models Using Smoothing and Sieve Methods

Keywords: Dynamic discrete choice; numerical solution; Monte Carlo; sieves.

\thispagestyle{empty}

\pagenumbering{arabic} \setcounter{page}{1}

Introduction

Discrete Decision Processes (DDPs) are widely used in economics to model forward-looking discrete decisions. For their implementation, researchers are required to solve the model which generally cannot be done in closed form. Instead, a number of methods have been proposed for solving the model numerically; see, e.g., rust2008dynamic for an overview. We propose two novel methods for approximating the solutions to a general class of Markovian DDP models in terms of either the so-called integrated or expected value function. These two functions are relevant for estimation of DDP's and for welfare analysis of policy experiments. Our framework allows for both continuous and discrete state variables, non-separable utility functions and unrestricted dynamics. As such, we cover most relevant models used in empirical work. The proposed implementation of model and estimators are found to be computationally very efficient, and at the same time providing precise results with small approximation errors due to the use of simulations and sieve methods.

Our first proposal proceeds in three steps: First, we develop smoothed simulated versions of the Bellman operators that returns the integrated and expected value functions as fixed points. Next, we approximate the unknown value function by a sieve, that is, a parametric function class, thereby turning the problem into a finite-dimensional one. Finally, we solve for the parameters entering the chosen sieve using projection-based methods. When the chosen sieve is linear in the parameters, the approximate solution can be computed using an iterative procedure where each step is on closed form.

As an alternative to the above sieve-based method, we also adapt and generalize the so-called “self-approximating” method proposed in Rust1997curse to our setting: We design the importance sampler used in the simulated Bellman operators so that the corresponding expected and integrated value functions can be solved for directly without the use of sieves. In comparison with the sieve approach, the self-approximating solution method has the advantage that it will not suffer from any biases due to function approximations. But at the same time, the importance sampler used in its implementation will generally have a larger variance compared to the class of samplers that can be used for the sieve method. This larger variance also translates into a larger simulation bias of the self-approximating solution due to the non-linear nature of the problem. Thus, neither method strictly dominates the other.

Our two procedures, the sieve-based and self-approximating one, differ from existing proposals in three important aspects: First, we solve for either the integrated or expected value function instead of the value function itself. This reduces the dimensionality of the problem since we integrate out any i.i.d. shocks appearing in the model before solving it. Moreover, while the value function is non-differentiable, the integrated and expected value functions are generally smooth which means that our sieve method performs better compared to existing ones that aim at approximating the value function. Second, we allow for a general class of importance samplers in the simulation of the Bellman operator; these can be designed to reduce variances and biases due to simulations. Third, we smooth the simulated Bellman operator by replacing the $\max$-function appearing in its expression by a smoothed version where the degree of smoothing is controlled by a parameter akin to the bandwidth in kernel smoothing methods. This is similar to the logit-smoothed accept-reject simulator of probit models as proposed by McFadden1989; see also Fermanian&Salanie2004, Kristensen&Shin20112 and Iskhakovetal2017. The smoothing turns the problem of solving for the integrated and expected value functions into differentiable ones. In particular, the exact solutions to the smoothed simulated Bellman equations become smooth as functions of state variables and any underlying structural parameters. This in turn means that standard sieves, such as polynomials, will approximate the exact solutions well and that we can control the error rate due to function approximation. Moreover, if used in estimation, standard numerical solvers can be employed in computing estimators of the structural parameters. The smoothing entails an additional bias but this can be controlled for by suitable choice of aforementioned smoothing parameter.

The smoothing device also facilitates the theoretical analysis of the approximate value functions since it allows us to use a functional Taylor expansion of it. This expansion is then used to analyze the leading numerical error terms of the approximate value functions due to simulations, smoothing and function approximations. In particular, under regularity conditions, we show that the approximate value function will converge weakly towards a Gaussian process which is the first result of its kind to our knowledge. These results allow researchers to, for example, build confidence intervals around the approximate value function and should be useful when analyzing the impact of value function approximation when used in welfare analysis and estimation of structural parameters. They may also be potentially helpful in designing selection rules for number of basis functions and the smoothing parameter.

A numerical study investigates the performance of the solution methods in practice. We implement the proposed methods for the engine replacement model of Rust1987 and investigate how smoothing, number of basis functions and number of simulations affect the approximation errors. We also investigate how the procedures are affected by the dimensionality of the problem and how derivative-based solvers affect computation times. We find that the sieve method generally performs best of the two methods: It is computationally faster and in most situations provides a better approximation in terms of bias and variance. Moreover, the sieve method is found to also work well in higher dimensions with its bias and variance being fairly stable as we increase the the number of state variables of the model. In contrast, variances of the self-approximating method increase dramatically as the number of state variables increases and so appears to be less robust. Finally, the errors due to simulations and function approximation behave according to theory and are found to vanish at the expected rates.

Our proposed methods share similarities with the ones developed in, amongst others, Arcidiaconoetal2013, Keane&Wolpin1994, Munos&Szepesvari2008, Norets2012, Pal&Stacurski2013 and Rust1997curse who also use simulations and/or sieve methods to solve DDP's. However, except for Keane&Wolpin1994, the methods proposed in these papers approximate the value function while ours target the integrated or expected value function which are more well-behaved (smooth) objects and therefore easier to approximate. Moreover, in contrast to the cited papers, we employ importance sampling and smoothing in our implementation which comes with the aforementioned computational advantages. From a theory perspective, we provide a more complete asymptotic analysis of the approximate integrated and expected value functions. On the other hand, Munos&Szepesvari2008 and Rust1997curse provide an analysis of the computational complexity of solving for the value function and so the theories of this paper and these studies complement each other.

The remains of the paper are organized as follows: Section (ref) introduces a general class of DDP's and their corresponding value functions. In Section (ref), we develop our smoothed simulated versions of the Bellman operators that the integrated and expected value functions are fixed points to. We then show how to (approximately) solve these simulated Bellman equations in Section (ref). An asymptotic theory of the approximate value function is presented in Section (ref), while the results of the numerical experiments are found in Section (ref). Appendix (ref) contains some general results for approximate solutions to fixed point problems, while proofs of the main results can be found in Appendix (ref).

Model

We consider the following DDP where a single agent at time $t\geq1$ solves

where $\mathcal{D}=\left\{ 1,...,D\right\} $ is the set of alternatives, $u(S_{t},d)$ is the per-period utility, $0<\beta<1$ is the discount factor, $S_{t}$ is a set of state variables that follows a controlled Markov process with transition kernel $F_{S}\left(S_{t}|S_{t-1},d_{t-1}\right)$ and the so-called value function $\nu$ solves the following fixed-point problem,

Following Rust1987 and many subsequent empirical specifications, we assume that $S_{t}=\left(Z_{t},\varepsilon_{t}\right)\in\mathcal{Z}\times\mathcal{E}\subseteq\mathbb{R}^{d_{Z}}\times\mathbb{R}^{d_{\varepsilon}}$ where $Z_{t}$ and $\varepsilon_{t}$ satisfy the following conditional independence condition, \[ F_{S}\left(Z_{t},\varepsilon_{t}|Z_{t-1},\varepsilon_{t-1},d_{t-1}\right)=F_{\varepsilon}\left(\varepsilon_{t}|Z_{t}\right)F_{Z}\left(Z_{t}|Z_{t-1},d_{t-1}\right). \] In many cases $F_{\varepsilon}\left(\varepsilon_{t}|Z_{t}\right)=F_{\varepsilon}\left(\varepsilon_{t}\right)$ in which case $\varepsilon_{t}$ is an i.i.d. sequence and so can be thought of as idiosyncratic shocks to utility. If no shocks are present in the model, we can always choose $\varepsilon_{t}=\emptyset$ to be an empty variable so that $S_{t}=Z_{t}$. Throughout, we will assume that the support $\mathcal{Z}$ is a compact set. This is done to simplify the theoretical analysis since it, for example, implies that value functions defined below will lie in the space of bounded functions on $\mathcal{Z}$, $\mathbb{B}\left(\mathcal{Z}\right)$, equipped with the sup-norm, $\left\Vert v\right\Vert _{\infty}=\sup_{z\in\mathcal{Z}}\left|v\left(z\right)\right|$. At the same time, we allow the support of the error term, $\mathcal{E}$, be unbounded and for both countable and continuously distributed state variables.

In the above formulation, the model is characterized by the value function $\nu(s)$. However, it is possible to rewrite the models in terms of either the so-called integrated value function or the expected value function and solve for these instead. These are defined as \[ v(Z_{t})=E\left[\nu(Z_{t},\varepsilon_{t})|Z_{t}\right]=\int_{\mathcal{E}}\nu\left(Z_{t},e\right)dF_{\varepsilon}\left(e|Z_{t}\right), \] and \[ V(Z_{t},d_{t})=E\left[\nu(Z_{t+1},\varepsilon_{t+1})|Z_{t},\varepsilon_{t},d_{t}\right]=E[v(Z_{t+1})|Z_{t},d_{t}]=\int_{\mathcal{E}}v(z^{\prime})dF_{Z}\left(z^{\prime}|Z_{t},d_{t}\right), \] respectively, where we have used the conditional independence assumption. Observe that given $v\left(z\right)$, we can recover $V\left(z,d\right)=E[v(Z_{t+1})|Z_{t},d_{t}]$ which in turn can be used to compute $\nu(S_{t})=\max_{d\in\mathcal{D}}\left\{ u(Z_{t},\varepsilon_{t},d)+\beta V\left(Z_{t},d\right)\right\} $. Thus, there is no loss in focusing on the integrated and expected value function, except that in the former case we need to compute $E[v(Z_{t+1})|Z_{t},d_{t}]$ numerically to obtain the value function. Furthermore, in many cases the expected value function itself is of interest. For example, the conditional choice probabilities, which are needed for counterfactuals and for estimation, take as input the relative expected value function, \[ P\left(d_{t}=d|Z_{t}=z\right)=M_{u,d}(\beta\Delta V(z)|z),\:M_{u,d}(r|z)=\frac{\partial M_{u}(r|z)}{\partial r\left(d\right)}, \] where $\Delta V(z,d)=V(z,d)-V(z,D)$, $d\in\mathcal{D}$, and$M_{u}(r|z)$ is a generalized version of the so-called social surplus function defined as, for any $r=(r(1),...,r(D)),$

It is also useful in welfare analysis of policy experiments where we wish to see how a policy change will affect the expected present value of lifetime utility (i.e., the expected value function).

Except for a few special cases, analytical expressions of $v$ and $V$ are not available and so numerical approximations have to be employed. We will here develop numerical methods for solving for either $v$ or $V$ instead of $\nu$ for the following reasons: First, $\nu$ is a function of $s=\left(z,\varepsilon\right)$ while $V$ and $v$ are functions of $z$ alone and therefore their approximations are lower-dimensional problems. Second, $\nu$ is non-differentiable due to the max-function in ((ref)); in contrast, $v(z)$ and $V(z,d)$ are both smooth functions of $z$ if $F_{\varepsilon}(e|z)$ and $F_{Z}\left(z^{\prime}|z,d\right)$ are. If there is no i.i.d. component in the model, $\varepsilon_{t}=\emptyset$, then $\nu\left(s\right)=\nu\left(z\right)=v(z)$ and so the integrated value function becomes non-smooth. In contrast, $V(z,d)$ remains smooth even in this case. The functions $v$ and $V$ each solves their own fixed-point problem: Taking conditional expectations on both sides of eq. ((ref)), $V$ can be expressed as the solution to

where, with $M_{u}$ defined in eq. ((ref)),

Here and in the following, we let $V\left(z\right)=\left(V(z,1),....,V(z,D)\right)^{\prime}$ denote the $D\times1$-vector of expected value function and similar for other objects. With this notation, we can represent the fixed-point problem in vector form, $V(z)=\Gamma(V)(z)$, where

Next, to derive the fixed-point problem that $v$ solves, again take conditional expectations on both sides of eq. ((ref)) but now only condition on $Z_{t}$ to obtain

Combining this with eq. ((ref)),

where \[ \bar{\Gamma}(v)(z)=M_{u}\left(\left.\beta\int_{\mathcal{Z}}v\left(z^{\prime}\right)dF_{Z}(dz^{\prime}|z)\right|z\right). \] Under regularity conditions provided below, $\Gamma$ and $\bar{\Gamma}$ are contraction mappings and so $V$ and $v$ are well-defined and unique. The above two transformations of the original problem into the ones for either the integrated or expected value function are particular cases of the general class of transformations analyzed in Ma&Stachurski.

Simulated Bellman operators

As a first step towards a computationally feasible method for solving for either $v$ or $V$, we here develop simulated versions of their two Bellman operators and then introduce the smoothing device. To allow for added flexibility and precision in the implementation and to cover as special case a modified version of Rust's self-approximating solution method, we employ importance sampling: Let $\Phi_{Z}\left(z^{\prime}|z,d\right)$ and $\Phi_{\varepsilon}\left(e|z\right)$ be conditional importance sampling distribution functions as chosen by the researcher. These have to be chosen such that $F_{Z}(\cdot|z,d)$ and $F_{\varepsilon}\left(\cdot|z\right)$ are absolutely continuous w.r.t. $\Phi_{Z}\left(\cdot|z,d\right)$ and $\Phi_{\varepsilon}\left(\cdot|z\right)$, respectively, with Radon-Nikodym derivatives $w_{Z}\left(\cdot|z,d\right)\geq0$ and $w_{\varepsilon}\left(\cdot|z\right)\geq0$ so that

We will throughout assume that eq. ((ref)) is satisfied. In the leading case $dF_{Z}=f_{Z}d\mu_{Z}$ and $d\Phi_{Z}=\phi_{Z}d\mu_{Z}$ for some measure $\mu_{Z}$ in which case $w_{Z}=f_{Z}/\phi_{Z}$ and similar for the sampling of $\varepsilon_{t}$. The above covers the case where $F_{Z}(z^{\prime}|z,d)$ is a continuous distribution (in which case $\mu_{Z}$ is the Lesbesque measure), a discrete distribution (in which case $\mu_{Z}$ is the counting measure) and the mixed case. With discrete finite support, we could in principle compute the exact Bellman equation and its corresponding solution and so would not need to resort to numerical methods. But if the discrete support is large this may still be computationally very demanding and so even in this case the numerical methods developed below may be computationally attractive, c.f. Arcidiaconoetal2013.

Given the chosen importance sampler, we can rewrite $\Gamma(V)(z,d)$ as \[ \Gamma(V)(z,d)=\int_{\mathcal{Z}}\int_{\mathcal{E}}\max_{d^{\prime}\in\mathcal{D}}\left\{ u(s^{\prime},d^{\prime})+\beta V(z^{\prime},d^{\prime})\right\} w\left(s^{\prime}|z,d\right)d\Phi(s^{\prime}|z,d) \] where $s^{\prime}=\left(z^{\prime},e^{\prime}\right)$ and \[ w\left(s^{\prime}|z,d\right)=w_{\varepsilon}\left(e^{\prime}|z^{\prime}\right)w_{Z}\left(z^{\prime}|z,d\right),\:\Phi\left(s^{\prime}|z,d\right)=\Phi_{\varepsilon}(e^{\prime}|z^{\prime})\Phi_{Z}\left(z^{\prime}|z,d\right). \] For any given candidate $V$, we can then approximate this integral by Monte Carlo methods: First generate $N\geq1$ i.i.d. draws, $Z_{i}\left(z,d\right)\sim\Phi_{Z}\left(\cdot|z,d\right)$ and $\varepsilon_{i}\left(z,d\right)\sim\Phi_{\varepsilon}(\cdot|Z_{i}\left(z,d\right))$, $i=1,...,N$, and then compute

where $S_{i}\left(z,d\right)=\left(Z_{i}\left(z,d\right),\varepsilon_{i}\left(z,d\right)\right)$ and

Note here that we normalize the importance weights so that $\sum_{i=1}^{N}w_{N,i}\left(z,d\right)=1$. This is done to ensure that $\Gamma_{N}$ is a contraction mapping on $\mathbb{B}\left(\mathcal{Z}\right)^{D}$. Similarly, we approximate $\bar{\Gamma}\left(v\right)$ by

where again we normalize the weights to ensure $\bar{\Gamma}_{N}$ is a contraction on $\mathbb{B}\left(\mathcal{Z}\right)$, \[ w_{Z,N,i}\left(z,d\right)=\frac{w_{Z}\left(Z_{i}\left(z,d\right)|z,d\right)}{\sum_{i=1}^{N}w_{Z}\left(Z_{i}\left(z,d\right)|z,d\right)},\:w_{\varepsilon,N,i}\left(z,d\right)=\frac{w_{\varepsilon}\left(\varepsilon_{i}\left(z,d\right)|z\right)}{\sum_{i=1}^{N}w_{\varepsilon}\left(\varepsilon_{i}\left(z\right)|z\right)}. \] When $\varepsilon_{t}=\emptyset$, the simulated Bellman operator $\bar{\Gamma}_{N}$ includes as special cases the ones considered in Rust1997curse (who chooses $\Phi_{Z}$ as the uniform distribution on $\mathcal{Z}$) and Pal&Stacurski2013 (who chooses $\Phi_{Z}=F_{Z}$).

Example 1 (continued). Suppose we can compute the integral w.r.t. $\varepsilon_{t}$ analytically. In this case, the following simplified version of the simulated Bellman operator can be employed,

where $G_{\lambda}$ was defined in eq. ((ref)). Importantly, the $\max$-function has been replaced by its smoothed version $G_{\lambda}\left(\cdot\right)$.

If $F_{\varepsilon}(e|z)$ and $F_{Z}(z'|z,d)$ are smooth functions w.r.t. $z$ then $\Gamma(V)(z)$ and $\bar{\Gamma}(v)(z)$ will be smooth functions of $z$ as well. In contrast, the general versions of $\Gamma_{N}(V)(z,d)$ and $\bar{\Gamma}_{N}(v)(z)$ are non-smooth due to the presence of the $\max$-function in their definitions which does not get smoothed for finite $N$. This in turn implies that their corresponding fixed points, $V_{N}\left(z\right)=\Gamma_{N}(V_{N})\left(z\right)$ and $v_{N}\left(z\right)=\bar{\Gamma}_{N}(v_{N})\left(z\right)$, will be non-differentiable w.r.t. the state variables, $z$, and w.r.t. any underlying structural parameters in the model. This is an unattractive feature for two reasons: First, estimation and counterfactuals will be non-smooth problems. Second, the theoretical analysis of $V_{N}$ and $v_{N}$ becomes more complicated.

To resolve this issue, we take inspiration from the additive model in Example 1 and propose to smooth the simulated Bellman operators by replacing the “hard” $\max$-function appearing in eqs. ((ref)) and ((ref)) by its smoothed version $G_{\lambda}(r)$ defined in eq. ((ref)). This yields the following smoothed simulated operators,

where $u(z,\varepsilon_{j}\left(z\right))=\left(u(z,\varepsilon_{j}\left(z,1\right),1),...,u(z,\varepsilon_{j}\left(z,D\right),D)\right)$ and other vector functions are defined similarly. Setting $\lambda=0$ in eqs. ((ref))-((ref)), we recover the original non-smooth versions defined in eqs. ((ref))-((ref)). Thus, the smoothed versions are generalized versions of the original ones. The use of $G_{\lambda}(r)$ in place of $\max_{d\in\mathcal{D}}r\left(d\right)$ generates an additional bias in the approximate solutions, but this can be controlled for by suitable choice of $\lambda$. We now interpret $\lambda>0$ as a smoothing parameter that plays a role similar to that of the bandwidth in kernel regression estimation. Elementary calculations show

so that $G_{\lambda}\left(r\right)\rightarrow\max_{d\in D}r\left(d\right)$, as $\lambda\rightarrow0$, uniformly in $r\in\mathbb{R}^{D}$. Thus, the smoothing entails a bias of order $O_{P}$$\left(\lambda\right)$. We discuss the choice of $\lambda$ in practice in the next section.

In some situations, $G_{\lambda}(r)$ appears in the Bellman operators as an inherent feature of the model specification in which case no smoothing bias will be present. We saw this in Example 1 and it extends to the following class of models: Suppose that $\varepsilon_{t}=\left(\varepsilon_{t}^{\text{\ensuremath{\left(1\right)}}},\varepsilon_{t}^{\left(2\right)}\right)$ with $\varepsilon_{t}^{\text{\ensuremath{\left(1\right)}}}=\left(\varepsilon_{t}^{\text{\ensuremath{\left(1\right)}}}\left(1\right),...,\varepsilon_{t}^{\text{\ensuremath{\left(1\right)}}}\left(D\right)\right)$ are mutually independent extreme value shocks that enter the per-period utility additively,

where as in Example 1 $\lambda>0$ is scale parameter that determines the impact of $\varepsilon_{t}^{\text{\ensuremath{\left(1\right)}}}\left(d\right)$ on the per-period utility. By the same arguments as in Example 1, we find that the expected value function in this case solves $\Gamma_{\lambda}\left(V\right)=V$ where \[ \Gamma_{\lambda}(V)(z)=\int_{\mathcal{Z}}\int_{\mathcal{E}}G_{\lambda}\left(u(z^{\prime},e^{\prime}))+\beta V(z^{\prime})\right)dF_{\varepsilon^{\left(-1\right)}}(e^{\prime}|z)dF_{Z}(z'|z,d), \] and $\Gamma_{N,\lambda}$ in eq. ((ref)) is clearly an unbiased simulated version of $\Gamma_{\lambda}$. Similarly, $\bar{\Gamma}_{N,\lambda}$ is an unbiased estimator of $\bar{\Gamma}_{\lambda}$. To summarize, if the original model of interest contains an additive extreme value term, which is the case in many empirical papers, $G_{\lambda}$ appears as part of the model and so no smoothing bias will be present in our proposed simulated Bellman operators.

The above shows that the smoothing device corresponds to adding structural shocks to the DDP of interest. In earlier work on solving DDPs, researchers have in some cases done the opposite and removed structural errors in order to facilitate the numerical solution of the model; see Lumsdaineetal1992 for one example of this. This was, however, done in the context of discrete state variables with a small number of support points in which case removing continuous structural errors meant that the Bellman operators could be evaluated analytically. Our method is aimed at models where the state variables are either continuous or have a very large discrete support in which case simulations are required in the first place to evaluate the Bellman operator. Once simulations are introduced, there is little computational gains from removing shocks from the model and instead the introduction of smoothing facilitates solving and analyzing the corresponding solution.

Approximate value functions

The smoothed simulated Bellman operators $\Gamma_{N}(V)\left(z,\lambda\right)$ and $\bar{\Gamma}_{N}(v)\left(z,\lambda\right)$ in eqs. ((ref))-((ref)) are functionals that, for given function $V$ and $v$, depend on $\left(z,\lambda\right)$. We will here and in the following treat them as functionals that take given function $V\left(z,\lambda\right)$ and $v\left(z,\lambda\right)$, respectively, and map them into functions of $\left(z,\lambda\right)\in\mathcal{Z}\times\left[0,\bar{\lambda}\right]$ for some $\bar{\lambda}>0$. This simplifies the analysis of the impact of smoothing. In particular, under suitable regularity conditions, $\Gamma_{N}$ and $\bar{\Gamma}_{N}$ are contraction mappings on $\mathbb{B}\left(\mathcal{Z}\times\left[0,\bar{\lambda}\right]\right)^{D}$ and $\mathbb{B}\left(\mathcal{Z}\times\left[0,\bar{\lambda}\right]\right)$, respectively, where $\mathbb{B}\left(\mathcal{Z}\times\left[0,\bar{\lambda}\right]\right)$ denotes the space of bounded functions with domain $\mathcal{Z}\times\left[0,\bar{\lambda}\right]$. Thus, they have unique fixed points $V_{N}\left(z,\lambda\right)$ and $v_{N}\left(z,\lambda\right)$ solving

In practice, we will only solve for the particular value of $\lambda$ as chosen by us, but for the theory it proves helpful to treat the solutions as mappings defined on $\left(z,\lambda\right)\in\mathcal{Z}\times\left(0,\bar{\lambda}\right)$. However, solving these two simulated Bellman equations are not generally feasible since these are infinite-dimensional problems. We here present two ways to reduce the problems to finite-dimensional ones. The first method is a generalized version of the so-called self-approximating method proposed in Rust1997curse while the second one uses projection-based methods as advocated by Pal&Stacurski2013.

Self-approximating method

Rust1997curse proposed to turn the infinite-dimensional problems in eq. ((ref)) into a finite-dimensional ones by choosing the importance sampling to be based on marginal, instead of conditional distributions. In our generalized version this corresponds to restricting $\Phi_{Z}\left(z^{\prime}|z,d\right)=\Phi_{Z}\left(z^{\prime}\right)$ for some marginal distribution $\Phi_{Z}\left(\cdot\right)$ so that the draws $Z_{i}\sim\Phi_{z}\left(\cdot\right)$ and $\varepsilon_{i}\sim\Phi_{\varepsilon}(\cdot|Z_{i})$, $i=1,...,N$ no longer depend on $\left(z,d\right)$. In this case, for a given value of $\lambda\in\left[0,\bar{\lambda}\right]$, the fixed-point problems in eq. ((ref)) reduce to the following two sets of $N$ nonlinear equations,

for $k=1,...,N$, that can be solved for w.r.t. $\left\{ V_{N,\lambda,k}:k=1,...,N\right\} $ and $\left\{ v_{N,\lambda,k}:k=1,...,N\right\} $, respectively. Here, $V_{N,k}=V_{N}\left(Z_{k},\lambda\right)$ and $v_{N,k}=v_{N,\lambda}\left(Z_{k},\lambda\right)$, $k=1,...,N$. Each of the two sets of equations have a unique solution due to the contracting property of $\Gamma_{N,\lambda}$ and $\bar{\Gamma}_{N,\lambda}$. Once, for example, eq. ((ref)) has been solved, the approximate expected value functions can be evaluated at any other value $z$ by \[ V_{N}\left(z,\lambda\right)=\sum_{i=1}^{N}G_{\lambda}\left(u\left(S_{i}\right)+\beta V_{N,i}\right)w_{N,i}\left(z\right). \] Note that $V_{N}\left(z,\lambda\right)$ is a smooth function even if $\lambda=0$ as long as $w_{N,i}\left(z\right)$ is smooth and so smoothing is not needed for this property to hold when marginal samplers are employed. However, without smoothing, the set of equations ((ref)) become non-smooth w.r.t. the variables $\left\{ V_{N,k}:k=1,...,N\right\} $ and so cannot be solved using derivative-based methods. Thus, the numerical implementation of the self-approximating method still benefits from smoothing.

In addition to smoothing, the above self-approximating method differs from Rust's original proposal in two other ways: First, while Rust1997curse solved for the value function $\nu\left(z,\varepsilon\right)$, we here solve for either $V\left(z,\lambda\right)$ or $v\left(z,\lambda\right)$ for a fixed value of $\lambda$. As explained in Section (ref), $V$ and $v$ convey the same information as $\nu$ and at the same time they are of lower dimension in terms of variables and are more smooth, features which facilitate their numerical approximation. Moreover, our formulation allows for the following generalized version of the simulated Bellman equations for $v_{N}$,

where we allow for different number of draws from $\Phi_{\varepsilon}$ and $\Phi_{Z}$. In particular, we can choose $\tilde{N}$ as large as we wish (thereby decreasing the variance of the problem) without increasing the number of variables that need to be solved for ($N$). A similar generalization of the simulated Bellman equations for $V_{N}$ is possible. Second, we here only require that the state dynamics together with the chosen importance sampler satisfy ((ref)); in contrast, Rust1997curse assumed that $S_{t}$ was continuously distributed with compact support and chose as importance sampler the uniform distribution with same support. Thus, our version allows for a broader class of models and samplers.

The self-approximating method may not always work well: First, finding a marginal distribution $\Phi_{Z}\left(\cdot\right)$ so that ((ref)) holds can be difficult in some models. For example, in many specifications with continuous dynamics, the transition density $f_{Z}\left(z^{\prime}|z,d\right)$ of $Z_{t}$ will have singularities, e.g., $\lim_{z^{\prime}\rightarrow z}f_{Z}\left(z^{\prime}|z,d\right)=+\infty$, in which case $w_{Z}\left(z^{\prime}|z,d\right)=f_{Z}\left(z^{\prime}|z,d\right)/\phi_{Z}\left(z^{\prime}\right)$ is not well-defined no matter how we choose $\phi_{Z}\left(z^{\prime}\right)$. And even if ((ref)) does hold, the use of marginal samplers instead of conditional ones will generally lead to a larger variance of the solutions since the “marginal” draws $Z_{1},...,Z_{N}$ do not adapt to the changing shape of $F_{Z}\left(\cdot|z,d\right)$ as a function of $z$. In particular, many of the draws may fall outside of the support of $F_{Z}\left(\cdot|z,d\right)$ and so are “wasted” in which case a large $N$ is required to achieve a reasonable approximation; see Section (ref) for an example of this. This issue tends to become more severe in higher dimensions ($d_{Z}$ is “large”) since the volume of the support shrinks, and so the self-approximating method will generally suffer from a built-in curse-of-dimensionality. This curse-of-dimensionality does not appear in the subclass of models that Rust1997curse focused on where it was assumed that $S_{t}|S_{t-1},d_{t-1}$ has support $\left[0,1\right]^{\dim\left(S_{t}\right)}$ for all values of $S_{t-1},d_{t-1}$.

Finally, given that $V_{N,\lambda}$ and $v_{N,\lambda}$ are solutions to non-linear equations, a large variance in the simulated Bellman operator translates into a large bias as is well-known from non-linear GMM estimators. This can be controlled for by choosing $N$ large. But large $N$ means that numerically solving either ((ref)) or ((ref)) becomes computationally very costly. These issues motivate us to pursue a sieve-based solution strategy.

Sieve-based method

We now return to the general versions of the simulated Bellman operators and so again allow for conditional importance samplers. Let $\mathcal{\bar{V}\subseteq\mathbb{B}\left(\mathcal{Z}\right)}$ be a suitable function space that $z\mapsto v_{N}\left(z;\lambda\right)$ defined in ((ref)) is known to lie in; see below for more details on this. We then choose a finite-dimensional function space (commonly called a sieve in the econometrics literature) $\mathcal{\bar{V}}_{K}=\left\{ v_{K}\left(\cdot;\alpha\right):\mathcal{Z}\mapsto\mathbb{R}|\alpha\in\mathcal{A}_{K}\right\} \subseteq\bar{\mathcal{V}}$, where $\mathcal{A}_{K}\subseteq\mathbb{R}^{K}$ is a parameter set with $K<\infty,$ that provides a good approximation to functions in $\bar{\mathcal{V}}$. Similarly, we let $\mathcal{V}\subseteq\mathbb{B}\left(\mathcal{Z}\right)^{D}$ be a space of $D$-dimensional vector functions that the solution $z\mapsto V_{N}\left(z;\lambda\right)$ to ((ref)) lie in and $\mathcal{V}_{K}=\left\{ V_{K}\left(\cdot;\alpha\right):\mathcal{Z}\mapsto\mathbb{R}^{D}|\alpha\in\mathcal{A}_{K}\right\} \subseteq\mathcal{V}$ be our sieve for this space. Let

be the corresponding projections for given (pseudo-) norms $\left\Vert \cdot\right\Vert _{\mathcal{\bar{V}}}$ and $\left\Vert \cdot\right\Vert _{\mathcal{V}}$ as chosen by us as well. We then approximate $V_{N,\lambda}$ and $v_{N,\lambda}$ by the solutions to the projected Bellman equations,

These are finite-dimensional problems of size $K$. When $K$ is small relative to $N$, which will generally be the case, the above problems are computationally much more tractable compared to the corresponding self-approximating ones. Note here that these projection-based approximations are different from the least-squares approximations that would solve $\min_{V\in\mathcal{V}_{K}}\left\Vert V-\Gamma_{N,\lambda}(V)\right\Vert _{\mathcal{V}}$ and $\min_{v\in\mathcal{\bar{V}}_{K}}\left\Vert v-\bar{\Gamma}_{N,\lambda}(v)\right\Vert _{\mathcal{\bar{V}}}$, respectively. In particular, by suitable choice of the projection operators, $\Pi_{K}\Gamma_{N,\lambda}$ and $\bar{\Pi}_{K}\bar{\Gamma}_{N,\lambda}$ will be contraction mappings w.r.t. $\left\Vert \cdot\right\Vert _{\infty}$ guaranteeing that $\hat{V}_{N,\lambda}$ and $\hat{v}_{N,\lambda}$ exist and are unique. The following discussion focuses on the integrated value function approximation since it carries over with only minor modifications to the one of the expected value function. We discuss their numerical implementation in further detail in the subsection below.

The projection operator $\bar{\Pi}_{K}$ can be thought of as a function approximator with the approximation error being $v-\bar{\Pi}_{K}\left(v\right)$ for a given function $v$. Roughly speaking, the projection-based method approximates $v_{N}$ in ((ref)) by $\hat{v}_{N}=\bar{\Pi}_{K}\left(v_{N}\right)$ which incurs an additional sieve approximation error, $v_{N}-\bar{\Pi}_{K}\left(v_{N}\right)$. The smoothness of $v_{N}$ here proves helpful since many well-known sieves are able to provide good approximations of smooth functions using a low-dimensional space (“small” $K$). Due to these features, our proposed projection-based solutions will generally suffer from quite small additional biases relative to the exact simulated solution. This is in contrast to existing projection-based solution methods, such as the one in Pal&Stacurski2013, that aim at approximating the value function $\nu\left(s\right)$ which is non-differentiable.

The smoothness of $v_{N}$ here help guiding us in choosing the sieve: It allows us to restrict $\mathcal{\bar{V}}$ to a suitable smoothness class and then import existing approximation methods for smooth functions as developed in the literature on numerical methods and nonparametric econometrics. A leading example is the class of linear function approximations where the finite-dimensional function space takes the form of $\mathcal{\bar{V}}_{K}=\left\{ \alpha^{\prime}B_{K}\left(z\right):\alpha\in\mathbb{R}^{K}\right\} $ for a set of basis functions $B_{K}\left(z\right)=\left(b_{1}\left(z\right),...,b_{K}\left(z\right)\right)^{\prime}$. The basis functions can be chosen as, for example, Chebyshev polynomials or B-splines that are able to approximate smooth functions well. However, other non-linear function space are possible such as wavelets, artificial neural networks and shrinkage-type function approximators such as LASSO, where the additional constraints are imposed on $\alpha$; we refer to Chen2007 for a general overview of different function approximators and constrained sieve estimators. We also allow for flexibility in terms of the chosen norms $\left\Vert \cdot\right\Vert _{\mathcal{\bar{V}}}$ with a leading example being $\left\Vert v\right\Vert _{\mathcal{\bar{V}}}=\sum_{i=1}^{M}v^{2}\left(z_{i}\right)$ for a set of design points $z_{1},...,z_{M}\in\mathcal{Z}.$ Very often the $M\geq1$ design points will be chosen in conjunction with the sieve.

The above procedure does not suffer from any of the above mentioned issues of the self-approximating method: We can use conditional importance samplers freely which can be designed to control the variance of the simulated Bellman operators; and the dimension of the problem remains $K$ irrespectively of the number of draws $N$. The main drawback is that unique solutions to eqs. ((ref)) do not necessarily exist for a given choice of $N$ and $K$. A sufficient condition for this to hold is that $\bar{\Pi}_{K}$ is a non-expansive operator w.r.t $\left\Vert \cdot\right\Vert _{\infty}$, that is, $\left\Vert \bar{\Pi}_{K}\left(v_{1}\right)-\bar{\Pi}_{K}\left(v_{2}\right)\right\Vert _{\infty}\leq\left\Vert v_{1}-v_{2}\right\Vert _{\infty}$ for any two functions $v_{1},v_{2}\in\mathcal{\bar{\mathcal{V}}}$, since this translates into $\bar{\Pi}_{K}\Gamma_{N}$ being a contraction mapping. However, while $\bar{\Pi}_{K}$ is non-expansive w.r.t. $\left\Vert \cdot\right\Vert _{\bar{\mathcal{V}}}$ by definition, it is not necessarily non-expansiveness w.r.t $\left\Vert \cdot\right\Vert _{\infty}$. Pal&Stacurski2013 provide some examples of projections that are non-expansive w.r.t. $\left\Vert \cdot\right\Vert _{\infty}$, but these are unfortunately computationally expensive to use in general. But $\bar{\Pi}_{K}$ will generally be close to non-expansive w.r.t $\left\Vert \cdot\right\Vert _{\infty}$ asymptotically as $K\rightarrow\infty$ for a wide range of sieves and pseudo-norms in the sense that \[ \left\Vert \bar{\Pi}_{K}\right\Vert _{op,\infty}:=\sup_{v\in\bar{\mathcal{V}},\left\Vert v\right\Vert =1}\left\Vert \bar{\Pi}_{K}\left(v\right)\right\Vert _{\infty}\leq\sup_{v\in\bar{\mathcal{V}},\left\Vert v\right\Vert =1}\left\Vert \bar{\Pi}_{K}\left(v\right)-v\right\Vert _{\infty}+1, \] where the first term in the last expression will go to zero in great generality as $K\rightarrow\infty$ for many popular sieves (see next subsection for details). Given that $\bar{\Gamma}_{N}$ is a contraction with Lipschitz coefficient $\beta<1$, this in turn implies that $\bar{\Pi}_{K}\bar{\Gamma}_{N}$ will be a contraction mapping for all $K$ large enough. This will be used in our asymptotic analysis of the algorithm. Unfortunately, it is generally not known how large $K$ should be chosen to ensure $\bar{\Pi}_{K}\bar{\Gamma}_{N}$ is a contraction. But in our numerical experiments we did not experience any convergence problems.

Numerical implementation of the two methods

We here discuss in more detail the numerical implementation of the self-approximating and projection-based methods. First, the researcher has to choose the importance sampling distributions, the smoothing parameter $\lambda$ and, in the case of the projection-based method, the function approximation method. Second, given these choices, either eq. ((ref)) or ((ref)) has to be solved for. As before, the discussion here focuses on solving for the integrated value function since most results and arguments for this case carries over with very minor modifications to the expected value function.

Importance sampler

The choices of $\Phi_{Z}$ and $\Phi_{\varepsilon}$ determine the variance of $\bar{\Gamma}_{N}$ and should ideally be tailored to minimize it. In the case of projection-based methods, where we can choose $\Phi_{Z}$ and $\Phi_{\varepsilon}$ as conditional distributions, we can rely on the already existing theory for efficient importance sampling for how to do so; see Chapter 3 in Robert&Casella2004 for an introduction. In our numerical experiments, we did not experiment with different choices and throughout set $\Phi_{Z}=F_{Z}$ and $\Phi_{\varepsilon}=F_{\varepsilon}$.

In the case of the self-approximating method, the choice of $\Phi_{Z}$ is restricted to the class of marginal distributions. Generally, this entails a large variance of the corresponding simulated Bellman operators. It will hold in great generality that $\left(Z_{t},d_{t}\right)$ has a stationary distribution, say, $F_{S}^{*}\left(z,d\right)$. In this case, a suitable choice would be the marginal of this, $\Phi_{Z}\left(z\right)=F_{S}^{*}\left(z,D\right)$. However, the stationary distribution depends on the value function and so is rarely available on closed form; so this strategy requires an initial exploration of the model and its solution. Alternatively, one can try to construct a good approximation of the stationary approximation through a mixture Markov model on the form $\Phi_{Z}\left(z^{\prime}\right)=\sum_{d\in\mathcal{D}}\int\omega_{d}\left(z\right)F_{Z}\left(z^{\prime}|z,d\right)d\mu_{Z}\left(z\right)$ for a set of pre-specified mixture weights $\omega_{d}\left(z\right)\geq0$. In the numerical experiments, we follow Rust1997curse and choose $\Phi_{Z}\left(z\right)$ as the uniform distribution on $\mathcal{Z}$ which we conjecture is far from optimal in many cases, and so more research in this direction is needed.

Smoothing

The use of $G_{\lambda}\left(r\right)$ in place of $\max_{d\in\mathcal{D}}r\left(d\right)$ generally generates an additional bias in the corresponding integrated value function of order $O\left(\lambda\right)$. At the same time, the variance of $v_{N}$ is an increasing function of $\lambda$. Thus, ideally we would like to choose $\lambda$ to balance these two effects. A natural criterion would be to minimize the so-called integrated mean-square-error, $\lambda^{*}=\arg\min_{\lambda\geq0}E\left[\int_{\mathcal{Z}}\left\Vert v_{N,\lambda}\left(z\right)-v\left(z\right)\right\Vert ^{2}dF_{Z}\left(z\right)\right]$, where $F_{Z}\left(z\right)$ is a suitably chosen distribution such as the stationary one of $Z_{t}$. Since $v\left(z\right)$ is unknown and we cannot evaluate the expectations, $\lambda^{*}$ cannot be solved for but cross-validation methods can be used instead. This could in principle be done along the same lines as bandwidth selection for smoothed empirical cdfs, see Bowmanetal1998. However, this is computationally somewhat burdensome. Moreover, in our numerical experiments we found that the quality of the approximate value function was quite insensitive to the choice of $\lambda$ and so in practice we recommend using very little smoothing such as $\lambda=0.01$.

Function approximation

As mentioned earlier, many approximation architectures are available in the literature. In our numerical experiments we focus on the class of linear function approximators where $\bar{\mathcal{V}}_{K}=\left\{ \alpha'B_{K}\left(z\right):\alpha\in\mathbb{R}^{K}\right\} $ for a set of pre-specified basis functions $B_{K}\left(z\right)\in\mathbb{R}^{K}$. For a given set of $M\geq1$ design points in $\mathcal{Z}$, $z_{1},...,z_{M}$, eq. ((ref)) then becomes

The design points may either be random or deterministic and can be chosen relative to the basis functions to ensure that $\bar{\Pi}_{K}$ is easy to compute and provides a good approximation for a broad class of functions. The performance of most function approximations will depend on the smoothness of the function of interest.

A standard smooth function class often considered in approximation theory is the following: For any vector $a=\left(a_{1},...,a_{d_{Z}}\right)\in\mathbb{N}_{0}^{d_{X}}$, let $D^{a}f\left(x\right)=\partial^{\left|a\right|}f\left(x\right)/\left(\partial x_{1}^{\alpha_{1}}\cdots\partial x_{d_{Z}}^{a_{d_{z}}}\right)$, where $\left|a\right|=a_{1}+\cdots+a_{d_{z}}$, be the corresponding partial derivative. For $\alpha>0$, let $\underline{\alpha}\geq0$ be the greatest integer smaller than $\alpha$. For any $\underline{\alpha}$ times differentiable function $f\left(x\right)$, we then define

and let $\boldsymbol{\mathbb{C}}_{r}^{\alpha}\left(\mathcal{X}\right)$ be the space of all $\underline{\alpha}\geq0$ times continuously differentiable functions $f:\mathcal{X}\mapsto\mathbb{R}$ with $\Vert f\Vert_{\alpha,\infty}<r$. Due to smoothing, $v_{N}\in\boldsymbol{\mathbb{C}}_{r}^{\alpha}\left(\mathcal{Z}\right)$, for some $r<\infty$, if $u$ and $F_{Z}$ are sufficiently smooth (see Theorem (ref)). We can therefore import existing results for approximation methods for functions in $\boldsymbol{\mathbb{C}}_{r}^{\alpha}\left(\mathcal{Z}\right)$:

As can be seen from the above example, standard polynomial tensor product approximations suffer from the well-known computational curse of dimensionality: To reach a given level of error tolerance, the total number of basis functions $K$ has to grow exponentially as $d_{Z}$ increases. This issue can be partially resolved by using more advanced function approximation methods:

As noted earlier, there is no guarantee that a given function approximator is non-expansive. But this can, in principle, be examined numerically for a given choice of $\bar{\Pi}_{K}$. For the least-squares projection, this amounts to solving, for a given choice of basis functions and grid points,

When $M$ and/or $\dim\mathcal{Z}$ is large this may be computationally demanding and instead one can obtain a lower bound by restricting $z$ to only take values on the chosen set of grid points: With $\mathbf{B}_{K,M}\in\mathbb{R}^{K\times M}$containing the basis functions evaluated at the grid points, we can represent $\bar{\Pi}_{K}$ when only evaluated at chosen grid points $z_{1},....,z_{M}$ in terms of \[ \mathbf{P}_{K,M}=\mathbf{B}_{K,M}'\left[\mathbf{B}_{K,M}\mathbf{B}_{K,M}'\right]^{-1}\mathbf{B}_{K,M}\in\mathbb{R}^{M\times M}. \] In particular, it is easily checked that with the supremum in ((ref)) being only taken over $z\in\left\{ z_{1},...,z_{M}\right\} $, $\left\Vert \bar{\Pi}_{K}\right\Vert _{op,\infty}=\left\Vert \mathbf{\mathbf{P}_{K,M}}\right\Vert _{op,\infty}$. Furthermore, $\left\Vert \mathbf{P}_{K,M}\right\Vert _{op,\infty}\leq1$ if and only if \[ \max_{i=1,...,M}\sum_{j=1}^{M}\left|p_{ij}\right|\leq1, \] where $p_{ij}$ is the $\left(i,j\right)$th element of $\mathbf{P}$, c.f. Lizotte2011.

Solving for the approximate value functions

Computing the simulated self-approximating solution or the projection-based one can be done using three different numerical algorithms: Successive approximation (SA), Newton-Kantorovich (NK), or a combination of the two. The latter corresponds to the hybrid solution method proposed in Rust1988Siam. We here discuss the implementation of these algorithms with focus on the sieve-based approximation of $v$; the implementations of the sieve approximation of $V$ and the self-approximating solutions of either of the two follow along the same lines. The main difference between solving for $V$ or $v$ is that the latter involves smaller computational burden since it is a scalar function while the former is a $D$-dimensional vector function.

SA utilizes that (for $K$ chosen large enough), $\bar{\Pi}_{K}\bar{\Gamma}_{N,\lambda}$, is a contraction mapping which guarantees that the following algorithm will converge towards the solution to ((ref)),

for $k=1,2,...$, given some initial guess $\hat{v}_{N}^{\left(0\right)}$. In the leading case of ((ref)), this can be expressed as a sequence of least-squares problems that are easily computed: $\hat{v}_{N}^{\left(k\right)}\left(z\right)=\hat{\alpha}_{k}^{\prime}B_{K}\left(z\right)$ where \[ \hat{\alpha}_{k}=\left[\sum_{i=1}^{M}B_{K}\left(z_{i}\right)B_{K}\left(z_{i}\right)^{\prime}\right]^{-1}\sum_{i=1}^{M}B_{K}\left(z_{i}\right)\bar{\Gamma}_{N,\lambda}(\hat{\alpha}_{k-1}^{\prime}B_{K})\left(z_{i}\right)^{\prime}\in\mathbb{R}^{K}, \] for $k=1,2,...$, given some initial guess $\hat{\alpha}_{0}$. In the case where $\varepsilon_{t}=\emptyset$ or when the model is on the form eq. ((ref)) with $\varepsilon_{t}^{\text{(1}}$ being extreme-valued distributed and $\varepsilon_{t}^{\left(-1\right)}=\emptyset$, \[ \bar{\Gamma}_{N}(\alpha^{\prime}B_{K})(z_{i};\lambda)=G_{\lambda}\left(u\left(z\right)+\beta\alpha^{\prime}\sum_{j=1}^{N}B_{K}\left(Z_{j}\left(z_{i}\right)\right)w_{Z,N,j}\left(z_{i}\right)\right), \] and so $\sum_{j=1}^{N}B_{K}\left(Z_{j}\left(z_{i}\right)\right)w_{z,N,j}\left(z_{i}\right)$, $i=1,...,M,$ only need to be computed once and then recycled in each iteration; in contrast, the simulated averages appearing in $\Gamma_{N}(\hat{\alpha}_{k-1}'B_{K})\left(z_{i}\right)$, $i=1,...,M$, have to be recomputed in each step of the SA algorithm. Thus, in this special case, it is faster to (approximately) solve for $v_{N,\lambda}$ instead of $V_{N,\lambda}$. While SA is guaranteed to converge globally when $\bar{\Pi}_{K}\bar{\Gamma}_{N}$ is a contraction, the rate of convergence will be slow with the error vanishing at rate $\beta^{k}$,

To speed up convergence, we therefore follow Rust1988Siam and combine SA with NK iterations since NK converges with a quadratic rate once a given guess of the value function is close enough to the fixed point. Moreover, in situations where $\bar{\Pi}_{K}\bar{\Gamma}_{N}$ is expansive, NK is still guaranteed to converge locally. Since both the self-approximating and sieve-based methods solve finite-dimensional problems, the NK algorithm for these are equivalent to Newton's method. First consider the sieve-based method where we focus on the least-squares projection as given in ((ref)). We are then seeking $\hat{\alpha}$ solving the following $K$ equations, \[ \bar{S}_{N,K}\left(\alpha;\lambda\right)=0, \] with \[ \bar{S}_{N,K}\left(\alpha;\lambda\right)=\alpha-\left[\sum_{i=1}^{M}B_{K}\left(z_{i}\right)B_{K}\left(z_{i}\right)^{\prime}\right]^{-1}\sum_{i=1}^{M}B_{K}\left(z_{i}\right)\bar{\Gamma}_{N}\left(\alpha^{\prime}B_{K}\right)\left(z_{i};\lambda\right). \] The corresponding derivatives of the left-hand side as a function w.r.t. $\alpha$ can be expressed in terms of the Hadamard differential of $\bar{\Gamma}_{N}$ w.r.t. $v$,

where $dv:\mathcal{Z}\times\left[0,\bar{\lambda}\right]\mapsto\mathbb{R}$ is the direction and

The partial derivatives of $\bar{S}_{N,K}\left(\alpha;\lambda\right)$ then becomes \[ \bar{H}_{N,K}\left(\alpha;\lambda\right)=I_{K}-\left[\sum_{i=1}^{M}B_{K}\left(z_{i}\right)B_{K}\left(z_{i}\right)^{\prime}\right]^{-1}\sum_{i=1}^{M}B_{K}\left(z_{i}\right)\nabla\bar{\Gamma}_{N}\left(\alpha^{\prime}B_{K}\right)\left[B_{K}\right]\left(z_{i};\lambda\right)^{\prime}\in\mathbb{R}^{K\times K}. \] With these definitions, the NK algorithm takes the form \[ \hat{\alpha}_{k}=\hat{\alpha}_{k-1}-\bar{H}_{N,K}^{-1}\left(\hat{\alpha}_{k-1};\lambda\right)\bar{S}_{N,K}\left(\hat{\alpha}_{k-1};\lambda\right). \]

The NK algorithm for the self-approximating method is on the same form, except that we now solve directly for the value function at the $N$ draws. With slight abuse of notation, let $v_{N}=\left\{ v_{N,\lambda}\left(Z_{k};\lambda\right):k=1,...,N\right\} $ be the vector of integrated values across the set of draws solving $\bar{S}_{N}\left(v_{N};\lambda\right)=0$ where

for $k=1,...,N$. The corresponding derivatives is $\bar{H}_{N}\left(\alpha;\lambda\right)=\left(\bar{H}_{N,1}\left(\alpha;\lambda\right),...,\bar{H}_{N,N}\left(\alpha;\lambda\right)\right)^{\prime}\in\mathbb{R}^{N\times N}$ where, with $\mathbf{1}_{N}=\left(1,...,1\right)^{\prime}\in\mathbb{R}^{N}$, \[ \bar{H}_{N,k}\left(\alpha;\lambda\right)=I_{N}-\nabla\bar{\Gamma}_{N}\left(v_{N}\right)\left[\mathbf{1}_{N}\right]\left(z_{k};\lambda\right)\in\mathbb{R}^{N}. \] Finally, we note that the NK algorithm for the expected value function takes a similar form with the functional differential of $\Gamma_{N}$ given by

where $dV\left(z\right)=\left(dV\left(z,1\right),...,dV\left(z,D\right)\right)^{\prime}.$

Comparing the NK algorithm for the self-approximating and the sieve-based method, we note that the former involves inverting a $N\times N$-matrix while the latter a $K\times K$-matrix. As pointed out earlier, the self-approximating method generally needs $N$ to be chosen quite large to achieve a precise simulated version of the Bellman operator, in particular in higher dimensions, and so the NK algorithm for this method may become numerically infeasible in some cases. While the projection-based method also suffers from a curse of dimensionality, since the number of basis functions, $K$, has to be quite large in higher dimensions to achieve a reasonable approximation, it is less severe and is implementable for higher-dimensional models. If more advanced function approximation methods are employed, even better performance can be achieved.

Theory

We here develop an asymptotic theory for the self-approximating and sieve-based methods. We first establish some important properties of the smoothed simulated Bellman operators and their exact solutions, $v_{N}$ and $V_{N}$ defined in ((ref)). These are then used in the asymptotic analysis of the self-approximating solution method and the sieve-based one. This analysis will rely on two general results for estimated solutions to fixed point problems as stated in Theorems (ref) and (ref) in the appendix. The asymptotic analysis will mostly focus on $V_{N}$ and $\hat{V}_{N}$ since our results for these easily translate into similar results for the approximate integrated value function. For example, $v_{N}(z;\lambda)=M_{N,u}\left(\beta V_{N}\left(z;\lambda\right)|z;\lambda\right)$, where \[ M_{N,u}\left(r|z;\lambda\right)=\sum_{j=1}^{N}G_{\lambda}\left(u(z,\varepsilon_{j}\left(z\right))+r\right)w_{\varepsilon,N,j}\left(z\right), \] and so the asymptotic results for $V_{N}$ in conjunction with the functional Delta method can be used to obtain similar results for $v_{N}$.

Without loss of generality, we assume that the draws can be written as

for some i.i.d. draws $U_{i}\sim P_{U}$, $i=1,...,N$, and some functions $\psi_{Z}$ and $\psi_{\varepsilon}$. We then define, with $\psi=\left(\psi_{Z},\psi_{\varepsilon}\right)$ and for any given function $V\left(z,\lambda\right)$,

so that

where expectations are taking over $U\sim P_{U}$, and

where $w_{\psi,N}\left(U_{i};z,d\right)=w_{\psi}\left(U_{i};z,d\right)/\sum_{j=1}^{N}w_{\psi}\left(U_{j};z,d\right).$ Here, and in the following, we let $V_{0}$$\left(z,\lambda\right)$ denote the exact solution to

As explained earlier, we here define the two operators to take a given function $V\left(z,\lambda\right)$, $\left(z,\lambda\right)\in\mathcal{Z}\times\left[0,\bar{\lambda}\right]$, for some given $\bar{\lambda}>0$, and map them into another function with domain $\mathcal{Z}\times\left[0,\bar{\lambda}\right]$. In particular, $V_{0}$$\left(z,0\right)$ and $V_{N}\left(z,0\right)$ are the non-smoothed ($\lambda=0$) exact and simulated solutions, respectively. We then impose the following regularity conditions on the model and chosen importance sampler, where we recall the function norm defined in eq. ((ref)) and the function set $\boldsymbol{\mathbb{C}}_{r}^{\alpha}\left(\mathcal{Z}\right)$ defined below this equation:

Assumption (ref) share some similarities with the regularity conditions found in Rust1988Siam who considered an additive version of our general model. Importantly, we only require that $Z_{t}$ has bounded support while $\varepsilon_{t}$ can have potentially unbounded support. This is in contrast to Pal&Stacurski2013 and Rust1997curse who require both components to be bounded. We conjecture that the subsequent results can be generalized to also hold in the case of $\mathcal{Z}$ unbounded but then our conditions and arguments would have to be changed. For example, the existence of unique fixed points would have to be verified in a function space equipped with a weighted $\sup$-norm and with additional moment conditions on $Z_{t}$, see, e.g., Norets2010differentiability. Similarly, our empirical process results would need to be established using bracketing conditions with weighted norms and additional moment conditions, see, e.g., Section 2.10.4 in VW. Similar to the results in Rust1988Siam, Assumption (ref) implies $\Gamma(V)\in\boldsymbol{\mathbb{B}}\left(\mathcal{Z}\times\left(0,\bar{\lambda}\right)\right)^{D}$ for all $V\in\mathcal{\boldsymbol{\mathbb{B}}}\left(\mathcal{Z}\times\left(0,\bar{\lambda}\right)\right)^{D}$, see below. This particular result actually holds under the weaker requirement that $E\left[\Vert u_{\psi}\left(U;\cdot\right)\Vert_{0,\infty}\Vert w_{\psi}\left(U;\cdot\right)\Vert_{,\infty}\right]<\infty$ but the existence of the second order moment is needed for the subsequent asymptotic analysis of $V_{N}\left(z;\lambda\right)$ and so we impose this restriction throughout.

Assumption (ref) impose smoothness conditions on the model and the chosen samplers in terms of the state variables $Z_{t}$. These conditions imply that the expected and integrated value functions will be smooth too. Note here that the degree of smoothness $\alpha$ is left unrestricted at this stage and so the functions are not required to be differentiable, merely Lipschitz. While the focus is on models with continuous state variables our theory also covers models with discrete state space. In this case, we can dispense of Assumption (ref) and instead rely on the more general Theorem (ref) which implies that the subsequent results still go through when $Z$$_{t}$ is discrete.

We first establish existence and uniqueness of the (generally) infeasible simulated solutions and show that they inherit the smoothness properties of $z\mapsto u_{\psi}\left(U;\cdot\right)$ and $z\mapsto w_{\psi}\left(U;z\right)$. This feature of the approximate solutions is important for two reasons: First, it allows us to show uniform convergence of certain functionals as part of our proof of weak convergence. Second, we can control the approximation error due to the use of sieves later on.

The above is a fixed $N$ result with the bound on $V_{N}$, $r_{N}$, being random since it depends on the particular set of draws. We derive a deterministic bound on $r_{N}$ as $N\rightarrow\infty$ below. The condition $\inf_{z\in\mathcal{Z}}\sum_{i=1}^{N}w_{\psi}\left(U_{i};z\right)>0$ will hold with probability approaching 1 (w.p.a.1) as $N\rightarrow\infty$ and so can be dropped in our asymptotic analysis. Next, we analyze the effect of smoothing on the exact and simulated value function:

This shows that the smoothing can be controlled for by suitable choice of $\lambda$ both asymptotically ($N=+\infty$) and for any finite number of simulations ($N<\infty$). Also note that the above result holds independently of the smoothness properties of the unsmoothed exact and simulated solutions.

Self-approximating method

In this section, we provide an analysis of the smoothed fixed point, $V_{N}$, to $\Gamma_{N}$ defined in eq. ((ref)) thereby allowing for general importance samplers. As a special case, we obtain an asymptotic theory for the self-approximating solution method (where $\Phi_{z}\left(z^{\prime}|z,d\right)$ is restricted to be marginal distribution). The results for $V_{N}$ will then in turn be used in the analysis of the corresponding sieve-based methods in the next section.

The first part of this theorem is similar to results found in Rust1997curse and Pal&Stacurski2013 who also show $\sqrt{N}$-convergence of their value function approximation. Importantly, the convergence result holds uniformly over the smoothing parameter $\lambda$ and so there is no first-order effect from smoothing if $\lambda$ vanishes sufficiently fast. Specifically, for any sequence $\lambda_{N}$ satisfying $\sqrt{N}\lambda_{N}\rightarrow0$, Theorems (ref) and (ref) yield $\sup_{z\in\mathcal{Z}}\left\Vert V_{N}\left(z,\lambda_{N}\right)-V\left(z,0\right)\right\Vert =O_{p}(1/\sqrt{N})$. This is similar to convergence of smoothed empirical cdf where the indicator function is replaced by a smoothed version; this also does not affect the convergence rate as long as the smoothing bias is controlled for. The second part of the theorem appears to be a new result and shows that if the problem is smooth enough, the first-order partial derivatives of $V_{N}\left(z,\lambda\right)$ also converge uniformly over $z$ with rate $\sqrt{N}/\lambda$. Since we need $\lambda\rightarrow0$ to kill the smoothing bias, this could seem to imply that the first-order derivatives converge with slower than $\sqrt{N}$-rate. However, we conjecture that the derived rate is not sharp and that $\sqrt{N}$-convergence does actually hold. The proof of this appears to require a more delicate and refined arguments, however, and so we leave this for future research.

The above result is then in turn used to derive the asymptotic distribution of $V_{N}\left(z,\lambda\right)$ uniformly in $\left(z,\lambda\right)\in\mathcal{Z}\times\left(0,\bar{\lambda}\right)$. Here, the smoothing proves important since it allows us to generalize the standard arguments used in the analysis of finite-dimensional extremum estimators to our setting: We first expand the “first-order condition”, $V_{N}-\Gamma_{N}(V_{N})=0$, around $V_{0}=\Gamma(V_{0})$ to obtain, with $\nabla\Gamma_{N}$ defined in ((ref)),

where the rate of the remainder term follows from Theorem (ref). Next, employing empirical process theory, we show that $\sqrt{N}\left\{ \Gamma(V_{0})-\Gamma_{N}(V_{0})\right\} \rightsquigarrow\mathbb{G}$ in $\boldsymbol{\mathbb{B}}\left(\mathcal{Z}\times\left(0,\bar{\lambda}\right)\right)^{D}$ for a Gaussian process $\mathbb{G}\left(z,\lambda\right)$ with covariance kernel

Finally, we show that $\nabla\Gamma_{N}(V_{0})\left[dV\right]\rightarrow^{P}\nabla\Gamma(V_{0})\left[dV\right]$ uniformly over $dV\in\boldsymbol{\mathbb{C}}_{2r}^{1}\left(\mathcal{Z}\times\left[0,\bar{\lambda}\right]\right)^{D}$ with $r<\infty$ given in Theorem (ref). Since $V_{N}-V_{0}\boldsymbol{\in\mathbb{C}}_{2r}^{1}\left(\mathcal{Z}\times\left[0,\bar{\lambda}\right]\right)^{D}$, we conclude that:

The above result implies, for example, that $\sqrt{N}\left\{ V_{N}\left(z,\lambda\right)-V_{0}\left(z,\lambda\right)\right\} \rightarrow^{d}N\left(,\Omega_{V}\left(z,\lambda,z,\lambda\right)/N\right)$ as $N\rightarrow\infty$ for any given $\left(z,\lambda\right)$, where \[ \Omega_{V}\left(z_{1},\lambda_{1},z_{2},\lambda_{2}\right)=\int\int r^{*}\left(z_{1}^{\prime},\lambda_{1}^{\prime}|z_{1},\lambda_{1}\right)\Omega\left(z_{1}^{\prime},\lambda_{1}^{\prime},z_{2},\lambda_{2}\right)r^{*}\left(z_{2}^{\prime},\lambda_{2}^{\prime}|z_{2},\lambda_{2}\right)^{\prime}d\left(z_{1}^{\prime},\lambda_{1}^{\prime}\right)d\left(z_{2}^{\prime},\lambda_{2}^{\prime}\right), \] where $r^{*}$ is the Riesz representer of $dV\mapsto\left\{ I-\nabla\Gamma(V_{0})\right\} ^{-1}\left[dV\right]\left(z,\lambda\right)$. Thus, it allows us to construct (pointwise or uniform) confidence bands for the expected value function. We expect that the result will also be useful in analyzing the impact of value function approximation when used in estimation. This could be done by combining the above weak convergence result with, e.g., the results for approximate estimators found in Kristensen&Salanie2017.

The proof of Theorem (ref) proceeds by verifying the two high-level conditions of the “master” Theorem (ref) where the same weak convergence result is obtained under more general conditions. Theorem (ref) allows us to replace the smoothness conditions in Assumption (ref) with some other conditions implying that $V_{N}-V_{0}$ is situated in a function set with finite entropy. One example would be to impose restrictions on $u$ and $w$ so that the value function and its estimator are both monotone functions, c.f. Pal&Stacurski2013, in which case we could then appeal to Theorem 2.7.5 in VW to obtain the results of Theorem (ref).

We conjecture that a similar weak convergence result will hold for the non-smoothed value function approximation ($\lambda=0$). However, the proof of such a result would require different arguments and seemingly stronger assumptions. In particular, the current proof only requires the empirical process $\left(z,\lambda\right)\mapsto\Gamma_{N}(V_{0})\left(z,\lambda\right)$ to converge weakly. To allow for non-smooth value function approximation, we conjecture that we would now need to show that the empirical process $\left(V,z\right)\mapsto\Gamma_{N}(V)\left(z,0\right)$ converges weakly over a suitable function set that the estimated non-smooth solution, $V_{N}$, would be situated in. For this to hold, the uniform entropy of the function set would need to be finite. Standard choices of function sets are smooth classes, but $V_{N}$ and its limit $V_{0}$ are both non-smooth now and so the proof appears to be rather delicate.

Finally, for a complete analysis that takes into account the smoothing bias, we state the following corollary to Theorem (ref): For any $\lambda_{N}\rightarrow0$ such that $\lambda_{N}\sqrt{N}\rightarrow0$, $\sqrt{N}\{V_{N}\left(\cdot;\lambda_{N}\right)-V_{0}\left(\cdot;0\right)\}\rightsquigarrow\mathbb{G}_{V}\left(\cdot,0\right).$

Sieve-based approximation of value functions

We now proceed to analyze the asymptotic properties of the sieve-based approximate value function, $\hat{V}_{N}$. To this end, we use the following decomposition of the over-all error,

where the second term converges weakly towards a Gaussian process, c.f. Theorem (ref). What remains is to control the first term which is due to the sieve approximation; this is done by imposing the following high-level assumption on the projection operator when applied to a function set $\mathcal{V}$ which is chosen so that $V_{N}\in\mathcal{V}$ w.p.a.1.:

This is a high-level condition that requires the chosen function approximation method to have a uniform error rate over the function class $\boldsymbol{\mathbb{C}}_{r}^{\alpha}\left(\mathcal{Z}\right)$ which we know $z\mapsto V_{0}\left(z,\lambda\right)$ belongs uniformly in $\lambda$ under Assumption (ref). As discussed earlier, one could replace Assumption (ref) with other regularity conditions that ensure $z\mapsto V_{0}\left(z,\lambda\right)$ is sufficiently regular (e.g., monotonic) in which case $\boldsymbol{\mathbb{C}}_{r}^{\alpha}\left(\mathcal{Z}\right)$ in Assumption (ref) should be modified accordingly. Assumption (ref) is satisfied for standard polynomial approximators with $\rho_{K}=\log\left(K\right)/K^{\left(s+1\right)/d}$, c.f. Section (ref). Compared to results on sieve approximations of value functions found elsewhere in the literature, our rate is better since we are here seeking to approximate the expected value function that is situated in $\boldsymbol{\mathbb{C}}_{r}^{\alpha}\left(\mathcal{Z}\right)$. In contrast, sieve-based approximations developed in other papers, such as Munos&Szepesvari2008 and Pal&Stacurski2013, try to approximate the value function which is at most Lipschitz and for such functions the approximation error will be larger in general. In the case of $\mathcal{Z}$ being finite, we have $\sup_{V\in\mathcal{V}}\left\Vert \Pi_{K}\left(V\right)-V\right\Vert _{\infty}=0$ for $K>\left|\mathcal{Z}\right|$ under great generality and so there will be no asymptotic bias component due to sieve approximations in this case.

The second part of Theorem (ref) together with the fact that $\Gamma_{N,\lambda}(V_{\lambda})-\Gamma_{\lambda}(V_{\lambda})=O_{P}\left(1/\sqrt{N}\right)$, c.f. Proof of Theorem (ref), now yield the following result:

The discussions following Theorems (ref) and (ref) carry over to the above result. In particular, the rate result still goes through when no smoothing is employed ($\lambda=0$) but the current proof of the asymptotic distribution result requires smoothing ($\lambda>0$). Compared to the rate results for $V_{N,\lambda}$, the projection-based method suffers from an additional error due to the sieve approximation, $O_{P}\left(\rho_{K}\right)$. This can be interpreted as a bias term, while $O_{p}(1/\sqrt{N})$ is its variance component which is shared with $V_{N,\lambda}$. The requirement that $\sqrt{N}\rho_{K}\rightarrow0$ is used to kill the sieve bias term so that $\hat{V}_{N}$ is centered around $V_{0}$.

The above result provides a refinement over existing results where a precise rate for the bias is not available; see, e.g., Lemma 5.2 in Pal&Stacurski2013. It shows that there is an inherent computational curse-of-dimensionality built into our projection-based value function approximation when polynomial interpolation is employed: In high-dimensional models, a large number of basis functions are needed which in turn increases the computational effort. In the case of polynomial approximations, the rate condition becomes $\sqrt{N}\log\left(K\right)/K^{\left(s+1\right)/d}\rightarrow0$ and so, as $d$ increases, we need $K$ to increase faster with $N$ to kill the sieve bias component. However, also note that $K$ as no first-order effect on the variance and so there is no bias-variance trade-off present. In particular, we can let $K$ increase with $N$ as fast as we wish and so our procedure should in principle also work for models with high-dimensional state space. That is, we can achieve $\sqrt{N}$-rate regardless of the dimension of the problem and so our method does not suffer from any statistical curse-of-dimensionality. However, this requires choosing $K$ large enough in order to control the sieve approximation bias which will increase computation time as the dimension grows. Thus, there is a potential computational curse-of-dimensionality.

Numerical results

In this section we examine the numerical performance of the proposed solution algorithms with focus on how the theoretical results derived in the previous sections translate into practice and how different features of model and implementation affect their performances.

We focus exclusively on approximating the integrated value function, $v\left(z\right)$, and measure the performance of a given approximate solution, say, $\tilde{v}\left(z\right)$ in terms of its pointwise bias, variance and mean-square error (MSE) defined as $Bias\left(z\right):=E[\tilde{v}\left(z\right)]-v\left(z\right)$, $Var\left(z\right):=Var\left(\tilde{v}\left(z\right)\right)=E\left[(\tilde{v}\left(z\right)-E[\tilde{v}\left(z\right)])^{2}\right]$ and $MSE\left(z\right)=Bias^{2}\left(z\right)+Var\left(z\right)$, respectively. As overall measures we use uniform bias, variance and MSE, $\left\Vert Bias\right\Vert $$_{\infty}=\sup_{z\in\mathcal{Z}}\left|Bias\left(z\right)\right|$, $\left\Vert Var\right\Vert $$_{\infty}=\sup_{z\in\mathcal{Z}}\left|Var\left(z\right)\right|$ and $\left\Vert MSE\right\Vert $$_{\infty}=\sup_{z\in\mathcal{Z}}\left|MSE\left(z\right)\right|$. Given that the exact solution $v\left(z\right)$ is unknown, we replace this by a very precise approximate solution computed in the following way: First, instead of using simulations in the computation of the Bellman operator, we utilize that the state transitions follow a Beta distribution in the chosen model (see below) and so we can use nodes and weights based on Jacobi polynomials to compute it using numerical integration. We then implement the sieve method using $K=60$ Chebyshev polynomials and $N=60$ sets of quadrature nodes and weights. The ``exact'' solution was computed by successive approximation until a contraction tolerance of machine precision was reached. We approximate the point-wise bias and variance of a given method through $S\geq1$ independent replications of it: Let $\tilde{v}_{1}\left(z\right),....,\tilde{v}_{S}\left(z\right)$ be the solutions obtained across the $S$ replications, where $S$ generally was chosen to 2,000. We then approximate the mean by $\hat{E}[\tilde{v}\left(z\right)]=\frac{1}{S}\sum_{s=1}^{S}\tilde{v}_{s}\left(z\right)$ which in turn is used to obtain the following pointwise bias and variance estimates, $\hat{Bias}\left(z\right)=\hat{E}[\tilde{v}\left(z\right)]-v_{0}\left(z\right)$, and $\hat{Var}\left(z\right)=\frac{1}{S}\sum_{s=1}^{S}(\tilde{v}_{s}\left(z\right)-\hat{E}[\tilde{v}\left(z\right)])^{2}$. Based on these, we approximate $\left\Vert Bias\right\Vert $$_{\infty}$ and $\left\Vert Var\right\Vert $$_{\infty}$ by the maximum pointwise biases and variances over a uniform grid over $[0,1000]$ of size 500 in the univariate case and a uniform grid over $[0,1000]^{2}$ of size 250.

To implement the sieve-based method, we need to choose the sieve space used in constructing $\Pi_{K}$. We here focus on Chebyshev basis functions and B-Splines as discussed in Section (ref)\footnote{For more details on their implementation, see appendix (ref) .}.

A model of optimal replacement

To provide a test bed for comparison of the sieve-based approximation method, we use the well-known engine replacement model by Rust1987. Rust's model has become the basic framework for modeling dynamic discrete-choice problems and has been extensively used in other studies to evaluate the performance of alternative solution algorithms and estimators. While the model and its solution is well described in many papers, for completeness we briefly describe our variation of it below.

We consider the optimal replacement of a durable asset (such as a bus engine) whose controlled state $Z_{t}\in\mathbb{R}_{+}$ is summarized by the accumulated utilization (mileage) since last replacement. In each period, the decision maker faces the binary decision $d_{t}\in\mathcal{D=\mathrm{\left\{ 0,1\right\} }}$ whether to keep ($d_{t}=0)$ or replace $(d_{t}=1$) the durable asset with a fixed replacement cost $RC>0$. If the asset is replaced, accumulated usage $Z_{t}$ regenerates to zero. The maintenance/operating costs are assumed to be linear in usage $Z_{t}$, $c(Z_{t})=\theta_{c}\cdot0.001\cdot Z_{t}$. The state and decision dependent per period utility is then given by $\bar{u}(Z_{t},d_{t})+\varepsilon_{t}(d_{t})$ where $\bar{u}(Z_{t},d_{t})=\left(RC+c(0)\right)I\left\{ d_{t}=0\right\} +c\left(Z_{t}\right)I\left\{ d_{t}=1\right\} $ and the utility shocks $\varepsilon_{t}=\left(\varepsilon_{t}(0),\varepsilon_{t}(1)\right)$ are i.i.d. extreme value and fully independent of $Z_{t}$. This specification is a special case of Example 1 with $\lambda=\sigma_{\varepsilon}$ and so the simulated Bellman operator takes the form ((ref)) where $G_{\sigma_{\varepsilon}}\left(\cdot\right)$ appears as part of the model. Thus, there is no smoothing bias present in the baseline model. In Section (ref), we investigate the effect of smoothing by pretending that we are not able to integrate out $\varepsilon_{t}$ analytically in the baseline model and instead we simulate both $Z_{t}$ and $\varepsilon_{t}$ and then include our smoothing device in the computation of the simulated Bellman operator.

We assume that $Z_{t}$ (in absence of the replacement decision) follows a mixture of a discrete distribution with a probability mass $\pi>0$ at zero and a linearly transformed Beta distribution with shape parameters $a$ and $b$ and scale parameter $\sigma_{\varepsilon}>0$. Thus,

where $F_{+}\left(z^{\prime}|z,d\right)$ has density $f_{+}\left(z^{\prime}|z,d\right)=f_{\beta}\left((z^{\prime}-z)/\sigma_{Z};a,b\right)/\sigma_{Z}$, $\pi>0$ is the probability of no usage and $f_{\beta}(x;a,b)$ is the probability density function of the Beta distribution with shape parameters $a,b$. Note here that it has bounded support $\left(0,\sigma_{Z}\right)$ so that $f_{+}\left(z^{\prime}|z,d\right)=0$ for $z^{\prime}<z$ or $z^{\prime}-z>\sigma_{Z}$. This is in line with the discretized model in the original formulation in Rust1987 where monthly mileage were only allowed to take a few discrete values and monthly mileage is naturally bounded above and below (busses never drives backwards and there are limits how far a bus can drive within a month). We introduce probability mass $\pi$ at $z^{\prime}=z$ to allow for the possibility that the asset is not used in a given period and thereby can end in the same state with positive probability when $\pi>0$. As explained below, this feature turns out to be quite important for the applicability of the self-approximating method of Rust (1997). Note that the support of $Z_{t}$ is unbounded (the positive half line) and therefore the theory does not apply directly, since we throughout assumed bounded support. However, we expect that the theory extends to the unbounded case after suitable modifications, c.f. discussion following Assumptions (ref)-(ref).

In the numerical illustrations below we use the following set of benchmark parameter values unless otherwise specified: We set replacement cost to $RC=10$ and the cost function parameter to $\theta_{c}=2$ so that $RC$ is 5 times as large as $c(1000)$. This implies a large variation in the probability of replacement over $Z_{t}$ compared to Rust1987 and a more curved value function. The parameters indexing the transition density $f_{+}\left(z^{\prime}|z,d\right)$ are set to $\sigma_{Z}=15$, $a=2$, $b=5$ and $\pi=10^{-10}$ as default. This implies a quite sparse transition density, which is similar to the fitted model in Rust1987. In Figure (ref) we plot the corresponding “exact” solution as described earlier. Importantly, since the transition density is an analytic function the value function is also analytic and so well-approximated by polynomial interpolation methods.

Numerical implementation of simulated Bellman operators

The simulated Bellman operators in ((ref)) and ((ref)) require the user to choose an importance sampling distribution. For the self-approximating solution method we need to choose a marginal sampler, $d\Phi_{Z}\left(z^{\prime}|z,d\right)=\phi_{Z}\left(z^{\prime}\right)dz^{\prime}$. We follow Rust1997curse and choose $\phi_{Z}\left(z^{\prime}\right)=I\left\{ 0<z^{\prime}<z^{\max}\right\} $ as a uniform density with support support $\left[0,z^{\max}\right]$ for some truncation point $0<z^{\max}<\infty$ chosen by us. First note that this entails that the simulated Bellman operator used for the self-approximating value function is biased since we do not sample from the full support $\mathcal{Z}=\mathbb{R}_{+}$; however, this bias can be controlled by choosing $z^{\max}$ large enough. We will explain below why we do not choose $\phi_{Z}\left(z^{\prime}\right)$ as a density with support $\mathbb{R}_{+}$. Using a uniform sampler, the corresponding Radon-Nikodym derivative takes the form $w_{Z}\left(z^{\prime}|z,d\right)=\pi\delta\left(z^{\prime}-z\right)+(1-\pi)f_{+}\left(z^{\prime}|z,d\right)$ where $\delta\left(\cdot\right)$ denotes Dirac's delta function. We approximate this by $\hat{w}_{Z}\left(z^{\prime}|z,d\right)=\pi I\left\{ z^{\prime}=z\right\} +(1-\pi)f_{+}\left(z^{\prime}|z,d\right)$ which entails another small approximation error. For the sieve-based version, we simply choose $\Phi_{Z}\left(z^{\prime}|z,d\right)=F_{Z}\left(z^{\prime}|z,d\right)$ and so $w_{Z}\left(z^{\prime}|z,d\right)=1$.

As explained in Section (ref), using a marginal importance sampler creates issues since it fails to adapt to the particular shape of the support of $F_{Z}\left(z^{\prime}|z,d\right)$. In particular, for a given choice of $z$, many of the draws from $\phi_{Z}\left(z^{\prime}\right)$ will tend to fall outside the support of $f_{Z}\left(z^{\prime}|z,d\right)$ and so will not contribute. In contrast, when $\phi_{Z}\left(z^{\prime}|z,d\right)=f_{Z}\left(z^{\prime}|z,d\right)$, the draws from $\phi_{Z}$ will by construction fall within the support of $f_{Z}\left(z^{\prime}|z,d\right)$. This can be seen in Figure (ref) where we have plotted the random draws obtained from the two different importance samplers used for the sieve-based and self-approximating solutions together with the actual support of $f_{Z}\left(z^{\prime}|z,d\right)$. In the left-hand side panel we have plotted pairs of the uniform draws, $\left(Z_{i},Z_{j}\right)$ for $i,j=1,...,N$, used for Rust's self-approximating method with $N=400$ and $z^{max}=1,000$, while in the right-hand side we have plotted $\left(z_{i},Z_{j}\left(z_{i},d\right)\right)$ where $z_{i}$ are uniform draws and $Z_{j}\left(z,d\right)\sim f_{Z}\left(\cdot|z,d\right)$. In both cases, we have marked the pairs for which the corresponding density, $f_{Z}\left(Z_{j}|Z_{i},d\right)$ and $f_{Z}\left(Z_{j}\left(z_{i},d\right)|z_{i},d\right)$, respectively, is positive. Clearly, the use of a marginal importance sampling density leads to very poor coverage of the actual support of $f_{Z}\left(z^{\prime}|z,d\right)$ as $z$ varies while by construction $\Phi_{Z}\left(z^{\prime}|z,d\right)=F_{Z}\left(z^{\prime}|z,d\right)$ does an excellent job. This translates into the former simulated Bellman operator exhibiting much larger variance compared to the latter.

This issue is further amplified when we introduce the normalization given in eq. ((ref)): Suppose that we had not included a discrete component $\pi I\left\{ z^{\prime}=z\right\} $ in the model. Then, with $Z_{i}\sim U\left[0,z^{\max}\right]$, $w_{N,Z,i}(Z_{j},d)=f_{+}\left(Z_{i}|Z_{j},d\right)/{\textstyle \sum}_{k=1}^{N}f_{+}\left(Z_{k}|Z_{j},d\right).$ Since $f_{+}(z^{\prime}|z,d)$ has bounded support, it often happens that ${\textstyle \sum}_{k=1}^{N}f_{+}\left(Z_{k}|Z_{j},d\right)=0$ for even large values of $N$ and so the simulated Bellman operator is not even well-defined. This issue will vanish as $N\rightarrow\infty$, but this on the other hand increases the computational burden since the self\textendash approximating method require us to solve for the value function at the $N$ draws. Introducing the discrete component in the model resolves this issue since now $w_{Z,N,i}(Z_{j},d)=\hat{w}_{Z}\left(Z_{i}|Z_{j},d\right)/{\textstyle \sum}_{k=1}^{N}\hat{w}_{Z}\left(Z_{k}|Z_{j},d\right),$ where ${\textstyle \sum}_{k=1}^{N}\hat{w}_{Z}\left(Z_{k}|Z_{j},d\right)>0$ for all $j=1,...,N$ by construction. Thus, $\pi>0$ functions as a regularization device.

Why not choose $\phi_{Z}\left(z\right)$ as a density with unbounded support in order to avoid the issue of truncation? In our initial experimentation, we did try out sampling from distributions with unbounded support, but the above numerical issues became even more severe in this case since the resulting draws are even more dispersed. Figure (ref) shows how the solution depends on $z^{\max}$. The effect of the truncation $z^{\max}$ will be model specific and in practice experimentation is required. If we, for example, simply set $z^{\max}=1,000,000$, the variance of the simulated Bellman operator becomes very large for a given $N$ due to the issue with undefined sample weights $w_{N,i}\left(z,d\right)$ mentioned above. At the same time, choosing $z^{\max}$ too small leads to a large bias. To balance the bias and variance, we ended up using $z^{\max}=1000$ which all subsequent numerical results for the self-approximating method is based on. Finally, we would like to stress that none of these issues appear for the sieve-based method.

Convergence properties and computation times

We first investigate the convergence properties of our solution methods for given choice of $K$ and $N$. Do they converge and if so how fast?

Global convergence properties of sieve method

As demonstrated in Theorem (ref), the simulated Bellman operators are always contraction mappings and so the self-approximating method is guaranteed to converge using successive approximations. In contrast, $\Pi_{K}\bar{\Gamma}_{N,\lambda}$ is not necessarily a contraction and so global convergence of the sieve method may fail, c.f. discussion in Section (ref). A sufficient condition for global convergence is $||\Pi_{K}||_{op,\infty}<1/\beta$ and we saw that $||\mathbf{P}_{K}||_{op,\infty}>1$ implies $||\Pi_{K}||_{op,\infty}>1$. However, even if $||\mathbf{P}_{K}||_{op,\infty}>1$, successive approximation may still converge: Across various parameter values of model, choices of sieve spaces and number of simulations, we did not encounter any failure of the sieve method to converge and the resulting approximate solution was well-behaved. This finding held across various initializations of the solution algorithms (initial choice of sieve coefficients). For example, we implemented the sieve method using $M=64$ evaluation points and using either $K=1$ or $K=4$ Chebyshev basis functions. We found that $||\mathbf{P}_{1}||_{op,\infty}=1$ while $||\mathbf{P}_{4}||_{op,\infty}>1.78$ and so the sieve method was guaranteed to converge for $K=1$ but not for $K=4$. Nevertheless, the method of successive approximations did in fact converge to a tolerance of $10^{-12}$ for both $K=1$ and $K=4$.

Successive approximation versus Newton-Kantorovich

In Section (ref) we advocated a hybrid of successive approximation (SA) and Newton-Kantorovich (NK) where we start with SA to ensure global convergence, and switch to NK iterations once the domain of attraction has been reached since NK generally converges faster. We illustrate this attractive feature of the NK algorighm in Figure (ref) where we have plotted the log residual error of the current value function approximation (relative to the “exact” solution) against the iteration count for the SA and NK algorithms, respectively, for four different values of $\beta$. As expected, the convergence of the SA algorithm requires a very large number of iterations ($>1000$) with computation time increasing in $\beta$, where as NK converges after less than 10 iterations and with the value of $\beta$ having little effect on its performance.

Figure (ref) is silent about the over-all computation time of SA relative to NK. Compared to SA, each NK iteration is more expensive since the former only requires computing the simulated Bellman operator evaluated at the value function obtained in the previous step while the latter, in addition, requires computing its functional derivative and inverting a $K\times K$ dimensional matrix for the integrated value function and a $KD\times KD$ dimensional matrix for the expected value function, c.f. Section (ref). With $K$ large, one could therefore fear that NK would become computationally too expensive.

In Figure (ref) we report best of 10 run-times for various levels of $K$ and $\beta$ and tolerance levels of SA and NK where we also include set-up time (time spent on initial computations before starting the actual algorithm). As expected, we find that NK is the faster of the two algorithms when $\beta$ is relatively large and $K$ is relatively small. With $K=5$ NK is faster across all levels of $\beta$ while for $K=100$ and $K=500$, SA is faster for moderate values of $\beta$. However, as we shall subsequently see, with $K=5$ the sieve method carries almost no bias and so choosing $K$ larger (such as 100 or 500) is actually unnecessary here and is only included here to illustrate potential issues with NK for models where a large number of sieve terms are needed to obtain a good approximation of the value function. Moreover, in most empirical applications, $\beta$ is chosen to be larger than 0.99 in which case NK still dominates SA even with $K=500$.

Approximation quality

We here investigate how the approximate value function is affected by the number of draws and the chosen projection basis. The goal is to demonstrate the rate results of the theoretical sections, and to compare the two types of basis functions spaces that we described above. We will take a partial approach and first fix $N$ to study the role of $K$, and then afterwards fix $K$ to study the role of $N$. All subsequent results are for the case of $\beta=0.95$. This is to save space. We implemented the methods for other values of $\beta$ and since the numerical results were qualitatively the same, we have left these out. The main difference in the numerical results is that higher values of $\beta$ tend to shift the overall level of the value function upwards and add more curvature to it. This in turn generally leads to an increase in the absoute bias and variance numbers. However, in terms of percentage bias and variance, the performance of the methods were very similar across different values of $\beta$.

Effect of varying $K$ for projection-based value function approximation

The theory for the projection-based value function approximation informs us that the choice of the basis functions will have a first-order effect on the bias while only a second-order effect on the variance. In particular, we expect $Bias\left(z\right)$, as defined in the beginning of this section, to satisfy $Bias\left(z\right)\cong\Pi_{K}\left(v\right)\left(z\right)-v\left(z\right)$, c.f. discussion following Theorem (ref), while $Var\left(z\right)$ should be much less affected by $K$. The actual size of the bias obviously depends on the curvature and smoothness of $v_{0}$ and the particular choice of basis functions. But we know that $v$ is an analytic function and with only moderate curvature, c.f. Figure (ref) and so expect it to be well-approximated by a small number of polynomial basis functions.

This is confirmed by the pointwise bias and standard deviation, $\sqrt{Var\left(z\right)}$, reported in Figure (ref): First, as can be seen in the left-hand side panel of Figure (ref), first-order B-splines lead to significantly larger point-wise bias compared to the other two sieve bases, namely second-order B-splines and Chebyshev polynomials. This is accordance with theory since we know that a smooth function is better approximated by higher-order polynomials, c.f. the error rates reported in Example 1 as a function of $s$. At the same time, second-order B-splines and Chebyshev polynomials exhibit very similar biases for a given choice of $K$.

The right-hand side panel of Figure (ref) shows the point-wise standard deviation across different choices of $K$ for the three different sieve bases. Consistent with the theory, the standard deviation of the value function approximation is not very sensitive to the particular choice of the sieve basis and the number of basis functions uses. That is, the sieve basis mostly affect the bias with only minor impact on the variance.

Finally, we examine how the bias behaves as we further increase $K$. Figure (ref) plots $\left\Vert Bias\right\Vert _{\infty}$ as a function of $K$. Similar to Figure (ref) we see much more rapid convergence when smooth basis functions are used, and with little improvement for $K$ greater than 9. This is not surprising given the reported shape of $v$. The second-order B-splines and Chebyshev basis functions produce very similar fits, even if they are evaluated on different grids and the B-splines have very different properties compared to Chebyshev polynomials. Indeed, the curves are practically overlapping. This is in accordance with the asymptotic theory that predicts that higher-order B-splines and Chebyshev polynomials should lead to similar biases. Moreover, the theory informs us that if $v$ is analytic, and this is the case in this particular implementation, we should expect the bias to vanish with rate $O(K^{-K})$ when using polynomial interpolation. The bias indeed does go to zero very quickly and so the numerical results support the theory.

Simulation errors, rates of convergence and asymptotic normality

We now compare the errors due to simulations and the rates with which these vanish for the two solution methods. For both methods, theory tells us that $N$ should have a first-order effect on the variance of the approximate value function which is supposed to vanish at rate $1/N$, c.f. Theorems (ref) and (ref). Our asymptotic theory is, on the other hand, silent about the size of simulation bias and the rate with which it should vanish with. However, we can think of both the sieve-based and self-approximating method as a nonlinear GMM-estimator where the simulated Bellman operator defines the sample moments. Importing results for GMM estimators, see, e.g., Newey&Smith2004, we should expect the simulation bias to be of order $1/N$.

In Figure (ref) we investigate this prediction by plotting $\left\Vert Bias\right\Vert $$_{\infty}$ and $\left\Vert \sqrt{Var}\right\Vert $$_{\infty}$ for the sieve-based method (left panels) and for the self-approximating method (right panels) for two different choices of $\sigma_{\varepsilon}$ and for across different values of $N$. To examine the rate with which the simulation bias and variance vanish we estimate the following an exponential regressions by NLS $||\sqrt{Var}||_{\infty}=\exp(\alpha_{SD}+\rho_{SD}\ln(N))$ and $||Bias||_{\infty}=\exp(\alpha_{Bias}+\rho_{Bias}\ln(N))$ where $\rho_{SD}$ and $\rho_{Bias}$ measures the rate that $||\sqrt{Var}||_{\infty}$ and $||Bias||_{\infty}$ vanishes with respectively. The resulting regression fit estimates are reported in both Figure (ref) as well as in Table (ref). In Table (ref) we present bias and standard deviation for $N=500$ as well as their rates of convergence both methods; with various values of $K$ for the sieve approximation method.

According to the theory, the variance should vanish with rate $1/N$ for both methods and we therefore expect $\rho_{SD}=-0.5$ so that $||\sqrt{Var}||_{\infty}$ vanish with $1/\sqrt{N}$. For the projection based method, we see that the rate with which the standard deviation shrinks to zero is indeed close to $-0.5$ for all values of $K>1$ and irrespectively of the value of $\sigma_{\varepsilon}$. For the self-approximating method we estimate the rate to $\rho_{SD}=-0.541$ when $\sigma_{\varepsilon}=100$, which is in line with the theory. However, $||\sqrt{Var}||_{\infty}$ is found to vanish with rate $1/N^{3.6}$ for $\sigma_{\varepsilon}=15$. This seems to indicate that the asymptotic theory developed in Theorems (ref) and (ref) do not provide a very accurate approximation of the performance of the self-approximating method for small and moderate choices of $N$ when the support of $Z_{t}|\left(Z_{t-1}=z,d_{t}=1\right)$ is small ($\sigma_{Z}=15$). We conjecture that the discrepancy between theoretical predictions and numerical results for the self-approximating method is due to the aforementioned issues with the marginal importance sampler discussed in Section (ref): Many of the draws are not used in the computation of the simulated Bellman operator because they fall outside the support of $Z_{t}|Z_{t-1}=z$ for a given choice of $z$. Thus, the effective number of draws is smaller than $N$ and changes as $z$ varies.

For the projection based method, the main source of bias is due to the sieve projection. From Figure (ref), we see that, with $N=200$ and $K=9$, the sieve-based methods using second-order B-splines or Chebyshev polynomials have virtually no bias, and both Figure (ref) as well as in Table (ref) also confirms that we practically eliminate by approximate the value function using Chebychev polynomials with $K=20$. However, there still remains a small bias that vanishes as $N$ grows. For small $K$, we see that the bias is roughly independent of $N$. As $K$ increases so does the dependence on $N$. However, even for $K=20$ where we estimate $\rho_{Bias}$ to be $0.21$ and $0.30$ for $\sigma_{Z}=15$ and $\sigma_{Z}=100$ respectively, the rate of convergence is far from $1/N$. This is probably due to the presence of higher-order bias components that our asymptotic theory does not account for.

For the self-approximating method, there is no sieve projection bias but a larger simulation induced bias that decreases with $N$. We obtain rate estimates of $1/N^{1.7}$ and $1/N^{1.4}$ for the bias when $\sigma_{Z}=15$ and $\sigma_{Z}=100$ respectively; these are slightly faster than expected but not too far from the theoretical predictions of $1/N$. For the self-approximating method, bias constitute more than half of RMSE when $N<600$ for $\text{\ensuremath{\sigma_{Z}=15}}$ (or $N<400$ for $\text{\ensuremath{\sigma_{Z}=100}}$), but since $||Bias||_{\infty}$ decays faster than $||\sqrt{Var}||_{\infty}$ , the simulation bias eventually becomes second order for large $N$.

Comparing $||MSE||_{\infty}$ for $N=500$ we find that the sieve-based method clearly dominates the self-approximating method when $\sigma_{Z}=15$, whereas the self-approximating method performs best when $\sigma_{Z}=100$. This is not entirely surprising since a large value of $\sigma_{\varepsilon}$ implies a large conditional support of $Z_{t}$ in which case the draws of the marginal sampler are more likely to fall within the support, c.f. the discussion in Subsection (ref). Thus, the over-all error of the self-approximating method will tend to be smaller when $\sigma_{Z}$ is large. The opposite is the case for the sieve based method which becomes more precise for smaller value of $\sigma_{Z}$ since the variance of the simulated Bellman operator used for this method gets smaller as $\sigma_{Z}$ gets smaller. This shows that there is considerably scope for improving the performance of the sieve-based method by more careful design of the sampling method.

Theorems (ref) and (ref) state that when $N$ is large, the approximate value functions should be normally distributed. We here investigate whether this asymptotic approximation is useful in practice by looking at the pointwise distribution of the approximate solutions obtained through both methods. In Figure (ref), we plot the distribution of $\left(\tilde{v}\left(z\right)-E\left[\tilde{v}\left(z\right)\right]\right)/\sqrt{Var\left(\tilde{v}\left(z\right)\right)}$ for $z=500$, where $\tilde{v}$ denotes a given approximation method, together with the standard normal distribution. It is here important to note we do not center the estimate around $v\left(z\right)$ but instead around $E[\tilde{v}(z)]$ ; this is due to the sizable bias of the self-approximating method. For the sieve-based method, we see that its normalized distribution is quite close to the standard normal irrespectively of the value of $\sigma_{Z}$. In contrast, the normal distribution is a poor approximation for the self-approximating method when $\sigma_{Z}=15$ when $N=500$; we expect this is due to the fact that the effective number of draws is quite small and so the asymptotic approximation is poor in this case. As expected the approximation gets better as $N$ and/or $\sigma_{Z}$ increases.

Effect of smoothing

The results reported above did not involve any smoothing bias. We now numerically study the effect of smoothing. This is done by, instead of integrating out the i.i.d. extreme value taste shocks $\varepsilon_{t}$ analytically as we have done so far, using Monte Carlo simulations to evaluate this part of the integral and then introducing the smoothing device to ensure that the simulated Bellman operator remains smooth. While this may appear somewhat artificial, the merit of doing this exercise is that we can use the same “exact” solution as benchmark as used above.

In Figure (ref) we plot the sup-norm of the mean squared error, $||\text{MSE}||_{\infty}$ as a function of $\lambda$ (the smoothing scale parameter) for the sieve-based method using $K=4$ or $8$ Chebyshev polynomials (similar results were obtained for the self-approximating method and so are left out). For $K=4$, the MSE increases monotonically as a function of $\lambda$ while for $K=15$ the bias due to smoothing is non-monotonic in $\lambda$. In both cases, at $\lambda=0$, any remaining biases are due to either sieve-approximation or simulations. Importantly, the bias due to smoothing is negiglible (relative to the other biases) for small and moderate values of $\lambda$ while the variance is largely unaffected. We have no theory or heuristics for choosing an optimal $\lambda$ to optimally balance bias and variance due to smoothing but the current numerical results indicate that choosing a quite small $\lambda$ value works well.

Performance in the bivariate case

We now examine how the solution methods perform in the bivariate case ($d_{Z}=\dim\left(Z_{t}\right)=2$) in order to see if there is any curse of dimensionality built into the two methods. We do this for two different models as described below.

An Additive DDP

We here follow the approach of Arcidiaconoetal2013 and rust1997comparison and build a $d_{Z}$-dimensional model by adding up $d_{Z}$ independent versions of the univariate model considered so far. That is, we choose the utilities and state dynamics as $\bar{u}(z,\varepsilon,d)=\sum_{i=1}^{d_{Z}}u(z_{i},\varepsilon_{i},d_{i})$ and $\bar{F}_{Z}(z^{\prime}|z,d)=\prod_{i=1}^{d_{Z}}F_{Z}(z_{i}^{\prime}|z_{i},d_{i})$, where $z=\left(z_{1},...,z_{d_{z}}\right)$ and $d=\left(d_{1},...,d_{d_{Z}}\right)$, with $F_{Z}(z_{i}^{\prime}|z_{i},d_{i})$ and $u(z_{i},\varepsilon_{i},d_{i})$ denoting the state transition and per-period utility in the univariate case as described in Section (ref). Note here that $\left(Z_{t,i},\varepsilon_{t,i}\right)$ and $\left(Z_{t,j},\varepsilon_{t,j}\right)$ are fully independent of each other, $i\neq j$ and the number of alternatives are $2^{d_{Z}}$, where $d_{Z}$ Thus, the model considers the joint replacement decision of $d_{Z}$ assets whose stochastic usages $\left(Z_{t,1},...,Z_{t,d_{Z}}\right)$ are mutually independent. Conveniently, the integrated value function of this multidimensional problem, $\bar{v}(z_{1},...,z_{d_{Z}})$, is simply the sum of the solutions to each of the underlying univariate models, $\bar{v}(z_{1},...,z_{d_{Z}})=\sum_{i=1}^{d_{Z}}v(z_{i})$, where $v(z_{i})$ is the solution to the univariate model in Section (ref). This is a rather simplistic multivariate model but it comes with the major advantage that we can obtain a very accurate approximation of the exact solution by simply adding up the “exact” solution found for the univariate case. With a more complicated multidimensional structure, the computational cost of finding the “exact” solution is much higher. However, when implementing our solution methods, we forgo forgo the knowledge of the additive structure of the solution and so treat the above model as a “proper” multivariate problem.

Simulation error

Given the issues with the self-approximating method for small values of $\sigma_{Z}=15$, we here focus exclusively on the case $\sigma_{Z}=100$. To get a sense of the pointwise performance of the self-approximating method, we plot the pointwise bias and standard deviation for this method with $N=3000$ in Figure (ref) together with the pointwise errors of the corresponding replacement (choice) probabilities. The overall shape and level of the integrated value function is quite well captured, and the same is true for the policy. However, the approximation errors tend to get larger out in the tails of the distribution and some of this comes from the fact that the issues with the marginal sampler used for the self-approximating method are amplified here. The problems are especially present in the off-grid evaluations, where we often have very few draws in a given region where we want to evaluate the value function or policies.

Next, we examine $||Bias||_{\infty}$ and $||\sqrt{Var}||_{\infty}$ for both methods as we increase $N$. These are plotted in Figure (ref) where it should be noted that the reported range of $N$ reported on the $x$-axis of the two figures differ substantially. This is due to the fact that the self\textendash approximating method became numerically unstable for $N$ smaller than 1,400 while no such issues were present for the sieve-based method. As in the univariate case, both bias and variance of the two methods vanish as $N$ increases. However, comparing Figures (ref) and (ref), while the errors of the sieve-based method in the bivariate case is of a similar magnitude as in the univariate case, the errors of the self-approximating method are much larger in the bivariate case. This seems to indicate a certain type of curse-of-dimensionality in this particular application of the self-approximating method. This is caused by the issues with the marginal importance sampler employed for this method.

Sieve approximation error

In the implementation of the sieve-based method we use as sieve basis the tensor product of univariate Chebyshev polynomials or B-splines. That is, given, say, $J$ univariate basis functions, say, $p_{1},...,p_{J}$, we construct our bivariate basis functions as $B_{i,j}\left(z_{1},z_{2}\right)=p_{i}\left(z_{1}\right)p_{j}\left(z_{2}\right)$ for $i,j=1,...,J$ yielding a total of $K=J^{2}$ bivariate basis functions. In particularly, we do not exploit the additive structure of the problem since we are interested in the practical contents of Theorems (ref) where no particular sparsity/special structure of the model is assumed to be known.

However, in practive, Chebyshev polynomials very easily pick up the additive structure and effectively sets the coefficients of the cross-product terms to zero. This is illustrated in Table (ref) in Appendix (ref), where we report the coefficients for one particular projection-based bivariate value function estimate using a tensor product of $J=5$ Chebyshev polynomials. However, this is due to the particular properties of the Chebyshev polynomials and is not enforced by us in the implementation. For example, if we instead use B-splines, the “estimated” coefficients of the cross-product terms were significantly different from zero, c.f. Table (ref) in Appendix (ref).

In the left-hand side panel (a) of Figure (ref), we report the uniform bias of the projection-based method with $N$ chosen very large for the additive bivariate model. We find that the bias vanishes as $K$ increases as in the one-dimensional model. However, convergence is now slower in $K$ relative to the univariate case and we require $K=50$ to obtain a sieve approximation bias of $10^{-2}$ while $K=7$ sufficed in the univariate case. This is consistent with theoretical error rates for polynomial interpolation where the rate slows down as the dimension of the problem increases, c.f. Section (ref).

A non-additive DDP

One concern with the numerical results reported for the bivariate additive model in the previous section is that they may understate the curse of dimensionality of the sieve method: The true value function is by construction additive in the two state variables and so interaction terms do not appear. This in turn implies that the computational complexity of solving this particular model is relatively low; in particular, the solution should be well-approximated by lower-dimensional sieves ($K$ small).

To investigate how the sieve method performs when applied to a more complex, non-additive model, we here consider a slightly more complicated bivariate model where we include a multiplicative interaction term so that maintenance and replacement costs of the two busses interact, $\bar{u}(z,d)=\sum_{i=1}^{2}u(z_{i},d_{i})-u(z_{1},d_{1})u(z_{2},d_{2})/20$. Such a structure could, for example, reflect that capacity constraints make it more costly to simultaneously replace the engines of both busses. The resulting value funtion will have a more complicated multidimensional structure and so we expect that the computational cost of our sieve method should be higher in this scenario.

The sieve approximation bias of our solution method for this model is reported in the right-hand side panel (b) in Figure (ref). Compared to panel (a) \textendash the additive case \textendash we see that more sieve terms are required in order to reach a specific absolute error level in the model with interactions. In Table (ref) in Appendix (ref) the coefficients on the first ten basis functions in each dimension and their interactions are reported. Compared to the Chebyshev-based solution earlier we see quite significant coefficients on the coefficients for the cross-terms. However, the coefficients on the basis functions tend to zero quite quickly as K increase. The sup-norm of the difference in the value function at $40.000$ evaluation grids is on the order of $10^{-5}$ when comparing the solutions with $K=50^{2}=2500$ and $K=30^{2}=900$ basis functions, and individual coefficients fall below $10^{-6}$ for univariate basis functions and cross products beyond the 22nd univariate basis functions, and below $10^{-8}$ around the 30th basis functions.However, it is important to stress that a large $K$ here only comes with a computational cost while a large $K$ has little effect on the variance of the sieve method. All together, we find that the sieve-based solution method works well also in higher dimensions, in particular when the model has a particular structure that can be utilized in the solution method.

Conclusion

We have proposed two novel methods for numerical computation of either the so-called integrated or expected value functions in a general class of dynamic discrete choice models.. Both methods rely on a smoothed simulated version of the Bellman operators definining the integrated and expected values functions. The smoothing facilitates both the practical implementation and the theoretical analysis of the approximate value functions. Under regularity conditions, we develop an asymptotic theory for the two methods as the number of simulations used to compute the simulated Bellman operators diverge. A set of numerical experiments show that our first method, the so-called self-approximating method can be somewhat unstable while the second one, which relies on sieve methods, apppears much more numerically robust. The next step is to develop methods for choosing the number of simulations, sieve basis functions and smoothing parameter $\lambda$ in a given setting so that the resulting approximate solution is of a good quality. Another area of research is to investigate how the proposed solution methods can be used for the estimation of dynamic discrete choice models.