Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
108,572 characters · 17 sections · 43 citation commands
Training Neural Networks Embedded in Dynamic Discrete Choice Models
\affil{Kellogg School of Management, Northwestern University}
Dynamic programs are notoriously difficult to estimate. To reduce the computational burden, classic estimators leverage two properties commonly found in dynamic discrete choice models: likelihood concavity and utility function linearity. When the empirical likelihood is concave, the optimizer finds the (pseudo) MLE on the first try---it doesn't get led astray by local optima or saddle points. And when the utility function is linear, the estimation problem decouples into two much easier problems: one that solves a set of Bellman equations, and the other that maximizes the empirical likelihood, given the corresponding fixed point solutions. Further, utility function linearity almost always implies likelihood function concavity, so we generally enjoy both benefits when the utility function is affine.
Unfortunately, affine functions have an Achilles heel: dimensionality. Splines and polynomial basis functions yield excellent approximations---when there's a single state variable. However, the number of “factors” or linear basis functions required to achieve a good approximation grows exponentially with the dimension of the state space.
Neural networks mitigate the problem of dimensionality, as they can recognize patterns in high-dimensional spaces with a parsimonious set of parameters. Breathtaking advances in neural networks mean that we must no longer model organic human preferences with flat hyperplanes. Indeed, a growing literature demonstrates the efficacy of neural networks to model utility functions in discrete choice models.\footnote{This structural approach is distinct from earlier work that applies machine learning primarily as a predictive classifier in discrete choice settings (VANCRANENBURGH2022100340), or from approaches that use neural networks to map moments to parameters Wei25.} HAN2022166 and SIFRINGER2020236 replace fixed coefficients with neural network layers, allowing preferences to vary flexibly with observed characteristics and capturing heterogeneity that linear specifications cannot accommodate. hsieh25 further formalize this approach by proposing neural utility representations that preserve economic structure while accommodating behavioral nonlinearities. aouad25 and ARKOUDI2023 extend these ideas to general random utility maximization frameworks and high-dimensional choice environments. However, all these recent works correspond to static models.
We bring neural networks to dynamic models by introducing them to the celebrated dynamic discrete choice problem.\footnote{A notable recent utilization of neural networks in dynamic settings appears in kallus for the estimation of choice probabilities and continuation values. In contrast, we embed the neural network directly into the economic model as a structural primitive to recover the underlying utility function.} This model was an instant classic, winning its creator, John Rust, the 1992 Ragnar Frisch Medal, and inspiring four survey articles: rust94, Aguirregabiria2009, Arcidiacono2011a, and Kasahara2018. We update this model by casting its utility function non-parametrically: rather than a pre-specified functional form, we allow utilities to arise from a neural network of arbitrary complexity (we make no other changes besides this one substitution: e.g., the same Bellman equation governs the system dynamics).\footnote{ This non-parametric reframing echoes \citesp{rust1987optimal}[1001] original paper, which cautioned against an over-reliance on “closed-form” specifications:
}
Unfortunately, embedding a neural network into Rust's framework sacrifices the two key properties that facilitate estimation: concavity and linearity. Losing concavity means running the optimizer from many starting points, as few gradient ascent paths will end at the global optimum. And losing linearity means nesting the dynamic programming problem within the estimation problem, so that the difficulty of the first multiplies, rather than adds to, the difficulty of the second (contrary to popular belief, even CCP estimators suffer from nested fixed points, in which additional likelihood gradient ascents necessitate additional Bellman solutions: see Section (ref)). Further, losing these two speed-ups is especially burdensome when training neural networks, given their high intrinsic training costs. Indeed, estimating a neural network utility function with traditional methods would require solving at least one Bellman equation for every neural network backpropagation.
We develop a new estimator that can efficiently train a neural network embedded in dynamic discrete choice problem. Rather than maximize a (pseudo) likelihood problem, as other estimators do, we minimize the norm of a random compression of the first-order conditions. Whereas the old maximum likelihood problem decouples from the underlying dynamic programming problem only when the utility function is linear, our new norm minimization problem decouples always. Put differently, for the maximum likelihood problem, we must solve a Bellman equation every time we evaluate the objective (when the utility function is non-linear), but for the norm minimization problem, we must solve upfront slightly more Bellman equations than there are utility parameters (e.g., neural network weights), after which we can evaluate the objective with ease. This change shortens the time between neural network backpropagations from hours to seconds (see Section (ref)). Our estimator thus frees us from a reliance on linearity.
Our estimator also frees us from a reliance on concavity (or convexity, in the norm minimization case) by enabling rapid and inexpensive exploration of the optimization landscape---precisely what is required when optimizing a non-concave objective, where exhaustive search becomes unavoidable. The curse of dimensionality ensures that almost all general purpose dynamic discrete choice estimators spend almost all of their time solving Bellman-style fixed points.\footnote{ We write `general purpose' estimators to exclude the special-purpose estimators designed for the rare cases in which the econometrician can circumvent Bellman's equation. See Section (ref). } Our estimator is no exception: the norm minimization time is negligible relative to the fixed point solution time. And since we can recycle our fixed point solutions across optimization runs, we can thus blanket the objective function with many inexpensive optimization runs. For example, initializing our optimization routine with 1000 different starting values scales the estimation time by only a factor of 7.7\ (see Section (ref)). In contrast, initializing a traditional estimator with 1000 different starting values would take 1000 times as long.
Further, our Bellman solutions are independent of the neural network configuration, which means that estimating many versions of the model is essentially no more costly than estimating a single version. Accordingly, our estimators enable the econometrician to iteratively craft the neural network to the data---assigning the most flexibility where the sample is the richest---rather than blindly specifying its configuration upfront.
We begin with some basic discrete choice theory. Consider an agent that chooses one action from a set of actions $\mathbb{A}\equiv\{1,\dots,A\}$. The utility of choosing action $a \in \mathbb{A}$ is $e_a + c_a$, where $c \in \mathbb{R}^{A}$ is a vector of deterministic values and $e \in \mathbb{R}^{A}$ is a vector of stochastic shocks that resolves from continuous density function $g$, which has full support over $\mathbb{R}^A$. McFadden's social surplus function denotes the expected utility of the top choice, unconditional on the random component:
The Williams-Daly-Zachary theorem establishes that the gradient of the social surplus function yields the corresponding choice probabilities:
Choice probability vector $p(c)$ resides in $(A-1)$-simplex $\mathbb{P}$.
Next, define $I^{+}$ as the matrix that maps $(v_{1}, \dots, v_{A-1})'$ to $(v_{1}, \dots, v_{A-1}, 0)'$, and define $\Delta$ as the matrix that maps $(v_{1}, \dots, v_{A})'$ to $(v_{1} - v_{A}, \dots, v_{A-1} - v_{A})'$. Composing these two matrices yields a normalizing operator: matrix $I^{+}\Delta$ maps $(v_{1}, \dots, v_{A})'$ to $(v_{1} - v_{A}, \dots, v_{A} - v_{A})'$. Note that normalized value vector $I^{+}\Delta c$ yields the same choice probabilities as initial value vector $c$:
Hence, we have
Next, for $v\in \mathbb{R}^{A-1}$ define
This function maps differenced value vector $\Delta c \in \mathbb{R}^{A-1}$ to choice probability vector $p(c) \in \mathbb{P}$:
In fact, this mapping is invertible Hotz1993, so there exists $\rho^{-1}$ that satisfies
The following result establishes that the random discrete problem admits an equivalent deterministic continuous problem.
For a concrete example, suppose the elements of $e$ are independent standard Gumbels. In this case, $\nu$ and $\epsilon$ respectively coincide with the log-sum-exp function and the information entropy function (both increased by Euler--Mascheroni constant $\gamma$), and $\rho$ and $\rho^{-1}$ respectively coincide with the softmax function and the log-odds function:
It is straightforward to verify that Proposition (ref) holds in this case.
We are now ready to define the dynamic discrete choice problem. We observe an agent confronted with a discrete-time, infinite-horizon Markov decision process. In each period, the agent chooses action variable $a \in \mathbb{A}\equiv\{1,\dots,A\}$ in response to state variables $e$ and $x$. The agent observes $e \in \mathbb{R}^{A}$ and $x \in \mathbb{X} \equiv \{1,\dots,X\}$, but we do not observe the former, which will serve as our statistical error term. Variable $e$ is a multivariate random vector with continuous density function $g$, which has full support over $\mathbb{R}^A$. Taking action $a$ in state $(e, x)$ sets the probability density of the next-period state, $(e', x')$, to $f_a(x'|x)g(e')$, and yields contemporaneous utility $e_a + u_a^{\theta}(x)$. Our objective is to estimate utility parameter $\theta \in \Theta \subset \mathbb{R}^{t}$, which will ultimately represent the weights of a neural network. In each period, the agent chooses the action that yields the largest expected discounted utility, under discount factor $\beta \in [0,1)$.
We could frame $\beta$ as a parameter to estimate, but doing so slightly complicates our procedure (see Appendix Section (ref)), so we will follow the usual practice, and take it as given. We similarly omit unobserved heterogeneity from our baseline model; however, accommodating it is straightforward as our estimators are compatible with standard two-stage procedures (see Appendix Section (ref)). Also, transition kernel $f$ will almost certainly be parameterized, but these parameters are trivial to pre-estimate, so we take the transition probabilities as given, without loss of generality.
The following Bellman equation specifies the agent's value function, which quantifies their expected discounted utility:
Our dynamic program has state space $\mathbb{R}^A \times \mathbb{X}$ and action space $\mathbb{A}$. But infinite state spaces are tricky, since they require interpolation, and discrete action spaces are a nuisance, since they do not yield first-order conditions, as rust explains:
Fortunately, Proposition (ref) enables us to translate our finite-action-space, infinite-state-space dynamic program to an analogous infinite-action-space, finite-state-space dynamic program. This transformation will yield first-order conditions, and will allow us to disregard the unobserved state variable. We drop this error term by integrating over it, defining a new value function that characterizes the agent's expected discounted utility, conditional on $x$, but not on $e$:
where $u^{\theta}(x)$ and $f(x'|x)$ are length-$A$ vectors with $a$th elements $u_{a}^{\theta}(x)$ and $f_{a}(x'|x)$.
Whereas Bellman equation (ref) characterizes a dynamic program with state $(e, x) \in \mathbb{R}^{A}\times \mathbb{X}$ and action $a \in \mathbb{A}$, Bellman equation (ref) characterizes a dynamic program with state $x \in \mathbb{X}$ and action $p \in \mathbb{P}$. Thus, rather than elements of $\mathbb{A}$, we can equivalently suppose that the agent chooses conditional choice probability (CCP) measures over $\mathbb{A}$. Under this framing, the policy functions that govern the agent's behavior are mappings from $\mathbb{X}$ to $\mathbb{P}$. We will use $\mathbb{X}^{\mathbb{P}}$ to denote the set of such mappings. Further, we will use $p$ to denote a generic element of $\mathbb{P}$ (i.e., a single CCP vector), and use $P$ to denote a generic element of $\mathbb{X}^{\mathbb{P}}$ (i.e., a collection of CCP vectors indexed by $\mathbb{X}$).
Adhering to the convention of using capital letters for objects that span over $\mathbb{X}$, define $E(P)$, $P_{a}$, $U_{a}^{\theta}$, and $V^{\theta}$ as the length-$X$ vectors with $x$th elements $\epsilon(P(x))$, $P_{a}(x)$, $u_{a}^{\theta}(x)$, and $v^{\theta}(x)$, and define $F_{a}$ as the $X\times X$ matrix with $xx'$th element $f_{a}(x'|x)$. Next, define $U_{P}^{\theta} \equiv E(P) + \sum_{a \in \mathbb{A}} \text{diag}(P_{a}) U_{a}^{\theta}$ and $F_{P} \equiv \sum_{a \in \mathbb{A}} \text{diag}(P_{a}) F_a$ as the utility function and transition kernel associated with policy $P$. And, finally, define $C_{V}^{\theta}$ as the $X \times A$ matrix with $a$th column $U_{a}^{\theta} + \beta F_{a}V$. Each element of $C_{V^{\theta}}^{\theta}$ is a choice specific value (CSV)---i.e., an expected discounted utility, net of the current error term.
Proposition (ref) establishes that these CSVs underlie the first-order conditions associated with Bellman equation (ref):
The expression above relies on the convention that $\rho^{-1}$ and $\rho$ are evaluated row-wise when applied to a matrix. For example, the $x$th row of system (ref) specifies that
and the $a$th column specifies that
System (ref) implies the following optimal policy:
In the expression above, $P_{V}^{\theta} \in \mathbb{X}^{\mathbb{P}}$ specifies the optimal policy to follow today when $V \in \mathbb{R}^{X}$ specifies the value tomorrow. And whereas $P_{V}^{\theta}$ maps value functions to policy functions, $V_{P}^{\theta}$ maps policy functions to value functions:
Specifically, $V_{P}^{\theta}$ characterizes the value of following policy $P$ forever. Composing ${P}_{V}^{\theta}$ and $V_{P}^{\theta}$ yields the policy iteration operator:
The expression above yields the optimal policy for the agent to follow today, given that they will follow policy $P$ forever thereafter. It underlies the policy iteration algorithm, which is second only to the celebrated value iteration algorithm, which iteratively applies the following operator:
If $V$ characterizes tomorrow's value, then $T_{V}^{\theta}$ characterizes today's value. Bellman equation (ref) establishes that $V^{\theta}$ is the unique solution to the following fixed point equation:
rust devised the first dynamic discrete choice estimator. His nested fixed point (NFXP) algorithm maximizes the empirical likelihood of the observed actions:
and $n_{a}$ is a length-$X$ vector whose $x$th element specifies the number of action-$a$, state-$x$ observations in our sample. Since this log likelihood depends on the optimal policy, which depends on the solution of (ref), Rust's algorithm solves a Bellman equation every time it evaluates the empirical likelihood function. In other words, it nests a dynamic programming fixed point problem within a maximum likelihood problem.
However, (ref) is not the only nested fixed point equation that bedevils the NFXP algorithm, as the gradient of the likelihood function depends on the following $t$ fixed points:\footnote{ When calculating this derivative, it is useful to note that $\tfrac{\partial}{\partial \theta}V^{\theta} = \tfrac{\partial}{\partial \theta}V_{P}^{\theta} \big|_{P = P^{\theta}}$, which follows from the envelope theorem. }
Specifically,
and $\tfrac{\partial}{\partial j} \rho_{a}(M)$ denotes the $X \times X$ Jacobian of $\rho_{a}(M)$ with respect to the $j$th column of matrix $M$. Further, the Hessian matrix of the likelihood function depends on the $t(t+1)/2$ distinct fixed points defined in the following:
Specifically,
Accordingly, applying a Newton step to the NFXP likelihood function has traditionally required solving the $t + t(t+1)/2$ fixed points defined in (ref) and (ref). nfxp_manual's nfxp_manual algorithm avoids such Newton steps, since
However, Rust's algorithm does solve the $t$ fixed points in (ref) for each gradient ascent step.
To avoid solving a full dynamic program for each likelihood evaluation, Hotz1993 developed the conditional choice probability (CCP) estimator, which “does not require econometricians to explicitly solve the valuation functions used to characterize optimal decision rules via backwards recursion methods." There are now many CCP estimators, but essentially all of them solve some version of the following problem: \footnote{ Although often expressed in slightly different forms, equation (ref) is a core component of all general-purpose, infinite-horizon CCP estimators. For example, it corresponds to \cites{Aguirregabiria2002a} equation (8), \cites{Aguirregabiria2007} equation (14), \cites{Pakes2007a} equation (9), \cites{Pesendorfer2008} equation (6), \cites{Bajari2007} first equation on page 1343, \cites{Hotz1994} equation (4.2), and to the infinite-horizon analog of \cites{Hotz1993} equation (3.12). Aguirregabiria2002a, Aguirregabiria2007, Pakes2007a and Pesendorfer2008 re-express fixed point (ref) in terms of $(I - \beta F_{\hat{P}})^{-1}$. But for large problems, it would generally be faster to repeatedly solve (ref) than explicitly compute $(I - \beta F_{\hat{P}})^{-1}$. Bajari2007 and Hotz1994 use forward simulation to approximate the equivalent sequential representation, $V_P^\theta = \sum_{t=0}^\infty \beta^t F_P^t U_P^\theta$. Finally, Hotz1993 works with the finite-horizon analog of (ref). To see this with their notation, note that
Denote the expected continuation value given $H_t$ with $V(H_t) \equiv \sum_{k = 1}^J p_{tk} v_k(H_{t})$, and the policy induced transition with $F_{pt}(H_{t+1}|H_t) \equiv \sum_{j=1}^{J} p_{tj} F_j(H_{t+1}|H_t)$. Then
Index feasible histories at time $t$ as as $\{H^1_t, \dots, H^{n_t}_t\}$. Let $V_t \in \mathbb{R}^{n_t}$ have $i$th entry $V(H^i_t)$, and $U_{t+1} \in \mathbb{R}^{n_{t+1}}$ have $j$th entry $U(p_{t+1}, H^j_{t+1})$. Define $F_{pt}$ as the $n_t \times n_{t+1}$ matrix with $ij$th entry $F_{pt}(H^j_{t+1}|H^i_t)$. Then $V_t = F_{pt}U_{t+1} + F_{pt}V_{t+1}$. In a stationary environment with generic states, including geometric discounting $\beta \in [0,1)$, letting $T \to \infty$ gives the infinite–horizon analog: $V = \beta F_p U + \beta F_p V$. }
It is generally understood that CCP estimators do not suffer from nested fixed points. This belief is incorrect. For example, evaluating the CCP likelihood requires evaluating $P_{\hat{P}}^{\theta}$, which in turn requires solving fixed point equation (ref), with $P = \hat{P}$. However, CCP fixed point (ref) is simpler than NFXP fixed point (ref): the former solves the value function that corresponds to policy $\hat{P}$, which is pre-estimated, whereas the latter solves the value function that corresponds to policy $P^{\theta}$, which is computed. Put differently, evaluating the CCP likelihood function requires implementing a single policy-iteration step, whereas solving the NFXP likelihood requires iterating policy-iteration steps to convergence.
The CCP likelihood function also has similar gradient and Hessian nested fixed points. Specifically, $\tfrac{\partial}{\partial \theta} \mathcal{L}_{\text{CCP}}^{\theta}$ and $\tfrac{\partial^{2}}{\partial \theta_{i}\partial \theta_{j}} \mathcal{L}_{\text{CCP}}^{\theta}$ depend on the solutions to
as $\tfrac{\partial}{\partial \theta} \mathcal{L}_{\text{NFXP}}^{\theta}$ and $\tfrac{\partial^{2}}{\partial \theta_{i}\partial \theta_{j}} \mathcal{L}_{\text{NFXP}}^{\theta}$ depend on the solutions to (ref) and (ref). In fact, (ref) and (ref) are no easier to solve than (ref) and (ref), conditional on a fixed policy. Thus, essentially all of the computational gains from the CCP estimator stem from avoiding the repeated computation of $P^{\theta}$.
\cites{Bray_strong_convergence} strong convergence (SC) estimator is the same as the CCP estimator, except it replaces nested fixed point (ref) with the following nested fixed point:
where $\Psi$ is the matrix that maps $(v_{1}, \cdots, v_{X})$ to $(v_{1} - \bar{v}, \cdots, v_{X} - \bar{v})$, for $\bar{v} \equiv \sum_{i = 1}^{X}v_{i}/X$. Specifically,
However, it is simple to show that $P_{\Psi V}^{\theta} = P_{V}^{\theta}$, which implies that $\mathcal{L}_{\text{SC}}^{\theta} = \mathcal{L}_{\text{CCP}}^{\theta}$. Thus, SC is simply a different implementation of CCP. But this implementation will be faster when the Markov chain under policy $\hat{P}$ is ergodic, in which case successive approximations of (ref) will converge strictly faster than successive approximations of (ref), a property known as strong convergence Bray_endogenous.
We can leverage strong convergence for computing gradients and Hessians as well, since $\tfrac{\partial}{\partial \theta} \mathcal{L}_{\text{SC}}^{\theta}$ and $\tfrac{\partial^{2}}{\partial \theta_{i}\partial \theta_{j}} \mathcal{L}_{\text{SC}}^{\theta}$ depend on the solutions to
which converge faster than the solution to (ref) and (ref) when the state variables follow an ergodic Markov chain.
So far we have made only incremental progress: the SC fixed points are easier to solve than the CCP fixed points, which in turn are easier to solve than the NFXP fixed points. However, a real breakthrough arises in one special case: when the utility function is linear. If $U^{\theta}_a = W_a\theta$, for some $X \times t$ matrix $W_a$, then the CCP and SC fixed points unnest, and these estimators become substantially more efficient (aguirregabiria1999dynamics; Aguirregabiria2009; Bray_strong_convergence). For the CCP case, note that $U^{\theta}_a = W_a\theta$ implies
And for the SC case, we apply $\Psi$ where appropriate. As you see, there are still fixed points in the system above, but these fixed points do not depend on $\theta$. Hence, we could solve $V_{\hat{P}}^{E}$ and $V_{\hat{P}}^{W}$ upfront, after which $V_{\hat{P}}^{\theta}$, $\tfrac{\partial}{\partial \theta} V_{\hat{P}}^{\theta}$, and $\tfrac{\partial^{2}}{\partial \theta \partial \theta}V_{\hat{P}}^{\theta}$ are available in closed form. Hence, after solving the two fixed point equations above, we could evaluate $\mathcal{L}_{\text{CCP}}^{\theta}$, $\tfrac{\partial}{\partial \theta} \mathcal{L}_{\text{CCP}}^{\theta}$, and $\tfrac{\partial^{2}}{\partial \theta_{i}\partial \theta_{j}} \mathcal{L}_{\text{CCP}}^{\theta}$ for any number of $\theta$ values, with negligible marginal effort. Computing $V_{\hat{P}}^{E}$ and $V_{\hat{P}}^{W}$ may be costly in large problems, but once obtained, maximizing the CCP likelihood is trivial.\footnote{ Another way to unnest the CCP fixed points would be to compute $(I - \beta F_{\hat{P}})^{-1}$ upfront. But this matrix inversion is not practical in general, as it requires $O(X^3) $ operations with standard algorithms (or $ O(X^{2.373}) $ operations with esoteric algorithms), whereas solving for $V_{\hat{P}}^{\theta}$ requires only $ O(X^2) $ operations Judd1998. }
Our contribution will be to extend this fixed point unnesting technique to dynamic programs with non-linear utility functions.
The following proposition will underpin our estimation approach.
We will use this proposition to create estimators whose fixed points unnest for all utility functions in the next section. But before that, we will illustrate the value of Proposition (ref) by showing how it accelerates the computation of $\tfrac{\partial}{\partial \theta_{i}} \mathcal{L}_{\text{NFXP}}^{\theta}$, $\tfrac{\partial^{2}}{\partial \theta_{i}\partial \theta_{j}} \mathcal{L}_{\text{NFXP}}^{\theta}$, $\tfrac{\partial}{\partial \theta_{i}} \mathcal{L}_{\text{CCP}}^{\theta}$, and $\tfrac{\partial^{2}}{\partial \theta_{i}\partial \theta_{j}} \mathcal{L}_{\text{SC}}^{\theta}$.
To begin, consider a dynamic program with transition function $\tilde{F}_{P} \equiv F_{P}$ and utility function $\tilde{U}_{P}^{\theta} \equiv \sum_{a \in \mathbb{A}} \text{diag}(P_{a}) \tfrac{\partial}{\partial \theta_{i}}U_{a}^{\theta}$. Line (ref) establishes that restricting this dynamic program to policy $\tilde{P} \equiv P^{\theta}$ yields value function $\tilde{V}_{\tilde{P}} \equiv \tfrac{\partial}{\partial \theta_{i}} V^{\theta}$. Next, define $\lambda$ as the solution to (ref) under policy $\tilde{P}$ and weighting vector
In this case, line (ref) establishes that
And, with this, we can express the NFXP likelihood gradient in a way that depends on the one fixed point in (ref), rather than $t$ fixed points in (ref):
We can speed up the computation of $\tfrac{\partial^{2}}{\partial \theta_{i}\partial \theta_{j}} \mathcal{L}_{\text{NFXP}}^{\theta}$ with the same trick. For this case, consider a dynamic program with state transition matrix $\breve{F}_{P} \equiv F_{P}$ and utility function
Line (ref) establishes that restricting this dynamic program to policy $\breve{P} \equiv P^{\theta}$ yields value function $\breve{V}_{\breve{P}} \equiv \tfrac{\partial^{2}}{\partial \theta_{i} \partial \theta_{j}}V^{\theta}$. Thus, for $w$ and $\lambda$ defined above, (ref) establishes that
And, with this, we can express the NFXP likelihood Hessian in a way that depends on the one fixed point in (ref) and the $t$ fixed points in (ref), rather than the $t(t+1)/2$ fixed points in (ref):
The $\tfrac{\partial}{\partial \theta_{i}} \mathcal{L}_{\text{CCP}}^{\theta}$, $\tfrac{\partial}{\partial \theta_{i}} \mathcal{L}_{\text{SC}}^{\theta}$, $\tfrac{\partial^{2}}{\partial \theta_{i}\partial \theta_{j}} \mathcal{L}_{\text{CCP}}^{\theta}$, and $\tfrac{\partial^{2}}{\partial \theta_{i}\partial \theta_{j}} \mathcal{L}_{\text{SC}}^{\theta}$ simplifications are analogous: Proposition (ref) decreases the number of fixed points required to evaluate the gradients from $t$ to 1, and decreases the number of fixed points required to evaluate the Hessians from $t(t+1)/2$ to $t + 1$. These reductions are meaningful when the utility function is parameterized by the weights of a neural network, in which case $t$ could be in the hundreds or thousands.
Proposition (ref) reduces the number of nested fixed points in NFXP and CCP, but it does not eliminate them. For example, the $\lambda_{i}$ used to streamline the computation of $\tfrac{\partial}{\partial \theta_{i}} \mathcal{L}_{\text{NFXP}}^{\theta}$ depends on $\theta$, and therefore must be recomputed at each gradient ascent step. The difficulty is that NFXP, CCP, and SC do not naturally align with Proposition (ref), since their objective functions involve many distinct linear combinations of the value function. For example, the CCP likelihood depends on $(F_a(x) - F_A(x)) V_{\hat{P}}^{\theta}$ for every $x \in \mathbb{X}$ and every $a \in \mathbb{A} \setminus \{A\}$, where $F_a(x)$ denotes the $x$th row of $F_a$. Consequently, unnesting the CCP fixed points via Proposition (ref) would require solving (ref) for $X(A-1)$ different weighting vectors $w$.
Rather than fit Proposition (ref) to pre-existing estimators, as we did in the previous section, we will now fit an estimator to Proposition (ref). Whereas NFXP and CCP operate in the CCP space---depending on $P^{\theta}$ or $P_{\hat{P}}^{\theta}$---our estimator operates in the CSV space---depending on $C_{V_{\hat{P}}^{\theta}}^{\theta}$. Put differently, NFXP and CCP exploit the optimal policy, line (ref), while we exploit the associated first-order conditions, line (ref). This shift is crucial: the optimal policy is nonlinear in the value function, whereas the first-order conditions are linear in value function. This linearity will enable us to construct an empirical objective that depends on only a small number of linear combinations of the value function.
The first-order conditions in (ref) imply that
Our unnested fixed point (UFXP) estimator minimizes the sum of squares of $m > t$ random linear functionals of these limit conditions (although $m = t + 1$ suffices, we usually set $m$ slightly larger):
and $Z_{i}$ is an $X \times (A-1)$ matrix comprised of independent standard normal random variables, for $i \in \{1, \cdots, m\}$.
Objective $\mathcal{Q}_{Z}^{\theta}$ fits seamlessly with Proposition (ref). To see this, note that
$\lambda_{i}$ is the corresponding fixed point of (ref), and $Z_{ia}$ is the $a$th column of $Z_{i}$. The key, of course, is replacing $w_{i}'V_{\hat P}^{\theta}$ with $\lambda_{i}'U_{\hat P}^{\theta}$, as doing so exchanges $\theta$-dependent fixed point (ref) for $\theta$-independent fixed point (ref). Accordingly, after solving $\lambda_{1}, \cdots, \lambda_{m}$, we can evaluate our estimator's objective for any number of $\theta$ values, with negligible marginal cost. These $m$ dual-variable vectors also yield the objective's gradient and Hessian matrix, which also do not require any fixed point computations to evaluate (see Appendix Section (ref) for the explicit expressions of $\tfrac{\partial}{\partial \theta} \mathcal{Q}_{Z}^{\theta}$ and $\tfrac{\partial^2}{\partial \theta_j \partial \theta_k} \mathcal{Q}_{\text{Z}}^{\theta}$).
The following theorem establishes the consistency and asymptotic normality of our UFXP estimator.
\cites{cameron} equation (6.24) suggests that we can improve our estimator's asymptotic performance by replacing random weights $z$ with the following “optimal” weights:
where $\hat \Sigma(x)$ and $\hat \Gamma(x)$ are analogs to $\Sigma(x)$ and $\Gamma(x)$ with $\hat P(x)$ replacing $P^{\theta}(x)$.
Unfortunately, the optimal weights depend on $\theta$, so using them requires a two-step estimator. In the first step, we calculate $\hat{\theta}_{\text{UFXP}}$, and in the second step we calculate
Note, if the utility function is linearly parameterized ($U^{\theta}_a = W_a\theta$) then the OUFXP estimator simplifies to:
where $\Lambda$ is the $t \times X$ matrix whose $i$th row is $\lambda_i$, and $W(x)$ is the $A \times t$ matrix whose $a$th row corresponds to the $x$th row of $W_a$.\footnote{ See MiessiSanches2016 and Dearing2019 for other closed-form estimators of dynamic programs with linearly parameterized utility functions. }
The following result justifies the name of our new estimator.
When interpreting Theorem (ref), keep in mind that we have treated the model's state transition matrices as given. Thus, this theorem indicates that OUFXP is as asymptotically efficient as any other dynamic discrete choice estimator that first pre-estimates state transition probabilities (which is essentially all of them).
We will compare our estimators to other prominent dynamic discrete choice estimators with a simplified version of the structural econometric model of Section (ref). Our estimators will tackle this toy model with ease, but traditional estimators will struggle. The difficulty traditional estimators face with this simplified model suggests that the more elaborate version we present in Section (ref) would be effectively unestimatable without our approach.
Our model concerns a grocery store's infinite-horizon inventory management problem. The key object of interest is the inventory holding cost. Most stylized inventory models---such as the newsvendor, linear-quadratic, and $(s, S)$ specifications---capture inventory costs with one or two parameters. But with over a million observations in our actual empirical sample (which we will introduce in the next section), we can afford to estimate the inventory holding cost function nonparametrically.
The store's holding cost function depends on three state variables: $r \in \{1, 2, 3, 4, 5, 6\}$ indexes the expectation of the current daily demand, $o \in \{1, 2, 3\}$ indexes the amount of other-product inventory held, and $i \in \{0, \cdots, 29\}$ records the number of units of the focal product held. The holding cost combines a baseline linear cost, a step penalty, and a weighted combination of concave and convex tails (see Appendix Section (ref) for further details).
Figure (ref) depicts holding cost $h(r, o ,i)$ (and our estimates of it, which we will discuss later). To be clear, this function is ad hoc: there is no economic rationale for this particular functional form. However, we believe this function represents roughly the correct level of complexity, as many operational factors could kink and twist a store's holding cost function: e.g., the focal-product holding cost could increase with the other-product inventory level, due to crowding in the backroom; the holding cost could increase with the demand level, as more popular products are allocated more valuable real estate in the back room; and the holding cost could be discontinuous, as overflow inventory is stored elsewhere, once the designated shelf is full.
Importantly, the holding cost function is not only irregular, but is irregular in three dimensions. This makes holding costs difficult to capture with a linear basis (e.g., splines perform poorly in higher dimensions). We believe most true, organic utilities would exhibit such multi-dimensional idiosyncrasies---a belief that motivates our embrace of neural networks.
The holding cost is just one part in a four-part utility function:
The second term in the expression above is the cost of unfulfilled demand, the third term is the cost of shipping, and the final term is a statistical error term. More specifically, $\eta$ is the loss associated with one unsatisfied customer, $\kappa$ is the cost of receiving one shipment of inventory, $d$ is the daily customer demand, $q$ is the amount of inventory the store orders from its supplier, and $e_{q}$ is an independent Gumbel shock associated with order quantity $q$. The utility function above is analogous to the one presented in Section (ref), which we develop in more detail.
The action variable is order quantity $q \in \mathbb{Q} \equiv \{0, 6, 12, 18\}$, which the store manager chooses with the aim of maximizing expected discounted utility, under discount factor $\beta = 0.9997$ (which corresponds to an annual discount factor of $0.9997^{365} = 0.896$).\footnote{Fixing this parameter is not necessary, as we could estimate it with the procedure outlined in Appendix Section (ref))} Finally, our model has the following transition kernel
In addition to $h$, we formally define $f^r$, $f^o$, and $f_{q}^i$ in Appendix Section (ref). The model we have described has $6 \cdot 3 \cdot 30 = 540$ distinct states. To gauge how estimation times vary with the size of the state space, we also consider a version of this model with $6 \cdot 3 \cdot 300 = 5400$ states (see Appendix Section (ref)).
For a given version of our model, we begin by solving the agent's optimal policy, and using that to simulate a sample with 1,073,679 observations (the size of the panel we will use in Section (ref)). Next, we use this random sample to pre-estimate CCPs and transition probabilities (see Appendix Section (ref) for estimation details). And finally, we estimate $\eta$, $\kappa$, and $h$ from our simulated sample, $\hat{P}$, $\hat{f}_{q}$, and $\beta$.
Note, we treat $h$ as a model primitive, rather than the parameters that define it. We estimate this function with a multi-layer perceptron that receives six features---$r$, $o$, $i$, $ro$, $ri$, $oi$---and returns a corresponding holding cost. Our neural network has one output layer of width 1 and $\gamma_{1}$ hidden layers of width $\gamma_{2}$. We try three configurations: a wide architecture with $\gamma = (2, 6)$, a balanced architecture with $\gamma = (4, 4)$, and a deep architecture with $\gamma = (6, 3)$. The wide network comprises 91 parameters (78 weights and 13 biases), the balanced network comprises 93 parameters (76 weights and 17 biases), and the deep network comprises 85 parameters (66 weights and 19 biases). For each architecture, we consider two possible activation functions, ReLU and softplus, for a total of six neural network specifications. Finally, we apply an output transformation that enforces nonnegativity and makes the function evaluate to zero when $i = 0$. More specifically, we set $h(r, o, i) = \text{ReLU}(n(r, o, i) - n(r, o, 0))$, where $n(r, o, i)$ is the output of the neural network.
Overall, specifying the neural network entails several design choices: the architecture, the activation function, the set of inputs (e.g., interactions $ro$, $ri$, and $oi$ rather than quadratics $r^{2}$, $o^{2}$, and $i^{2}$), and the final transformation function (e.g., the imposition of economic restrictions $h(r,o,0) = 0$ and $h(r,o,i) \ge 0$). Neural networks are highly customizable.
We estimate $\eta$, $\kappa$, and $h$ with the UFXP, OUFXP, NFXP, CCP, and SC estimators presented in Section (ref), as well as the mathematical program with equilibrium constraints (MPEC) estimator of \cites{mpec}.\footnote{The MPEC estimator avoids nested fixed points by recasting Bellman's problem from a subroutine required to evaluate the likelihood function into a constraint that need only hold at the optimum:
Whereas NFXP enforces (ref) each time the likelihood is evaluated, MPEC enforces it only for the final parameters. Since it computes only one Bellman equation, this formulation seems easier to solve, but that is not necessarily the case, as comment_mpec demonstrated.} NFXP, CCP, SC, and MPEC are general-purpose estimators, capable of estimating any model that conforms to the framework presented in Section (ref). These four are therefore natural benchmarks for our estimators, which are also general purpose. We do not compare UFXP and OUFXP to specialized estimators that leverage structural properties of the dynamic program.\footnote{ Examples of such special-purpose estimators are those that leverage finite dependence, which a dynamic program exhibits if there exists a policy under which the distribution of the state variables after a set number of transitions is independent of the current action Hotz1994, Arcidiacono2011. Another example is the Euler equation approach of Aguirregabiria2013, which replaces the inversion of $I-\beta F_{\hat P}$ with local first-order conditions that offer computational advantages only when the transition kernel exhibits strong local structure (such as sparse or low-dimensional action-dependent dynamics) so that the set of states that are reachable within two steps is small. Similarly, the pseudo-value function estimator of Dearing2019 is computationally advantageous only if some policy $P$ arranges matrix $I - \beta F_{P}$ in a form that is easy to invert, such as triangular, banded, block-diagonal, or arrow. Other estimators, such as MPEC, exploit sparsity in the state transition matrix. } The model we estimate does not exhibit any exploitable special structure (while common in single-dimensional problems, such structure is rare in higher dimensions).
We maximize $\mathcal{L}_{\text{NFXP}}^{\theta}$, $\mathcal{L}_{\text{CCP}}^{\theta}$, $\mathcal{L}_{\text{SC}}^{\theta}$, $-\mathcal{Q}_{Z}^{\theta}$, and $-\mathcal{Q}_{Z^{\hat{\theta}_{\text{UFXP}}}}^{\theta}$ with four optimizers: Knitro Interior-Point Direct barrier method knitro, Adam (Adaptive Moment Estimation) adaptive gradient descent method adam, L-BFGS (Limited-memory Broyden–Fletcher–Goldfarb–Shanno) quasi-Newton method lbfgs, and Trust-region Newton method Nocedal2006Numerical. (Following mpec, we maximize $\mathcal{L}_{\text{MPEC}}^{\theta, V}$ with Knitro only, as the other three optimizers have difficulty with MPEC's non-linear constraint).
We provide all four optimizers with analytical objectives and gradients, and additionally provide Trust and Knitro with analytical Hessians. We compute the NFXP, CCP, and SC gradients and Hessians with the dual fixed point representations facilitated by Proposition (ref) (see Section (ref)). Doing so reduces the number of fixed points solved by a factor of at least $t/2 \ge 85/2$. Consequently, we implement these benchmark estimators more efficiently than usual.
Since the estimators' objective functions are not concave, multiple starts will be required to achieve global optimality (or close to it). For all estimations of the 540-state model, we use 100 random starting points. We also use 100 random starting points for the 5400-state model; however, since this model is much more difficult to estimate, we limit our focus in this case to the L-BFGS and Trust optimizers, and to the balanced-architecture, softplus-activation network.\footnote{We exclude MPEC from the 5400-state experiments because it is incompatible with these optimizers, and casual experimentation suggests it performs as poorly with the 5400-state model as it does with the 540-state model.}
Following standard neural network training practices, we draw our initial neural network weights from a Kaiming uniform distribution kaiming15, draw our neural network biases from a uniform distribution with bounds defined by the inverse square root of the number of inputs of each respective layer, and draw $\eta$ and $\kappa$ from a Uniform(0, 10) distribution. Finally, for MPEC, we draw the initial elements of $V$ independently from the standard normal distribution.
For the UFXP estimator's weighting matrices, $Z_1, \dots, Z_m$, we set $m = 100$. To give more weight to the more precisely estimated elements of $\hat{P}$, we draw the $(x,q)$th element of matrix $Z_i$ from an independent normal distribution with mean 0 and variance $\tfrac{n_{x,q}n_{x,0}}{n_{x,q} + n_{x,0}}$, where $n_{x,q}$ denotes the number of times action $q$ is observed in state $x$. For the OUFXP estimator, we construct weighting matrix $Z^{\hat{\theta}_{\text{UFXP}}}$ from the best first-stage UFXP estimate for each neural network configuration. Further, we warm-start each OUFXP optimization from its respective UFXP estimate rather than from random starting points.
We execute each estimation as an independent, single-threaded job allocated to a single CPU core on Northwestern University's Quest HPC cluster, using nodes equipped with Intel Xeon processors. We allocated 4 GB of RAM for the 540-state model and 64 GB for the 5400-state model. To ensure a fair comparison of computational time, we submitted the jobs in parallel and executed them in a randomized order across all combinations of estimators, neural network configurations, optimizers, and random seeds for initializations. Finally, the computer cluster subjected every estimation run to a seven-day time limit, from the start of the job.
We execute 12600 optimization runs for the 540-state model: 2400 runs (100 starting values times six neural network configurations times four optimizers) for each of NFXP, CCP, SC, UFXP, and OUFXP, and 600 runs (100 starting values times six neural network configurations times one optimizer) for MPEC. And we execute 1000 runs for the 5400-state model: 200 runs (100 starting values times one neural network configuration times two optimizers) for each of NFXP, CCP, SC, UFXP, and OUFXP. Only UFXP and OUFXP successfully return estimates for all optimization runs across both models. NFXP, CCP, and SC failed to converge within the seven-day time limit for 3.6%, 5.3%, and 0.4%\ of the 540-state runs, respectively, and for 62%, 60%, and 13%\ of the 5400-state runs. While all 600 MPEC runs for the 540-state model finished within the time limit, 7.2%\ of them failed to converge to a feasible solution satisfying (ref).
Figure (ref) demonstrates that our neural network approximations can accurately recover the retailer's complex holding cost function. The plots compare the true $h$ for the 540-state model with the best estimate, $\hat{h}$, from all successful runs of each estimator. We define “best" as the estimate that minimizes ${\vert\kern-0.25ex\vert \hat{h} - h \vert\kern-0.25ex\vert}_{1}$, where ${\vert\kern-0.25ex\vert \cdot \vert\kern-0.25ex\vert}_{1}$ denotes the $\ell_1$ norm. The functional $R^{2}$ that treats each state as one observation indicates that the best NFXP, CCP, SC, UFXP, and OUFXP neural network approximations explain 98.2%, 98.4%, 98.9%, 97.3%, and 97%\ of the variation in $h$ for the 540-state model, and 95.4%, 97.4%, 98.6%, 97.7%, and 97.4%\ for the 5400-state model, respectively.\footnote{ The seven-day computational limit may have lowered the $R^2$ values for the NFXP and CCP estimators of the 5400-state model, as optimization runs that explored the parameter space more thoroughly were more likely to time out. } In contrast, the best MPEC approximation for the 540-state model explains only 57.8%\ of the variation in $h$. Hence, out of 600 MPEC runs, not a single one yielded a reasonable estimate of $h$.
Figure (ref) plots the empirical CDF of the 540-state model's holding cost error across the 100 $\hat{\theta}$ initializations, for each neural network, optimizer, and estimator. These empirical CDFs demonstrate that we cannot attribute MPEC's poor performance to Knitro, as that optimizer performs well for the other five estimators. The MPEC objective---with its 540 additional parameters and constraints---simply seems harder to maximize. However, the MPEC estimator is not the only one that yields poor estimates: 49%, 48%, 50%, 66%, and 67%\ of NFXP, CCP, SC, UFXP, and OUFXP estimates are worse than the poor MPEC estimate depicted in Figure (ref).
This performance inconsistency underscores the need for multiple optimization starts. To make this point clear, we will define an estimate as “inadequate” if its holding cost error, ${\vert\kern-0.25ex\vert \hat{h} - h \vert\kern-0.25ex\vert}_{1}$, is 50% larger than the smallest error for the given (neural network, optimizer, estimator) triple. For example, the minimum (Softplus Balanced, L-BFGS, NFXP) error is 0.0609, so the (Softplus Balanced, L-BFGS, NFXP) estimates for which ${\vert\kern-0.25ex\vert \hat{h} - h \vert\kern-0.25ex\vert}_{1} \ge 1.5 \cdot 0.0609 = 0.0913$ are inadequate. Overall, 81%, 80%, 81%, 85%\, and 86%\ of NFXP, CCP, SC, UFXP, and OUFXP estimates are inadequate. Accordingly, ensuring less than a 1% probability of reporting an inadequate estimate would require running the optimizer from at least $\lceil\log(0.01)/\log(0.80)\rceil =21$ starting points.
Our unnested estimators are relatively more competitive when we initialize the optimizer with $R = 21$ starting points than when we initialize them with $R = 1$ starting point, as Table (ref) demonstrates. This table reports the ratios of the average estimation times of NFXP, CCP, SC, and MPEC divided by the average estimation times of UFXP\footnote{See Table (ref) in Appendix for an OUFXP analog.} (recall that we truncate all estimation times to seven days, a transformation that works against us, since no UFXP run lasted longer than a week). Increasing $R$ makes UFXP relatively more efficient, as it enables the estimator to amortize the time spent solving dual-variable fixed points across more optimization starts. For example, for the (Softplus Balanced, L-BFGS, 5400-state) case, computing $\lambda_1, \dots , \lambda_m$ requires 872\ seconds, whereas maximizing $-\mathcal{Q}_{Z}^{\theta}$ requires an average of 20\ seconds. Accordingly, the average estimation time when we initialize UFXP with $R$ starts is $872 + 20 R$. In contrast, the corresponding average estimation times when we initialize NFXP, CCP, and SC with $R$ starts are $433586 R$, $443627 R$, and $273470 R$. Note that whereas the average UFXP time is affine in $R$, the average NFXP, CCP, SC, and MPEC times are linear in $R$, since each optimization run solves an independent set of fixed points.
But even without this multi-start advantage, the unnested estimator is still markedly more efficient: across all $R = 1$ cases, UFXP is between 435\ and 7575\ times faster than NFXP, between 419\ and 5669\ times faster than CCP, between 59\ and 826\ times faster than SC, and between 7\ and 37\ times faster than MPEC. And increasing the number of optimization starts accentuates this advantage: across all $R = 21$ cases, UFXP is between 861\ and 13234\ times faster than NFXP, between 903\ and 9903\ times faster than CCP, between 116\ and 4414\ times faster than SC, and between 9\ and 64\ times faster than MPEC. (The MPEC times are moot, since that estimator did not yield sensible estimates.)
Table (ref) compares the estimators with two other performance metrics: the workload and the span. The workload is the average number of fixed points solved per optimization run. For example, UFXP has a workload of $m = 100$, since it solves $\lambda_{1}, \cdots, \lambda_{m}$ upfront and recycles these values across all $R$ jobs. The OUFXP workload is between 270 and 286, since the number of fixed points associated with $\tfrac{\partial}{\partial \theta}V^{\theta}_{\hat P}$ and with the second-stage dual variables varies with the neural network configuration.\footnote{ The OUFXP workload calculation is predicated on solving the second-stage fixed points for only the best of the first-stage estimates. } In contrast, the NFXP workload is between 10946\ and 300925,\footnote{ This actually understates the computational burden of NFXP. Following \cites{nfxp_manual} manual, we implement a “polygorithm” that combines value iterations (contraction mappings) and policy iterations (Newton-Kantorovich steps), only the latter of which comprise fixed point equations that contribute to our tables. For example, the average (Softplus Balanced, L-BFGS) NFXP run not only sequentially solves 77709\ fixed point equations, but also sequentially implements an additional 21764\ Bellman contractions. } and the CCP workload is between 1624\ and 188398. (And these workloads would have been much larger had we not used the dual-fixed point technique outlined in Section (ref) to streamline the corresponding gradient and Hessian computations.)
The second metric reported in Table (ref) reflects the fact that equations that can be solved in parallel are less burdensome than those that must be solved in sequence. Specifically, the span records the average length of an optimization run's computational critical path, under maximal parallelization---that is, the expected number of batches of fixed point equations that must be solved in sequence. For example, UFXP always has a span of 1, since $\lambda_{1},\cdots,\lambda_{m}$ can be solved in parallel. In contrast, OUFXP always has a span of 3, because it solves the fixed points that define (i) the first-stage dual variables, (ii) the Jacobian matrix $\tfrac{\partial}{\partial \theta}V^{\theta}_{\hat P}$ (which influences $Z^{\hat \theta_{\text{UFXP}}}$), and (iii) the second-stage dual variables. The fixed points in each of these three stages can be parallelized, but one stage can't begin until the preceding stage has ended, so switching from UFXP to OUFXP triples the span. However, switching from OUFXP to NFXP or CCP increases the span by more than three orders of magnitude.
We now present an empirical model that we believe would be intractable to estimate without our approach. Our model is an extension of the toy example studied in Section (ref). Whereas we previously estimated a single-stage inventory model, we will now estimate a multi-echelon supply chain model. And whereas we previously had three state variables, we will now have six (the three from before, plus in-transit inventory, upstream inventory, and upstream demand).
We will estimate our model with a 1,073,679-observation panel of daily data, spanning from February 2, 2011 to December 31, 2014. The sample describes the movement of 38 laundry detergent products through a supply chain that comprises one upstream distribution center (DC) and 67 downstream stores. (We focus on detergents, since their bulk makes them particularly difficult to store.) For each store--product--day, we observe: (i) retail prices, wholesale prices, and discounts, (ii) customer sales, (iii) the amount of inventory at the store, at the DC, and in the pipeline in between, and (iv) the store's inventory order, and the size and timing of the corresponding shipment from the DC.
A large Shanghai retailer owns all stores in our sample, but each store is overseen by a manager whose performance bonuses and promotion prospects depend on the functioning of their particular store. The managers, therefore, prioritize local store performance over global chain performance. For example, managers will generally squabble over scarce upstream inventory ration_gaming.
As before, the action variable is order quantity, $q \in \mathbb{Q} \equiv \{0, 6, 12, 20\}$, which the store manager chooses with the aim of maximizing expected discounted utility, under discount factor $\beta = 0.9997$. Each possible order quantity corresponds to an idiosyncratic utility shock, which we collect in a vector of independent Gumbels, $e \in \mathbb{R}^{|\mathbb{Q}|}$. Unobserved state vector $e$ is accompanied by observed state vector $x\equiv (r,o,i,j,k,\ell)$. The first three observed state variables are analogous to their Section (ref) counterparts: $r \in \{1, 2, 3, 4, 5\}$ indexes the mean demand for the focal detergent at the given store, $o \in \{1, 2, 3\}$ indexes the aggregate inventory of all but the focal detergent at the given store, and $i \in \mathbb{I}$ denotes the inventory level of the focal detergent at the given store, where $|\mathbb{I}| = 100$ and $\max \mathbb{I} = 149$.\footnote{ Capping $i$ at 149 truncates this variable’s top 10%. We perform this truncation because the right tail of $i$ spans an extremely large range. Indeed, we believe modeling the tail behavior of $i$ would likely require an entirely different approach. Thus, rather than estimate $h$ over the full domain of $i$, we more modestly estimate this function for non-tail values of $i$. To mitigate distortions associated with this truncation, we set the rows of $Z_1, \dots, Z_m$ that correspond to $i > \max \mathbb{I} - \max \mathbb{Q} = 149 - 20 = 129$ to zero. In other words, we assign zero empirical weight to observations for which the agent could order inventory beyond 149 units. }
The final three state variables describe the supply chain for the focal product: $j \in \mathbb{Q}$ specifies the amount of inventory in transit between the DC and the given store, $k \in \{1, 2, 3\}$ indexes the amount of inventory held at the distribution center, and $\ell \in \{1,2,3\}$ indexes the mean demand faced by the DC (i.e., the expected sum of store orders). Upstream demand $\ell$ is relevant because it is predictive of future upstream inventory: the DC roughly follows an $(s, S)$ inventory model, with demand state $\ell$ determining the slope of the line from $S$ to $s$. Upstream inventory $k$ is relevant because the DC is less likely to fulfill an order when its stocks run low ration_gaming: e.g., the DC fulfills 91% of orders when its inventory level is above the first decile but only 34% when its inventory level is below the first decile. And pipeline inventory $j$ is relevant because there is a lag between an order's placement and its fulfillment (unlike in the toy model): conditional on fulfillment, orders arrive after one day with probability 0.13, after two days with probability 0.67, and after three or more days with probability 0.20. To simplify our model dynamics, we suppose that all fulfillments take two days, so that today's $j$ reports the inventory the store will receive tomorrow in response to the order it placed yesterday.
The total utility from the focal product in the given day is
where $d$ is the daily customer demand, $e_{q}$ is an independent Gumbel shock associated with order quantity $q$, $\eta$ is the value of a sale, $h(r, o ,i)$ is the inventory holding cost, and $\kappa$ is the cost of receiving a shipment of inventory. Unfortunately, this utility does not fit into the framework presented in Section (ref), as it depends on demand $d$. However, we can integrate over this random variable to create the following equivalent utility function:
We can replace $\tilde{u}$ with $u$, because maximizing the expected discounted sum of $\tilde{u}_q(r, o, i, j, d) + e_{q}$ yields the same policy as maximizing the expected discounted sum of $u_q(r, o, i, j) + e_{q}$.
Utility function $u$ is the same as that used in Section (ref), except the shipping cost depends on $j$ rather than $q$, since not all orders translate into shipments. As before, the object of interest is holding cost function $h$, which we will model with neural networks. Note that state variables $k$ and $\ell$ do not influence current utilities, but they do influence future utilities, via the following transition kernel:
As you see, today's $\ell$ influences tomorrow's $k$, today's $k$ influences tomorrow's $j$, and today's $j$ influences tomorrow's $i$.
Transition kernels $f^r$, $f^o$, $f^k$, $f^\ell$ are general, but transition kernels $f^i$ and $f^j$ reflect the following laws of motion:
where random demand $d$ resolves from general distribution $f^{d}(d | r)$ and order fulfillment indicator variable $b$ resolves from a Bernoulli distribution $f^{b}$ whose mean depends on $k$.
Our sample does not include mean demands, so we begin by pre-estimating the expected inventory draw-downs of the DC and each store, by product and day.
We estimate the mean downstream demand using a neural network: First, we input a rich set of demand features ($y$) into a multilayer perceptron with two hidden layers of widths 64 and 32, and an output layer of width 1. Second, we exponentiate the neural network's output and treat the result as the mean of a Poisson distribution to evaluate the empirical likelihood of the observed sales. Specifically, we model the reduced-form probability of observing sales quantity $s$ given features $y$ as $\exp(m(y))^s \exp(-\exp(m(y))) / s!$, where $m(y)$ is the network's output for feature vector $y$. We train this model by maximizing the corresponding Poisson log-likelihood, $\sum_{y} \sum_{s} n_{ys} \log(\exp(m(y))^s \exp(-\exp(m(y))) / s!)$, where $n_{ys}$ is the number of observations with feature set $y$ and sales $s$. Finally, we set the expected mean demands to the corresponding fitted values, $\exp(m(y))$, and then map each mean demand to a demand index $r \in \{1, 2, 3, 4, 5\}$ (see Appendix Section (ref) for details).
We estimate the DC's mean demand analogously, but for this case we start with the upstream demand, which is the sum of store orders, and end with upstream demand state index $\ell \in \{1,2,3\}$.
In addition to pre-estimating demand means, to construct $r$ and $\ell$, we must also pre-estimate CCPs and state transition probabilities. For CCPs, we have too many states to use a simple frequency estimator, so we again use a neural network: First, we input the state vectors $x \in \mathbb{X}$ into a multilayer perceptron with two hidden layers of widths 64 and 32, and an output layer of width 4 (since $|\mathbb{Q}| = 4$). Second, we apply a softmax function to the neural network's output, and use the result to evaluate the empirical log-likelihood of the observed order quantities. Specifically, we model the reduced-form probability of ordering $q$ in state $x$ as $\exp(m_q(x)) / \sum_{q' \in \mathbb{Q}} \exp(m_{q'}(x))$, where $m_q(x)$ is the $q$th element of the network's output vector at state $x$. We train this model by maximizing the corresponding log-likelihood, $\sum_{x \in \mathbb{X}} \sum_{q \in \mathbb{Q}} n_{xq} \log(\exp(m_q(x)) / \sum_{q' \in \mathbb{Q}} \exp(m_{q'}(x)))$, where $n_{xq}$ is the number of state-$x$ and order-$q$ observations. Finally, we set $\hat{P}$ to the corresponding fitted values.
To be clear: we do not need our unnested estimators to accommodate these incidental networks, which are trained in an offline pre-estimation stage.
After pre-estimating CCPs, we then pre-estimate state transition probabilities. We estimate Poisson arrival rates for $f^d$ and Bernoulli probabilities for $f^{b}$, then combine these with the corresponding laws of motion to estimate $f^{i}$ and $f^{j}$ on the discretized state space (see Appendix Section (ref) for details). We estimate the other transition kernels, $f^r$, $f^o$, $f^k$, and $f^{\ell}$ nonparametrically with their empirical distributions. For example, upstream inventory $k$ gradually decreases, as the DC sends shipments to stores, and then jumps up, when the DC receives a shipment from its supplier; rather than separately estimate the DC's inbound and outbound shipments, we estimate the aggregate inventory adjustment---i.e., the final next-period inventory level, $k'$, given the current-period inventory level, $k$, and demand rate, $\ell$. As before, we combine estimates $\hat{f}^r$, $\hat{f}^o$, $\hat{f}^k$, $\hat{f}^\ell$, $\hat{f}^d$, and $\hat{f}^{b}$ into a corresponding $\hat{f}_{q}$ estimate.
After estimating the CCPs and the transition kernel, we are ready to estimate $\theta$. We estimate these structural parameters with UFXP. We decompose $h$ into a base inventory cost and a multiplicative congestion penalty: $h(r, o, i) = h^{i}(r, i) (1 + h^{o}(r, o))$, where functions $h^{i}(\cdot)$ and $h^{o}(\cdot)$ stem from separate neural networks, which together comprise 77 parameters (60 weights and 17 biases).
We set $h^{i}(r, i) \equiv \sum_{j = 1}^{i} w_j n(r, j)^2$, where $n(r, j)$ is the output of a multilayer perceptron that uses a single hidden layer of width 11 to map $r$, $i$, and interaction $ri$ to a scalar output, and $w_j$ is the fraction of discrete inventory levels that correspond to inventory state $j$. For example, $w_{80} = 2 / 150$: the denominator reflects the 150 possible inventory levels (0–149 units) and the numerator reflect the two levels (80 and 81 units) that map to state $i = 80$. By design $h^{i}(r, i)$ is non-decreasing in $i$, and is zero when $i$ is zero.
We set $h^{o}(r, o) \equiv \sum_{k = 2}^{o} m(r, k)^2$, where $m(r, k)$ is the output of another multilayer perceptron that uses a single hidden layer of width 4 to map $r$, $o$, and interaction $ro$ to a scalar output. Similarly, by design $h^{o}(r, o)$ is non-decreasing in $o$, and is zero when $o$ is one.
We initialize our neural networks with 1000 random starting values, following the methodology outlined in Section (ref). We use the Trust optimizer and recycle the same $\lambda_{1},\cdots, \lambda_{m}$ dual fixed points for each of the $R = 1000$ minimizations of $\mathcal{Q}_{Z}^{\theta}$. We construct the $m = 300$ weighting matrices, $Z_1, \dots, Z_m$, by drawing their elements from a standard normal distribution (with the proviso that the rows that correspond to $i > \max\mathbb{I} - \max \mathbb{Q} = 149 - 20 = 129$ are zeroed out, as explained in footnote (ref)).
We use the computer setup described in Section (ref), to make the estimation times roughly comparable to those reported in Section (ref).
We run the estimation with 1000 random starts. Our unnesting approach makes this multi-start strategy computationally feasible: it takes 208.4\ minutes to solve dual fixed points $\lambda_{1}, \cdots, \lambda_{m}$, but only 1.4\ minutes, on average, to convert them into a set of estimates, $\hat{\theta}_{\text{UFXP}}$.
Our estimates illustrate the intrinsic strangeness of neural networks. In most traditional settings, the estimator's criterion function will feature a unique global optimum and a smattering of clearly inferior local optima (if any). This is not the case for neural networks. Rather than a single point, neural networks often have a high-dimensional valley of nearly equivalent solutions. And this is precisely what we find: out of our 1000 optimization runs, 383\ yield a value that is within 0.5% of overall minimum $\mathcal{Q}_{Z}^{\hat{\theta}}$. Because optimality encompasses a region rather than a point, an honest assessment of the utility function requires an ensemble of estimates, which reinforces the need for many optimization starts.
Across the ensemble, the estimated shortage cost $\hat{\eta}$ averages 0.41\ with a standard deviation of 0.008, and the fixed ordering cost $\hat{\kappa}$ averages 6.79\ with a standard deviation of 0.003. We plot the holding cost functions that correspond with our ensemble of 383\ near-optimal estimates in Figure (ref). Despite some variation across the ensemble, an overall pattern clearly emerges: congestion state $o$ has a negligible impact on costs, whereas the demand state $r$ significantly alters the shape of $h$. The influence of $r$ is unsurprising because this variable captures structural heterogeneity across stores, with larger stores exhibiting consistently higher $r$ values.
Overall, the estimated holding cost function proves to be neither convex nor concave, let alone linear or quadratic. This last fact contradicts most of the literature on inventories, which generally imposes linear or quadratic forms.
The key to this work is the switch from the primal linear program
to the dual linear program
In the primal problem, utility parameter $\theta$ enters the constraint; in the dual, it enters the objective. Consequently, the primal solution varies with $\theta$, whereas the dual solution does not. By strong duality, both programs have the same value, implying
Thus, after solving dual variable vector $\lambda$, we can evaluate $w'V_{\hat{P}}^{\theta}$ for many values of $\theta$ with ease.
This observation motivates a dynamic discrete choice estimator that depends on value function $V_{\hat{P}}^{\theta}$ only through inner products $w_i' V_{\hat{P}}^{\theta}$, for $i = 1, \dots, m$, where $m$ is not much more than the length of the utility parameter vector, $\theta$. After computing the corresponding dual variables, $\lambda_1, \dots, \lambda_m$, such an estimator would be trivial to evaluate.
Restricting attention to linear functionals of $V_{\hat{P}}^{\theta}$ precludes likelihood-based estimators, which are intrinsically nonlinear. Thus, we base our estimator on the following first-order conditions:
First, we approximate these first-order conditions by replacing $P^{\theta}$ and $V^{\theta}$ with $\hat{P}$ and $V_{\hat{P}}^{\theta}$. And then we minimize the squared Euclidean norm of a random low-dimensional compression of this approximation:
In the expression above, it suffices to set $Z_{ia}$ to a length-$X$ vector of independent standard normal random variables, for $i \in \{1, \cdots, m\}$, and to set $m$ as small as $t + 1$, where $t$ is the length of $\theta$. Expressed this way, our objective does not appear any easier to evaluate than the CCP estimator's objective. However, unlike the CCP objective, we can entirely remove the value function from ours, since
As you see, our estimator still depends on the $m$ fixed points that define $\lambda_{1}, \cdots, \lambda_{m}$, but it does not depend on any nested fixed points (i.e., those whose solutions change with $\theta$). Of course, UFXP could not be efficient, as its random compression is ad hoc. However, the inefficient UFXP estimates give way to efficient OUFXP estimates.
Our unnested estimators are motivated by a related unnesting of the CCP fixed point that arises when the utility function is linear. In that case, the UFXP and CCP unnesting schemes are comparable. Our approach, therefore, meaningfully improves upon the state of the art only when the utility function is nonlinear. Fortunately, a promising class of nonlinear utility functions is drawing increasing attention: those arising from neural networks.
Embedding a neural network in a dynamic discrete choice model sacrifices not just linearity, but all semblance of smoothness---needless to say, the objective function will not be concave. This means the estimation problem will require brute force---initializing the optimizer with many different starting points. Our technique makes this practical, as each start reuses the same dual fixed point solutions, $\lambda_1, \dots, \lambda_m$. In contrast, traditional estimators must start from scratch for each set of initial values. To make this concrete, note that solving the dual fixed points took 208.4\ minutes for our 54000-state model, whereas converting these fixed points into utility parameter estimates took only 1.4\ minutes, on average. Hence, initializing our estimator with 1000 starting points took only $(1000 \cdot 1.4 + 208.4) / (1.4 + 208.4) = 7.7$\ times longer than initializing it with a single starting point. Compare this increase with the 1000-fold increase that would beset NFXP, CCP, SC, or MPEC if we started them with 1000 starting points. Our approach enables us to swarm the objective function with hundreds of cheap gradient ascent paths. UFXP likewise amortizes fixed point computations across optimizers and neural network specifications, making large-scale optimization ensembles computationally tractable.
And our estimators not only enable faster ensembles---they also enable faster individual runs. First, the UFXP fixed points can be computed in parallel, whereas the NFXP and CCP fixed points must be computed in sequence. Second, UFXP solves little more than $t$ fixed points, whereas traditional NFXP and CCP implementations solve $t$ fixed points for each gradient ascent, of which there could be hundreds. However, our dual fixed point technique reduces the number of fixed points required to evaluate the NFXP and CCP gradients from $t$ to one.
Now, for convenience, we will gather our our empirical results: