Robust Structural Estimation under Misspecified Latent-State Dynamics

Robust Structural Estimation under Misspecified Latent-State Dynamics

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

\onehalfspacing

Introduction

Dynamic structural models are useful tools for counterfactual policy analysis in various fields of economics. The dynamic process of latent variables is a key feature of these models, as it captures the persistence of unobserved factors that affect agents' decisions over time. Examples of potentially serially correlated latent variables include product characteristics in demand estimation (nair2007intertemporal,schiraldi2011automobile,gowrisankaran2012dynamics), search costs in consumer search (koulayev2014search), firm productivity in trade (piveteau2021empirical), patent profitability in optimal stopping (pakes1984patents), quality in technology adoption (de2019subsidies), health shocks in insurance (fang2021life), and beliefs about ability in labor economics (miller1984job,arcidiacono2025college).

Assumptions about the dynamic process governing the serial dependence of latent variables are central to the estimation and counterfactual analysis in dynamic structural models. These assumptions capture agents' uncertainty about the future, such as a consumer's uncertainty about future product characteristics in demand estimation. The misspecification of these assumptions can lead to biased estimates of future continuation value, which in turn bias counterfactual predictions (e.g., welfare and elasticity). However, the direction and magnitude of this bias are unclear, because the models are dynamic and nonlinear. This raises the need for sensitivity analysis of empirical results to these distributional assumptions.

In this paper, we propose a framework to quantify this sensitivity by perturbing a reference dynamic process to compute bounds on a scalar parameter of interest. The scalar parameter (e.g., welfare and elasticity) is a function of model primitives, such as model parameters, the distribution of latent variables, and the value function in dynamic structural models. In our proposed framework, the bounds on this scalar parameter are the solutions to constrained optimization problems whose feasible region (identified set) is defined by moment conditions for estimation, structural constraints (e.g., the Bellman equation), and the perturbation set around the reference transition distribution. The distribution that achieves the bound is called the worst-case distribution.

A central challenge is to define the perturbation set in a way that simplifies computation while maintaining key structural features of the distribution of latent variables. For time-homogeneous models, the structural feature is the stationarity condition, i.e., the perturbed dynamic process must be stationary. For finite-horizon, time-inhomogeneous models, the perturbed trajectory of latent variables must be Markovian. In addition, the terminal distribution is fixed because it can be nonparametrically point identified (lewbel2000semiparametric,matzkin2007nonparametric). Because the stationarity and Markov conditions are functional constraints imposed directly on the distribution being optimized, they are computationally difficult to impose. We contribute to the sensitivity analysis literature (e.g., christensen2023counterfactual) by providing a computationally tractable framework to deal with these constraints.

There are further practical challenges. First, the properties (e.g., closed-form, smoothness) of the worst-case distribution are typically unclear, complicating the choice of approximation methods. Second, the expectations that define the scalar parameter, the model-implied moments, and the structural constraints are all calculated with respect to the perturbed distribution. Because this distribution is itself an optimization variable that changes during the optimization, the numerical integration can be difficult to implement. Third, the value function in dynamic structural models is infinite-dimensional due to the serial dependence of latent variables, which complicates the optimization problem further.

To address these challenges, we define the perturbation set as a Kullback–Leibler (KL) divergence ball around the reference distribution. For the time-homogeneous case, the reference distribution is a joint distribution of current and future latent variables. For the time-inhomogeneous case, it is the joint distribution of the entire trajectory of latent variables. The KL radius controls the size of the perturbation set. We then employ the Optimal Transport (OT) framework to impose structural constraints on the distribution of latent variables. In the time-homogeneous case, the stationarity condition requires that the marginal distributions of current and future latent variables of the perturbed joint distribution coincide, which can be imposed using OT. This marginal constraint, together with the KL divergence penalty from the problem’s Lagrangian, yields a computationally tractable Entropic Optimal Transport (EOT) problem. In the time-inhomogeneous case, the perturbed trajectory must be Markovian, and the terminal distribution is fixed. We show how to formulate this as an EOT problem.

The EOT problem can be solved using its dual formulation, which has three key advantages over the primal formulation. First, it provides a closed-form expression for the worst-case distribution and characterizes its smoothness. During our proposed optimization algorithm, this closed-form expression is used to update the value function in dynamic structural models. Second, the Sinkhorn algorithm (sinkhorn1967concerning,cuturi2013sinkhorn) allows us to solve the EOT problem and compute the worst-case distribution efficiently. Finally, because the expectations in the dual problem are taken with respect to the reference distribution, an appropriate numerical integration method can be chosen in advance.

Then we consider three complementary sensitivity measures to interpret the results. First, the global sensitivity approach computes the largest deviation from the reference value. It progressively increases the KL radius until the bounds flatten. We show that it provides a tractable approximation to the nonparametric bounds on the scalar parameter when the KL divergence constraint is removed. Moreover, we derive an explicit upper bound on the approximation error, which enables us to control the error within a desired level. Second, the local sensitivity approach analyzes the effect of small perturbations. It computes the right derivatives of the bounds with respect to the KL radius at zero, which serves as our local sensitivity measure. Finally, the robustness metric approach, inspired by spini2021robustness, computes the smallest deviation from the reference distribution required to produce sensitive results. It is the smallest KL divergence from the reference distribution required for the scalar parameter's value to fall below a user-specified threshold (e.g., 5% below the reference value). In addition to our three sensitivity measures, we can also estimate an alternative model and set the radius as the KL divergence between the alternative and reference models.

For large sample properties, we propose an estimator for the bound, establishing its consistency and convergence rate. To this end, we first establish the consistency and convergence rate of the estimator of the identified set, following chernozhukov2007estimation. We then derive the asymptotic distribution of the plug-in estimator by proving the Hadamard directional differentiability of the bound.

We apply our framework to an infinite-horizon dynamic demand model for new cars in the UK, Germany, and France. We consider the sensitivity of the price elasticity and consumer surplus from an additional \$3,000 electric vehicle subsidy. In the model, the indirect utility of purchasing is a latent variable due to unobserved product characteristics, and its transition is typically modeled as an AR(1) process (e.g., schiraldi2011automobile,gowrisankaran2012dynamics). For the price elasticity, we find that the French market is the least sensitive to the distributional assumption (at most 6.20% deviation from the reference elasticity), while the German market is the most sensitive (at most 15.24% deviation). We also find that this sensitivity is relatively stable over time for all three markets. For the consumer surplus, the German market is also the most sensitive (at most 102.75% deviation from the reference consumer surplus), while the UK and French markets are less sensitive (at most 25.17% and 24.73% deviation, respectively). Importantly, the results remain economically meaningful even under the worst-case distribution, with the consumer surplus remaining at least \$309 million for the German market, \$1,243 million for the French market, and \$2,584 million for the UK market.

We also apply our framework to a finite-horizon dynamic labor supply model for taxi drivers in New York City. We consider the sensitivity of the elasticity of stopping work and the Frisch elasticity of labor supply. In the model, the market-level supply shock is a latent variable, and its transition is also modeled as an AR(1) process. For the elasticity of stopping work, we find that weekday drivers' elasticity is more sensitive in the morning, while weekend drivers' elasticity is more sensitive in the afternoon. For the Frisch elasticity, both weekday and weekend drivers' elasticities are sensitive to the distributional assumption, with at most 76.83% and 42.84% deviation from the reference elasticity, respectively.

Related Literature

Identification and estimation. Many papers have focused on identification in a range of dynamic structural models including finite mixture models (kasahara2009nonparametric,luo2022identification,higgins2023identification,higgins2025learning), unobservable Markov processes (hu2012nonparametric), and counterfactual conditional choice probabilities in dynamic binary choice models (norets2014semiparametric). berry2023instrumental uses the generalized instrumental variable approach for unobserved state variables in dynamic discrete choice (DDC) models. In fixed effects DDC models, aguirregabiria2021sufficient considers identification of structural parameters using sufficient statistics, while aguirregabiria2024identification studies identification of average marginal effects. hwang2024identification employs proxy variables for the latent variables. kalouptsidi2021linear proposes the Euler Equations in Conditional Choice Probabilities (ECCP) estimator. By leveraging finite-dependence properties and cross-sectional data, it identifies structural parameters in the presence of serially correlated market-level unobserved variables without distributional assumptions. arcidiacono2011conditional adapts the Expectation-Maximization algorithm to estimate DDC models with discrete latent types. norets2009inference extends the Bayesian estimation of DDC models by imai2009bayesian to allow for serially correlated latent variables, while blevins2016sequential proposes a sequential Monte Carlo method. chiong2016duality shows identification in DDC models is an optimal transport problem under serial independence of utility shocks. In addition, chen2011sensitivity, schennach2014entropic, and fan2023partial,fan2025partial consider inference of finite-dimensional parameters in the presence of an infinite-dimensional parameter, namely the distribution of latent variables.

Our framework contributes to this literature in three ways. First, we focus directly on a scalar parameter (e.g., welfare, elasticity) rather than on the identified set of model primitives. Second, to analyze sensitivity, we consider a perturbation set around the reference distribution rather than the set of all distributions. Third, we complement identification strategies that do not rely on distributional assumptions. In the labor supply application, we use the ECCP estimator of kalouptsidi2021linear to point identify the utility parameters and then conduct sensitivity analysis of the labor supply elasticity with respect to assumptions about the dynamic process of the market-level supply shock.

