The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
68,537 characters
Semiparametric Bayesian Estimation of Dynamic Discrete Choice Models
\title{Semiparametric Bayesian Estimation of Dynamic Discrete Choice Models}
\author{Andriy Norets \thanks{Department of Economics, Brown University; andriy$\_$[email removed]} and Kenichi Shimizu \thanks{Department of Economics, University of Alberta; [email removed] } \thanks{This version: \today}}
\maketitle
\begin{abstract}
We propose a tractable semiparametric estimation method for structural dynamic discrete choice models.
The distribution of additive utility shocks in the proposed framework is modeled by location-scale mixtures of extreme value distributions with varying numbers of mixture components. Our approach exploits the analytical tractability of extreme value distributions in the multinomial choice settings and the flexibility of the location-scale mixtures.
We implement the Bayesian approach to inference using Hamiltonian Monte Carlo
and an approximately optimal reversible jump algorithm.
In our simulation experiments, 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.
We develop theoretical results on approximations by location-scale mixtures in an appropriate distance and
posterior concentration of the set identified utility parameters and the distribution of shocks in the model.
\end{abstract}
\begin{keyword}
Dynamic Discrete choice, Bayesian nonparametrics, set identification, location-scale mixtures, MCMC, Hamiltonian Monte Carlo, reversible jump
\end{keyword}
\section{Introduction}
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, \citet{Rust_handbook:94} and \citet{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 \citet{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 \cite{Rust:87} and \citet{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, \citet{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. \citet{Aguirregabiria:10} shows the nonparametric identification of the shock distribution under particular assumptions on the per-period payoffs. \citet{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. \cite{BuchholzShumXu:20} provide identification results for the per-period payoffs when the observed state is continuous.
The framework of \cite{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, \citet{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, \citet{Norets_ddc_mult:11} shows that it is in principle possible to extend the method from \citet{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 \cite{Norets2017mcmc}.
Similarly to \citet{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 \cite{Rust:87}, our approach delivers estimation results that are consistent with the previous literature on
semiparametric estimation (\citet{NoretsTang2013}).
For the multinomial choice model of medical care use and work absence from \citet{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, \cite{Rust_handbook:94}, \cite{MagnacThesmar:02},
\citet{NoretsTang2013}, \cite{AbbringDaljord2020}, and \cite{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{section:setup} describes the general model setup. In Section \ref{section:semi}, we introduce our semiparametric framework. In Section \ref{section:bayesian}, we describe the Bayesian estimation method. Section \ref{section:theory} presents theoretical results. Sections \ref{section:rust} and \ref{section:gilleskie} contain the applications.
Possible framework generalizations are discussed in Section \ref{sec:extensions}.
Derivations, proofs, and implementation details are given in appendices.
\section{General Model Setup}
\label{section:setup}
In the infinite-horizon version of the model, the decision maker maximizes the expected discounted sum of the per-period payoffs
\begin{align}
\label{eq:sequence_formulation}
\max_{d_t,d_{t+1},\ldots} E_t \left(\sum_{j=0}^\infty \beta^j u\left( x_{t+j},d_{t+j},\epsilon_{t+j} \right) \right),
\end{align}
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 \cite{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 \cite{Rust:87}, \citet{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{sec:extensions}.
Under mild regularity conditions (\cite{BhattacharyaMajumdar:89}), the decision problem in \eqref{eq:sequence_formulation} admits the following Bellman representation
\begin{align}
\label{eq:emax_def}
Q(x) = \int \max_{j=0,1,\ldots,J} \bigg[ u(x,j) + \beta G_x^j Q + \epsilon_j \bigg] dF(\epsilon),
\end{align}
where $Q$ is called the \textit{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
\begin{align}
\label{eq:ccp_def}
p(d | x) = \int 1\bigg\{ u(x,d) + \beta G_x^d Q + \epsilon_d \geq u(x,j) + \beta G_x^j Q + \epsilon_j, \forall j \bigg\} dF(\epsilon).
\end{align}
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
\begin{align}
\label{eq:likelihood}
\log L(D^n) = \sum_{i,t} \log p(d_{it} | x_{it}).
\end{align}
\cite{Rust:87} proposed to solve the dynamic optimization problem by first
iterating on the Bellman equation \eqref{eq:emax_def} 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 \eqref{eq:ccp_def} and evaluate the likelihood function at a given $\left(u;\beta;G \right)$.
Alternatively, \cite{JuddSu:08} proposed to use constrained optimization to maximize the likelihood function subject to \eqref{eq:emax_def}.
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 \eqref{eq:emax_def} and \eqref{eq:ccp_def},
\begin{align}
\label{eq:dyn_logit}
p(d| x) = \frac{e^{u(x,d) + \beta G_x^d Q}}{\sum_{j=0}^J e^{u(x,j) + \beta G_x^j Q}}, \;
Q(x)= \log \sum_{d=0}^J e^{u(x,d) + \beta G_x^d Q}.
\end{align}
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, \citet{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.
\section{Semiparametric Model}
\label{section:semi}
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
\begin{equation}
\label{eq:phi_gumbel}
\phi(z; \mu,\sigma) = \prod_{j=1}^J \frac{1}{\sigma} \phi\left( \frac{z_j-\mu_j}{\sigma} \right),\, \mbox{where } \phi(z_j) = e^{-z_j-\gamma -e^{-z_j-\gamma}}
\end{equation}
is the univariate Gumbel density
and $\gamma$ is the Euler-Mascheroni constant. Some relevant properties of the Gumbel distribution are outlined in Appendix \ref{EV}.
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
\begin{align}
\label{eq:mixture}
\epsilon_t \sim \sum_{k=1}^m \omega_k \phi(\cdot ; \mu_k,\sigma_k ),
\end{align}
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
\cite{McLachlanPeel:00} and \cite{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
(\cite{Rousseau:10}, \cite{ShenTokdarGhosal2013}, and \cite{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 \eqref{eq:mixture},
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 \eqref{eq:mixture}.
Thus, it is reasonable to expect that the specification \eqref{eq:mixture} is very flexible. Indeed, in Section \ref{section:theory},
we show that it can approximate smooth multivariate densities arbitrarily well in an appropriate distance so that the
conditional choice probabilities and the \textit{Emax} function implied by the model with \eqref{eq:mixture} approximate those from the model with an arbitrary smooth density for $\epsilon_t$.
The model specification with \eqref{eq:mixture} also possesses attractive analytical properties.
If a normalization $\epsilon_{t0} = 0$ is not imposed and $(J+1)$-dimensional version of \eqref{eq:mixture} is used, then
$Q$ and $p$ could be expressed as mixtures of the appropriately recentered and rescaled expressions from the dynamic logit model \eqref{eq:dyn_logit}.
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.
\begin{lemma}\label{lm:ccp_emax}
Suppose $\epsilon \sim \sum_{k=1}^m \omega_k \phi(\cdot ; \mu_k,\sigma_k )$. Then,
\begin{equation}\label{eq:ccp}
p(d|x)=
\begin{cases}
\sum_{k=1}^m \omega_k \exp\left[-e^{-a_{kx}} \right], & \text{if }\ d=0 \\
\sum_{k=1}^m \omega_k \exp\left[\frac{u(x,d) + \beta G_x^d Q +\mu_{dk}}{\sigma_k} - A_{kx} \right] \left\{1- \exp\left[-e^{-a_{kx}} \right] \right\}
, & \text{if }\ d=1,\ldots,J;
\end{cases}
\end{equation}
\begin{equation}\label{eq:emax}
Q(x) = \sum_{k=1}^m \omega_k \sigma_k \left[A_{kx} + E1( e^{-a_{kx}} ) \right], \mbox{ where } E1(z) = \int_{z}^{\infty} e^{-t} / t dt,
\end{equation}
\begin{align*}
A_{kx} &= \log \sum_{j=1}^J \exp \left( \frac{u(x,j) + \beta G_x^j Q +\mu_{jk}}{\sigma_k} \right) \mbox{ and }
a_{kx} = \frac{u(x,0) + \beta G_x^0 Q}{\sigma_k} + \gamma - A_{kx}.
\end{align*}
\end{lemma}
The derivations of \eqref{eq:ccp} and \eqref{eq:emax} can be found in Appendix \ref{sec:ccp_emax}.
The derivatives of \eqref{eq:ccp} and \eqref{eq:emax} that are useful for the model solution and estimation are given in
Appendix \ref{sec:derivatives}.
Similarly to \cite{Rust:87}, we obtain the solution of the Bellman equation \eqref{eq:emax} by a Newton-Kantorovich method
described in Appendix \ref{sec:NK}.
\section{Inference}
\label{section:bayesian}
\subsection{Motivation and Overview of the Bayesian Approach}
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.
\subsection{Normalizations}
\label{subsection:normalizations}
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
\eqref{eq:mixture} 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 \eqref{eq:mixture}.
For comparisons of the estimation results with the dynamic logit estimates and the identified sets in
\cite{NoretsTang2013},
the parameter draws are renormalized for reporting as follows
\begin{equation}
\label{eq:renorm_params}
s \cdot \bigg ( \hat{\theta}_{1:J}^{dl} + \sum_{k=1}^m \omega_k \mu_k, \hat{\theta}_{J+1}^{dl}, \theta_{J+2:d_\theta}
\bigg),
\end{equation}
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 \cite{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.
\subsection{Priors}
Let us introduce the prior distributions for the parameters of the mixture in \eqref{eq:mixture}.
We use the following prior distributions on the number of mixture components and the mixing weights,
\begin{align}
\label{eq:m_prior}
\Pi(m) & \propto e^{-\underbar{A}_m m (\log m)^\tau}, \\
\label{eq:omega_prior}
\Pi(\omega_1,\ldots,\omega_m|m) & = \mbox{Dirichlet}(\underbar{a}/m,\ldots,\underbar{a}/m),
\end{align}
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 \eqref{eq:m_prior} and \eqref{eq:omega_prior} perform well in applications and deliver optimal posterior contraction rates in nonparametric multivariate density estimation by mixtures of normal distributions, see, for example, \cite{ShenTokdarGhosal2013} and \cite{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, \cite{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}$.
\subsection{MCMC Algorithm}
\label{sec:mcmc_details}
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 \citet{Norets2017mcmc} for simulating the number of mixture components.
HMC is a very popular and efficient MCMC algorithm; see, for example, \citet{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{sec:auxiliary_res_det}.
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 \cite{NoretsPelenis_Eca_2022}, see
\cite{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
\begin{align}
\alpha(m^\ast,m) &= \frac{p(D^n|m^\ast,\psi_{1m^\ast}) \Pi(\psi_{1m^\ast}|m^\ast)\Pi(m^\ast)}
{p(D^n|m,\psi_{1m}) \Pi(\psi_{1m}|m)\Pi(m) } \nonumber
\\
&
\label{eq:alpha_ar}
\cdot
\left(
\frac{1\{m^\ast=m+1\}}{\tilde{\pi}_m(\psi_{m+1}|\psi_{1m},Y)} + 1\{m^\ast=m-1\}\tilde{\pi}_{m-1}(\psi_{m}|\psi_{1m-1},D^n)
\right).
\end{align}
Innocuous random relabeling of mixture components increases the acceptance probability for attempts to delete $m$-th mixture component ($m^\ast=m-1$).
\cite{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.
\begin{algorithm}[h]
\setstretch{1.00}
\caption{Estimation Algorithm}
Step 1: Pre-estimate/fix the observed state transition probabilities $G$.\\
Step 2: Tune HMC hyperparameters for each number of mixture components $m=1,...,\bar{m}$, where $\bar{m}$ is a pre-specified positive integer.\\
Step 3: Initialize the parameters $m^{(0)}$ and $\psi_{1m^{(0)}}^{(0)}$.\\
Step 4: Run MCMC by iterating the following two steps.\\
\KwOutput{The posterior draws: $\left\{m^{(\ell)}, \psi_{1m^{(\ell)}}^{(\ell)} \right\}$, $\ell=1,...,L$.}
\For{$\ell \in \{0,1,...,L-1\}$}{
1) Simulate $m^{(\ell+1)}|\ldots$ using optimal reversible jump.
\begin{itemize}
\item Propose $m^\ast$ with $Pr(m^\ast=m^{(\ell)}+1|m^{(\ell)})=Pr(m^\ast=m^{(\ell)}-1|m^{(\ell)})=1/2$.
\item If $m^\ast=m^{(\ell)}+1$, then also simulate $\psi_{m^{(\ell)}+1}^{(\ell)}$ from
$\tilde{\pi}_{m^{(\ell)}+1}(\cdot|D^n,\psi_{1m^{(\ell)}}^{(\ell)})$.
\item If $m^\ast=m^{(\ell)}-1$, choose $k$ randomly from $\{1,\ldots,m^{(\ell)}\}$ and exchange values of $\psi_k^{(\ell)}$ and $\psi_{m^{(\ell)}}^{(\ell)}$.
\item With probability $\min\{1, \alpha(m^\ast,m^{(\ell)})\}$ accept $m^{(\ell+1)}= m^\ast$; otherwise, $m^{(\ell+1)}=m^{(\ell)}$.
\end{itemize}
2) Simulate $\psi_{1m^{(\ell+1)}}^{(\ell+1)}|\ldots$ using HMC.
\begin{itemize}
\item Initialize HMC algorithm by the current value $\psi_{1m^{(\ell+1)}}^{(\ell)}$
and the hyperparameters specific to $m=m^{(\ell+1)}$.
If HMC hyper parameters have not yet been tuned/obtained for $m=m^{(\ell+1)}$, then obtain and store them.
\item Perform one or several iterations of the HMC algorithm to get $\psi_{1m^{(\ell+1)}}^{(\ell+1)}$.
\end{itemize}
}
Step 5: Use MCMC draws $\left\{m^{(\ell)}, \psi_{1m^{(\ell)}}^{(\ell)} \right\}$ to conduct inference on parameters or functions of interest.
\label{alg:mcmc}
\end{algorithm}
\FloatBarrier
The Matlab code for the MCMC algorithm and replication instructions for the estimation results in the applications in Sections \ref{section:rust} and \ref{section:gilleskie} are publicly available.\footnote{\url{https://anorets.github.io/papers/mix_ddcm_code.zip}}
\pagebreak
\section{Approximation Results and Asymptotics}
\label{section:theory}
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 \cite{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.
\subsection{Approximation Results}
\label{sec:approx_results}
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 \textit{Emax} function and the conditional choice probabilities are continuous in that distance as shown in the following lemma.
\begin{lemma}
\label{lm:LipschitzCont_Q_p}
Suppose (i) $|u(x,j)|\leq \bar{u}<\infty$ for all $x \in X$ and $j=0,1,\ldots,J$;
(ii) under $F$, the density for $\epsilon_j -\epsilon_d$ is bounded for all $j \ne d$;
(iii) under $F$, $E(|\epsilon_j|)$ is finite for all $j$.
Then, the \textit{Emax} function and the conditional choice probabilities are locally Lipschitz continuous in $F$,
\[
\sup_x |Q(x;F)-Q(x;\tilde{F})| \leq C \cdot \rho(F,\tilde{F}),
\]
\[
\sup_{d,x} |p(d|x;F)-p(d|x;\tilde{F})| \leq C^\prime \cdot \rho(F,\tilde{F}),
\]
where constants $C$ and $C^\prime$ depend on $\beta$, $\bar{u}$ and the bounds on the densities and moments in conditions (ii)-(iii).
\end{lemma}
The lemma holds irrespective of whether the innocuous normalization $\epsilon_{t0}=0$ is imposed.
Its proof is given in Appendix \ref{sec:proofs}.
The following lemma shows that densities satisfying smoothness and finite moment conditions can be approximated by mixtures of Gumbels in distance $\rho$.
\begin{lemma}
\label{lm:approx_by_mgumbles}
Let $f$ be a density on $\mathbb{R}^I$ satisfying a moment existence condition
\[\int ||\mu||_2 f(\mu) d\mu < \infty,
\]
and a smoothness condition
\begin{equation}
\label{eq:f_smooth}
|f(z+h) - f(z) | \leq ||h||_2 L_f(z) e^{\tau ||h||_2},
\end{equation}
for some $\tau>0$ and an envelope function $L_f(\cdot)$ such that
\begin{equation}
\label{eq:L_inf_smooth}
\int (1+ |z_i| ) L_f(z) dz < \infty , i=1,\ldots,I.
\end{equation}
Then, for any $\delta>0$, there exist $(m,\omega,\mu,\sigma)$
where
$m \in \mathbb{Z}^+$,
$\omega_j \in [0,1]$ with $\sum_{j=1}^m \omega_j=1$,
$\mu_j \in \mathbb{R}^I$, $j=1,\ldots,m$, and
$\sigma>0$
such that
\[
\rho\left(f( \cdot ), \sum_{j=1}^m \omega_j \phi(\ \cdot \ ; \mu_j, \sigma) \right) < \delta.
\]
\end{lemma}
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{sec:proofs}. The proof uses only smoothness and tail conditions on $\phi$ that are shown to hold for Gumbel densities
in Lemmas \ref{lm:lemma_bounded_integral} and \ref{lm:Lipschitz_Gumbel_location} in Appendix \ref{sec:proofs}.
Thus, Lemma \ref{lm:approx_by_mgumbles} 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.
\begin{lemma}
\label{lm:Fcont_params}
Let $F^1$ and $F^2$ denote two mixtures of Gumbel densities on $\mathbb{R}^{J}$ with densities
$
f^i(\epsilon) = \sum_{k=1}^m \omega_k^i \phi\left( \epsilon; \mu_k^i, \sigma^i \right)
$. Then, for a given $\delta>0$ and $F^1$, there exists $\tilde{\delta}>0$ such that
for any $F^2$ with parameters satisfying:
$|\sigma^1-\sigma^2|< \tilde{\delta}$, $|\omega_k^1-\omega_k^2|< \tilde{\delta}$, and
$|\mu_k^1-\mu_k^2|< \tilde{\delta}$, $k=1,\ldots,m$, we have
$\rho ( F^1, F^2 ) < \delta$.
\end{lemma}
\subsection{Posterior Consistency}
\label{sec:post_cons}
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$.
\begin{theorem}
\label{th:post_cons}
Let $(\theta_0,F_0)$ be the data generating values of parameters.
Suppose
(i) The observed state space is finite, $X=\{1,\ldots,K\}$;
(ii) $G^d$, $d=0,\ldots,J$ and the distribution of the initial observed state $x_{i1}$ are known and fixed;
(iii) $\forall x \in X$, $\exists t \in \{1,\ldots,T\}$, such that $Pr(x_{it}=x)>0$;
(iv) $u(x,d; \theta)$ is continuous in $\theta$;
(v) $F_0$ satisfies the conditions of Lemma \ref{lm:approx_by_mgumbles};
(vi) For any $\delta>0$, $\Pi(B_\delta(\theta_0))>0$, where $B_\delta(\theta_0)$ is a ball with radius $\delta$ and center $\theta_0$;
(vii) For any $\delta>0$, positive integer $m$, $\mu_k \in \mathbb{R}^J$, $\sigma_k>0$, $w_k \geq 0$, $k=1,\ldots,m$, $\sum_{k=1}^m w_k=1$,
$\Pi(B_\delta(\mu_1,\sigma_1, \ldots, \mu_m,\sigma_m, w_1, \ldots,w_{k-1})|m)>0$.
Then, for any $\delta>0$,
\[
\Pi \big( \theta,F: ||P(\theta_0,F_0)-P(\theta,F) ||> \delta \big |D^n \big) \rightarrow 0 \mbox{ almost surely.}
\]
\end{theorem}
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 \cite{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, \cite{Norets_ddcm_diff_cont:09}; and, thus, Lipschitz continuity of $P(\theta,F)$ in $F$ from Lemma \ref{lm:LipschitzCont_Q_p} delivers continuity of $P(\theta,F)$ in $(\theta,F)$.
The finite mixture approximation result in Lemma \ref{lm:approx_by_mgumbles},
the continuity of $P(\theta,F)$ in $(\theta,F)$, the continuity of
finite mixtures in parameters in Lemma \ref{lm:Fcont_params},
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{sec:extensions}.
Theorem \ref{th:post_cons} 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,
\cite{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{th:post_cons},
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
\begin{align*}
&\Pi\bigg(\theta \in B_{\tilde{\delta}}(\theta_0), F \in B_{\tilde{\delta}}(F_0) \bigg| P \in B_{\delta}(P(\theta_0,F_0))\bigg) =
\frac{\Pi(\theta \in B_{\tilde{\delta}}(\theta_0), F \in B_{\tilde{\delta}}(F_0))}{\Pi\big (P \in B_{\delta}(P(\theta_0,F_0))\big)}
\\
&\geq
\Pi\big(\theta \in B_{\tilde{\delta}}(\theta_0), F \in B_{\tilde{\delta}}(F_0)\big)>0,
\end{align*}
which suggests that the conditional prior would not rule out the data generating parameter values.
\subsection{Exact Matching of CCPs}
\label{sec:exact_match_ccps}
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.
\begin{lemma}
\label{lm:exact_match_ccps}
Suppose
(i) The observed state space is finite, $X=\{1,\ldots,K\}$;
(ii) $(\theta_0,F_0)$ are the data generating values of parameters;
(iii) $F_0$ has finite first moments and a density that is positive on $\mathbb{R}^J$.
Then there exists a finite mixture of Gumbels $F$ such that
$P(\theta_0,F_0)=P(\theta_0,F)$. An upper bound on the number of mixture components in $F$ depends only on $K$ and $J$.
\end{lemma}
The result in Lemma \ref{lm:exact_match_ccps} 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{sec:proofs}.
Lemma \ref{lm:exact_match_ccps} can be used to relax the smoothness assumptions on $F_0$ in
the posterior consistency results of Theorem \ref{th:post_cons}. Specifically,
conditions (v) in Theorem \ref{th:post_cons} can be replaced by conditions (iii) in Lemma \ref{lm:exact_match_ccps}; in the proof,
the approximation results in Lemma \ref{lm:approx_by_mgumbles} can be replaced
by the exact CCPs matching results in Lemma \ref{lm:exact_match_ccps}.
Nevertheless, the approximation results in Lemma \ref{lm:approx_by_mgumbles} 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{sec:post_cons}.
Finally, while asymptotically the number of mixture components is bounded for the exact matching in Lemma \ref{lm:exact_match_ccps} and has to increase to infinity for the approximation results in Lemma \ref{lm:approx_by_mgumbles}, 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).
\subsection{Inference for Identified Sets}
\label{sec:inference_IS}
\citet{NoretsTang2013} and generalizations of their results to the multinomial case in \citet{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.
\cite{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. \citet{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.
\cite{KlineTamer2016} further study this approach and \cite{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 \citet{NoretsTang2013} and \cite{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
\begin{equation}
\label{eq:B_I_eta}
B_{1-\alpha}^{I_\eta} = \bigcup_{P \in B_{1-\alpha}^P} I_\eta(P).
\end{equation}
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 \eqref{eq:B_I_eta} 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{section:gilleskie}.
\section{Application I: Rust's Binary Choice Model}
\label{section:rust}
\citet{NoretsTang2013} propose a method for computing identified sets for parameters in dynamic binary choice models and apply their method to the \citet{Rust:87}'s model. In this section, we show that our semiparametric model can also recover the identified set for that model.
\subsection{\citet{Rust:87}'s Optimal Bus Engine Replacement Problem}
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$.
\citet{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 \citet{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 \cite{NoretsTang2013} for a given fixed vector of CCPs.
Estimation results for simulated data are presented for the multinational choice application in Section \ref{section:gilleskie}.
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 \eqref{eq:renorm_params}.
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.
\begin{align*}
& \pi(m=k) \propto e^{-\underbar{A}_m k(\log k)^\tau}, \, \underbar{A}_m=0.05, \, \tau=5, \\
& (\omega_1,\ldots,\omega_m) \sim \mbox{Dirichlet}(\underbar{a}/m,\ldots,\underbar{a}/m), \, \underbar{a}=10, \\
&\mu_{k} \sim 0.5N(2.5,1^2)+0.5N(-3,7^2),\\
& \log \sigma_{k} \sim 0.4N(0,1^2)+0.6N(-6,1^2), \; \log \sigma \sim N(0,0.01^2).
\end{align*}
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{fig_Rust_id} 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{fig_Rust_id} 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.
\begin{figure}[ht]
\begin{subfigure}[b]{0.3\linewidth}
\centering
\includegraphics[width=\linewidth]{Rust_id_prior}
\caption{prior}
\end{subfigure}
\begin{subfigure}[b]{0.3\linewidth}
\centering
\includegraphics[width=\linewidth]{Rust_id_N3}
\caption{$N=3$}
\end{subfigure}
\begin{subfigure}[b]{0.3\linewidth}
\centering
\includegraphics[width=\linewidth]{Rust_id_N10}
\caption{$N=10$}
\end{subfigure}
\caption{\small{
Utility parameter draws.
Prior (left). Posterior draws for $N=3$ (center) and $N=10$ (right).
500,000 MCMC iterations. The first 100,000 draws discarded as burn-in.
Every 10th draw is shown.
The true identified set is shown by solid black lines.
The red dot corresponds to the point-identified utility parameter values under the logistic distribution for $\epsilon$.
}}
\label{fig_Rust_id}
\end{figure}
To assess the convergence of the MCMC algorithm, consider Figure \ref{fig_Rust_id3D} 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{fig_Rust_trace} in Appendix \ref{sec:rust_implement_details} shows additional evidence of MCMC convergence including trace plots of $m$, $\sum_{k=1}\omega_k\mu_k$, and other parameters.
\begin{figure}[ht]
\begin{subfigure}[b]{0.4\linewidth}
\centering
\includegraphics[width=\linewidth]{Rust_id3D_N3}
\caption{$N=3$}
\end{subfigure}
\begin{subfigure}[b]{0.4\linewidth}
\centering
\includegraphics[width=\linewidth]{Rust_id3D_N10}
\caption{$N=10$}
\end{subfigure}
\caption{\small{
Utility parameter draws. The z-axis represents the MCMC iterations.
The true identified set is shown in black.
}}
\label{fig_Rust_id3D}
\end{figure}
\FloatBarrier
\section{Application II: Gilleskie's Multinomial Model of Medical Care Use and Work Absence}
\label{section:gilleskie}
In this section, we illustrate our methodology using a multinomial choice model of medical care use and work absence from \citet{Gilleskie_Eca:98}.
An extension of \citet{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.
\subsection{The model}\label{model}
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 \citet{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$.
\subsubsection{Alternatives}
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$.
\subsubsection{State variable transitions}
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.
\subsubsection{Utility}
\label{sec:utility}
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.
\begin{alignat}{4}
u(d_t=1,x_t,\epsilon_t)
&= \theta_1 +&&\theta_4 1(t=0) &&+\theta_6 C(x_t,1)1(t>0) &&+ \epsilon_{t1}, \nonumber \\
u(d_t=2,x_t,\epsilon_t)
&= \theta_2 +&&\theta_4 1(t=0) &&+\theta_6 C(x_t,2)1(t>0) &&+\epsilon_{t2},\nonumber \\
u(d_t=3,x_t,\epsilon_t)
&= \theta_3 +&&\theta_4 1(t=0) &&+\theta_6 C(x_t,3)1(t>0) &&+\epsilon_{t3}, \nonumber \\
u(d_t=0,x_t,\epsilon_t)
&= &&\theta_5 1(t=0) &&+ \theta_6 C(x_t,0)1(t>0) &&+\epsilon_{t0}, \nonumber
\end{alignat}
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{sec:utility_derivation}.
\subsection{Estimation}
\label{sec:gill_estimation}
For data generation we use parameter values based on estimates in \citet{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{sec:giil_data_gen} 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{sensitivity} show that the obtained estimation results are not substantively affected by moderate changes in the prior.
Appendix \ref{gilleskie_logit} shows results when extreme value distribution of shocks is used for the data generation.
\subsubsection{Estimation Results and Counterfactuals}\label{results}
We use 20,000 MCMC iterations to explore the posterior distribution.
Figure \ref{fig_Gilleskie_m} 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
\begin{figure}[ht]
\begin{subfigure}[b]{0.4\linewidth}
\centering
\includegraphics[width=0.8\linewidth]{Gilleskie_trace_m}
\caption{\small{Trace plot of $m$ }}
\end{subfigure}
\begin{subfigure}[b]{0.4\linewidth}
\centering
\includegraphics[width=0.8\linewidth]{Gilleskie_pmf_m}
\caption{\small{p.m.f. of $m$ }}
\end{subfigure}
\caption{\small{Trace plot and p.m.f. of $m$ }}
\label{fig_Gilleskie_m}
\end{figure}
\FloatBarrier
Figure \ref{fig_Gilleskie_posterior_theta} shows the posterior densities of the utility function parameters in the
location and scale normalization corresponding to the original model in \cite{Gilleskie_Eca:98} described
in Section \ref{subsection:normalizations}. The corresponding trace plots presented in
Figure \ref{fig_Gilleskie_trace} in Appendix \ref{appsec:traceplots} provide evidence that the MCMC algorithm converged.
\begin{figure}[ht]
\begin{subfigure}[b]{0.7\linewidth}
\centering
\includegraphics[width=\linewidth]{Gilleskie_posterior_theta}
\end{subfigure}
\caption{\small{Posterior densities of utility parameters (solid) and the data generating values (dashed). }}
\label{fig_Gilleskie_posterior_theta}
\end{figure}
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{fig_Gilleskie_EvEa} and Table \ref{table_estimation}, 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{table_estimation} and panel (b) in Figure \ref{fig_Gilleskie_EvEa} 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{fig_Gilleskie_EvEa} and Table \ref{table_estimation}
present credible intervals for the identified sets of $E(v)$, $\hat{B}_{0.95}^{I_{E(v)}}$, that are introduced in Section \ref{sec:inference_IS}.
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.
\begin{figure}[ht]
\begin{subfigure}[b]{0.5\linewidth}
\centering
\includegraphics[width=\linewidth]{Gilleskie_Ev}
\caption{\small{Mean doctor visits }}
\end{subfigure}
\begin{subfigure}[b]{0.5\linewidth}
\centering
\includegraphics[width=\linewidth]{Gilleskie_Ev_cf}
\caption{\small{ Mean c.f. doctor visits }}
\end{subfigure}
\caption{\small{Actual and counterfactual (c.f.) expected number of doctor visits.
Posterior density (blue solid),
95\% confidence interval (red dashed),
95\% HPD interval (blue dotted),
$\hat{B}^{I_{E(v)}}_{0.95}$ (light blue dash-dotted),
and data generating value (black solid). }}
\label{fig_Gilleskie_EvEa}
\end{figure}
\FloatBarrier
\begin{table}[h]
\centering
\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}
\begin{tabular}{l*{6}{c}}
\toprule
& \multicolumn{3}{c}{MLE under logit} & \multicolumn{3}{c}{Semiparametric Bayes} \\
\cmidrule(lr){3-4}\cmidrule(lr){5-7}
& \multicolumn{1}{c}{True}
& \multicolumn{1}{c}{Est.} & \multicolumn{1}{c}{95\%CI (length)}
& \multicolumn{1}{c}{Est.} & \multicolumn{1}{c}{95\%HPD (length)}
& \multicolumn{1}{c}{$\hat{B}^{I_{E(v)}}_{0.95}$ (length)} \\
\midrule
\addlinespace
$E(v)$ & 1.49 & 1.35 & [1.15, 1.55] (0.40) & 1.37 & [1.17, 1.58] (0.41) & [0.99, 1.84] (0.85) \\
\addlinespace
$c.f.E(v)$ & 2.31 & 1.70 & [1.46, 1.94] (0.48) & 1.72 & [1.28, 2.18] (0.90) & [1.10, 3.99] (2.88) \\
\addlinespace
\bottomrule
\end{tabular}
\caption{Estimated (Est.) $E(v)$ and $c.f.E(v)$:
the MLE with its 95\% confidence interval and
the posterior mean with the 95\% HPD interval and $\hat{B}^{I_{E(v)}}_{0.95}$. }
\label{table_estimation}
\end{table}
\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 \cite{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{sec:inference_IS} 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.
\section{Extensions}
\label{sec:extensions}
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 \eqref{eq:likelihood} 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{th:post_cons} 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 \cite{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{sec:inference_IS} would use $G$ as point identified reduced form parameters along with the CCPs.
As demonstrated by \cite{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{lm:ccp_emax} (analytical expressions for CCPs and \textit{Emax}), Lemma \ref{lm:LipschitzCont_Q_p} (continuity of CCPs and \textit{Emax} in $F$), Lemma \ref{lm:approx_by_mgumbles} (approximations by mixtures of Gumbels), and Lemma \ref{lm:Fcont_params} (continuity of mixtures in parameters) are proved for arbitrary/continuous $X$. A version of the posterior consistency result in Theorem \ref{th:post_cons} 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{lm:LipschitzCont_Q_p} 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{lm:exact_match_ccps} 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 (\cite{KristensenMogensenMoonSchjerning2021}) or random grids (\cite{Rust:97}, \cite{ImaiJainChing:09}, \cite{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 \cite{AguirregabiriaMiraGames:07} and \cite{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 \cite{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.
\section{Conclusion}
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.
\bibliographystyle{latex_common/ecta}
\bibliography{latex_common/allreferences}