Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
68,550 characters · 24 sections · 79 citation commands
Semiparametric Bayesian Estimation of Dynamic Discrete Choice Models
A dynamic discrete choice model is a dynamic program with discrete controls. These models have been used widely in various fields of economics, including labour economics, health economics, and industrial organization. See, for example, Rust_handbook:94 and Aguirregabiria_Mira:10 for literature surveys. In such models, a forward-looking decision-maker chooses an action from a finite set in each time period. The actions affect decision-maker's per-period payoff and the evolution of state variables. The decision-maker maximizes the expected sum of current and discounted future per-period payoffs.
Some state variables in these models are usually assumed to be unobserved by the econometrician (see, for example, page 1008 of Rust:87 for further discussion). Most of the previous work on estimation of dynamic discrete choice models imposes specific parametric assumptions on the distribution of the unobserved states or utility shocks. The most commonly used parametric assumption is that the unobserved states are extreme value independently identically distributed (i.i.d.). As shown in Rust:87 and Rust_handbook:94, under this assumption, the integrals over the unobserved states in the likelihood and the Bellman equations have closed form expressions, which considerably alleviates the computational burden of the model solution and estimation. At the same time, it is well known in the literature that imposing parametric distributional assumptions can be problematic, see, for example, Manski_book:99. Thus, it is desirable to relax these assumptions if possible.
There are several previous papers that treat the unobserved state distribution nonparametrically for the binary choice case. Aguirregabiria:10 shows the nonparametric identification of the shock distribution under particular assumptions on the per-period payoffs. NoretsTang2013 show that under an unknown distribution of the unobserved state and discrete observed states, the utility parameters and the unobserved state distribution are only set-identified. They also show how to compute the identified sets. BuchholzShumXu:20 provide identification results for the per-period payoffs when the observed state is continuous. The framework of ChristensenConnault2023 for structural models expressed through a finite number of moment equalities and inequalities can be used to check the sensitivity of counterfactuals to the variation of the utility shocks distribution within a neighborhood of extreme value distribution in dynamic discrete choice models with a finite observed state space; Rust's binary choice model of bus engine replacement is used in that paper for illustration.
For the multinomial choice case, chen_2017 uses exclusion restrictions (a subset of the state variables affects only current utility, but not state transition probabilities) to obtain identification and estimation results. In settings without exclusion restrictions, Norets_ddc_mult:11 shows that it is in principle possible to extend the method from NoretsTang2013 to compute the identified set in multinomial case, but it is computationally very difficult.
In this paper, we propose a tractable semiparametric estimation method applicable to the general multinomial choice case. It is based on modeling the unknown distribution of shocks by a finite mixture of extreme value distributions with a varying number of mixture components. Our approach exploits the analytical tractability of extreme value distributions and the flexibility of the location-scale mixtures. The unobserved utility shocks can be integrated out analytically in the likelihood function and the expected value functions, similarly to the case with extreme value distributed shocks. At the same time, we show that the location-scale mixtures can approximate densities from a large nonparametric class in an appropriate distance and that for any given distribution of utility shocks, a finite mixture of extreme value distributions can deliver exactly the same conditional choice probabilities. Posterior concentration on the identified sets of utility parameters and the distribution of shocks is an implication of these results. We implement the Bayesian approach to inference for the model using Hamiltonian Monte Carlo and an approximately optimal reversible jump Markov chain Monte Carlo (MCMC) algorithm from Norets2017mcmc. Similarly to NoretsTang2013, frequentist confidence sets for identified sets can also be computed from the MCMC output.
We apply our framework to binary and multinomial choice models. For the binary dynamic choice model from Rust:87, our approach delivers estimation results that are consistent with the previous literature on semiparametric estimation (NoretsTang2013). For the multinomial choice model of medical care use and work absence from Gilleskie_Eca:98, we demonstrate how uncertainty about model parameters and counterfactuals increases when the distributional assumptions on the shocks are relaxed. Moreover, we show that the standard dynamic logit model can deliver misleading results, especially about counterfactuals, when the shocks are not extreme value distributed. Our semiparametric approach delivers reliable inference in these settings.
Even when the distribution of the utility shocks is assumed to be known, parameters and counterfactuals in dynamic discrete choice models could still be only set identified under a variety of scenarios such as very flexible specifications of utility functions, lack of exclusion restrictions or variation in transitions for the observed state variables, and unknown time discount factors; see, for example, Rust_handbook:94, MagnacThesmar:02, NoretsTang2013, AbbringDaljord2020, and KalouptsidiScottSouzaRodrigues2021. Our estimation framework does not require any special adjustments to accommodate such scenarios since parameters and counterfactuals are already set identified under nonparametric specification of the distribution of shocks.
The rest of the paper is organized as follows. Section (ref) describes the general model setup. In Section (ref), we introduce our semiparametric framework. In Section (ref), we describe the Bayesian estimation method. Section (ref) presents theoretical results. Sections (ref) and (ref) contain the applications. Possible framework generalizations are discussed in Section (ref). Derivations, proofs, and implementation details are given in appendices.
In the infinite-horizon version of the model, the decision maker maximizes the expected discounted sum of the per-period payoffs
where $d_t \in \{0,1,\ldots,J \}$ is the control variable, $x_t \in X=\{1,\ldots,K\}$ is the state variable observed by the econometrician, $\epsilon_t = (\epsilon_{t0},\epsilon_{t1},\ldots,\epsilon_{tJ})^T \in R^{J+1}$ is the state variable unobserved by the econometrician, $\beta$ is the time discount factor, and $u(x_t,d_t,\epsilon_t)$ is the per-period payoff. The decision-maker observes both $x_t$ and $\epsilon_t$ at time $t$ before making the decision.
Following Rust:87 and the subsequent literature, we assume that (i) the per-period payoffs are additively separable in $\epsilon_t$, $u(x_t,d_t,\epsilon_t)= u(x_t,d_t) + \epsilon_{td_t}$; (ii) $\epsilon_t$'s are independent of other variables and independently identically distributed (i.i.d.) over time according to a distribution $F$ with zero mean; (iii) the observed states evolve according to a controlled Markov chain $G$ with transition probabilities $G_{x}^j =\{Pr(x_{t+1}|x_t=x,d_t=j), \, x_{t+1} \in X\}$ and an initial distribution $\{Pr(x_{1}), \, x_1 \in X\}$. The utility functions are assumed to depend on a vector of unknown parameters, $\theta \in \mathbb{R}^{d_\theta}$, that are estimated. Below, we often omit $\theta$ in $u(x_t,d_t; \theta)$ and related objects such as value functions for notation brevity. As in Rust:87, Gilleskie_Eca:98, and most of the literature, the time discount factor is assumed to be known; it is usually calibrated to imply a reasonable value of a risk free annual interest rate. In most applications, the observed state transition probabilities are estimated in a first stage prior to estimation of the preference parameters $\theta$ since it can be done directly from the observed state transitions with a relatively high precision and without solving the dynamic program. In line with that and for notation brevity we treat $G$ as known and fixed. Extensions of our results and methodology to continuous $X$ and unknown $G$ are discussed in Section (ref).
Under mild regularity conditions (BhattacharyaMajumdar:89), the decision problem in (ref) admits the following Bellman representation
where $Q$ is called the Emax function and $G_x^j Q$ denotes $E(Q(x_{t+1})|x_t=x,d_t=j)$. The conditional choice probability (CCP) can be expressed as
For a panel of observations, $D^n=\{x_{it}, d_{it},\, i=1,\ldots,n,\, t=1,\ldots,T\}$, of $n$ decision makers over $T$ time periods, the partial likelihood function (with the fixed $G_x^j$ pre-estimated from the observed transitions as is commonly done in practice) can be expressed as
Rust:87 proposed to solve the dynamic optimization problem by first iterating on the Bellman equation (ref) to get close to the fixed point $Q$ and then using a Newton method that quickly converges to the fixed point from a close starting point. With $Q$ at hand, one can compute the CCPs in (ref) and evaluate the likelihood function at a given $\left(u;\beta;G \right)$. Alternatively, JuddSu:08 proposed to use constrained optimization to maximize the likelihood function subject to (ref). In either scenario, assuming that $\epsilon_{tj}$'s are i.i.d. Gumbel (or extreme value type I) delivers analytical expressions for the integrals in (ref) and (ref),
In the resulting dynamic logit specification, the computational burden of the model solution and estimation is considerably alleviated. Hence, the dynamic logit is predominantly used in applications of the estimable dynamic discrete choice models. At the same time, the econometrics literature suggests that the distributional assumptions could be problematic in general, see, for example, Manski_book:99. In the following section, we specify a non-parametric model for the distribution of shocks for the general multinomial choice case that provides analytical simplifications comparable to those of the dynamic logit.
Rather than making a particular parametric assumption, we model the distribution of unobserved states using a flexible mixture specification. In order to reduce the number of parameters, we use an innocuous normalization $\epsilon_{t0} = 0$ (the agent's decisions and value functions do not change if $\epsilon_{t0}$ is subtracted from the per-period payoff $u(x_t,d_t,\epsilon_t)$ for all $d_t \in \{0,1,\ldots,J\}$).
For $\mu \in \mathbb{R}^J$ and $\sigma>0$, let us define a multivariate Gumbel density by
is the univariate Gumbel density and $\gamma$ is the Euler-Mascheroni constant. Some relevant properties of the Gumbel distribution are outlined in Appendix (ref).
For $\mu_k \in \mathbb{R}^{J}$, $\sigma_k \in \mathbb{R}_+$, $\omega_k \in [0,1]$, $k=1,\ldots, m$, and $\sum_{k=1}^m \omega_k =1$, we model the unknown density by a location-scale mixture of Gumbel densities
with a variable number of mixture components $m$ for which a prior distribution on the set of positive integers is specified. Mixture models are extensively used in econometrics and statistics literature, see monographs by McLachlanPeel:00 and FruhwirthSchnatter:06 for references. It is well known that location-scale mixtures with a variable or infinite number of components can approximate any continuous or smooth density arbitrarily well. For example, Bayesian models based on normal mixtures deliver optimal up to a log factor posterior contraction rates in adaptive estimation of smooth densities (Rousseau:10, ShenTokdarGhosal2013, and NoretsPelenis_Eca_2022). To develop intuition for this type of results note that the standard nonparametric density estimator based on kernel $\phi$ is a special case of (ref), or, alternatively and more in line with the actual proofs, the expectation of the standard kernel density estimator is a continuous mixture that can be discretized into a special case of (ref). Thus, it is reasonable to expect that the specification (ref) is very flexible. Indeed, in Section (ref), we show that it can approximate smooth multivariate densities arbitrarily well in an appropriate distance so that the conditional choice probabilities and the Emax function implied by the model with (ref) approximate those from the model with an arbitrary smooth density for $\epsilon_t$.
The model specification with (ref) also possesses attractive analytical properties. If a normalization $\epsilon_{t0} = 0$ is not imposed and $(J+1)$-dimensional version of (ref) is used, then $Q$ and $p$ could be expressed as mixtures of the appropriately recentered and rescaled expressions from the dynamic logit model (ref). However, even if the normalization $\epsilon_{t0} = 0$ is imposed, which is preferred as it reduces the dimension of the distribution we model nonparametrically, closed form expressions for $Q$ and $p$ are still available. They are presented in the following lemma.
The derivations of (ref) and (ref) can be found in Appendix (ref). The derivatives of (ref) and (ref) that are useful for the model solution and estimation are given in Appendix (ref). Similarly to Rust:87, we obtain the solution of the Bellman equation (ref) by a Newton-Kantorovich method described in Appendix (ref).
In estimation of models based on location-scale mixtures with a variable number of components, the econometrician faces several problems. First, the scale parameters need to be bounded away from zero; otherwise, the likelihood function is unbounded. Second, the likelihood function is a rather complex function of parameters with multiple modes. Third, the number of mixture components needs to be selected in the estimation procedure. Finally, there is usually considerable uncertainty about the estimated parameter values and it should be taken into account in model predictions and counterfactual analysis.
The Bayesian approach to inference and the associated simulation methods are well suited for solving these problems. Prior distributions can provide soft constraints for the scale parameters and an appropriate penalization for the number of mixture components or model complexity. MCMC methods can successfully explore very complex posterior or likelihood surfaces. Posterior predictive distributions for objects of interest automatically incorporate the uncertainty about parameter values including the number of mixture components.
In addition to the normalization $\epsilon_{t0} = 0$, the scale of $\epsilon_t$ can be innocuously normalized. Instead, to simplify the MCMC algorithm we impose a location and scale normalization on the parameters of the per-period payoffs and keep the location and scale of (ref) unrestricted. Specifically, consider a linear in parameters utility specification \[ u(x_t,j,\epsilon_t;\theta)=\theta_j + z_j(x_t)^\prime \theta_{J+1:d_\theta} + \epsilon_{tj}, \] where $z_j(x_t)$ are known functions of the observed state variables. In this specification, the intercepts $\theta_j$, $j=0,\ldots,J$, can be fixed to arbitrary values as long as the locations of $\epsilon_{tj}$, $j=1,\ldots,J$, are unrestricted.
In applications, we set $\theta_0=0$ and $\theta_j$, $j=1,\ldots,J$ to $\hat{\theta}_j^{dl}$, the estimates obtained under the dynamic logit specification. To normalize the scale of $\epsilon_t$, we assume that the sign of one of the coefficients, say $\theta_{J+1}$, is known and we keep this coefficient fixed (to the corresponding dynamic logit estimate, $\hat{\theta}_{J+1}^{dl}$). Thus, the MCMC algorithm produces draws of $\theta_{J+2:d_\theta}$ and the mixture parameters $(\psi_{1m},m)$ in (ref). For comparisons of the estimation results with the dynamic logit estimates and the identified sets in NoretsTang2013, the parameter draws are renormalized for reporting as follows
where the addition of $\sum_{k=1}^m \omega_k \mu_k$ to the intercepts corresponds to the zero mean for shocks and the scale factor is defined by the mixture parameters \[ s = \log 2\big/ E[\tilde{\epsilon}_{t1} 1(\tilde{\epsilon}_{t1}\geq M_{\tilde{\epsilon}_{t1}})], \] where $\tilde{\epsilon}_{t1} = \epsilon_{t1} - \sum_{k=1}^m \omega_k \mu_{1k}$, $M_{\tilde{\epsilon}_{t1}}$ denotes the median of $\tilde{\epsilon}_{t1}$, and $\epsilon_{t1} \sim \sum_{k=1}^m \omega_k \phi(\cdot ; \mu_{1k},\sigma_k )$. There are many possible scale normalizations. The particular scale normalization we use here reduces to the one introduced by NoretsTang2013 for the binary choice case. Let us emphasize that the normalizations discussed above are innocuous for estimation and counterfactual analysis as long as the assumed sign of $\theta_{J+1}$ is correct.
Let us introduce the prior distributions for the parameters of the mixture in (ref). We use the following prior distributions on the number of mixture components and the mixing weights,
where the hyperparameters $\underbar{a}$, $\underbar{A}_m$, and $\tau$ are specified in the applications below. For the theoretical results obtained in the present paper, we only need $\Pi(m)>0, \forall m$ and full support on the simplex for $\Pi(\omega_1,\ldots,\omega_m|m)$. Nevertheless, the functional forms in (ref) and (ref) perform well in applications and deliver optimal posterior contraction rates in nonparametric multivariate density estimation by mixtures of normal distributions, see, for example, ShenTokdarGhosal2013 and NoretsPelenis_Eca_2022. We allow the scale parameter $\sigma_k$ to have a multiplicative part $\sigma$ that is common across the mixture components: $\sigma_k=\tilde{\sigma}_k \cdot \sigma$. This multiplicative specification performs well in a variety of applications of location-scale mixture models (see, for example, Geweke:05) and is also important for the aforementioned optimal posterior concentration results for mixtures of normals. In the applications, we use finite mixtures of normals as flexible priors for $\log \sigma$, $\log \tilde{\sigma}_k$ and the location parameters $\mu_{kj}$.
Our MCMC algorithm for simulating from the model posterior distribution combines Hamiltonian Monte Carlo (HMC) for simulating parameters conditional on the number of mixture components and an approximately optimal reversible jump algorithm from Norets2017mcmc for simulating the number of mixture components. HMC is a very popular and efficient MCMC algorithm; see, for example, neal2012hmc for an introduction. HMC requires only evaluation of the likelihood and the prior and their derivatives. The proposals in HMC are obtained following the Hamiltonian dynamics on the parameter space that describe the movement of a puck on a friction-less surface with some initial random momentum. For implementing the HMC step of the algorithm we utilize the HMC sampler from the Matlab Statistics and Machine Learning toolbox. The package can choose HMC's parameters, such as a step size, automatically, and we perform this automatic initialization once for each value of $m$ that we encounter in the MCMC run. The package works only with unbounded parameters. Hence, we transform the bounded parameters, such as mixing weights and scales, for the HMC step. The form of the prior for the transformed parameters and the derivatives of the likelihood used in the algorithm are reported in Appendix (ref).
For the reversible jump algorithm, we need to transform the mixing weights into unnormalized weights $\gamma_k$, $k=1,\ldots,m$, so that they have interpretation under different values of $m$. Specifically, conditional on $m$, $\omega_k = \gamma_k/\sum_{l=1}^m \gamma_l$ and the Dirichlet prior on $(\omega_1,\ldots,\omega_m)$ corresponds to a gamma prior for the unnormalized weights: $\gamma_k|m \sim Gamma(\underline{a}/m,1)$, $k=1,\ldots,m$. Let $\psi_k = (\mu_k,\tilde{\sigma}_k,\gamma_k)$ and $\psi_{1m}=(\theta, \sigma,\psi_1,\ldots,\psi_m)$, where $\theta$ includes model parameters such as coefficients in the utility functions. With this notation, the likelihood function is denoted by $p(D^n|m,\psi_{1m})$.
The following short description of the reversible jump algorithm is adapted from NoretsPelenis_Eca_2022, see Norets2017mcmc for more details. Denote a proposal distribution for the parameter of a new mixture component $m+1$ by $\tilde{\pi}_{m+1}(\psi_{m+1}|D^n,\psi_{1m})$. The algorithm works as follows. Simulate proposal $m^\ast$ from $Pr(m^\ast=m+1|m)=Pr(m^\ast=m-1|m)=1/2$. If $m^\ast=m+1$, then also simulate $\psi_{m+1} \sim \tilde{\pi}_{m+1}(\psi_{m+1}|D^n,\psi_{1m})$. Accept the proposal with probability $\min\{1, \alpha(m^\ast,m)\}$, where
Innocuous random relabeling of mixture components increases the acceptance probability for attempts to delete $m$-th mixture component ($m^\ast=m-1$). Norets2017mcmc shows that an optimal choice of the proposal distribution $\tilde{\pi}_m$ is the conditional posterior $p(\psi_{m+1}|D^n,m+1,\psi_{1m})$. The conditional posterior can be evaluated up to a normalization constant; however, it seems hard to directly simulate from it and compute the required normalization constant. Hence, we use a Gaussian approximation to $p(\psi_{m+1}|D^n,m+1,\psi_{1m})$ as the proposal (with the mean equal to the conditional posterior mode, obtained by a Newton method, and the variance equal to the inverse of the negative of the Hessian evaluated at the mode).
The algorithm pseudo code is presented below.
\FloatBarrier
The Matlab code for the MCMC algorithm and replication instructions for the estimation results in the applications in Sections (ref) and (ref) are publicly available.\footnote{\url{https://anorets.github.io/papers/mix_ddcm_code.zip}}
In this section, we show that location-scale mixtures of Gumbel densities can arbitrarily well approximate densities from a large nonparametric class. These approximation results combined with the Schwartz:65's theorem imply a posterior consistency result for the set identified model parameters. We also show that a model with a finite mixture of Gumbels can exactly match the CCPs from a model with an arbitrary distribution of shocks.
Let us first define a distance for distributions of utility shocks: for $F_i$ with density $f_i$, $i=1,2$, \[ \rho(F_1,F_2) = \int (1+ \sum_{j=0}^J |\epsilon_j| ) | f_1(\epsilon) -f_2(\epsilon) | d\epsilon. \] This distance is appropriate for our purposes as the Emax function and the conditional choice probabilities are continuous in that distance as shown in the following lemma.
The lemma holds irrespective of whether the innocuous normalization $\epsilon_{t0}=0$ is imposed. Its proof is given in Appendix (ref).
The following lemma shows that densities satisfying smoothness and finite moment conditions can be approximated by mixtures of Gumbels in distance $\rho$.
We conjecture that the smoothness and tail conditions on $f$ in the lemma can be weakened at the expense of the proof simplicity. The lemma is proved in Appendix (ref). The proof uses only smoothness and tail conditions on $\phi$ that are shown to hold for Gumbel densities in Lemmas (ref) and (ref) in Appendix (ref). Thus, Lemma (ref) holds for more general location-scale mixtures. These generalizations do not seem essential and we do not elaborate on them here for brevity.
The final intermediate result that we need for establishing posterior consistency is the continuity of finite Gumbel mixtures in parameters in distance $\rho$, which we present in the following lemma.
Let us denote the short panel dataset by $D^n = \{d_{it},x_{it}, t=1,\ldots,T, \, i=1,\ldots,n\}$; the observations are assumed to be independently identically distributed over $i$, with a small $T$ and a large $n$. The utility function is parameterized by a vector $\theta \in \mathbb{R}^{d_\theta}$, $u(x,d; \theta)$. Let $P(\theta,F)=\{p(d|x;\theta,F), \, x \in X,\, d=0,\dots,J\}$ denote the collection of the CCPs for the distribution of shocks $F$ and parameters $\theta$.
The theorem shows that the posterior concentrates on the set of parameters and distributions of shocks $(\theta,F)$ such that their implied CCPs $P(\theta,F)$ are arbitrarily close to the data generating CCPs $P(\theta_0,F_0)$. To prove this result we use Schwartz:65 posterior consistency theorem: if the prior puts positive mass on any Kullback-Leibler neighborhood of the data generating distribution then the posterior puts probability converging to 1 on any weak neighborhood of the data generating distribution. Since $X$ is finite, the convergence in weak topology and Kullback-Leibler divergence for distributions on $\{d_{it},x_{it}, t=1,\ldots,T\}$ are equivalent to convergence for vectors $\{p(d|x), \, x \in X,\, d=0,\dots,J\}$ in a euclidean metric when $G^d$ and the distribution of the initial $x_{i1}$ are fixed and satisfy our theorem condition (iii). Thus, to obtain the conclusion of the theorem we only need to establish that the prior puts positive probability on any euclidean neighborhood of $P(\theta_0,F_0)$. First, note that when $u(x,d; \theta)$ is continuous in $\theta$, $P(\theta,F)$ is also continuous in $\theta$ in our settings, see, for example, Norets_ddcm_diff_cont:09; and, thus, Lipschitz continuity of $P(\theta,F)$ in $F$ from Lemma (ref) delivers continuity of $P(\theta,F)$ in $(\theta,F)$. The finite mixture approximation result in Lemma (ref), the continuity of $P(\theta,F)$ in $(\theta,F)$, the continuity of finite mixtures in parameters in Lemma (ref), and the theorem conditions (vi) and (vii) on the priors, imply a positive prior probability for any neighborhood of $P(\theta_0,F_0)$, and thus, the theorem conclusion. Possible extensions of the theorem (and the lemmas above) to continuous $X$ and unknown $G^d$ are discussed in Section (ref).
Theorem (ref) characterizes the support of the posterior in the limit but not its shape, which can also be of interest. Note that the data depend on $(\theta,F)$ only through CCPs $P(\theta,F)$ and the posterior for CCPs concentrates at $P(\theta_0,F_0)$. Therefore, the posterior for $(\theta,F)$ converges to the conditional prior $\Pi(\theta,F|P)$ at $P=P(\theta_0,F_0)$ under continuity conditions on $\Pi(\theta,F|P)$, see, for example, PlagborgMoller2019. As the distribution of shocks is an infinite dimensional object and the solution to the dynamic program does not have a simple explicit form, it appears difficult to characterize the conditional prior $\Pi(\theta,F|P)$, which is implied by the map $P(\theta,F)$ and the prior on $(\theta,F)$. Nevertheless, we can deduce from our approximation and continuity results that under the conditions of Theorem (ref), for $\delta > 0$ there exists $\tilde{\delta}>0$ such that $\theta \in B_{\tilde{\delta}}(\theta_0)$ and $F \in B_{\tilde{\delta}}(F_0)$ imply $P(\theta,F) \in B_{\delta}(P(\theta_0,F_0))$ and
which suggests that the conditional prior would not rule out the data generating parameter values.
In this subsection, we show that for a finite observed state space, our model formulation based on finite mixtures can exactly match the CCPs from a model with an arbitrary distribution of shocks.
The result in Lemma (ref) holds not only for mixtures of Gumbels but more generally for location-scale mixtures of distributions with finite first moments, which is evident from the proof presented in Appendix (ref). Lemma (ref) can be used to relax the smoothness assumptions on $F_0$ in the posterior consistency results of Theorem (ref). Specifically, conditions (v) in Theorem (ref) can be replaced by conditions (iii) in Lemma (ref); in the proof, the approximation results in Lemma (ref) can be replaced by the exact CCPs matching results in Lemma (ref). Nevertheless, the approximation results in Lemma (ref) have independent value. First, they hold for infinite and continuous $X$. Furthermore, they imply that the prior on the distribution of shocks is flexible in a sense that it puts positive probability on any metric $\rho$ neighborhood in a large nonparametric class of distributions, which suggests that the conditional prior for the distribution of shocks and parameters given CCPs, $\Pi(\theta, F | P)$, is also flexible as discussed at the end of Section (ref). Finally, while asymptotically the number of mixture components is bounded for the exact matching in Lemma (ref) and has to increase to infinity for the approximation results in Lemma (ref), in practice, a small number of mixture components delivers sufficiently good approximations and the exact matching requires a very large number of mixture components ($m=182$ for Rust's bus engine replacement model).
NoretsTang2013 and generalizations of their results to the multinomial case in Norets_ddc_mult:11 show that the utility parameters $\theta$ and the distribution of shocks $F$, and, thus, functions of $(\theta,F)$ such as results of counterfactual experiments, are set identified in the present settings. MoonSchorfheide:12 show that in contrast to the point identified regular settings, the Bayesian credible sets for set identified parameters do not have frequentist coverage properties and are too small from the classical perspective. NoretsTang2013 point out in their Section 3.2 that Bayesian and classical inference results can be reconciled if inference is performed on the identified sets. KlineTamer2016 further study this approach and Kitagawa:11 obtain related results under multiple priors for set identified parameters. In this subsection, we describe how credible and confidence sets for identified sets can be defined and computed from the output of our MCMC algorithm.
Suppose the data generating values of parameters are $(\theta_0,F_0)$ and we are interested in $\eta_0=g(\theta_0,F_0)$. The data generating values of CCPs, $P_0=P(\theta_0,F_0)$, can be consistently estimated from the observed data, and, thus, are considered known in the identification analysis. The identified set for $\eta_0$ is defined by \[ I_\eta(P_0) = \{\eta=g(\theta,F), \; \forall (\theta,F) \;s.t.\; P(\theta,F)=P_0\}. \] Following NoretsTang2013 and KlineTamer2016, we can define a posterior distribution on the space of identified sets $I_\eta(P)$ using the marginal posterior distribution on the CCPs $P$. Then, a $1-\alpha$-credible set for $I_\eta(P)$ can be defined from a $1-\alpha$-credible set for CCPs, $B_{1-\alpha}^P$, by
When $B_{1-\alpha}^P$ is a $1-\alpha$ highest posterior density set and the Bernstein - von Mises theorem holds for the point identified CCPs $P$, $B_{1-\alpha}^P$ asymptotically has a $1-\alpha$ frequentist coverage probability for $P(\theta_0,F_0)$ and, thus, $B_{1-\alpha}^{I_\eta}$ has at least a $1-\alpha$ frequentist coverage probability for the identified set $I_\eta(P_0)$ and $\eta_0$.
The sets in (ref) might be conservative; nevertheless, it would be prudent to report them in applications as in the limit they do not depend on the shape of the prior and possess both frequentist and Bayesian properties.
For a scalar $\eta_0$, an approximation to $B_{1-\alpha}^{I_\eta}$ can be computed from a sample of MCMC posterior draws $\{\theta^{(l)},F^{(l)},P^{(l)}=P(\theta^{(l)},F^{(l)}), \; l=1,\ldots,L\}$ as follows. First, we obtain an approximation to $B_{1-\alpha}^P$ \[ \hat{B}_{1-\alpha}^P = \left\{P:\; (P-\bar{P})^\prime \hat{\Sigma}^{-1} (P-\bar{P}) \leq \chi^2_{1-\alpha}(JK)\right\}, \] where $\bar{P}=\sum_{l=1}^L P^{(l)}/L$, $\hat{\Sigma} = \sum_{l=1}^L (P-\bar{P})(P-\bar{P})^\prime/L$, and $\chi^2_{1-\alpha}(JK)$ is the $1-\alpha$ quantile of the $\chi^2$ distribution with $JK$ degrees of freedom. Then, $B_{1-\alpha}^{I_\eta}$ is approximated by \[ \hat{B}_{1-\alpha}^{I_\eta}=\left[\min_{l:\: P^{(l)} \in \hat{B}_{1-\alpha}^P} \eta\left(\theta^{(l)},F^{(l)}\right), \max_{l:\: P^{(l)} \in \hat{B}_{1-\alpha}^P} \eta\left(\theta^{(l)},F^{(l)}\right)\right]. \] We report these sets along with the standard HPD sets in the application in Section (ref).
NoretsTang2013 propose a method for computing identified sets for parameters in dynamic binary choice models and apply their method to the Rust:87's model. In this section, we show that our semiparametric model can also recover the identified set for that model.
In each time period $t$, the agent decides whether to replace the bus engine ($d_t=1$) or not ($d_t=0$) given the current mileage $x_t$ of the bus. Replacing an engine costs $\theta_0$. If $d_t=0$, then the agent conducts a regular maintenance which costs $-\theta_1x$. The utility function of the agent is $u(x,0) = \theta_0 + \theta_1 x$ and $u(x,1) =\epsilon$. Rust:87 assumes that $\epsilon$ follows the logistic distribution. The mileage $x_t\in \{1, \ldots ,K=90 \}$ evolves over time following the transition probabilities: $Pr(x_{t+1}|x_t,d_t=0)=\pi_0$ for $x_{t+1}-x_t = 0$; $Pr(x_{t+1}|x_t,d_t=0)=\pi_1$ for $x_{t+1}-x_t = 1$; $Pr(x_{t+1}|x_t,d_t=0)=1-\pi_0-\pi_1$ for $x_{t+1}-x_t = 2$; and $Pr(x_{t+1}|x_t,d_t=0)=0$ otherwise. When the engine is replaced ($d_t=1$), the mileage restarts at $x_t=1$.
As in NoretsTang2013, we use the following data generating process: logistic distribution for $\epsilon$, $\theta_0 = 5.0727, \theta_1= -0.002293,\pi_0 = 0.3919, \pi_1 = 0.5953$ and the discount factor $\beta = 0.999$.
At the data generating parameters, we solve the dynamic program to obtain the vector of CCPs, $(p(d|1),\ldots,p(d|K))$ for $d=0,1$. Rather than using simulated observations in this exercise, we use the true CCPs as the sample frequencies and report the results for different sample sizes. Specifically, for a given $N \in Z^+$, we set $n_{dx}$, the number of times $d$ was chosen at each state $x$ as follows, $n_{0x} = p(0|x) \times N$ and $n_{1x}=p(1|x) \times N$ for $x=1,\ldots,K$. In this way, we can check if our MCMC algorithm for the semiparametric model specification can recover the identified set computed by the algorithm from NoretsTang2013 for a given fixed vector of CCPs. Estimation results for simulated data are presented for the multinational choice application in Section (ref).
Since we do not impose a location and scale normalization on the mixture specification for the distribution of $\epsilon$ and the model has only two utility parameters that are defined by the location and the scale, the values of $(\theta_0,\theta_1)$ corresponding to the location and scale of the logistic distribution are computed from $(\sigma,\omega_k,\mu_k,\sigma_{k}, \, k=1,\ldots,m)$ by formula (ref). We use the following flexible prior distributions that are tuned to spread the prior probability over a large region for $(\theta_0,\theta_1)$ that includes the identified set.
Below, we report the draws of $(\theta_0,\theta_1)$ obtained from the draws of $(\sigma,\omega_k,\mu_k,\sigma_{k}, \, k=1,\ldots,m)$ for the prior and the posterior for $N \in \{ 3, 10 \}$ and compare them to the true identified set. Panel (a) of Figure (ref) shows the prior draws of the utility parameters. First note that the prior draws are not uniformly distributed on the utility parameter space. In practice, it is difficult to come up with a prior for the distribution parameters that implies a uniform prior in the $\theta$ space. Second, many prior draws are outside of the identified set. The other two panels in Figure (ref) show posterior draws of utility parameters with different number of observations $N \in \{3, 10 \}$. The posterior concentrates more on the identified set as $N$ increases.
To assess the convergence of the MCMC algorithm, consider Figure (ref) showing the draws of utility parameters for 500,000 MCMC iterations. As can be seen from the figure, the chain sweeps through the identified set repeatedly during the MCMC run. Thus, we conclude that our approach can be used to recover the identified sets of utility parameters. Figure (ref) in Appendix (ref) shows additional evidence of MCMC convergence including trace plots of $m$, $\sum_{k=1}\omega_k\mu_k$, and other parameters.
\FloatBarrier
In this section, we illustrate our methodology using a multinomial choice model of medical care use and work absence from Gilleskie_Eca:98. An extension of NoretsTang2013 to the multinomial case seems computationally infeasible and, hence, in this application we provide comparisons of our method only with a dynamic logit specification.
In the model, individuals occupy one of $2$ distinct health states: well, $k=0$, or sick, $k=1$. An individual receives the utility associated with being well until contracting an illness of a specific type (we make the simplifying assumption that there is only one illness type, although Gilleskie_Eca:98 works with two illness types). An illness episode can last up to $T$ periods enumerated by $t=1,\ldots,T$; $t=0$ corresponds to the state of being well, $k=0$.
An individual who became sick makes decisions about doctor visits and and work absences. In each period $t$ of an illness, alternatives available to an employed individual who is sick are: $d_t=0$ - work and don't visit a doctor, $d_t=1$ - work and visit a doctor, $d_t=2$ - don't work and don't visit a doctor, and $d_t=3$ - don't work and visit a doctor. The utility of the agent depends on the elapsed length of the current illness $t$, the accumulated number of physician visits $v_t$, and the accumulated number of work absences $a_t$. The state variables observed by the econometrician and the agent at $t$ are $x_t=\left(t, v_t,a_t \right)$. Note that $k=1$ if and only if $t>0$, so $k$ does not appear in the definition of $x_t$.
The state variables evolve in the following way. An individual always starts with the state of being well, $x_0=(0,0,0)$. The individual contracts an illness and moves to the state $x=(1,0,0)$ with probability $\pi^S(H)$, where $H$ is a vector of exogenous indicators for health status and being between 45-64 years of age.
The accumulated number of physician visits $v_t$ and the accumulated number of illness-related absence from work $a_t$ both take values in $ \{0,1,\ldots, T-1\} $. They start from $v_1=a_1=0$ and evolve in the following way: $v_{t+1} = v_t +1(d_t=1 \text{ or } 3)$ and $a_{t+1} = a_t +1(d_t=2 \text{ or } 3)$.
In each illness period $t \in \{1,\ldots,T\}$, the individual recovers and returns to the state of being well with probability $\pi^W(x_t,d_t)$. Gilleskie parameterizes and estimates $\pi^W(x_t,d_t)$ and $\pi^S(H)$ prior to estimating the preference parameters. We use those estimates in our application.
The per-period consumption is defined as $C(x_t,d_t) = Y -\big[ PC 1(d_t=1 \text{ or } 3) + Y\big(1-L\Phi(x_t,d_t) \big)1(d_t=2 \text{ or } 3) \big]1(t>0)$, where $Y$ is the per-period labor income, $PC$ is the cost of a doctor visit, and $\Phi(x_t,d_t)=\exp(\phi_1+\phi_2a'(x_t,d_t))/[1+\exp(\phi_1+\phi_2a'(x_t,d_t))]$ is the portion of income that the sick leave coverage replaces, where $a'(x_t,d_t)$ is the value of $a_{t+1}$ given $(x_t,d_t)$. $L\in (0,1)$ is the sick leave coverage rate.
The per-period utilities can be expressed in the following form.
where $\theta_4=-\infty$ so that when $t=0$ the decision $d_t=0$ is always chosen. Since we do not restrict the location and scale of $\epsilon_t$, the values $ \theta_p$, $p=1,2,3$ can be set to arbitrary values and $\theta_5$ can be set to an arbitrary positive value in our semiparametric estimation procedure. More details on per-period utilities are given in Appendix (ref).
For data generation we use parameter values based on estimates in Gilleskie_Eca:98 for Type 2 illness with some adjustments so that the expected number of doctor visits and work absences roughly match with Gilleskie's sample. We let the data generating distribution of the utility shocks to be a two-component mixture of extreme value distributions: $ \sum_{k=1}^2 \omega_k \phi \left( \epsilon; \mu_k,\sigma_k \right) $. See Appendix (ref) for the data-generating parameter values. The panel data $\{x_{it}, d_{it},\, i=1,\ldots,n,\, t=1,\ldots,T\}$ is sequentially simulated for $n=100$ individuals and $T=8$ periods.
The priors are specified as follows, $\underbar{a}=10$, $A_m=0.05$, and $\tau=5$, $\mu_{jk} \sim N(0, 2^2)$, $\log \sigma_{k} \sim N(0, 1)$, $\log \sigma \sim N(0,0.01^2)$, and $\theta_6 \sim N(0,4^2)$. This gives normal prior on $\mu_{jk}$'s and $\theta_6$ with large variances. The log-normal prior on the component specific scale parameters also implies sufficiently large prior probabilities for large values of $\sigma_k$'s. Prior sensitivity checks presented in Appendix (ref) show that the obtained estimation results are not substantively affected by moderate changes in the prior. Appendix (ref) shows results when extreme value distribution of shocks is used for the data generation.
We use 20,000 MCMC iterations to explore the posterior distribution. Figure (ref) shows a trace plot and a p.m.f. of the number of mixture components $m$. The posterior has its peak at $m=2$ and the posterior probability of $m=1$ is small relative to its prior.
\FloatBarrier
\FloatBarrier
Figure (ref) shows the posterior densities of the utility function parameters in the location and scale normalization corresponding to the original model in Gilleskie_Eca:98 described in Section (ref). The corresponding trace plots presented in Figure (ref) in Appendix (ref) provide evidence that the MCMC algorithm converged.
The standard parametric approach that Gilleskie takes is to assume that the shocks are extreme value i.i.d. and to estimate the model by the maximum likelihood method. In Figure (ref) and Table (ref), we compare estimation results from our semiparametric method to those from the MLE. For the comparison, we use the expected number of doctor visits $E(v)$ implied by the model. It is a function of the model parameters and is of interest in the application. One of the main advantages of structural estimation is that it provides an attractive framework for counterfactual experiments. We consider a counterfactual experiment presented in Section 6.2 of Gilleskie's paper (Experiment 1). In this experiment, we are interested in the behavior of individuals when the coinsurance rate paid out of pocket is set to zero. Thus, we examine the counterfactual model solution when $PC=0$ and the transition probabilities and $(\theta,F)$ are unchanged. The last row of Table (ref) and panel (b) in Figure (ref) display the estimation results for $E(v)$ in the counterfactual environment.
In addition to the posteriors, the true values, the HPD credible intervals, and the MLE confidence intervals computed by the Delta method, Figure (ref) and Table (ref) present credible intervals for the identified sets of $E(v)$, $\hat{B}_{0.95}^{I_{E(v)}}$, that are introduced in Section (ref). While the point estimates of the increase in $E(v)$ are similar for both approaches, the 95% Bayesian credible interval for the identified set of the counterfactual $E(v)$ in the semiparametric model is up to 6 times wider than the 95% confidence interval for the MLE in the parametric setting. Importantly, the confidence interval by far misses the true counterfactual value of $E(v)$, while the credible interval for the identified set includes the true value and the HPD credible interval gets very close to it.
\FloatBarrier
\FloatBarrier
These results illustrate the following general observations. First, the dynamic logit MLE can deliver misleading point and set estimates, especially for counterfactuals, when the shocks are not extreme value distributed. Second, as is well known from MoonSchorfheide:12, the HPD credible intervals for set identified parameters might be too short at least from the classical perspective. Third, the credible intervals for the identified sets introduced in Section (ref) are more conservative than the standard HPD intervals and might be preferred under set identification. Finally, and most importantly, the whole posterior distribution (and not just point and set estimators) should be reported and taken into account in making decisions or policy recommendations.
In this section, we briefly discuss extensions of our results and methodology to models with unknown $G$, continuous $X$, multiple agents, and finite horizon.
The assumption of known transition and initial probabilities for the observed states, $G$, can be relaxed with the following notation changes. The expression for the likelihood function in (ref) needs to be replaced by \[ \sum_{i=1}^n \left( \sum_{t=1}^T \log p(d_{it} | x_{it}) + \sum_{t=2}^T \log Pr(x_{it}|x_{it-1},d_{it-1}) + \log Pr(x_{i1})\right), \] where $Pr(x_{it}|x_{it-1},d_{it-1})$ and $Pr(x_{i1})$ are the corresponding elements of $G$. Assume that a prior distribution on the elements of $G$ has a positive and continuous density at the DGP value $G_0$. To generalize the posterior consistency result in Theorem (ref) and its proof it suffices to note that the collection of CCPs $P(\theta,F,G)$ (with the explicit dependence on $G$ now) is continuous in $G$ in the present settings, see, for example Norets_ddcm_diff_cont:09. The conclusion of the theorem would be $\Pi \big( \theta,F,G: ||(G_0,P(\theta_0,F_0,G_0))-(G,P(\theta,F,G)) ||> \delta \big |D^n \big) \rightarrow 0$ almost surely. The construction of the credible sets for the identified sets in Section (ref) would use $G$ as point identified reduced form parameters along with the CCPs.
As demonstrated by SrisumaLinton2012, discretizations of continuous $X$ could nontrivially affect estimation results; hence, it would be desirable to extend our method to continuous $X$. The results in Lemma (ref) (analytical expressions for CCPs and Emax), Lemma (ref) (continuity of CCPs and Emax in $F$), Lemma (ref) (approximations by mixtures of Gumbels), and Lemma (ref) (continuity of mixtures in parameters) are proved for arbitrary/continuous $X$. A version of the posterior consistency result in Theorem (ref) would also hold for a continuous and bounded $X$ under additional assumptions that the CCPs and the transition densities for the observed states are bounded away from zero. In this case the continuity of CCPs in $F$ in the sup norm established in Lemma (ref) would deliver positive prior probability for the Kullback-Leibler neighbourhoods of the DGP, which is a sufficient condition for the Schwartz posterior consistency theorem. The norm on the CCPs in the conclusion of the theorem would be one implied by the weak topology, and it would be harder to interpret. Unknown $G$ parameterized by a finite dimensional parameter can be handled similarly to the extension discussed in the previous paragraph; treating $G$ nonparametrically for continuous $X$ would be more involved. The result on the exact matching of CCPs in Lemma (ref) would not hold for continuous $X$, but it is not essential for the proposed method. Solutions of the Bellman equations for continuous $X$ can be obtained using sieve approximations (KristensenMogensenMoonSchjerning2021) or random grids (Rust:97, ImaiJainChing:09, Norets_eca:09). However, the main challenge for implementing our method for models with continuous $X$ appears to be a fast and accurate computation of the derivatives of the CCPs and the likelihood, which we leave to future work.
Single agent models with extreme value distributed utility shocks were extended to dynamic discrete games with multiple agents in AguirregabiriaMiraGames:07 and PesendorferSchmidtDengler:08, where in the first stage choice probabilities are estimated nonparametrically and then structural parameters are estimated from equilibrium conditions and the estimated choice probabilities. A direct extension of our likelihood based method to dynamic games would have to handle the difficulties arising from possible multiple equilibria. The likelihood function would not be defined in this case (same values of shocks can be consistent with different choice probabilities). Flexible modelling of an equilibrium selection mechanism is a non-trivial problem, hence we do not pursue it in the present paper. Nevertheless, one way to exploit the ideas from our paper in the estimation of dynamic discrete games is to use the method of PesendorferSchmidtDengler:08 with several fixed location-scale mixtures of extreme value distributions instead of just one extreme value distribution to check the sensitivity of the estimation results to the distributional assumptions on the shocks.
To apply our method to models with a finite horizon, one would just need to use a backward induction method instead of a Newton method for solving the Bellman equations.
In this paper, we propose and implement a semiparametric Bayesian estimation method for dynamic discrete choice models that uses flexible mixture specifications for modeling the distribution of unobserved state variables. We establish approximation and posterior consistency results that provide frequentist asymptotic guarantees for the method. Our approach is shown to perform well in practice for binary and multinomial choice models. The computational costs of solving the dynamic program with our mixture specification are comparable to those of the dynamic logit. Even though the proposed MCMC algorithm requires many more iterations than a standard maximum likelihood method, the proposed framework is a robust and computationally tractable semiparametric alternative to the standard dynamic logit model; it provides more reliable inference, especially for counterfactuals.