Sensitivity analysis and robustness. In DDC models, kalouptsidi2021counterfactual,kalouptsidi2021identification relax common normalizations on the utility function. bugni2019inference considers the local misspecification of the transition density of observable variables, assuming the transition density is correctly specified in the limit (as the sample size goes to infinity) in DDC models. andrews2017measuring considers a setting in which the moments are locally misspecified under the reference distribution. Subsequent work armstrong2021sensitivity constructs near-optimal confidence intervals in such models. kitamura2013robustness considers the robust estimation under moment restrictions. bonhomme2022minimizing perturbs the reference model and assumes the size of the perturbation shrinks to zero as the sample size goes to infinity. chen2024robust relaxes the rational expectation assumption. gu2023dual considers the identification of scalar counterfactual parameters using optimal transport. spini2021robustness studies the robustness of policy effects to changes in the distribution of covariates. armstrong2025misspecification provides a selective review of misspecification in econometrics. Most closely related is christensen2023counterfactual, who conducts sensitivity analysis with respect to parametric assumptions about the distribution of latent variables in structural models. However, they focus on relaxing the marginal distribution assumption while maintaining serial independence.

We contribute to this literature in three ways. First, our focus on misspecified dynamic processes complements previous analysis of misspecification in dynamic structural estimation (e.g., bugni2019inference,kalouptsidi2021counterfactual,kalouptsidi2021identification,christensen2023counterfactual). Second, the size of our perturbation set is fixed and does not shrink with the sample size. Third, we provide a computationally tractable dual formulation for the optimization problem. While prior work (schennach2014entropic,gu2023dual,christensen2023counterfactual) uses duality to convert the infinite-dimensional problem to a finite-dimensional one, our dual problem is still infinite-dimensional. This is because the stationarity condition for time-homogeneous models and fixed terminal distribution for time-inhomogeneous models are infinite-dimensional, and serial dependence leads to an infinite-dimensional value function in dynamic structural models. However, by leveraging EOT duality, we provide a tractable implementation and demonstrate our approach through two empirical applications.

Distributionally robust optimization (DRO). The literature on DRO (kuhn2019wasserstein,rahimian2019distributionally,blanchet2022confidence,gao2023distributionally,wang2021sinkhorn) usually studies the uncertainty due to limited observability of data, noisy measurements, or estimation errors.

Our work is distinct in that we study the misspecification due to assumptions about the serial dependence of latent variables—a problem of model specification rather than data limitation. Moreover, because our framework can be treated as a moment-constrained DRO problem, we employ the minimax theorem (see fan1953minimax and ricceri1998minimax, Theorem 1.3) to exchange the order of the supremum over Lagrangian multipliers (for the moment conditions and structural constraints) and the infimum over distributions in the perturbation set. This exchange allows a direct application of duality results from the DRO literature, which simplifies our proof of the dual formulation significantly.

Applied work. Building on the seminal work on the estimation of DDC models (rust1987optimal,hotz1993conditional,aguirregabiria2002swapping,aguirregabiria2007sequential,pesendorfer2008asymptotic,arcidiacono2011conditional), most applied work assumes the serial independence of utility shocks. In the presence of serially correlated latent variables, parametric models are often used, as seen in the works of schiraldi2011automobile,gowrisankaran2012dynamics,blevins2018firm,piveteau2021empirical and others mentioned in the introduction.

We contribute to this literature by developing a computationally tractable framework for sensitivity analysis of scalar parameters of interest to these distributional assumptions. We also demonstrate our framework through two empirical applications: an infinite-horizon dynamic demand model, and a finite-horizon dynamic labor supply model.

Outline: Sections (ref) and (ref) present our framework for time-homogeneous and time-inhomogeneous models. (ref) establishes the large sample properties of our estimator. (ref) introduces three sensitivity measures. (ref) discusses practical implementation. Sections (ref) and (ref) present two empirical applications. (ref) concludes. Appendix (ref) presents additional examples. All proofs are in Appendix (ref).

Notation: Let $U \in \mathcal{U} \subseteq \mathbb{R}^{d}$ be a vector of latent variables with support $\mathcal{U}$ where $\mathcal{U}$ is assumed to be Polish. Let $\mathcal{P}(\mathcal{U})$ be the space of Borel probability measures on $\mathcal{U}$ and $\mathcal{B}(\mathcal{U})$ be the Borel $\sigma$-algebra on $\mathcal{U}$. In this paper, all measures are assumed to be absolutely continuous with respect to the Lebesgue measure. Denote by $\mathbb{E}_{F}\left[\cdot\right]$, $\mathbb{E}_{x}\left[\cdot\right]$ the expectations with respect to $F \in \mathcal{P}(\mathcal{U})$ and the probability distribution of the random variable $x$, respectively. For variables in a stationary dynamic context, a prime (e.g., $x'$) denotes the variable's value in the next period. Let $F_{1}, F_{2} \in \mathcal{P}(\mathcal{U})$, we write $F_{1} \ll F_{2}$ if $F_{1}$ is absolutely continuous with respect to $F_{2}$, and $F_{1} \otimes F_{2}$ as the product measure. For a finite-dimensional vector, denote by $\|\cdot\|_{p}$ the $p$-norm. Denote by $L^{p}(F)$ the space of functions for which $\int |f|^{p} dF < \infty$. For a set $\mathcal{A}$, let $\text{int}(\mathcal{A})$ denote its interior. Denote by $\mathbb{R}_{+}$ the set of non-negative real numbers, and $\mathbb{N}$ the set of natural numbers. Finally, let $\mathcal{J} := \left\{1,\cdots,J\right\}$.

Methodology for Time-Homogeneous Models

This section presents our methodology for time-homogeneous models. (ref) defines the perturbation set. (ref) gives two examples. (ref) presents the general framework and derives duality results. (ref) discusses how to perturb the stationary distribution.

Definition of Perturbation Set

We partition a vector of latent variables $U \in \mathcal{U} \subseteq \mathbb{R}^{d}$ into $2 \leq k \leq d $ subvectors, i.e., $U = (U_{1},\cdots,U_{k})$.\footnote{The vector $U$ can also contain observable variables.} Each subvector $U_{i} \in \mathcal{U}_{i} \subseteq \mathbb{R}^{d_{i}}$ has a marginal distribution $\nu_{i} \in \mathcal{P}(\mathcal{U}_{i})$ for $i = 1,\cdots,k$. Let $F_{0}$ denote the reference distribution for $U$. The perturbation set around this reference distribution is defined as:

where $\Pi(\nu_{1},\cdots,\nu_{k})$ is the set of joint distributions on $\mathcal{U}$ with marginals $\{\nu_{i}\}_{i=1}^{k}$, and $\delta \geq 0$ measures the “size” of the perturbation set, defined by the Kullback-Leibler (KL) divergence:

In time-homogeneous models, the marginal distribution constraints are used to impose the stationarity. Consider an unobserved stationary first-order Markov process $\{\xi_t\}_{t \in \mathbb{Z}}$ with state space $\Xi \subseteq \mathbb{R}^{d_\xi}$. The process is stationary with respect to $\nu_{0} \in \mathcal{P}(\Xi)$ if $\int F_{0}(d\xi'|\xi) \nu_{0}(d\xi) = \nu_{0}(d\xi')$ for all $\xi' \in \Xi$, where $F_{0}(d\xi'|\xi)$ is the reference transition kernel (e.g., conditional Gaussian distribution for an AR(1) process). This is equivalent to requiring that the marginal distributions of the joint distribution $dF_{0}(\xi,\xi') := F_{0}(d\xi'|\xi)\nu_{0}(d\xi)$ are both $\nu_{0}$, i.e., $F_{0} \in \Pi(\nu_{0},\nu_{0})$. Moreover, for any $F \in \Pi(\nu_{0},\nu_{0})$, its conditional density $F(d\xi'|\xi)$ preserves $\nu_{0}$ as a stationary distribution. For this example, let $U := (\xi,\xi')$. Then, the perturbation set is:

This definition allows us to perturb the transition kernel of the Markov process while keeping its stationary distribution unchanged. It can introduce non-linear dynamics into the latent variable process. For example, if the reference model is an AR(1) process, the perturbation set includes any nonlinear first-order Markov process with the same stationary distribution as the AR(1) process. (ref) discusses how to perturb the stationary distribution. In this case, we replace the condition $F \in \Pi(\nu_{0},\nu_{0})$ with $F \in \Pi(\nu,\nu)$, where $\nu$ is the perturbed stationary distribution that is in a neighborhood of $\nu_{0}$.

The marginal constraints with the KL constraint form a computationally tractable EOT problem whose implementation depends on its dual formulation (see Section (ref)).\footnote{For duality results of general divergence constrained OT problem, see bayraktar2025stability.}

Examples

Our framework applies to a variety of latent variables, including utility shocks, productivity characteristics, labor supply shocks, etc. This section focuses on a parametric model for latent variables beyond utility shocks in infinite and finite-horizon DDC models. Appendix (ref) considers: (i) serial independence of utility shocks in DDC models, and (ii) serial independence of consumption shocks in dynamic discrete-continuous choice models. The bounds on a scalar parameter are solutions to constrained optimization problems over the perturbation set subject to structural constraints (e.g., Bellman equation) and moment conditions.

Framework and Duality

We now present a general framework that nests the above examples. In general, the model is not point-identified when $\mathcal{F}$ is not a singleton.\footnote{For example, see schennach2014entropic,molinari2020microeconometrics.} Therefore, we propose to compute upper and lower bounds on the outcome of interest. Let the scalar parameter of interest, $s : \mathcal{U} \times \Theta \times \mathcal{V} \rightarrow \mathbb{R}$, be a function of the latent variable $U$, and the model primitives $(\theta,v) \in \Theta \times \mathcal{V}$. The lower bound is the solution to the following optimization problem:

where $ \mathcal{F} := \{ F \in \mathcal{P}(\mathcal{U}) \mid F \in \Pi(\nu_{1},\cdots,\nu_{k}), D_{KL}(F\|F_{0}) \leq \delta \}$.

The first constraint is a moment condition where the moment function $m : \mathcal{U} \times \Theta \times \mathcal{V} \rightarrow \mathbb{R}^{d_{P}}$ is finite-dimensional as we assume the observable variable $X \in \mathcal{X}$ has discrete support and stack the moment functions for each $x \in \mathcal{X}$. The second constraint is a structural constraint defined by $\psi : \mathcal{U} \times \Theta \times \mathcal{V} \times \mathcal{G} \rightarrow \mathbb{R}$ that is linear in the Lagrange multiplier function $g \in \mathcal{G}$ for a given $(\theta,v)$. Finally, $v \in \mathcal{V}$ is the solution to the structural constraint.

The \hyperref[Primal]{\textcolor{blue}{\ref*{Primal}}} problem can be intractable due to the optimization over $\mathcal{F}$. First, the properties (e.g., closed-form and smoothness) of the optimal $F^{*}$ are typically unknown, complicating the choice of approximation methods. Second, the marginal distribution conditions are functional constraints imposed directly on the distribution being optimized, which are computationally difficult to impose. Third, expectations are taken with respect to the perturbed distribution, making numerical integration difficult.

To overcome these issues, we derive the \hyperref[Dual problem]{\textcolor{blue}{\ref*{Dual problem}}} problem corresponding to the \hyperref[Primal]{\textcolor{blue}{\ref*{Primal}}}. The \hyperref[Dual problem]{\textcolor{blue}{\ref*{Dual problem}}} provides the closed-form of the optimal $F^{*}$, and characterizes its smoothness. Moreover, in the dual, the expectation is taken with respect to the reference distribution $F_{0}$. (ref) proposes a computationally tractable algorithm that utilizes the optimal $F^{*}$. To motivate the duality, consider the Lagrangian of the \hyperref[Primal]{\textcolor{blue}{\ref*{Primal}}}:

where $c(U;\theta,v,g,\lambda) := s(U;\theta,v) + \lambda^{T}m(U;\theta,v) + \psi(U;\theta,v,g)$, $\lambda \in \mathbb{R}^{d_{P}}$ is the Lagrange multiplier for the moment condition and $\lambda_{KL}$ is the Lagrange multiplier for the KL divergence constraint. For given $(\theta,v)$, under regularity conditions, we can swap the order of the infimum over $F$ and the supremum over $(\lambda,\lambda_{KL},g)$. Then, we can rewrite (ref) as:

The inner infimum is the Entropic Optimal Transport (EOT) problem (for $\lambda_{KL} > 0$)\footnote{See nutz2021introduction for a comprehensive introduction to the EOT problem. The EOT problem has close connection to the static Schr{\"o}dinger Bridge problem. It is called the Optimal Transport (OT) problem when $\lambda_{KL} = 0$ (see villani2009optimal). The optimal value is called the optimal OT value.}:

where $c(U;\theta,v,g,\lambda)$ is the cost function, and $\mathcal{C}(\theta,v,g,\lambda,\lambda_{KL})$ is called the optimal EOT value. The EOT problem is computationally fast to solve using the Sinkhorn algorithm, which relies on the duality of the EOT problem. Moreover, the closed-form and smoothness of the unique solution $F^{*}$ to the EOT problem can be derived from its duality (see (ref) and (ref)). (ref) reviews the EOT duality and the Sinkhorn algorithm. Next, we impose assumptions for the minimax theorem to swap the order of infimum and supremum, and the EOT duality to hold:

Assumptions (ref)-(ref) are mild. (ref) holds for indicator functions. There is no particular necessity to write $\rho(U)$ in log-density form in (ref). Our notation is chosen to simplify the expression in Assumptions (ref) and (ref). (ref) imposes the growth rate condition. It is satisfied for all Examples in (ref) if we assume $u$, $g$, and $v$ satisfy the growth rate condition. It ensures that $c(U;\theta,v,g,\lambda) \in L^{1}(F)$ for $\forall \ F \in \mathcal{F}$, and can also be used to show the convergence of the Sinkhorn algorithm and the convergence of optimal EOT value to optimal OT value as $\lambda_{KL} \downarrow 0$ (eckstein2022quantitative,eckstein2024convergence).\footnote{For the convergence of optimal EOT value to optimal OT value, lower semicontinuity alone is not sufficient (nutz2021introduction Example 5.1). A sufficient condition is the continuity condition on the cost function.}

(ref) establishes the duality and provides the closed-form of $F^{*}$. $\phi_{i}$ is the test function for the marginal distribution condition $\nu_{i}$ for $i = 1,\cdots,k$. The optimal $\left\{\phi_{i}^{*}\right\}_{i=1}^{k}$ are the optimal EOT potentials, which can be efficiently obtained using the Sinkhorn algorithm.

The \hyperref[Dual problem]{\textcolor{blue}{\ref*{Dual problem}}} is computationally tractable. First, it provides the closed-form of the worst-case distribution $F^{*}$, which can be used in the \hyperref[Primal]{\textcolor{blue}{\ref*{Primal}}} problem even if we want to solve it directly. In (ref), we propose an iterative algorithm that alternates between solving the EOT to obtain the worst-case distribution and updating $v$ by solving the structural constraint with the worst-case distribution. Second, for given $(\theta,v,g,\lambda,\lambda_{KL})$, the optimal $\left\{\phi_{i}^{*}\right\}_{i=1}^{k}$ (and thus $F^{*}$) can be efficiently computed using the Sinkhorn algorithm, which is computationally very fast. Third, the expectations in the \hyperref[Dual problem]{\textcolor{blue}{\ref*{Dual problem}}} are taken with respect to the marginal distributions and the reference distribution. Therefore, numerical integration methods can be determined in advance. Finally, we have:

(ref) shows the smoothness of $\left\{\phi_{i}^{*}\right\}_{i=1}^{k}$ for $\lambda_{KL} > 0$ (the EOT case). For $\lambda_{KL} = 0$ (the OT case), it is not straightforward to obtain the smoothness of the worst-case distribution (see villani2009optimal Chapter 12).

Perturbation of Stationary Distribution

This section discusses how to perturb the stationary distribution in the time-homogeneous setting. We consider the serial dependence of the latent variables, as detailed in the example in Section (ref). Let $U:=(\xi,\xi')$ be the vector of current and future latent variables. We define the perturbation set for the stationary distribution as $\mathcal{N} = \{\nu \in \mathcal{P}(\Xi) \mid D_{KL}(\nu\|\nu_{0}) \leq \delta_{1}\}$ where $\delta_{1} \geq 0$. The perturbation set for the joint distribution is:

Let $\kappa_{stationary}(\delta_{1},\delta,P)$ denote the lower bound on the scalar parameter of interest under this perturbation set. Under regularity conditions, we can swap the order of infimum over $F$ and the supremum over $(\lambda,\lambda_{KL},g)$:

where $\mathcal{C}(\theta,v,g,\lambda,\lambda_{KL},\nu)$ is the EOT problem with respect to the perturbed stationary distribution $\nu$: $\mathcal{C}(\theta,v,g,\lambda,\lambda_{KL},\nu) := \inf_{F \in \Pi(\nu,\nu)} \mathbb{E}_{F} \left[c(\xi,\xi';\theta,v,g,\lambda)\right] + \lambda_{KL} D_{KL}(F\|F_{0})$. Its dual formulation is:

Under regularity conditions, we can swap the order of infimum over $\nu$ and the supremum over $(\phi_{1},\phi_{2})$. Then, the inner infimum over $\nu$ is:

which is a KL-divergence distributionally robust optimization problem (see hu2013kullback,rahimian2019distributionally) whose dual formulation is:

where $\eta$ is the Lagrange multiplier for the KL divergence constraint for $\nu$. To summarize:

(ref) shows the dual formulation when the stationary distribution is perturbed. Although the Sinkhorn algorithm does not apply directly, we can sequentially update $(\phi_{1},\phi_{2})$ using the first-order optimality conditions like the Sinkhorn algorithm (see (ref) for a review of the Sinkhorn algorithm).

Methodology for Time-Inhomogeneous Models

This section extends our framework to the time-inhomogeneous setting. We begin by defining the perturbation set for these models in (ref). (ref) provides an example of a finite-horizon dynamic discrete choice model. (ref) presents the duality result. Finally, (ref) discusses how to perturb the initial distribution.

Definition of Perturbation Set

Consider a sequence of latent variables over a finite horizon, $U:=(\xi_{1},\xi_{2},\cdots,\xi_{T})$, that follows a first-order Markov chain. The reference joint distribution is the product of an initial distribution and transition kernels:

where $\nu_{1}(d\xi_{1})$ is the initial distribution, and $F_{t}(d\xi_{t+1}|\xi_{t})$ is the transition kernel from period $t$ to $t+1$. Let $\nu_{T}(d\xi_{T})$ be the terminal distribution implied by this process.

We consider perturbing the reference distribution while holding its initial and terminal distributions fixed, i.e.,

where $\Pi_{\text{Markov}}(\nu_{1},\nu_{T})$ is the set of all joint distributions over $\mathcal{U}$ that satisfy the first-order Markov property\footnote{That is, for any $t \leq T-1$, $F(d\xi_{t+1}|\xi_{1},\cdots,\xi_{t}) = F(d\xi_{t+1}|\xi_{t})$ almost surely under $F$.} and have $\nu_1$ and $\nu_T$ as their initial and terminal marginal distributions. Our perturbation set allows for any transition dynamics of the latent variables between the initial and terminal periods, as long as the overall process remains Markovian and the initial and terminal distributions are fixed.

We will discuss how to perturb the initial distribution in (ref). The terminal distribution can often be nonparametrically identified; therefore, we fix it. For instance, in a finite horizon DDC model, the terminal period is a static discrete choice problem where the distribution of the latent variable can be identified (lewbel2000semiparametric,matzkin2007nonparametric). Moreover, if $\xi_{t}$ is the market-level latent variable, and the model has the finite dependence property (see arcidiacono2011conditional)\footnote{For example, a model with a terminating action has the finite dependence property.}, then the utility parameters can be identified without a distributional assumption for the latent variables (see kalouptsidi2021linear), which in turn identifies the terminal distribution.

Example

Framework and Duality

The bound for the time-inhomogeneous case, $\kappa_{\text{TI}}(\delta, P)$, is defined similarly to the time-homogeneous case, but with the optimization performed over $\mathcal{F}_{\text{Markov}}$:

where “TI" stands for time-inhomogeneous. To solve \hyperref[eq: nonstationary]{\textcolor{blue}{\ref*{eq: nonstationary}}}, we follow the procedure in (ref). However, the set $\mathcal{F}_{\text{Markov}}$ is not necessarily convex, which prevents the proof strategy of the minimax duality in (ref) from being applied directly. Therefore, we propose to solve a relaxed problem where the Markov property condition is removed. We then show that, under certain reasonable conditions, the solution to the relaxed problem is Markovian, thereby also solving the original problem. The perturbation set for the relaxed problem is defined as:

where $\Pi(\nu_{1},\nu_{T})$ is the set of joint distributions with initial distribution $\nu_{1}$ and terminal distribution $\nu_{T}$. The relaxed problem is given by:

whose Lagrangian is:

where $c(U;\theta,v,g,\lambda) := s(U;\theta,v) + \lambda^{T}m(U;\theta,v) + \psi(U;\theta,v,g)$, $\lambda \in \mathbb{R}^{d_{P}}$ is the Lagrange multiplier for the moment condition and $\lambda_{KL}$ is the Lagrange multiplier for the KL divergence constraint. For given $(\theta,v)$, under regularity conditions, we can swap the order of the infimum over $F$ and the supremum over $(\lambda,\lambda_{KL},g)$. Then, we can rewrite (ref) as:

The inner infimum is:

which can be rewritten as the discrete-time dynamic Schr{\"o}dinger Bridge (SB) problem (see leonard2013survey,de2021diffusion). Because it only restricts the initial and terminal distributions, we can decompose it into two parts: the two-period marginal distribution part (the first and last period) and the conditional distribution part (the intermediate variables conditional on the first and last period). The second part is unconstrained, thereby having a closed-form solution. The first part is the static SB (or EOT) problem whose duality is similar to (ref). We impose the following assumptions for the minimax duality, decomposition, static SB duality, and the Markov property:

The boundedness condition in (ref) is stronger than (ref), which does not guarantee that $c(U;\theta,v,g,\lambda) \in L^{1}(F)$ for any $F \in \mathcal{F}_{\text{relaxed}}$. (ref) is a sufficient condition for the SB duality to hold. Finally, (ref) is the key to the Markov property of the solution to the \hyperref[eq: relaxed nonstationary]{\textcolor{blue}{\ref*{eq: relaxed nonstationary}}} problem, and is satisfied in (ref). It does not hold if the moment function depends on the entire path of the latent variables.

(ref) shows the duality for the \hyperref[eq: relaxed nonstationary]{\textcolor{blue}{\ref*{eq: relaxed nonstationary}}} problem. The difference between (ref) and (ref) is that (ref) does not fix the intermediate marginal distributions. Therefore, we can decompose (ref) into the sum of two-period marginal distribution ($F^{1,T}_{0}$) part, and the conditional distribution (the distribution of $(\xi_{2},\cdots,\xi_{T-1})$ given $\xi_{1}$ and $\xi_{T}$) part. The latter part is an unconstrained optimization problem, thereby having a closed-form solution. The first part is the static SB problem, whose duality is given by (ref).

(ref) is crucial for the solution to have the Markov property. Under this assumption, the cost function $c(U;\theta,v,g,\lambda)$ is also pairwise additive. Therefore, the density ratio in (ref) has the Markov property.\footnote{It can be treated as the (unnormalized) pairwise Markov random field wainwright2008graphical.} If there exists one optimal $\lambda_{KL}^{*} > 0$, then the \hyperref[eq: nonstationary]{\textcolor{blue}{\ref*{eq: nonstationary}}} problem is equivalent to the \hyperref[eq: relaxed nonstationary]{\textcolor{blue}{\ref*{eq: relaxed nonstationary}}} problem.

The \hyperref[eq: relaxed nonstationary]{\textcolor{blue}{\ref*{eq: relaxed nonstationary}}} problem can also be solved using the iterative algorithm proposed in (ref). There is an additional step to obtain the two-period auxiliary reference distribution $R_{1,T}$. The \hyperref[eq: nonstationary EOT]{\textcolor{blue}{\ref*{eq: nonstationary EOT}}} can also be solved using the Sinkhorn algorithm.

Perturbation of Initial Distribution

Let $\mathcal{N}$ be a convex closed set around the initial distribution $\nu_{1}$, e.g., $\mathcal{N} = \{\nu \in \mathcal{P}(\Xi) \mid D_{KL}(\nu\|\nu_{1}) \leq \delta_{1}\}$ where $\delta_{1} \geq 0$. Then, the perturbation set is defined as:

Let $\tilde{\kappa}_{\text{TI,Initial}}(\delta,P)$ be the lower bound on the scalar parameter for the \hyperref[eq: relaxed nonstationary]{\textcolor{blue}{\ref*{eq: relaxed nonstationary}}} problem with the perturbation set $\mathcal{F}_{\mathcal{N},\text{Relaxed}}$. The following minimax duality similar to (ref) holds:

To solve (ref), as shown in the proof of (ref) and (ref), we can decompose the problem into two parts: the two-period marginal distribution part, and the conditional distribution part. The first part requires solving:

where $\tilde{\Pi}(\nu,\nu_{T})$ is the set of distributions of $(\xi_{1},\xi_{T})$ whose marginal distributions are $\nu$ and $\nu_{T}$, respectively. It is equivalent to solving:

The inner infimum is an EOT problem. Let $\text{EOT}(\nu,\nu_{T})$ be its optimal value.

(ref) shows that (ref) is a convex optimization problem with respect to $\nu$. Moreover, $\phi^{*}$ is a result of the Sinkhorn algorithm, which can be used to search the optimal $\nu$ efficiently.

Large Sample Properties

This section establishes the large sample properties of the estimator for the bound. (ref) proposes a consistent estimator and shows its convergence rate. (ref) establishes the asymptotic distribution of the plug-in estimator for the bound.

Consistency and Convergence Rate

The bound on the scalar parameter is the projection of the identified set defined by the moment conditions and structural constraints onto the scalar parameter. We follow chernozhukov2007estimation to propose an estimator for the identified set and show its consistency and convergence rate. Let $P_{n}$ be an estimator for $P_{0}$ where $n$ is the sample size. Denote by $\epsilon_{n} \in \mathbb{R}_{+}$ the tolerance level for the moment conditions that goes to zero at a suitable rate as $n \to \infty$. Our estimators $\kappa(\delta,P_{n},\epsilon_{n})$, and $\tilde{\kappa}_{\text{TI}}(\delta,P_{n},\epsilon_{n})$ for the bounds replace the moment conditions by the approximate moment conditions.\footnote{To compute the estimator, the number of moment conditions is doubled due to the use of approximate moment conditions. The duality is similar to Theorems (ref) and (ref), thus is omitted for brevity.}

(ref) is mild. (ref) holds in single-agent DDC models. It rules out dynamic games with multiple equilibria. (ref) is also mild as the identified set for $\theta$ is usually nonempty under the reference distribution $F_{0}$. Moreover, if the identified set for $\delta = +\infty$\footnote{In this case, the KL divergence constraint is replaced by the absolute continuity constraint, i.e., $F \ll F_{0}$.} is nonempty, then (ref) implicitly assumes the radius $\delta$ is large enough. The smallest radius such that the identified set is nonempty can be estimated (see (ref)). (ref) implies that the identified set and its estimator are compact as $\mathcal{F}$ is compact (see (ref)) and thus $\mathcal{A}$ is compact (see (ref)).

Under (ref), the estimator for $\mathcal{A}_{I}$ is defined as:

The analysis of the consistency and convergence rate uses the Hausdorff Distance:

where $d(\alpha_{1},\mathcal{A}_{1}) := \inf_{\alpha_{2} \in \mathcal{A}_{2}} d(\alpha_{1},\alpha_{2})$ and $d(\alpha_{1},\alpha_{2})$ is a metric on $\mathcal{A}$.

(ref) is mild as we assumed the observable variable has discrete support, e.g., $P_{n}$ can be the frequency estimator. In practice, we can set $c_{n} \propto \log n$. Then, the convergence rate in (ref) is $\sqrt{n}$-consistent up to a logarithmic factor. We show some properties of the identified set and its estimator:

By the extreme value theorem, the infimum is achieved if the scalar parameter is continuous on $\mathcal{A}$, i.e., (ref) holds. Therefore, the optimization problem has a solution. Next, we impose the polynomial minorant condition as in chernozhukov2007estimation for the convergence rate of the estimator:

(ref) establishes the $\sqrt{n}$-consistency up to a logarithmic factor (if $c_{n} \propto \log n$). Then, we impose the following continuity assumption on the scalar parameter of interest:

Asymptotic Distribution

This section establishes the asymptotic distribution of $\kappa(\delta,P_{n})$ and $\tilde{\kappa}_{\text{TI}}(\delta,P_{n})$. To this end, we first show the Hadamard directional differentiability of $\kappa(\delta,P)$ and $\tilde{\kappa}_{\text{TI}}(\delta,P)$ with respect to $P$ at $P_{0}$ similar to christensen2023counterfactual. We begin with the definition of Hadamard directional differentiability:

Under (ref), we can restate the optimization problem as:

where $P(\alpha) := \mathbb{E}_{F} \left[m(U;\theta,v(\alpha))\right]$. Moreover, the identified set $\mathcal{A}_{I}$ is nonempty, which means the feasible set for the optimization problem is nonempty. By (ref), $\mathcal{A}_{I}$ is compact. Under (ref), the Extreme Value Theorem (see rudin1976principles Theorem 4.16.) implies that the infimum is attained. Denote by $\mathcal{A}_{I,\text{TH}}^{*}$, $\mathcal{A}_{I,\text{TI}}^{*}$ the nonempty sets of optimizers for the problems $\kappa(\delta,P_{0})$ and $\tilde{\kappa}_{\text{TI}}(\delta,P_{0})$, respectively.

To establish Hadamard directional differentiability of $\kappa(\delta,P)$ and $\tilde{\kappa}_{\text{TI}}(\delta,P)$ at $P_{0}$, we impose assumptions similar to those in bonnans2013perturbation Theorem 4.25.\footnote{We work on the primal problem to show the Hadamard directional differentiability, while christensen2023counterfactual works on the dual problem (see their Theorem 6.2).}

(ref) shows the asymptotic distribution of the bound's estimator. To conduct inference, we may follow the procedure in fang2019inference. In addition, the numerical delta method hong2018numerical combined with our practical implementation in (ref) can be used to overcome the computational challenge.

Interpreting the Results

In practice, we can estimate an alternative (parametric) model and set the radius to be the KL divergence between the alternative and the reference distribution. In addition, this section considers three complementary sensitivity measures to interpret the results: global sensitivity, local sensitivity, and robustness metric.

Global Sensitivity

The global sensitivity\footnote{See christensen2023counterfactual Theorem 2.1 for similar results. However, their results are silent about how large the radius should be so that the bounds are close to the nonparametric bounds.} approach progressively increases the radius until the bounds flatten. We show that it provides a computationally tractable approximation to the nonparametric bounds when the KL divergence constraint is removed. Moreover, we provide an explicit upper bound on the approximation error. We focus on the time-homogeneous case, for which the “nonparametric" perturbation set is:

After applying the minimax duality, we need to solve the following problem:

where the inner problem is an OT problem, which is computationally challenging in high-dimensional settings. The EOT is a computationally tractable approximation to the OT problem. Recall $\mathcal{C}(\theta,v,g,\lambda,\lambda_{KL}) := \inf_{F \in \Pi(\nu_{1},\cdots,\nu_{k})} \mathbb{E}_{F} \left[c(U;\theta,v,g,\lambda)\right] + \lambda_{KL} D_{KL}(F\|F_{0})$.

(ref) provides an explicit upper bound on the approximation error of $\mathcal{C}(\theta,v,g,\lambda,\lambda_{KL})$ to $\mathcal{C}(\theta,v,g,\lambda,0)$. For DDC models, the constants $L$ and $C$ can be explicitly characterized under additional conditions (see eckstein2022quantitative Lemma 3.5, and eckstein2024convergence Remark 2.1.) The upper bound strictly decreases to zero as $\lambda_{KL} \downarrow 0$. Therefore, we can choose a sufficiently small $\lambda_{KL}$ (or sufficiently large $\delta$) to achieve a desired accuracy for the approximation. Our framework thus approximates the nonparametric bounds in a computationally tractable way with an explicitly quantifiable approximation error.

Local Sensitivity

The local sensitivity\footnote{See bartl2021sensitivity Theorem 2.2 and christensen2023counterfactual Page 276 for similar results.} approach computes the right derivative of the bounds at $\delta = 0$, which measures the effect of a small perturbation of the reference distribution on the bounds. We show the right differentiability of the bounds with respect to $\delta$. Define:

where $\|\cdot\|_{Lip}$ is the Lipschitz constant, and $C_{3},C_{4},L$ are positive constants. We assume:

(ref) shows the right differentiability. We can also compute the derivative of the length of the bounds. In practice, we may need to compute it numerically due to the optimization over the set of optimizers and the Lagrange multipliers.

The Robustness Metric

The robustness metric is the smallest deviation from the reference distribution that can lead to sensitive results (spini2021robustness). In practice, we begin by estimating a reference scalar parameter, $\hat{s}_{F_{0}}$, under the reference distribution. If the perturbed scalar parameter $s_{F}$ is below a certain threshold, e.g., $\bar{s} = 0.95 \cdot \hat{s}_{F_{0}}$, then we may be concerned about the robustness of the results. The robustness metric is defined as:

where $\bar{s} \in \mathbb{R}$ is a user-specified threshold. The optimization problem searches for a distribution $F$ in the identified set that results in $\mathbb{E}_{F} \left[s(U;\theta,v)\right] \leq \bar{s}$ and is the closest to the reference distribution $F_{0}$ in terms of KL divergence.

We can plot the bounds against $\delta$ and then find the radius corresponding to $\bar{s}$. Alternatively, we can compute it directly by solving (ref) whose duality results are given by:

(ref) provides the dual formulation for computing the smallest radius in the time-inhomogeneous case. The equivalence holds as the regularization parameter is 1.

Practical Implementation

This section presents the practical implementation of the proposed framework. (ref) reviews the entropic optimal transport problem and the Sinkhorn algorithm. (ref) proposes a computationally feasible algorithm.

Entropic Optimal Transport and Sinkhorn Algorithm

This section reviews the Sinkhorn algorithm for the entropic optimal transport problem.\footnote{See sinkhorn1967concerning, cuturi2013sinkhorn and nutz2021introduction.} Let $(\mathcal{U}_{i}, \nu_{i})$ for $i=1, \dots, k$ be probability spaces, where $\mathcal{U}_{i}$ is the support for the random variable $U_{i}$. For a cost function $c: \mathcal{U}_{1} \times \dots \times \mathcal{U}_{k} \to \mathbb{R}$, the entropic optimal transport problem\footnote{If $F_{0} \neq F_{\otimes}$, then (ref) reformulates the problem as the (ref) problem.} with regularization parameter $\lambda_{KL} > 0$ is defined as:

whose dual is given by:

where $\phi_{i}$ is the test function for the marginal distribution constraint $\nu_{i}$. The dual problem is a concave maximization problem over $\{\phi_{i}\}_{i=1}^{k}$. The worst-case distribution is given by:

where the optimizers $(\phi_{1}^{*}, \dots, \phi_{k}^{*})$ are known as the optimal EOT potentials (also called Schr\"odinger potentials), which are the solutions to the Schr\"odinger equation (SE):

where $F_{\otimes, -i}$ is the product measure of all marginals except for the $i$-th marginal. The Schr\"odinger equations (ref) to (ref) can be interpreted as the variational first-order conditions for optimality (see nutz2021introduction Remark 3.4). Moreover, they also characterize the marginal constraints. To see this, define:

The marginal density can be obtained by integrating out the other marginals; therefore:

The Sinkhorn algorithm can be interpreted as a coordinate ascent scheme for the optimization problem (ref). It is a computationally fast\footnote{For its convergence rate, see peyre2019computational, carlier2022linear, and eckstein2022quantitative.}, iterative method for solving (ref)-(ref).

Proposed Algorithm

For given $(\theta, v, g, \lambda, \lambda_{KL})$, the EOT problem provides the worst-case distribution, $F^{*}$. This allows us to update $v$ by solving the structural constraint with $F^{*}$. We therefore propose an iterative algorithm to solve the minimax problem. The algorithm proceeds by alternating between updating $(\theta, v)$ and the dual (Lagrange multiplier) variables, $(g, \lambda, \lambda_{KL})$. After initializing all parameters, each iteration $t$ involves the following steps:

The computational cost per iteration mainly comes from solving the EOT problem and the structural constraint, which are both computationally fast. However, the number of iterations required for convergence can be much larger, as our optimization problem is a minimax problem that is potentially non-differentiable.

Empirical Application: Infinite Horizon DDC

This section applies our framework to an infinite-horizon dynamic demand for new cars in the UK, France, and Germany. Due to the unobserved product characteristics, the indirect utility of purchasing is the latent variable. To estimate the price elasticity and conduct welfare analysis of electric vehicles (EV) subsidy, we require a distributional assumption for the latent variable to solve the Bellman equation. Existing literature often uses an AR(1) process (e.g., schiraldi2011automobile,gowrisankaran2012dynamics), which may be misspecified. For example, the indirect utility may exhibit nonlinear dynamics. Therefore, we conduct a sensitivity analysis with respect to this reference distribution.

Data

We use the trim-level\footnote{In the automobile industry, a trim-level refers to a specific version of a vehicle model that comes with a particular set of features, options, and styling elements.} data from IHS Markit during the period from 2014 to 2023. The monthly level dataset contains sales, list price, and characteristics of car models in the UK, France and Germany, which are treated as three independent markets in our analysis. To construct the final dataset, we first aggregate data from the trim-level to the model-level. Then, we aggregate fuel types into: petrol, diesel, electric, and hybrid. We remove models whose total sales during the data period are less than 20,000.\footnote{In addition, we exclude car-month observations with sales below 150 units in the German market and below 5 units in the French market.} Finally, we adjust list prices by subtracting EV subsidies. The initial market size for January 2014 is calculated by subtracting the number of registered cars from the total population of each country. The market size is then updated each subsequent month by subtracting the total number of cars sold in the preceding period.

Table (ref) presents the summary statistics. The three markets offer around 141–215 products from around 23 to 30 brands per month. In terms of average sales per model, Germany has 81,401 units, closely followed by France with 81,357 units, and the UK with 70,682 units. The average price is around \$33,444 in the UK, \$27,830 in France, and \$36,117 in Germany.

The Model

The model is infinite horizon. At each month $t$, a consumer $i$ chooses $j \in \mathcal{J}_{t} \bigcup \{0\}$ where $\mathcal{J}_{t}$ is the set of available cars at $t$, and 0 is the outside option of not purchasing. Each car $j \in \mathcal{J}_{t}$ is characterized by a vector of observable characteristics $x_{jt}$, price $p_{jt}$, and an unobserved characteristic $\xi_{jt}$. The period utility of choosing $j$ is given by:

where $\varepsilon_{it}$ is a vector of i.i.d. type I extreme value utility shocks, and $x_{t},p_{t},\xi_{t}$ are the vectors of observable characteristics, prices, and unobserved characteristics for all cars in $\mathcal{J}_{t}$.

We assume a purchase is a terminating decision, i.e., consumers exit the market after the purchase. The conditional value function of purchasing car $j$ can be written as the sum of the current period utility and the flow utility after purchase:

where $\beta = 0.975$ is the discount factor. The inclusive value of purchasing is defined as:

Following schiraldi2011automobile and gowrisankaran2012dynamics, we assume:

Under the IVS assumption, $\omega_{t}$ is the only state variable, and the value function $V(\omega)$ is the solution to the smoothed Bellman equation:

where $v_{0}(\omega) = \beta \mathbb{E}\left[V(\omega')|\omega\right]$ and $v_{1}(\omega) = \omega$ are the conditional value functions of not purchasing and purchasing, respectively. The market share of car $j$ at time $t$ is given by:

First-Stage Estimation

In each market, the car with the highest total sales is set as the reference product, denoted by $r$.\footnote{The reference cars are Volkswagen Golf (Petrol) in Germany, Peugeot 208 (Petrol) in France, and Ford Fiesta (Petrol) in the UK. The reference car for each country is always available in that country's market.} Taking the log-odds ratio for cars $j$ and $r$ at time $t$ yields:

where $\Delta p_{jt} := p_{jt} - p_{rt}$, $\Delta x_{jt} := x_{jt} - x_{rt}$, and $\Delta \xi_{jt} := \xi_{jt} - \xi_{rt}$. The parameters $(\alpha,\theta)$ are identified\footnote{The intercept and reference fixed effects are not identified from (ref); they are absorbed into $\xi_{rt}$.} by the BLP instruments berry1993automobile. Moreover, $\Delta \xi_{jt}$ can be recovered by fitting the relative market share $\frac{s_{jt}}{s_{rt}}$, while the unobserved characteristic of the reference car, $\xi_{rt}$, cannot.

The exogenous characteristics, $x_{jt}$, include vehicle log-weight, log-horsepower, brand fixed effects, SUV fixed effect, and fuel type fixed effects.\footnote{The instruments are the exogenous product characteristics, average log-weight and log-horsepower of competitors' products, the proportion of competitors' products, the proportion of hybrid cars squared, and the number of brands. A competitor product is defined as a car whose brand is not that of $r$ or $j$.} Table (ref) presents the regression results. Price coefficients are negative and significant in all markets, ranging from -0.158 in France to -0.192 in Germany. Relative to petrol vehicles, EVs are preferred in the UK (0.019) and France (0.004), but not in Germany (0.000). Hybrid vehicles are valued positively across all three markets, with coefficients ranging from 0.029 in Germany to 0.037 in France. In contrast, diesel has a negative coefficient in the UK (-0.015) but positive effects in Germany (0.008) and France (0.006).

The Reference Distribution and Scalar Parameters of Interest

We first define the reference distribution and then introduce the scalar parameters of interest. After the first stage estimation, we can calculate $\omega_{t}$ up to $\xi_{rt}$ as:

As we have identified the utility parameters, the potential sensitivity of the empirical results solely arises from the distributional assumption on $\omega_{t}$.

The reference transition of $\omega_{t}$ to solve the Bellman equation is an AR(1) process:

where $\eta_{t}$ follows an i.i.d. normal distribution with mean 0 and variance $\sigma^{2}$.

The parameters $(\gamma_{0}, \gamma_{1}, \sigma^{2})$ are estimated using an iterative procedure. We begin with an initial guess of $(\gamma_{0}, \gamma_{1}, \sigma^{2})$ and circulate between: (i) solving the Bellman equation (ref), (ii) recovering $\{\omega_{t}\}_{t=1}^{T}$ from the market share of purchasing (the first part of (ref)), (iii) updating $(\gamma_{0}, \gamma_{1}, \sigma^{2})$ by refitting an AR(1) process (ref) until we find a fixed point. The reference distribution $F_{0}$ for $(\omega,\omega')$ is the product of the transition kernel of the estimated AR(1) process and its stationary distribution $\nu_{0}$. The perturbation set is defined as:

We consider two scalar parameters: (i) the industrywide price elasticity of demand, (ii) the welfare analysis of an additional EV subsidy. For both cases, the transition of $\omega_{t}$ is unchanged, i.e., we assume consumers' beliefs about the transition of $\omega_{t}$ stay the same.

For industrywide price elasticity at period $t_{1}$, we consider a 1% increase in the price of all cars. The future $\omega_{t_{1}+1}$ is conditional on:

and the industrywide price elasticity at time $t_{1}$ is:

where $s_{0}(\omega'_{t_{1}})$ is the model-implied market share of not purchasing.

For EV subsidy at $t_{1}$, we consider an additional 3,000 USD subsidy. Denote by $\mathcal{J}_{t_{1},\text{EV}}$ the set of EVs at $t_{1}$. The future $\omega_{t_{1}+1}$ is conditional on:

and the consumer surplus from the subsidy is given by:

where $M_{t_{1}}$ is the market size at time $t_{1}$.\footnote{The cost of the subsidy is $\text{Cost} = 3000 \cdot M_{t_{1}} \cdot \sum_{j \in \mathcal{J}_{t_{1},\text{EV}}} \frac{\exp\left(v_{j}(x_{jt_{1}},p_{jt_{1}}-3000,\xi_{jt_{1}})\right)}{\exp(V(\omega^{\text{EV}}_{t_{1}}))}$. }

Sensitivity Analysis

Our framework requires the constraints to be linear in $F$, while the Bellman equation (ref) is not. We first reformulate it to fit into our framework. By the Hotz-Miller Inversion Lemma (hotz1993conditional), we have:

Taking the log-odds ratio of purchasing and not purchasing, and using (ref), we have:

The above constraint is the fixed point problem on the market share space (see aguirregabiria2002swapping). The right-hand side is linear in the conditional distribution. The following lemma establishes the relationship between the fixed point problems in (ref) and (ref).

(ref) shows that solving the Bellman equation (ref) is equivalent to solving (ref). We further convert (ref) into an unconditional moment constraint by assuming that $s_{0}(\omega)$ is the solution to (ref) if and only if:

(ref) allows us to profile out $\{\omega_{t}\}_{t=1}^{T}$, which is useful for the implementation. A sufficient condition is that $s_{0}(\omega)$ is continuous and strictly decreasing, and its smallest and largest values are small and large enough. Under (ref), we have $s_{0}(\omega) \in C(\Omega)$. Moreover, we can expect that a higher inclusive value $\omega$ corresponds to a lower market share of not purchasing, i.e., $s_{0}(\omega)$ is decreasing in $\omega$. Therefore, (ref) is mild.

The last condition is the fixed point constraint similar to the procedure to estimate the AR(1) process. Suppose the distribution $F$ is used in (ref), and $\{\omega_{t}\}_{t=1}^{T}$ is the sequence of recovered inclusive values. Denote by $\hat{F}$ the estimator of the joint distribution for the pairs $\{(\omega_{t},\omega_{t+1})\}_{t=1}^{T-1}$. Then, our fixed point constraint is:

where $\epsilon_{T}$ is the tolerance level. To choose $\epsilon_{T}$, we estimate the joint distribution of inclusive values recovered from the AR(1) process by the kernel density estimator with Gaussian kernel and bandwidth selected by the 5-fold cross-validation. Then, we set $\epsilon_{T}$ to be the KL divergence between the kernel density estimator and the reference distribution.

For EV subsidy, by (ref), the consumer surplus (CS) is given by:

where $\omega^{\text{EV}}_{t_{1}} - \omega_{t_{1}}$ does not depend on $F$. Therefore, to bound CS, it is equivalent to bounding the change in the market share of purchase.

Putting everything together, the lower bound on the elasticity at $t_{1}$ is given by:

where $\omega'_{t_{1}}$ is replaced by $\omega^{\text{EV}}_{t_{1}}$ for the EV subsidy case. For $\delta = 0$, the reference distribution is the unique solution to the above problem. The corresponding elasticity is the reference industrywide elasticity.

Implementation

We adapt (ref) proposed in (ref). We will use numerical integration to discretize the support of the AR(1) process. Therefore, the recovered inclusive values $\{\omega_{t}\}_{t=1}^{T}$ are not differentiable with respect to the discretized market share function $s_{0}(\omega)$. To address this issue, we employ the MCMC optimization method. We alternate between solving the EOT problem to obtain the worst-case distribution, solving the Bellman equation to update $s_{0}(\omega)$, recovering the inclusive values, and checking the fixed point constraint. To derive a tractable dual formulation, we handle the fixed point constraint in a specific way. Because this constraint depends on an estimator, $\hat{F}$, that changes during optimization—potentially causing numerical instability—we first derive the dual formulation without it.\footnote{In principle, we can choose other constraints like the integral probability metrics and adapt the minimax theorem in (ref). However, it can lead to more optimization parameters.} Then, the fixed point constraint determines the acceptance/rejection of the candidate parameters in the Metropolis-Hastings step. Consider the following optimization problem:

Applying (ref), its dual is:

where $\mathcal{C}(s_{0},g,\lambda_{KL})$ is the EOT problem: $\mathcal{C}(s_{0},g,\lambda_{KL}) := \sup_{F \in \Pi(\nu_{0},\nu_{0})} \mathbb{E}_{F} \left[c(\omega,\omega';s_{0},g)\right] + \lambda_{KL} D_{KL}(F\|F_{0})$ whose cost function is $c(\omega,\omega';s_{0},g) = g(\omega)\left(\log(\frac{1 - s_{0}(\omega)}{s_{0}(\omega)}) - \omega + \beta \omega' - \beta \log(1-s_{0}(\omega'))\right)$ and the worst-case conditional distribution $F^{*}$ is given by:

where $\phi_{1}^{*}(\omega)$ and $\phi_{2}^{*}(\omega')$ are the optimal EOT potentials. During the optimization process, $dF^{*}(\omega'|\omega)$ is used to update $s_{0}(\omega)$ by solving (ref) using fixed point iteration.

As shown in (ref), the expectation in the dual is taken with respect to the reference distribution. Therefore, we discretize the estimated AR(1) process, which results in three approximation errors: (i) the Bellman equation is solved on the discretized support, (ii) $\{\omega_{t}\}_{t=1}^{T}$ are recovered approximately, and (iii) the elasticity is computed approximately. That is, we approximate $s_{0}(\omega)$ by the market share of the nearest grid point to $\omega$. Therefore, there is a trade-off between approximation error and computational cost. A finer discretization reduces the approximation error, while increasing the number of optimization parameters.

However, our dual formulation significantly improves computational efficiency. Suppose we discretize the AR(1) process into $N$ grid points. If we directly solve (ref), the number of optimization parameters is $O(N^{2})$ due to the transition matrix. In contrast, the number of optimization parameters is $O(N)$ in the dual formulation.

Algorithm (ref) summarizes the simulated annealing MCMC optimization algorithm (kirkpatrick1983optimization). It starts with the reference market share $s^{(0)}$, and alternate between proposing new parameters $(g', \lambda'_{KL})$, solving the EOT problem, solving the Bellman equation using the worst-case distribution, and accepting or rejecting the proposed parameters using the Metropolis-Hastings step based on the change in the elasticity penalized by the violation of the market share and fixed-point constraints. At each improvement step, it pools previous results across all radii for initialization. We choose 51 grid points, 5,000 MCMC steps, 5 optimization steps, 14 radii (the last is $10^{10}$), and 100 as the simulated annealing multiplier.

\setcounter{algocf}{2}

Results

We estimate two alternative transition densities for each market. The first assumes that the inclusive values are i.i.d. normally distributed. The second estimates a nonlinear AR(1) process using a cubic spline\footnote{The nonlinear AR(1) process is specified as: $g(\omega) = \sum_{k=1}^{N+3} \rho_k \Phi\left(\frac{\omega - a}{h} - (k-2)\right)$ where $a$ is the minimum of the discretized support of $\omega$, $h$ is the distance between two adjacent grid points, $\rho_{k}$ are parameters to be estimated, and $ \Phi(t) =

$. For this model, we set $N=4$ to avoid overfitting. Then, we discretize the reference AR(1) process into 4 grid points, and compute the KL divergence between the nonlinear AR(1) process and the reference AR(1) process.}. In the following figures, we plot the KL divergence between the alternative models and the reference model. The independent model is closer to the reference model with KL divergence between 0.04 to 0.67, while the nonlinear AR(1) process is farther away, with KL divergence between 3.94 to 10.09.

Figures (ref)-(ref) plot the bounds on the industrywide elasticities\footnote{schiraldi2011automobile finds a average long-run price elasticities ranging from -3.54 to -4.34 across different car segments for the Italian market. d2019automobile reports average elasticities between -3.94 and -6.40 across consumer groups for the French market. reynaert2021benefits finds a mean own-price elasticity of -5.45 for the European market. grieco2024evolution estimates an average elasticity of -5.36 for the U.S. market in 2015. Remmy2025 reports a mean price elasticity of -4.043 for the German market.} for the UK, Germany, and France in December 2023. The French market is the least elastic (reference elasticity: -4.048), while the Germany market the most elastic (reference elasticity: -6.073). The UK market's reference elasticity is -5.336. Based on our three sensitivity measures, we define the local (global) sensitivity as the ratio of the local (global) interval length to the reference value. For local sensitivity, we set $\delta=0.001$, while for global sensitivity, we set $\delta=10^{10}$. The robustness metric is defined as the smallest deviation from the reference distribution such that the elasticity can deviate by, for example, 2.5% from the reference elasticity.

The French market is the least sensitive in terms of local and global sensitivity, with 1.16% local deviation and 6.20% global deviation from the reference elasticity. The UK market is less sensitive locally (1.66%) than the German market (3.52%). They are both more sensitive globally (15.16% for the UK vs. 15.24% for Germany) than the French market. The bounds of UK market flatten around 0.178, while the French and German market flatten around 0.056. For the robustness metric, we consider 2.5% deviation from the reference elasticity. The UK market's robustness metric is around 0.018 for the upper bound and 0.008 for the lower bound. The French market's robustness metric is around 0.025 for the upper bound and 0.018 for the lower bound. The German market's robustness metric is around 0.002 for the lower bound and 0.003 for the upper bound. Therefore, in terms of robustness metric, the French market is also the most robust, while the German market is the least robust.

Figures (ref)-(ref) also plot the bounds on the consumer surplus from an additional \$3,000 EV subsidy. The subsidy is implemented between July and December 2023, when the reference consumer surplus is maximized. They are November for the UK (reference CS: \$2,880 million), October for France (reference CS: \$1,432 million), and September for Germany (reference CS: \$856 million). Overall, the EV subsidy is beneficial as the lower bounds are \$2,584 million for the UK, \$1,243 million for France, and \$309 million for Germany. The corresponding costs for subsidy are around \$12 million for the UK, \$23 million for France, and \$41 million for Germany. The costs are insensitive to the misspecification, as they only depend on the absolute change in the market share of purchases, and the conditional market share of EVs, instead of percentage change used for consumer surplus.

In terms of local sensitivity, the UK market is the least sensitive (2.52 % local deviation), the French market exhibits similar local sensitivity (4.72% local deviation), while the German market is the most sensitive (24.79% local deviation). In terms of global sensitivity, the UK market and French markets share similar sensitivity (25.17 % global deviation for the UK, and 24.73 % for France), while the German market is the most sensitive (102.75 % global deviation). For the robustness metric, if we consider 10% deviation from the reference CS, the UK and French markets' robustness metric are more than 0.018, while the German market's robustness metric is around 0.001. Therefore, the German market is also the least robust in terms of robustness metric.

Figure (ref) plots the time series of bounds on the industrywide elasticities. We set $\delta=0.001$ for the local deviation and $\delta = 1$ for the global deviation, as the bounds flatten at a maximum of 0.178 in December 2023. Large points in the figure indicate that the KL divergence constraint is binding—specifically, when the KL divergence between the worst-case distribution and the reference distribution exceeds $0.95 \cdot \delta$. When two consecutive points align horizontally, this indicates that increasing the radius does not affect the bounds, as exemplified by the UK market in November 2023. In terms of both local and global sensitivity, the UK and French markets are less sensitive than the German market. All three markets exhibit some sensitive periods. For the UK market, the upper bound's global deviation in February 2022 is around 30% away from the reference elasticity. For the French market, the upper bounds' global deviations in April 2021, and July 2022 are around 50% away from the reference elasticity. For the German market, the lower bound's global deviation in March 2022 and April 2023 is around 50% away from the reference elasticity. In terms of local sensitivity, the lower bounds' of German market in April 2021 is around 50% away from the reference elasticity.

Empirical Application: Finite Horizon DDC

This section applies our framework to a finite-horizon dynamic labor supply model for taxi drivers in New York City (NYC). In the model, the market-level supply shock is the latent variable. Our reference dynamic process is an AR(1) process. We consider the sensitivity analysis of the labor supply elasticity with respect to such distributional assumption.

The Data

We use data from New York City’s Taxi and Limousine Commission's (TLC) Taxi Passenger Enhancement Project (TPEP). The TPEP data contain a complete record of all trips operated by licensed drivers. The day shift starts at 5 AM and ends at 5 PM, and the night shift starts at 5 PM and ends at 5 AM. We choose a sample of 10,500 drivers that were active in 2013 as in kalouptsidi2021linear. We restrict the sample to day shift drivers who were only working during the day shifts. We aggregate the transaction-level data to the driver-hour level. In addition, we create 10 uniformly divided bins for weekdays (Monday-Thursday) and 4 bins\footnote{The number of bins is chosen so that a Gaussian distribution approximates the stationary distribution of the market-level supply shock recovered from the last period (see (ref)).} for weekends (Friday-Sunday) between the lowest and highest hourly earnings and calculate the average hourly earnings in each bin. Then, we remove important days (i.e., Memorial Day, the Fourth of July, and New Year's Eve). Finally, we restrict the sample to shifts that started between 5 AM and 8 AM, which accounts for 86.84% shifts for weekdays, and 71.84% for weekends. The final sample contains 3,562 drivers and 206 days for weekdays, and 3,322 drivers and 156 days for weekends.

(ref) presents the hourly summary statistics. The average hourly earnings range from \$24.22 at 4:00 PM to \$37.00 at 8:00 AM. The share of drivers who continue working is high in the early morning, with 100% of drivers working at 6:00 AM and 7:00 AM. This share starts to decline after 2 PM, and drops substantially to 52.24% and 53.59% by 4:00 PM. Therefore, we assume that drivers can only choose to stop working between 8 AM and 4 PM.

The Model

At the beginning of hour $t$ of day $m$, a taxi driver $i$ decides whether to continue working ($a=1$) or not ($a=0$). The decision to stop working is a terminating action, meaning the driver exits the market upon stopping. The period utility of working is given by:

where $k_{imt}$ is the number of hours worked, $w_{mt}$ is the average hourly earnings, $\varepsilon_{imt} := (\varepsilon_{i1mt},\varepsilon_{i0mt})$ is i.i.d. type I extreme value utility shocks, and $\xi_{mt}$ is an exogenously evolved stationary unobserved market-level supply shock. It captures the market-level time-variant unobserved heterogeneity such as weather, congestion, or city events.

We assume markets are i.i.d, and suppress the subscript $(i,m)$ for brevity. Let $u(k_{t},w_{t};\theta)$ be the deterministic part of period of utility of working up to $(\xi_{t},\varepsilon_{t})$. Let $\beta := 0.999999$ be the discount factor. The smoothed Bellman equation at time $t$ is given by:

where the conditional value functions of working and not working are given by:

First-Stage Estimation

With only one terminating action, the utility parameters cannot be identified using (ref). Therefore, we estimate the utility parameters using the Euler Equations in Conditional Choice Probabilities (ECCP) estimator introduced in kalouptsidi2021linear. By the Hotz-Miller Inversion Lemma (hotz1993conditional), we have:

where $p_{t}(k_{t},w_{t})$ is the CCP of not working. Combining (ref) and (ref) gives:

Without a distributional assumption for $\xi_{t}$, we cannot calculate the conditional expectation in (ref). However, the cross-sectional data allows us to estimate CCPs, denoted as $\hat{p}_{t}(k,w)$. We estimate $\hat{p}_{t}(k,w)$ by a flexible logit for each $t$. Let the expectational error be:

where $k_{t+1} = k_{t} + 1$. Then, we can rewrite (ref) as:

Therefore $\theta$ can be identified using an instrument for $\xi_{t} + \hat{e}(k_{t},w_{t},k_{t+1},w_{t+1},\xi_{t})$. Denote by $\mathcal{K}_{t}$ the set of possible hours worked at $t$.\footnote{Note that $|\mathcal{K}_{t}| = 4$ for $t \geq 9$ AM and $|\mathcal{K}_{t}| = 3$ for $t = 8$ AM.} The ECCP estimator stacks all $k \in \mathcal{K}_{t}$. A unit of observation is defined by day-hour.

Following kalouptsidi2021linear, we use the previous day's average hourly earnings for the same hour as the IV. (ref) shows the estimation results. The implied marginal value of time defined by $-\frac{\theta_{1} + 2\theta_{2}k}{\theta_{3}}$ ranges from \$0 at around $k=7$ hours to \$5.35 at $k=11$ hours for weekdays, and from \$0 at around $k=8$ hours to \$9.95 at $k=11$ hours for weekends.

The Reference Distribution and Scalar Parameters of Interest

We define the reference distribution and introduce the scalar parameters of interest. Let $N_{tk}$ be the number of drivers who has worked $k$ hours at hour $t$. The market-level $\xi_{t}$ equates the model-implied weighted average CCP of not working with the observed weighted average:

As we assume $\xi_{t}$ is stationary, its marginal distribution at $T$ is its stationary distribution. At hour $T$, drivers solve a static problem, and the CCP of not working is:

Therefore, the stationary distribution is identified by recovering $\xi_{T}$ to satisfy (ref) at $T$. (ref) plots the kernel density estimator of $\xi_{T}$. We fit a Gaussian distribution to the last period $\xi_{T}$. Denote its mean and standard deviation as $\mu_{\xi}$ and $\sigma_{\xi}$, respectively.

The reference model for $\xi_{t}$ is an AR(1) process:

where $\eta_{t}$ follows i.i.d normal distribution with mean 0 and variance $\sigma^{2}$. As we have identified its stationary distribution, we only need to solve the fixed point problem for $\rho$: we begin with an initial guess of $\rho$ and circulate between: (i) setting $\mu = \mu_{\xi}\cdot(1-\rho)$ and $\sigma = \sigma_{\xi} \cdot \sqrt{1-\rho^{2}}$, (ii) solving the Bellman equation using backward induction, (iii) recovering $\xi_{t}$ using (ref), and (iv) updating $\rho$ by refitting the AR(1) process to the recovered $\xi_{t}$ until convergence. The reference distribution $F_{0}$ for $(\xi,\xi')$ is the product of the transition kernel of the estimated AR(1) process and its stationary distribution $\nu_{0}$. The perturbation set is defined as:

We consider two scalar parameters of interest: the elasticity of stopping working, and the Frisch elasticity of labor supply. Both elasticities are at the individual level, meaning the demand side remains unchanged. Therefore, we keep the transitions of $w_{t}$ and $\xi_{t}$ fixed.

For the elasticity of stopping working, we increase the average hourly earnings from the current bin, $w_{mt}$, to the next bin, $w'_{mt}$. The weighted average of the elasticity at hour $t$ is:

As shown in (ref), the share of drivers who continue working begins to decline around 11 AM.Moreover, our model does not endogenize the initial entry decision (i.e., the choice of when to start a shift). To compute the Frisch elasticity, we assume 11 AM is the first hour drivers can choose to stop working, and consider a 1% increase in average hourly earnings beginning at 11 AM. The total expected hours worked at day $m$ is:

where $N_{mk}$ is the number of drivers whose hours worked is $k$ at 11 AM of day $m$, and $\bm{\xi}_{m},\bm{w}_{m}$ are the vectors of $\xi_{mt},w_{mt}$ for $t=11,\ldots,16$. Then, the Frisch elasticity is:

where $\bm{w}'_{m} = 1.01 \cdot \bm{w}_{m}$, and $p'_{}$ is derived from the Bellman equation with $\bm{w}'_{m}$.

Sensitivity Analysis

We convert (ref) into an unconditional moment constraint by assuming that $p_{t}(k,w,\xi)$ is the solution to (ref) for given $p_{t+1}(k+1,w,\xi)$ if and only if:

Let the term inside the expectation in (ref) be $\psi_{t}(k,w_{t},w_{t+1},\xi_{t},\xi_{t+1};u,p_{t},p_{t+1},g_{tk})$.

To profile out $\xi_{mt}$ from recovering the probability of stop working, we assume:

The final constraint is a fixed point constraint similar to the procedure used to estimate the AR(1) process. Suppose $F$ is used in (ref). Denote by $\hat{F}$ an estimator of the distribution of the pair $(\xi,\xi')$ using the recovered $\{\xi_{mt}\}_{m=1,t=1}^{M,T-1}$ from (ref). Then, our fixed point constraint is:

where $\epsilon_{M}$ is the tolerance level. As the sample size is large, we use Scott's Rule to initialize the bandwidth and then use 5-fold cross-validation to select bandwidth candidates around Scott's estimate that maximizes the log-likelihood. To choose $\epsilon_{M}$, we estimate the joint distribution of supply shocks recovered from the AR(1) process by the kernel density estimator with Gaussian kernel. Then, we set $\epsilon_{M}$ to be the KL divergence between the kernel density estimator and the reference distribution.

Then, the lower bound on the elasticity of stopping working at $t$ is:

The last period $T$ is a static problem, therefore $p_{T}$ is not an optimization parameter. For $\delta = 0$, the reference distribution is the unique solution to the above problem. The corresponding elasticity is the reference elasticity.

For the Frisch elasticity, we increase the earnings coefficient by 1% from 11 AM. This allows us to leave the discretized state space and its transition unchanged at the cost of solving additional Bellman equations. The lower bound on the Frisch elasticity is:

where $u'$ is the utility function with the earning coefficient increased by 1% from 11 AM.

Implementation

We adjust the procedure in (ref) to account for the finite horizon model. The main difference is that we solve the Bellman equation (ref) by backward induction. The number of Bellman equations (ref) is 31.\footnote{Note that $\mathcal{K}_{8} = \{1,2,3\}, \mathcal{K}_{t} = \{t-8,t-7,t-6,t-5\}$ for $t=9,\cdots,15$.} The number of optimization parameters is $31 \cdot N_{w} N_{\xi} + 1$ for the elasticity of stopping working, and $62 \cdot N_{w} N_{\xi} + 1$ for the Frisch elasticity, where $N_{w}$ is the number of bins for average hourly earnings, $N_{\xi}$ is the number of grid points for $\xi$, and 1 is for the KL divergence constraint. We choose $N_{\xi} = 99$, 5,000 MCMC steps for Frisch elasticities, 2,500 MCMC steps for elasticities of stopping working, 5 optimization steps, 14 radii (the last is $10^{10}$), and 100 as the simulated annealing multiplier. The covariance matrix for the random walk in (ref) is restricted to be diagonal.

Results

The alternative model is an independent model where the supply shocks are i.i.d. and follow the stationary distribution identified at $T$. The independent model is closer to the reference model for weekends than for weekdays. The KL divergence between the independent model and the reference model is 0.115 for weekends and 3.673 for weekdays.

(ref) plots the elasticity of stopping working from 9 AM to 3 PM for weekdays and weekends. We set $\delta=0.001$ for local deviation and $\delta = 1$ for global deviation. Large points indicate that the KL divergence constraint is binding. When two consecutive points align horizontally, this indicates that increasing the radius does not affect the bounds. The elasticity of stopping working is negative and decreases over time. In particular, the labor supply is inelastic before 11 AM on weekdays and 12 PM on weekends, and elastic after that. At 3 PM, the elasticity of stopping working is around -2 for weekdays and -2.5 for weekends. Overall, both weekday and weekend elasticities are not sensitive to the distributional assumption. For weekdays, elasticity in the morning is more sensitive than in the afternoon in terms of both local and global sensitivity, while it is the opposite for weekends.

(ref) plots the Frisch elasticity bounds for weekdays and weekends. The reference Frisch elasticity is 0.472 for weekends and 0.698 for weekdays. Our reference estimates are consistent with labor supply literature. For example, buchholz2023rethinking reports a Frisch elasticity ranging from 0.47-0.54 for NYC taxi drivers. For both weekdays and weekends, the bounds flatten around 0.032-0.056. In terms of local sensitivity, the results appear sensitive to the distributional assumption, with a deviation of 28.37% for weekdays and 21.76% for weekends. In terms of global sensitivity, the deviation is larger, with 76.83% for weekdays and 42.84% for weekends. For the robustness metric approach, we consider a 15% deviation from the reference Frisch elasticity. The weekday's robustness metric is around 0.001 for the upper bound, while the lower bound never reaches the 15% deviation. The weekend's robustness metric is around 0.008 for the upper bound, and 0.01 for the lower bound. Therefore, the weekday's lower bound is more robust than the weekend's lower bound, while the weekday's upper bound is less robust than the weekend's upper bound.

Conclusion

We propose a computationally tractable framework to quantify the sensitivity of outcomes of interest to misspecified latent-state dynamics in structural models. We derive bounds on a scalar parameter of interest by perturbing a reference dynamic process, while imposing a stationarity condition for time-homogeneous models or a Markovian condition for time-inhomogeneous models. We apply the approach to an infinite-horizon dynamic demand model for new cars in the UK, Germany, and France, and a finite-horizon dynamic labor supply model for taxi drivers in New York City.