The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
77,677 characters
Scalable likelihood-based inference for limited dependent variable models
\maketitle
\begin{abstract}
\noindent \footnotesize
Limited dependent variable models are central to empirical economics, but likelihood-based inference is infeasible when likelihoods involve high-dimensional integration over latent variables. This paper proposes Stochastically Estimated Gradient Ascent (SEGA), a scalable estimation approach for limited dependent variable models. Using Fisher's identity, SEGA replaces the intractable likelihood score with an unbiased augmented-data score evaluated at a single conditional draw of the latent variables, and embeds this score in a stochastic gradient ascent algorithm. With sufficiently many iterations, we show that SEGA is asymptotically equivalent to the infeasible maximum likelihood estimator. A variance estimator based on Fisher's and Louis' identities is proposed that allows inference to proceed in the usual manner. Applications to brand choice and household demand demonstrate the usefulness of SEGA for conducting inference in large-scale discrete-choice and censored-demand models.
\end{abstract}
{\bf Keywords:} Limited Dependent Variable Models, Maximum Likelihood Estimation, Stochastic Gradient Ascent
\\
{\bf JEL Classification:} C13, C15, C24, C25, C63
\thispagestyle{empty}
\clearpage
\setcounter{page}{1}
\section{Introduction}
Limited dependent variable models (LDVMs) are central to empirical economics because they allow researchers to model outcomes that are discrete, censored, truncated, or otherwise restricted by their support \citep{tobin1958estimation}. Such restrictions arise naturally in settings where researchers study a range of decisions made by economic agents, such as choices, participation, purchase incidence, sales, or program uptake. Recent applications of LDVMs appear in labor economics \citep{bonhomme2022discretizing}, development economics \citep{bhattacharya2024demand}, health economics \citep{pakes2021unobserved}, education economics \citep{kline2016evaluating}, marketing \citep{danaher2020advertising}, and digital economics \citep{kanazawa2026ai}.
While LDVMs can be used to model data arising from many economic settings, inference within this class of models remains challenging. On the one hand, the limited observability of the dependent variable means that observations only partially reveal the economic behavior of interest. As a consequence, large datasets are required to infer unknown parameters of interest with enough precision to be economically informative. On the other hand, estimation of such models is cumbersome even in small datasets since the resulting likelihood functions for many LDVMs involve high-dimensional integrals that are unavailable in closed-form. Consequently, for the large datasets necessary to produce economically useful inferences, applying maximum likelihood estimation (MLE) to conduct inference using LDVMs is often practically infeasible. To circumvent these issues, researchers are often forced to impose restrictive structures and to use approximate methods without the inferential guarantees associated with MLE.
This paper removes the need to calculate the high-dimensional integral required to evaluate the likelihood function, enabling inference in economically important LDVMs. To this end, we propose Stochastically Estimated Gradient Ascent (SEGA), a scalable estimation method for LDVMs that admit a latent-variable representation. The key idea is to use Fisher's identity to express the score of the likelihood as a conditional expectation of the latent variable augmented score \citep{Cappe}. We then use \textit{a single draw} from the conditional distribution of the latent variables therefore provides an unbiased estimate of the score, irrespective of the sample size or dimension of the parameter space. We embed this noisy score in a stochastic gradient ascent (SGA) algorithm, and smooth the iterates from this algorithm to obtain a low variance estimate. For inference, we combine Fisher's identity with Louis' identity to construct a sandwich variance estimator that avoids direct evaluation of the observed score and Hessian \citep{louis1982finding}. SEGA is widely applicable because many LDVMs used in economics have the three features needed for fast implementation: a latent-variable representation, an augmented likelihood with an analytically tractable score, and latent variables that can be generated efficiently from the conditional distribution.
We illustrate the empirical payoff of SEGA with two large-scale applications. The first is a multinomial probit (MNP) model of pasta brand choice using more than one million purchase observations. In this setting, a brand purchase identifies the chosen alternative but not the full ranking of alternatives, so large samples are needed to recover substitution patterns. SEGA allows us to estimate the unrestricted model, construct confidence intervals for own- and cross-price elasticities, and formally test covariance restrictions that are commonly imposed for tractability. The second application is a random effects Tobit model of household-level pasta demand using more than 200,000 households observed over 62 weeks. In this setting, observed sales record whether demand is positive but not latent desired demand when purchases are censored at zero. SEGA allows us to study heterogeneity in baseline demand, price sensitivity, and state dependence, and to test whether these sources of heterogeneity are independent. In both applications, standard likelihood-based approaches are infeasible at the required scale. SEGA therefore expands the class of empirically relevant LDVMs for which researchers can conduct frequentist likelihood-based inference.
We make two main theoretical contributions. First, for fixed sample sizes, we prove that the SEGA iterates concentrate around the infeasible MLE as the number of iterations increases. Second, when the number of iterations grows sufficiently fast with the sample size, we show that SEGA is asymptotically equivalent at first-order to the infeasible MLE, despite using only a single simulated latent variable draw at each iteration. The novelty in this contribution comes from the fact that, in contrast to standard SGA algorithms, simulation noise in SEGA is generated by latent variable draws from a conditional distribution that changes with the current parameter value, and not from subsampling observations \citep{polyak1992acceleration,moulines2011non,chen2020statistical}. While the noisy score we use to implement SEGA is path-dependent, and must be handled appropriately, Fisher's identity ensures that the estimated score is an unbiased estimate of the infeasible score that we leverage to obtain our theoretical results. To our knowledge,
this is the first result showing that simulated latent variable scores can be used inside a stochastic optimization algorithm while remaining asymptotically equivalent to the infeasible MLE.
SEGA is closest in spirit to the method of simulated scores (MSS) of \citet{hajivassiliou1998method}. Despite its attractive theoretical properties, MSS has long been known to be difficult to implement in practice \citep{cameron2005microeconometrics,train2009discrete}. MSS approximates the likelihood score using simulated latent variables and then directly solves the resulting simulated score equations. This requires the simulated score function to remain stable across candidate parameter values, which is challenging in models where the latent variables are constrained by the underlying economic model, such as in high-dimensional choice models.
Our numerical experiments directly demonstrate the sensitivity of MSS to the number of simulation draws, the accuracy of the conditional simulator, and the initialization of the algorithm, generating substantial simulation noise and numerical instability even in moderate-dimensional LDVMs.
SEGA provides an alternative approach: rather than eliminating simulation noise before solving the score equations, it leverages stochastic approximation and iterative averaging to curtail the simulation noise across the optimization path.
In practice, applied researchers often use simulated maximum likelihood (SML) or numerical quadrature to approximate the integrals in LDVM likelihoods. SML replaces likelihood contributions, such as the MNP choice probabilities, with Monte Carlo approximations \citep{borsch1993smooth}. The simulated score is biased for a fixed number of simulation draws. Consistency therefore requires the number of simulation draws to increase with the sample size. This requirement becomes more demanding in MNP models as the number of alternatives grows, because the choice probabilities involve higher-dimensional integrals and more simulation effort is needed to maintain accuracy. Numerical quadrature is commonly used in random effects LDVMs \citep{rabe2005maximum}. Its cost grows often exponentially with the dimension of the random effects, making it impractical with several sources of unobserved heterogeneity. As a result, SML and quadrature methods are feasible only in low-dimensional LDVMs, or in settings where strong model-specific restrictions make the integrals manageable. SEGA removes this limitation by avoiding accurate likelihood approximation at each parameter value.
Although less commonly used with LDVMs, expectation-maximization (EM) is another possible method for obtaining point estimators \citep{dempster1977maximum}. Rather than approximating the likelihood directly, EM algorithms alternate between computing conditional expectations of the latent variables given the observed data and current parameter values, and maximizing the resulting expected augmented likelihood. In LDVMs, the E-step is often the computational bottleneck because the required conditional moments are not available in closed form. Monte Carlo EM approximates these moments by simulation \citep{wei1990monte}, but accurate E-steps can require many conditional draws at each iteration. In addition, each M-step may still involve solving a nontrivial optimization problem. Recent work has developed more specialized EM-type methods for probit models. For example, \citet{ding2024computationally} propose a deterministic approximation to the required conditional moments using expectation propagation, focusing on estimation rather than inference. SEGA instead targets large-scale inference for a general class of LDVMs.
The remainder of the paper is organized as follows. Section~\ref{sec:method} introduces the model class, the leading examples, and the SEGA estimator and variance estimator. Section~\ref{sec:application} presents two large-scale empirical applications. Section~\ref{sec:theory} establishes the finite-sample concentration and asymptotic equivalence results. Section~\ref{sec:experiments} evaluates the finite-sample performance of SEGA through Monte Carlo experiments. Section~\ref{sec:conclusion} concludes.
\section{Setup, examples and a new algorithm}\label{sec:method}
\subsection{Setting and model class}
We observe a sequence of outcomes $y_i$ and covariates $x_i$, with $d_i=(y_i,x_i)$ and $i=1,\dots,n$, generated independently from some unknown data generating process $P_0$. In LDVMs, the observed outcome variable $y_i$ represents a restricted transformation of some unobservable random variable $z_i$, which we refer to as a latent variable throughout. We write this observation rule as $y_i=\mathcal{T}(z_i)$, where $\mathcal{T}$ is known. The latent variable $z_i$ is unrestricted, but the transformation $\mathcal{T}$ induces discreteness, censoring, truncation, or other restrictions in the observed outcome. The distribution of the latent variable is governed by the unknown parameters $\theta\in\Theta\subset\mathbb{R}^{d_\theta}$.
We consider LDVMs where, given $(z_i,x_i)$, the distribution of $y_i$ is modeled using a class of probabilistic models which we represent through the following conditional latent variable representation: for each $i=1,\dots,n$, and independently across $i$,
\begin{align}\label{eq:loglikelihood}
y_{i}&\mid z_i,x_i \;{\sim}\; p_\theta(y_{i}\mid z_i,x_i) ,\\
{z_{i}}&{\mid x_{i}\;{\sim}\; p_\theta(z_{i}|x_{i}).}
\end{align}
In many LDVMs, $p_\theta(y_i\mid z_i,x_i)$ is degenerate and simply encodes the observation rule $y_i=\mathcal{T}(z_i)$. The key feature of LDVMs is that the latent variables $z_{i}$ are required to construct the model, but are typically not themselves the object of inference. Thus, the main task in this framework is inference on $\theta$.
The key difficulty with maximum likelihood estimation in LDVMs is that the latent variables must be integrated out to conduct inference on $\theta$. The likelihood function is
\begin{equation*}
p_\theta(y_{1:n}|x_{1:n}) = \int p_\theta(y_{1:n}|z_{1:n},x_{1:n})p_\theta(z_{1:n}|x_{1:n})\mathrm{d} z_{1:n},
\end{equation*}
where $p_\theta(y_{1:n}|z_{1:n},x_{1:n}) = \prod_{i=1}^np_\theta(y_{i}|z_i,x_i)$, $p_\theta(z_{1:n}|x_{1:n}) = \prod_{i=1}^np_\theta(z_{i}|x_{i})$, and $r_{1:n} =(r_1^\top,\dots,r_n^\top)^\top$ for some random variable $r_i$. Denoting $\ell_n(\theta) = \log p_\theta(y_{1:n}|x_{1:n})$, the maximum likelihood estimator is
\begin{equation*}
\hat{\theta}_n = \underset{\theta\in \Theta}{\text{argmax} } \ \ell_n(\theta).
\end{equation*}
The normalized score function of $\ell_n(\theta)$ can be represented as
\begin{align}\label{Eq:Score}
m_n(\theta):=\frac{1}{n}\nabla_\theta \ell_n(\theta)= \frac{1}{n}\nabla_\theta \log\left[ \int p_\theta(y_{1:n}|z_{1:n},x_{1:n})p_\theta(z_{1:n}|x_{1:n})dz_{1:n}\right].
\end{align}This score function depends on a high-dimensional integral, and has no analytical solution except in simple models. To account for this issue, various estimators that seek to approximate the MLE are often used when conducting inference on $\theta$. We review the most common approaches to these settings in the confines of two commonplace examples.
\subsection{Examples and existing likelihood-based estimators}\label{sec:examples}
\subsubsection{Multinomial probit model}\label{sec:examplesMNP}
Multinomial probit models are widely used to analyze discrete choice behavior, as they allow for flexible substitution patterns among choice alternatives; see for instance, \citet{geweke2003bayesian}, \citet{natenzon2019random}, and \citet{khan2021inference}. To introduce the model, denote $y_i$ to be a multinomial choice for individual $i=1,\dots,n$, where $y_{i}=j$ if individual $i$ chooses choice alternative $j=0,1,\dots,J$. Let $z_i=(z_{i1},\dots,z_{iJ})^\top$ be a $J$-dimensional vector of continuous random variables. {The elements in $z_i$ can be interpreted as utility differences with respect to the base category $j=0$, where the differencing identifies the location of the utilities \citep{bunch1991estimability}.}
The multinomial outcome $y_{i}$ is determined by the maximum value of $z_i$:
\begin{align}\label{eq:Y_i}
y_i =
\begin{cases}
0, & \text{if } \max_{1\leq k\leq J} z_{ik}<0,\\
j, & \text{if } z_{ij}=\max_{1\leq k\leq J} z_{ik}>0.
\end{cases}
\end{align}
The latent utilities are modeled as
\begin{align}\label{eq:Z_i}
z_{ij} = x_{ij}^\top\beta +\varepsilon_{ij}, \quad \varepsilon_{i}=(\varepsilon_{i1},\dots,\varepsilon_{iJ})^\top\sim N(0_J,\Sigma),
\end{align}
where $x_{ij}$ is an $r$-dimensional vector, $\beta$ is an $r$-dimensional vector of coefficients, $\varepsilon_{i}$ is a $J$-dimensional normally distributed disturbance vector with mean zero and covariance matrix $\Sigma$. We use $x_i$ to denote the collection of alternative-specific regressors, with $x_i\beta=(x_{i1}^\top\beta,\ldots,x_{iJ}^\top\beta)^\top$. The parameter vector is $\theta=(\beta^\top,\mathrm{vech}({\Sigma})^\top)^\top$.
Taking $d_i=(y_i,x_{i})$ as the observed data, the likelihood for the MNP can be stated through the choice probabilities associated with choice $j\in\{0,\dots,J\}$. Denoting $P_{ij}(\theta) = \Pr(y_i=j\mid x_i)$, we have that $P_{i0}(\theta)=\Phi_J(0_J;x_i\beta,\Sigma)$ indicates the probability of the base category, where $\Phi_J(\cdot;\mu,\Sigma)$ is the CDF of a $J$-variate normal with mean $\mu$ and covariance matrix $\Sigma$. For the remaining categories, we have that
\begin{flalign}
P_{ij}(\theta)&=\int \left[\mathbb{I}\left\{x_{ij}^\top \beta+\varepsilon_{ij}\ge0\right\}\prod_{k\ne j}\mathbb{I}\left\{x_{ij}^\top \beta+\varepsilon_{ij}>x_{ik}^\top \beta+\varepsilon_{ik}\right\}\right]\phi_J(\varepsilon_i;0_J,\Sigma)\mathrm{d}\varepsilon_i,\label{Eq:intprobMNP}
\end{flalign}
where $\phi_J(\cdot;\mu,\Sigma)$ is the density of a $J$-variate normal with mean $\mu$ and covariance matrix $\Sigma$. If $P_{ij}(\theta)$ could be evaluated, the log-likelihood would take the simple multinomial form $\ell_n(\theta)=\sum_{i=1}^{n}\sum_{j=0}^{J}\mathbb{I}(y_i=j)\log P_{ij}(\theta)$. However, the choice probabilities \(P_{ij}(\theta)\) involve integration of a multivariate normal density over a truncated region of the latent utility space. For \(J\ge 3\), these integrals are not available in closed form \citep{botev2017normal}.
The simulated maximum likelihood (SML) estimator approximates the integrals in \eqref{Eq:intprobMNP} by simulation methods:
$$
\hat\theta_{SML}:=\operatorname*{arg\,max}_{\theta\in\Theta}\frac{1}{n}\sum_{i=1}^{n}\sum_{j=0}^{J}\mathbb{I}(y_i=j)\log \widehat{P}_{ij}(\theta),
$$
where $\widehat{P}_{ij}(\theta)$ denotes the estimate for ${P}_{ij}(\theta)$ \citep{bolduc1999practical}. However, even with efficient simulation and/or approximation methods, conducting exact estimations of the MNP model in modern large-scale data sets is challenging \citep{ding2024computationally}. Indeed, the approximate nature of SML results in biased estimators of the score, with this bias often exacerbated in larger choice sets and or large datasets \citep{borsch1993smooth}. An alternative to simulating the likelihood is to simulate the scores directly, ala MSS \citep{hajivassiliou1998method}, we defer a discussion of this approach to Section \ref{sec:MSS_discuss}.
\subsubsection{Random effects Tobit model}\label{sec:examplesTobit}
Random effects Tobit models are widely used to analyze censored outcomes in the presence of unobserved heterogeneity, with recent applications in marketing \citep{danaher2020advertising} and banking \citep{liu2023forecasting}. To introduce the model, let $n$ denote the total number of individuals and $T$ the total number of observations available for each individual. We denote $y_{it}\ge 0$ to be the limited dependent variable with $i =1,\dots,n$, and $t=1,\dots,T$.
The variable $y_{it}$ is modeled to be a function of the latent variable $y_{it}^*$, such that
\begin{equation}\label{Eq:linktobit}
p(y_{it}|y_{it}^*)=\left\{\begin{matrix} \mathbb{I}(y_{it}=0) & y_{it}^*\le 0,\\
\mathbb{I}(y_{it}=y_{it}^*)&y_{it}^*>0.\\\end{matrix}\right.
\end{equation}
This latent variable follows the normal distribution
\begin{equation}\label{Eq:yystar}
p_\theta(y_{it}^*|\alpha_i,x_{it}) = \phi_1\left(y_{it}^*;h_{it}^\top\beta+w_{it}^\top\alpha_i,\sigma^2\right),
\end{equation}
where the vector $h_{it}$ includes covariates with fixed effects $\beta$, and the vector $w_{it}$ includes covariates with random effects $\alpha_i\sim N(0_r,\Omega)$. The parameter vector is $\theta = (\beta^\top,\sigma^2,\text{vech}(\Omega)^\top)^\top$, where $\text{vech}(A)$ extracts the lower triangular elements of matrix $A$.
Taking $y_{i} = \left(y_{i1},\dots,y_{iT}\right)^\top$, $x_{i} = (x_{i1}^\top,\dots,x_{iT}^\top)^\top$, with $x_{it} = (h_{it}^\top,w_{it}^\top)^\top$, and $d_i = (y_i,x_i)$, the likelihood increment can be stated as
\begin{equation}\label{Eq:intprobTobit}
p_\theta(y_i|x_i) = \int g_\theta(\alpha_i,d_i)\phi_r(\alpha_i,0_r,\Omega)d\alpha_i.
\end{equation} where
\begin{equation*}
g_\theta(\alpha_i,d_i) = \prod_{\{t:y_{it}=0\}}\Phi_1\left(0;h_{it}^\top\beta+w_{it}^\top\alpha_i,\sigma^2\right) \prod_{\{t:y_{it}>0\}}\phi_1\left(y_{it};h_{it}^\top\beta+w_{it}^\top\alpha_i,\sigma^2\right),
\end{equation*}
Because $g_\theta(\alpha_i,d_i)$ contains products of normal CDFs and must be integrated over the random effects, $p_\theta(y_i\mid x_i)$ is not available in closed form except in special cases.
To tackle this problem, Gaussian-Hermite quadrature methods can be used, which yield the approximation
\begin{equation}\label{quadrature}
\widehat{p_\theta(y_i|x_i)} =
\frac{1}{\pi^{r/2}}
\sum_{q \in \{1,\dots,Q\}^r}
\left(
\prod_{j=1}^{r} \omega_{q_j}
\right)
g_\theta\!\left(\sqrt{2}\,\tilde{D}\zeta(q),d_i\right),
\end{equation}
where $\Omega = \tilde{D}\tilde{D}^\top$, $\zeta(q) = (\zeta_{q_1},\dots,\zeta_{q_r})$ and where $\{\zeta_s,\omega_s\}_{s=1}^Q$ denote the root, weight pairs of the Hermite polynomial of order $Q$ \citep{skrondal2004generalized}. A value $Q=30$ is often considered reasonable to accomplish high accuracy. However, note that the function $g_\theta(\cdot,d_i)$ must be evaluated a total $Q^r$ times, rendering quadrature methods practically infeasible beyond three dimensional random effects problems.
\subsection{Stochastically estimated gradient ascent}
As is clear from Sections \ref{sec:examplesMNP} and \ref{sec:examplesTobit}, in particular equations \eqref{Eq:intprobMNP} and \eqref{Eq:intprobTobit}, parameter estimation in LDVMs is made difficult through their dependence on the latent variables that must be integrated out before the likelihood function can be evaluated. Consequently, the likelihood score equations in \eqref{Eq:Score} depend on integrals which must generally be approximated.
If we are willing to conduct inference on the unknown parameters via an approximation to the MLE, then we can circumvent these computational difficulties and apply LDVMs to high-dimensional settings with millions of observations and hundreds of parameters. To do so, our starting point is to re-examine the score equations in \eqref{Eq:Score}, and apply integration by parts to see that (see, e.g., \citealp{Cappe})
\begin{flalign}
m_n(\theta)&=\frac1n\nabla_\theta
\log\left[\int p_\theta(y_{1:n}\mid z_{1:n},x_{1:n})p_\theta(z_{1:n}\mid x_{1:n})dz_{1:n}\right]\nonumber\\
&=\frac1n\int_{}\nabla_\theta \log p_\theta(y_{1:n},z_{1:n}\mid x_{1:n})p_\theta(z_{1:n}\mid d_{1:n})\mathrm{d} z_{1:n}.\label{Eq:FisherScore}
\end{flalign}
The gradient in \eqref{Eq:FisherScore} makes clear that if we can obtain draws from $p_\theta(z_{1:n}\mid d_{1:n})$, we could attempt to conduct inference on $\theta$ by estimating the score equations and applying modern gradient-based optimization methods.
While we cannot evaluate $m_n(\theta)$ in general, a feasible estimator exists:
\begin{align}\label{Eq:FisherScoreest}
\widehat{m}_n(\theta)= \frac{1}{n}\nabla_\theta \log p_\theta(y_{1:n},\tilde{z}_{1:n}\mid x_{1:n}),\quad \tilde{z}_{1:n}{\sim} p_\theta(\cdot\mid d_{1:n}).
\end{align}
While such an estimator is of course high-variance, it is nonetheless unbiased. Since $\widehat{m}_n(\theta)$ is an unbiased estimator of $m_n(\theta)$, we use $\widehat{m}_n(\theta)$ within a stochastic gradient ascent search algorithm to produce stochastic iterates of $\theta$. SGA is the backbone of modern machine learning algorithms due to its ability to scale to vast amounts of data, and millions of parameters. However, our use of the estimator $\widehat{m}_n(\theta)$ within an SGA algorithm is non-standard: SGA algorithms are ``stochastic'' not because they are based on gradients that depend on simulated data, but because they estimate the gradient through subsets (called mini-batches) of the data that are drawn randomly from the full data set \citep{robbins1951stochastic}.
In contrast to standard SGA, our implementation uses the entire dataset, but is based on estimating the gradient stochastically using a subset of the possible data that could have been simulated from $p_\theta(z_{1:n}\mid d_{1:n})$. To clarify this distinction, we refer to SGA algorithms based on the estimated gradient $\widehat{m}_n(\theta)$ as stochastically estimated gradient ascent (SEGA) algorithms.
Pseudo-code for the implementation of SEGA is given in Algorithm \ref{alg:sega}. The $k$-th parameter iterate in the SEGA path is denoted as $\theta^{(k)}_n$. The iterate $\theta^{(k)}_n$ is noisy, due to stochasticity induced by the latent variables. To account for this, more precise estimates can be obtained by smoothing the iterates to obtain
\begin{align}\label{eq:smooth}
\overline\theta_n=\frac{1}{k-k_0}\sum_{j=k_0+1}^{k}\theta^{(j)}_n,
\end{align}
where $1\le k_0\ll k$ denotes some initial ``burn-in'' period after which we smooth the iterates. The averaging step in \eqref{eq:smooth} is known as Polyak–Ruppert averaging \citep{polyak1992acceleration,ruppert1988efficient}. We refer to $\overline\theta_n$ as the SEGA estimator.
\begin{algorithm}[tb!]
\begin{algorithmic}[1]
\State{Initialize $\theta^{(0)}_n$, $k_0$, and $k$.}
\For{$j=0,\dots,k$}
\State{Draw $\widetilde{z}^{(j)}_{1:n}{\sim} p_{\theta^{(j)}_n}(\cdot\mid d_{1:n})$.}
\State{Calculate $\widehat{m}^{}_{n}({\theta}_n^{(j)})=\frac{1}{n}\nabla_\theta \log p_{ {\theta}_n^{(j)}}(y_{1:n},\widetilde{z}^{(j)}_{1:n})$.}
\State{Update
$\theta_n^{(j+1)}= \theta_n^{(j)}+\eta^{(j+1)}_{n}\odot\widehat{m}^{}_n(\theta^{(j)}_n)$.}
\EndFor
\State \Return{$\overline\theta_n=\frac{1}{k-k_0}\sum_{j=k_0+1}^{k}\theta^{(j)}_n$.}
\end{algorithmic}
\caption{Stochastically estimated gradient ascent algorithm}
\label{alg:sega}
\end{algorithm}
Step 5 in Algorithm \ref{alg:sega} requires a learning rate vector $\eta_n^{(j+1)}$. We use the ADADELTA method in \citet{zeiler2012adadelta}, which conducts exponential smoothing operations on the gradient. For $[\eta^{(j)}_{n}]_i$ denoting the $i$-th dimension of $\eta^{(j)}_{n}$, we have
\begin{equation}
[\eta^{(j)}_{n}]_i = \frac{\sqrt{\Delta_i^{(j-1)}+\epsilon}}{\sqrt{S_i^{(j)}+\epsilon}},
\end{equation}
where $\Delta^{(j)}=\rho \Delta^{(j-1)}+(1-\rho)\left\{\eta^{(j)}_{n}\odot\widehat{m}_n(\theta^{(j)}_n)\right\}^2$ and $S^{(j)} = \rho S^{(j-1)}+(1-\rho)\widehat{m}_n(\theta^{(j)}_n)^2$,
with $\rho=0.95$, $\epsilon=10^{-6}$, and $\Delta_i^{(0)}=S_i^{(0)}=0$. This learning-rate is self-adaptive, and the magnitude of $\Delta^{(j)}$ decreases on average with $j$, which in theory allows for a stopping rule on the algorithm based on $||\Delta^{(j)}||$. In practice, such statistics are noisy, and setting a large number of iterations after which the change in the smoothed parameters is checked for convergence works best.
\subsection{SEGA implementation in the two examples}\label{sec:implementation_examples}
SEGA requires evaluation of the augmented log gradient $\nabla_\theta \log p_\theta(y_{1:n},z_{1:n}\mid x_{1:n})$ in closed-form and simulation from the distribution $p_\theta(z_{1:n}\mid d_{1:n})$. This feature allows for the direct application of SEGA in both classes of examples considered in Section \ref{sec:examples}.
\subsubsection{SEGA for the multinomial probit model}
Since the scale of the covariance matrix of the latent utilities is not identified in a
multinomial probit model, we follow \cite{ding2024computationally} and fix the trace of the precision matrix so that $\text{Trace}(\Sigma^{-1}) = J$, and where $\Sigma^{-1} = CC^\top$. Denote the vector $\psi = \text{vech}(C)$, which is a vector of dimension $d_C=J(J+1)/2$. The trace restriction imposes $\text{Trace}{(\Sigma^{-1})} = \sum_{i=l}^{d_C} \psi_l^2= J$. Similar to \citet{loaiza2021scalable} we parametrize $\psi$ in terms of a spherical coordinate system with radius $\sqrt{J}$ and angles $\kappa = (\kappa_1,\dots,\kappa_{n-1})^\top$, where
\begin{equation}\label{eq:transformation}
\psi_{l}(\kappa) = \begin{cases}
\sqrt{J}\cos\kappa_{1} & \text{for $l=1$},\\
\sqrt{J}\cos\kappa_{l}\prod_{j = 1}^{l-1}\sin\kappa_{j} & \text{for $1<l<d_C$},\\
\sqrt{J}\prod_{j = 1}^{l-1}\sin\kappa_{j} & \text{for $l=d_C$}.
\end{cases}
\end{equation}
To ensure positive diagonal elements of $C$, set $0<\kappa_{\frac{J(J+1)}{2}-\frac{(J-i+1)(J-i+2)}{2}+1}\le \frac{\pi}{2}$ for $i = 1,\dots,J-1$, and constrain the remaining angles to $[0,\pi]$. The angles are transformed into the real line by $\kappa_l = \Phi(\xi_l)\text{UB}_l$, where $\text{UB}_l$ denotes the upper bound of the $l^{\text{th}}$ angle. Note that $\Sigma$ is now implicitly a function of $\xi$, and $\theta = \left(\beta^\top,\xi^\top\right)^\top$.
Under this parametrization, the MNP model in \eqref{eq:Y_i} and \eqref{eq:Z_i} delivers a closed-form estimator of the score in \eqref{Eq:FisherScoreest}. The required gradient and distribution are
\begin{flalign}
\nabla_\theta \log p_\theta(y_{1:n},z_{1:n}\mid x_{1:n})&=\sum_{i=1}^n (\nabla_\beta \log \phi_J(z_i;x_i\beta,\Sigma)^\top,\nabla_{\xi} \log \phi_J(z_i;x_i\beta,\Sigma)^\top)^\top, \label{eq:grad_augpost_mnp}\\
p_\theta(z_{1:n}\mid d_{1:n})&=\frac{\prod_{i=1}^n\phi_J(z_i;x_i\beta,\Sigma)I[A(y_i)z_i<0]}
{\int \prod_{i=1}^n\phi_J(z_i;x_i\beta,\Sigma)I[A(y_i)z_i<0]\mathrm{d} z_{1:n}}, \label{eq:pz_mnp}
\end{flalign}
Thus, $p_\theta(z_{i}\mid d_{i})$ is multivariate truncated normal with truncation matrix $A(y_i)=I_J$ if $y_i=0$ and $A(y_i)=I_J+(-e_j-1)e_j'$ otherwise, where $e_j$ is a zero vector with a one as the $j$th element. This distribution can be simulated from via Gibbs methods. Further details are provided in Supplemental Appendix~\ref{A:mnp_fisher}.
\subsubsection{SEGA for the random effects Tobit model}
In the random effects Tobit model, the latent variables $z_{1:n}$ are represented by the random effects $\alpha_i$ and the latent variables $y_{it}^\ast$. To facilitate optimization, all parameters are transformed to the real line. The variance parameter is transformed as $c = \log(\sigma^2)$. The covariance matrix of the random effects is written down in terms of its Cholesky decomposition $\Omega = \tilde{D}\tilde{D}^\top$, where $\tilde{D}$ is a lower triangular matrix with diagonal elements $\tilde{d}_{i,i} = \exp(\delta_{i,i})$, and off diagonal elements $\tilde{d}_{i,j} = \delta_{i,j}$.
Denote $D$ to be the lower triangular matrix comprised of the elements $\delta_{i,j}$, and define
$\delta= \text{vech}(D)$ to be the vector of unique parameter vectors that characterize $\Omega$. The parameter vector is $\theta = (\beta^\top,c,\delta^\top)^\top$.
While of a different nature to the MNP model, the random effects Tobit model in \eqref{Eq:linktobit} and \eqref{Eq:yystar} also delivers the gradient and distribution required by the estimator of the score in~\eqref{Eq:FisherScoreest}:
\begin{flalign}\label{eq:tobit_sega}
\nabla_\theta \log p_\theta(y_{1:n},z_{1:n}|x_{1:n})&=\sum_{i=1}^n\left(\sum_{t=1}^T\nabla_{\beta,c}\log\phi_1(y_{it}^*;\mu_{it},\sigma^2)^\top,\nabla_{\delta}\log\phi_r(\alpha_i;0_r,\Omega)^\top\right)^\top,\notag \\
p_\theta(y_{1:n}^*\mid d_{1:n},\alpha_{1:n})&=\prod_{\{(i,t):y_{it}=0\}}\frac{\phi_1(y_{it}^*;\mu_{it},\sigma^2)I(y_{it}^*<0)}
{\Phi_1(0;\mu_{it},\sigma^2)}\prod_{\{(i,t):y_{it}>0\}} I(y_{it}^*=y_{it}),\notag \\
p_\theta(\alpha_{1:n}\mid d_{1:n},y^*_{1:n})&=\prod_{i=1}^n\phi_r(\alpha_i;\bar{\alpha}_i,V_{i}),
\end{flalign}
\sloppy with $\mu_{it}=h_{it}^\top\beta+w_{it}^\top\alpha_i$, $\bar{{\alpha}}_i =\frac{1}{\sigma^2}V_i\left(\sum_{t=1}^{T}w_{it}{y}_{it}^*-\left[\sum_{t=1}^{T}w_{it}h_{it}^\top\right]\beta\right)$, and $V_i = \left[\Omega^{-1}+\frac{1}{\sigma^2}\sum_{t=1}^{T} w_{it}w_{it}^\top\right]^{-1}$. The latent variables $y^*_{1:n}$ can be simulated from a univariate truncated normal and $\alpha_{1:n}$ from a multivariate normal via Gibbs methods. Further details are provided in Supplemental Appendix~\ref{A:tobit_fisher}.
\subsection{Inference after SEGA}\label{sec:inference}
Although the SEGA estimator is obtained from stochastic gradients rather than direct likelihood maximization, in Section~\ref{sec:theory} we show that $\bar\theta_n$ in \eqref{eq:smooth} is asymptotically equivalent to the infeasible MLE under regularity conditions. Hence, inference can be based on the same sandwich covariance matrix as the infeasible MLE; namely, $V=\mathcal{H}^{-1}\mathcal{I}\mathcal{H}^{-1}$, where $\mathcal{H}=-\mathbb{E}\nabla_{\theta}^2\log p_\theta (d_i)$ is the expected Hessian of the log likelihood and $\mathcal{I}=\text{Var}\{\nabla_\theta \log p_\theta (d_i)\}$ the variance of the score.
\subsubsection{A variance estimator for SEGA}
Estimating $V$ requires estimates of $\mathcal{H}$ and $\mathcal{I}$, which are not directly available because the observed likelihood score and Hessian involve integrals over latent variables. We therefore estimate these quantities using Louis' identity \citep{louis1982finding} and Fisher's identity, respectively, where conditional expectations are approximated by simulation.
Given simulated data, $z_{i,s}\stackrel{iid}{\sim} p_\theta(z_i|d_i)$, for $\ell(\theta
,z_i) = \log p_\theta(y_i,z_i|x_i)$, and $\widehat{\nabla_\theta\ell_i(\theta)} = \frac{1}{S}\sum_{s=1}^S\nabla_\theta\ell(\theta,z_{i,s})$, our proposed estimator is
\begin{align}\label{eq:Vhat}
\widehat{V} = \widehat{\mathcal{H}(\bar{\theta}_n)}^{-1}
\widehat{\mathcal{I}(\bar{\theta}_n)}
\widehat{\mathcal{H}(\bar{\theta}_n)}^{-1},
\end{align}
where
\begin{align*}
\widehat{\mathcal{I}(\theta)} =& \frac{1}{n}\sum_{i=1}^n\widehat{\nabla_\theta\ell_i(\theta)}\widehat{\nabla_\theta\ell_i(\theta)}^\top-\left(\frac{1}{n}\sum_{i=1}^n\widehat{\nabla_\theta\ell_i(\theta)}\right)\left(\frac{1}{n}\sum_{i=1}^n\widehat{\nabla_\theta\ell_i(\theta)}\right)^\top,\\
\widehat{\mathcal{H}(\theta)} =& \frac{1}{n}\sum_{i=1}^n \frac{1}{S}\sum_{s=1}^S\nabla_\theta^2 \ell(\theta,z_{i,s})+
\nabla_\theta \ell(\theta,z_{i,s})\nabla_\theta \ell(\theta,z_{i,s})^\top-\widehat{\nabla_\theta\ell_i(\theta)}\widehat{\nabla_\theta\ell_i(\theta)}^\top.
\end{align*}
Each of the terms can be estimated separately using the simulated data $z_{i,s}$. Fisher's identity allows the observed score for each unit to be recovered by averaging the augmented score $\nabla_\theta \log p_\theta(y_i,z_i|x_i)$ over independent draws from the conditional distribution of the latent variables $p_\theta(z_i|d_i)$. The estimator $\widehat{\mathcal{I}}(\theta)$ is therefore the sample variance of these estimated observed scores. To estimate $\mathcal{H}$, we use Louis' identity, which expresses the observed Hessian as the conditional mean of the augmented Hessian $\nabla^2_\theta \log p_\theta(y_i,z_i|x_i)$ plus a correction for the conditional variation of the augmented score.
Tests and confidence intervals can then be based on $\widehat V$, while for nonlinear economic quantities, such as elasticities and marginal effects, we use a \citet{krinsky1986approximating} parametric bootstrap: we draw parameters from the estimated asymptotic distribution of $\overline\theta_n$ and transform each draw into the corresponding quantity of interest.
\subsubsection{Evaluating the variance estimator in the SEGA examples}
The evaluation of $\widehat V$ requires $S$ independent draws from $p_\theta(z_i|d_i)$. Since the Gibbs samplers used in Section~\ref{sec:examples} for SEGA point estimation produce autocorrelated draws, they cannot be used for estimation of $\widehat{V}$ without modifications.
For the multinomial probit model, we use the accept-reject sampling method proposed by \citet{botev2017normal} for generating i.i.d.\ draws from the truncated multivariate normals implied by $p_\theta(z_{1:n}\mid d_{1:n})$. This approach is computationally more demanding than the Gibbs sampler used in SEGA, and is therefore solely for the construction of $\widehat V$.
In the random effects Tobit model, we employ a two step generation process based on $p_\theta(z_i|d_i) = p_\theta(y_i^*|d_i)p_\theta(\alpha_i|d_i,y_i^*)$. The density $p_\theta(y_i^*|d_i)$ can be written as a truncated multivariate normal, for which we can again use the \citet{botev2017normal} method for generating i.i.d.\ draws. The density $p_\theta(\alpha_i|d_i,y_i^*)$ is multivariate normal and also part of the Gibbs sampling scheme in SEGA.
Given the simulated latent variables, both models reduce to Gaussian augmented likelihoods, so the augmented score $\nabla_\theta \log p_\theta(y_i,z_i|x_i)$ and augmented Hessian $\nabla^2_\theta \log p_\theta(y_i,z_i|x_i)$ entering Fisher's and Louis' identities are available in closed form. Details are provided in Supplemental Appendix~\ref{A:MNPse} and~\ref{A:Tobitse}.
\subsection{SEGA and the method of simulated scores}\label{sec:MSS_discuss}
In contrast to standard simulation-based estimators, SEGA attempts to solve the score equation $m_n(\theta) = 0$ stochastically, rather than deterministically. Thus, SEGA iterates will not satisfy the equation $m_n(\theta^{(k)}_n)=0$ for any finite $k$, and in general neither will $\overline\theta_n$. As such, SEGA iterates cannot be interpreted as Z-estimators (\citealp{van2000asymptotic}, Ch. 5) or simulated method of moments estimators \citep{pakes1989simulation}.
Nevertheless, SEGA estimators share a similar philosophy to simulation-based econometric estimators, such as the MSS estimators presented in \cite{hajivassiliou1996simulation}. {Indeed, the replacement of the intractable score $m_n(\theta)$ with an average over a large number of simulated draws is the key idea in MSS. In particular,} the MSS estimator is also based on an approximation to the scores of LDVMs, and can be thought of as trying to estimate $\theta$ by solving a simulated version of the scores
\begin{equation}
\widehat\theta_{\mathrm{MSS}}^{(R)}\in\left\{\theta\in\Theta:\frac{1}{n}\sum_{i=1}^n\frac{1}{R}\sum_{s=1}^R\widetilde m_i\{\theta,z_i(\theta;\epsilon_{i,s})\}=0\right\},\label{eq:MSS}
\end{equation}
where $z_i(\theta;\epsilon_{i,s})$ is drawn from $p_\theta(z_{i}\mid d_{i})$ based on the fixed set of random numbers $\epsilon_{i,s}$, for $i\in\{1,\dots,n\}$ and $s\in\{1,\dots,R\}$, that do not change across the algorithm, and where $\tilde{m}_i\{\theta,z_{i}(\theta;\epsilon_{i,s})\}$ denotes a per-unit estimator of the corresponding per-unit infeasible score equations $m_i(\theta)$, where $m_n(\theta)=\frac{1}{n}\sum_{i=1}^{n}m_i(\theta)$.
For a given $\theta$ value, implementation of MSS requires generating $z_i(\theta;\epsilon_{i,s})$ using a set of fixed random numbers $\epsilon_{i,s}$ to ensure that $z_i(\theta;\epsilon_{i,s})$ cannot change in a stochastic manner. This requirement is necessary for consistent estimation in MSS. As can be seen by \eqref{eq:MSS}, the MSS estimator is a Z-estimator based on numerically solving an approximation to the infeasible score equations $0=m_n(\theta)$. Consequently, as has been known since at least \cite{pakes1989simulation}, if MSS allows the random draws $\epsilon_{i,s}$ to change across evaluations of the criterion in \eqref{eq:MSS}, the MSS estimator $\widehat\theta_{\mathrm{MSS}}^{(R)}$ will not be consistent.
However, simulating draws $z_i(\theta;\epsilon_{i,s})$ that only depend on the fixed set of random numbers $\epsilon_{i,s}$ is no simple task in LDVMs. In general, such an approach is infeasible in any setting where generating random variables requires accept/reject sampling or other techniques that do not allow for the random numbers to be fixed in their implementation.
In cases where Gibbs sampling is used to simulate draws of $z_i(\theta;\epsilon_{i,s})$, the accuracy of the resulting parameter estimates depend on the interaction between the sample size, $n$, the number of simulations used to estimate the gradient, $R$, and the number of steps used in the Gibbs sampler to draw each $z_i(\theta;\epsilon_{i,s})$, as well as its initialization.
Using a Gibbs sampler with MSS that is poorly initialized, or that uses an insufficient number of steps, results in simulated variables $z_i(\theta;\epsilon_{i,s})$ that will not be drawn from the correct distribution and the resulting $\widehat\theta_{\text{MSS}}^{(R)}$ can be a poor approximation of the MLE. Furthermore, all else equal, the larger the dimension of $z_i(\theta;\epsilon_{i,s})$, the larger the number of steps we are required to use within the Gibbs sampler to ensure that we are drawing from the correct distribution; initialization of the Gibbs sampler also becomes more costly as the dimension increases. In our numerical experiments conducted in Section \ref{sec:MSS_Emp}, we demonstrate that even in moderate dimensions the MSS estimator requires a large number of steps within the Gibbs sampler to accurately approximate the MLE, which becomes computationally prohibitive even in moderate dimensions.
The requirement that the latent variables be simulated as $z_i(\theta;\epsilon_{i,s})$, with $\epsilon_{i,s}$ fixed across evaluations, is completely absent in SEGA: each $\widetilde{z}^{(k)}_{i}\stackrel{iid}{\sim} p_{\theta^{(k)}_n}(\cdot\mid d_{i})$
is generated based on new random numbers at each iteration. SEGA bypasses this requirement since its iterations are, by construction, random and since SEGA only seeks to guarantee that the limit, as $k$ diverges, of the random solution path solves $m_n(\theta)=0$.
\section{Empirical applications}\label{sec:application}
To illustrate the empirical relevance of SEGA, we apply the method to two large-scale LDVMs for which likelihood-based inference is infeasible with existing approaches: a multinomial probit model for pasta brand choice and a random effects Tobit model for household-level pasta demand. We implement SEGA with 100,000 iterations, use the final 25{,}000 iterations to construct the smoothed estimates, and use $S=20{,}000$ draws to construct the variance estimator. Details are deferred to Supplemental Appendix~\ref{A:mnp_application} and \ref{A:tobit_application}.
\subsection{Brand choice and substitution in the pasta market}\label{sec:app_mnp}
We apply the MNP model to a data set on pasta brand purchases using SEGA. This consumer choice data set, collected at a leading U.S. grocer, includes more than one million purchases and is made available by the Dunnhumby data platform\footnote{https://www.dunnhumby.com/source-files/} as ``Carbo-Loading: A Relational Database". The final sample consists of purchases without coupons of the ten top-selling pasta brands, excluding private labels, amounting to 1,070,436 observations. We include an intercept and the log price for each brand in the model. \citet{loaiza2023fast} and \citet{loaiza2024hybrid} fitted this MNP model to this data using approximate Bayesian estimation methods.
Application of the MNP model to data problems of this magnitude -- large choice set and more than a million observations -- with the theoretical guarantees associated with maximum likelihood estimation are only feasible through the use of SEGA.
Due to the differenced utility representation in the MNP model, the model parameters do not provide direct economic interpretation. Instead, empirical interest typically focuses on the price elasticities and the substitution patterns implied by the model parameters.
Figure~\ref{fig:ped} shows the estimated own-price elasticity of the pasta brand `Barilla' in Panel (a) and the cross-price elasticity of `San Giorgo' with respect to the price of `Barilla' in Panel (b). The prices of the other pasta brands are fixed at their mean. The solid lines represent the point estimates for the price elasticities, and the dashed lines the 95\% confidence intervals. As expected, the estimated own-price elasticity is negative, and the confidence intervals are narrow due to the large sample size. The positive cross-price elasticity indicates that increases in the price of Barilla raise the probability of purchasing San Giorgio, consistent with the two brands being substitutes.
\begin{figure}[tb!]
\caption{Price elasticities for two pasta brands}
\centering
\includegraphics*[width=\textwidth]{figures_application/ped.eps}
\begin{flushleft}
\footnotesize
This figure shows the price elasticity of the pasta brands `Barilla' (Panel (a)) and `San Giorgo' (Panel (b)) as a function of the price of Barilla, with the prices of the other pasta brands fixed at their mean. The solid lines represent the price elasticities and the dashed lines the 95\% confidence intervals.
\end{flushleft}
\label{fig:ped}
\end{figure}
Empirical applications of the MNP model often impose covariance restrictions. \citet{paetz2018utility} discuss the implied substitution patterns of an identity covariance matrix, \citet{rossi2012bayesian} justify a diagonal matrix by arguing that choice behavior is often characterized by large differences in relative variance between choice alternatives, and \citet{geweke1994alternative} show that assuming an identity covariance matrix for the undifferenced utilities implies an equicorrelated matrix on the differenced utilities. SEGA enables formal testing of these specifications via Wald tests. The restrictions $H_0:\Sigma=I_J$, $H_0:\sigma_{jk}=0$ for all $j\neq k$, and $H_0:\Sigma=\frac{1}{2}(I_J+\iota_J\iota_J^\top)$ are all decisively rejected, with test statistics of $27{,}322$, $24{,}094$, and $9{,}451.6$, respectively, far exceeding their corresponding 1\% critical values of $68.709$, $58.619$, and $68.709$.
These results indicate that the unrestricted covariance structure captures substitution patterns that are strongly rejected by commonly imposed restricted specifications.
\subsection{Household heterogeneity in pasta demand}\label{sec:app_tobit}
We fit the random effects Tobit model to a data set constructed from the same consumer choice data used in Section~\ref{sec:app_mnp}. We select the purchases of private label pasta products that do not involve coupons from $n=203,965$ households across $T=62$ weeks. We model in each week $t$ the total dollar sales per ounce per household $i$, conditional on a set of promotion indicators, a price index, lagged sales, and weekly dummies. The 16 promotion variables indicate the number of shopping trips during the week with a private label pasta product in the weekly mailer or in a temporary in-store display, with different indicators for different locations in the mailer or display. The standardized log price index is constructed as the arithmetic mean across the price of all products in a particular week, where the weekly price for each product is approximated as the average price over all transactions for that product in that week. The lagged dollar sales for each household in each week is measured as an average across the past four weeks. The model includes a random intercept and random coefficients for the price index and lagged sales.
Figure~\ref{fig:tobit} shows the estimated probability of positive pasta sales in Panel (a) and the estimated marginal effect of the log price index on expected sales in Panel (b), both as a function of the household-specific lag-sales coefficient. Promotions are set to zero, and Panel (a) conditions on the average price index and zero lagged sales, and Panel (b) on the median price index across observations with positive sales and the median positive lagged sales. Panel (a) shows that households with stronger lag-sales effects have lower baseline purchase inclination. Panel (b) shows that the price marginal effect becomes more negative for households with stronger lag-sales effects, while the uncertainty band widens at high lag-sales coefficients. Together, these patterns suggest that households with low baseline demand are more state dependent and, when they have recently purchased, exhibit stronger price responsiveness.
\begin{figure}[tb!]
\caption{Household heterogeneity and state dependence in pasta demand}
\centering
\includegraphics*[width=\textwidth]{figures_application/tobit.eps}
\begin{flushleft}
\footnotesize
This figure shows the probability of positive pasta sales with average price index and zero lagged sales (Panel (a)) and the marginal effect of the log price index on expected sales with median price index and median positive lagged sales (Panel (b)), as a function of the household-specific lag-sales coefficient. Promotions are set to zero. The solid lines represent the point estimates and the dashed lines the 95\% confidence intervals.
\end{flushleft}
\label{fig:tobit}
\end{figure}
The random effects Tobit model allows for household heterogeneity in baseline purchase inclination, price sensitivity, and state dependence. Restricted specifications often assume that these sources of heterogeneity are independent, or that state dependence is homogeneous across households. SEGA enables formal testing of such restrictions via Wald tests. We first test whether the random effects are independent, $H_0:\Omega_{0p}=\Omega_{0\ell}=\Omega_{p\ell}=0$, where the subscripts denote the random intercept, price-index coefficient, and lag-sales coefficient. We also test the economically salient restriction that baseline purchase inclination and state dependence are unrelated, $H_0:\Omega_{0\ell}=0$. Both restrictions are decisively rejected, with both test statistics larger than 2900, exceeding their corresponding 1\% critical values of 11.345 and 6.635.
\section{Theoretical properties of SEGA}\label{sec:theory}
In this section, we show that SEGA produces inferences that are first-order asymptotically equivalent to the infeasible MLE when the number of iterations in the algorithm scale appropriately with $n$. Critically, unlike standard simulation-based estimators, such as SML and MSS, we do not require that the number of simulations used to approximate the score diverges: first-order asymptotic equivalence follows using a single simulated draw of the latent variables at each iteration of the algorithm.
\subsection{Finite-sample concentration}
In this section, we analyze the randomness in SEGA due to the simulated latent variables generated along the iteration path. To make the algorithms dependence on the simulated variables $z_{1:n}$ explicit, in what follows we write $\widehat{m}_n(\theta;z_{1:n})$.
The iterates in Algorithm~\ref{alg:sega} can be written as
\begin{align}
\theta^{(k+1)}_n&=\theta^{(k)}_n+\eta^{(k+1)}_nm_n(\theta^{(k)}_n)-\eta_n^{(k+1)}\left\{m_n(\theta^{(k)}_n)-\widehat{m}_n(\theta^{(k)}_n,\widetilde{z}^{(k)}_{1:n})\right\},
\label{eq:new}
\end{align}
where the notation $\widetilde{z}^{(k)}_{1:n}$ clarifies that the latent variables used to estimate the score are drawn at the $k$-th step of the algorithm, and are therefore conditional on $\theta_n^{(k)}$. The first term in \eqref{eq:new} is the infeasible gradient-ascent update based on the exact normalized score \(m_n(\theta)\). The second term can be viewed as simulation noise: under Fisher's identity, this noise is a martingale difference sequence (MDS) conditional on the past iterates for any $k\ge1$. Hence, the theoretical behavior of $\theta_n^{(k+1)}$ will be determined by how we can control the MDS term. The formulation in \eqref{eq:new} also shows that the randomness in $\theta_n^{(k+1)}$ is due to the conditional simulation of $\widetilde{z}_{1:n}^{(k)}$ - generated conditionally on $\theta_n^{(k)}$ and the observed data $d_{1:n}$ - clarifying that this randomness depends on the observed path of iterates $\{\theta_n^{(0)},\theta_n^{(1)},\dots,\theta_n^{(k)}\}$.
These observations imply that SEGA estimators have different theoretical behavior to that of other simulation-based estimators, such as MSS. In particular, the conditional simulation of $\widetilde{z}_{1:n}^{(k)}$ introduces feedback between the simulated datasets and the updated parameters values: a draw of $\widetilde{z}_{1:n}^{(k)}$ in the tail of the distribution $p_\theta(\widetilde{z}_{1:n}^{}\mid d_{1:n})$ could deliver an iterate $\theta_n^{(k+1)}$ that lies quite far away from the MLE. Hence, the analysis of any individual parameter iterate may not be that informative, as such draws of $\widetilde{z}_{1:n}^{}$ will occur in practice. Consequently, the most useful way to understand the behavior of the sequence $\theta_n^{(k+1)}$ is by analyzing its average behavior, where we average over $\widetilde{z}_{1:n}$ and $d_{1:n}$.
To investigate the average behavior of $\theta_n^{(k)}$, we maintain a set of regularity conditions for the conditional expectation of $\widehat{m}_n(\theta;\widetilde{z}_{1:n})$ that are similar to those encountered in the literature on SGA methods.
Let $\mathbb{E}$ denote the expectation with respect to the observed and simulated data. Let $(\mathcal F_k)_{k\ge0}$ denote the sigma-field generated by the observed data, the initial value $\theta_n^{(0)}$, and all simulated latent variables used up to iteration $k-1$. We maintain the following assumptions.
\begin{assumption}\label{ass:mds}
The following conditions hold with probability one (wp1).
\noindent (i) For each $k\ge1$:
$$
\mathbb{E}[\widehat{m}_n(\theta^{(k)}_n;z^{(k)}_{1:n})\mid \mathcal{F}_{k}]=m_n(\theta^{(k)}_n);
$$
\noindent(ii) There exists a $\nu\in\mathbb{R}_+$ such that, for each $k\ge1$:
$$
\mathbb{E}\left[\|\widehat{m}_n(\hat\theta_n;z^{(k)}_{1:n})\|_2^2\mid\mathcal{F}_{k}\right]\le\nu^2;
$$
\noindent(iii) For any $\mathcal{F}_{k}$ measurable $\theta_2,\theta_1\in\Theta$, there exists an $L>0$ such that
$$
\mathbb{E}[\|\widehat{m}_n(\theta_2;z^{(k)}_{1:n})-\widehat{m}_n(\theta_1;z^{(k)}_{1:n})\|^2_2\mid\mathcal{F}_{k}]\le L^2\|\theta_2-\theta_1\|^2_2.
$$
\end{assumption}
Assumption~\ref{ass:mds} collects the regularity conditions required for the simulated score. Part~(i) is the conditional unbiasedness implied by Fisher's identity. This is automatically satisfied by SEGA if $z^{(k)}_{1:n}$ is an exact draw from its conditional distribution. If draws are produced by Gibbs sampling, this condition assumes the Gibbs chain to be at stationarity. Part~(ii) requires the stochastic score to have a uniformly bounded conditional second moment at the MLE $\hat\theta_n$. Part~(iii) is a mean-square Lipschitz condition in the parameter, evaluated using the same simulated latent variables.
In order to provide non-asymptotic, in $n$ and $k$, convergence results for $\theta_n^{(k)}$ in Algorithm~\ref{alg:sega}, we require one additional assumption.
\begin{assumption}\label{ass:conv} The objective function $\ell_n(\theta)$ is continuously differentiable and $\mu$-strongly concave: for some $\mu>0$, and any $\theta_2,\theta_1\in\Theta$,
$$\ell_n\left(\theta_2\right) \leq \ell_n\left(\theta_1\right)+\left\langle\nabla_\theta \ell_n\left(\theta_1\right), \theta_2-\theta_1\right\rangle-\frac{n\mu}{2}\left\|\theta_1-\theta_2\right\|_2^2 .$$
\end{assumption}
Assumption~\ref{ass:conv} is a global strong-concavity condition on the sample log likelihood; this condition is equivalent to the requirement that the normalized objective $n^{-1}\ell_n(\theta)$ is $\mu$-strongly concave, so the normalized score $m_n(\theta)=n^{-1}\nabla_\theta\ell_n(\theta)$ is strongly monotone with curvature $\mu$. This condition is stronger than what is required for asymptotic likelihood theory, where local curvature around the MLE is sufficient. We impose it here to obtain a finite-$n$ and finite-$k$ concentration bound for the SEGA iterates and $\overline\theta_n$. An asymptotic version can be obtained under weaker local conditions.
While Assumptions~\ref{ass:mds}-\ref{ass:conv} may appear strict, {related conditions have been maintained by, e.g., \cite{moulines2011non} and \cite{chen2020statistical} when providing non-asymptotic convergence results for different SGA algorithms than SEGA.} Moreover, with effort, these conditions can be verified in the MNP example under low-level sufficient conditions, see Supplemental Appendix~\ref{A:lemmas} for details.
\begin{lemma}\label{lemma:MNP}
Under the regularity conditions in Supplemental Appendix~\ref{A:lemmas}, Assumptions~\ref{ass:mds}-\ref{ass:conv} are satisfied for the MNP model.
\end{lemma}
To obtain our key results, we require additional structure on the stepsize which we impose through the following assumption.
\begin{assumption}\label{ass:iters}
The step-size satisfies $\eta_n^{(k)}\propto\eta_n k^{-\alpha}$, with $\alpha\in(1/2,1)$, and $0<\underline{c}\le \eta_n\le \overline{c}$ (wp1), with $\eta_n$ being $\mathcal{F}_0$ measurable.
\end{assumption}
Assumption \ref{ass:iters} restricts how fast $\eta_n^{(k)}$ decreases to zero. The condition $\alpha\in(1/2,1)$ is the usual Robbins--Monro range, which ensures that $\sum_k \eta_n^{(k)}=\infty$ and $\sum_k \{\eta_n^{(k)}\}^2<\infty$.
Although the original formulation of ADADELTA given in Algorithm \ref{alg:sega} performs well in practice, its update need not satisfy Assumption \ref{ass:iters}. One may therefore use ADADELTA during the initial iterations and, once the smoothed iterates plateau, switch to the scaled update in Assumption \ref{ass:iters}.
We are now ready to state our first main result, which bounds the deviation between the MLE, $\hat\theta_n$, and the SEGA iterates. Proofs of all subsequent results are deferred to Appendix~\ref{A:proofs}.
\begin{theorem}\label{thm:main}
Under Assumptions \ref{ass:mds}-\ref{ass:iters},
$$
\mathbb{E}\|\theta^{(k)}_n-\hat\theta\|^p\lesssim k^{-\frac{p\alpha}{2}}(1+\mathbb{E}\|\theta^{(0)}_n-\hat\theta_n\|^2)^{p/2},\quad p\in\{1,2\}.
$$
\end{theorem}
The individual iterates $\theta^{(k)}_n$ are noisy due to the simulation of the latent variables $z_{1:n}$. Given this, we may obtain more precise inferences by instead considering the smoothed version of the iterates $\overline\theta_n$. It is fairly direct to show that $\overline\theta_n$ satisfies a version of Theorem \ref{thm:main}.
\begin{corollary}\label{corr:smoothed}
Under Assumptions \ref{ass:mds}-\ref{ass:iters},
$$\mathbb{E}\|\overline\theta_n-\hat\theta\|^p\lesssim k^{-\frac{p\alpha}{2}}(1+\mathbb{E}\|\theta^{(0)}_n-\hat\theta_n\|^2)^{p/2},\quad p\in\{1,2\}.
$$
\end{corollary}
\subsection{Asymptotic equivalence to the infeasible MLE}
We now study the asymptotic behavior of $\bar\theta_n$ as both the sample size \(n\) and the number of iterations \(k_n=k(n)\) diverge.
The result in Theorem \ref{thm:main} states that the average deviation between $\theta^{(k)}_n$ and $\hat\theta_n$ can be bounded by a constant multiple of $k^{-\alpha/2}$ and the difference between the starting value of the algorithm, $\theta^{(0)}_n$, and the MLE $\hat\theta_n$. So long as $\mathbb{E}\|\theta_n^{(0)}-\hat\theta_n\|\le C_n$, which includes the case where $\theta_n^{(0)}$ is an arbitrary non-random starting point restricted to lie in $\Theta$, with $\Theta$ compact, then we can immediately prove convergence in distribution of $\overline{\theta}^{}_n$. To state such a result, we make the relationship between $k$ and $n$ formal via the following assumption.
\begin{assumption}\label{ass:iters2}
The number of iterations, $k_n=k(n)$, satisfies $\sqrt{n} / k_n^{\alpha / 2}=o(1)$ as $n \rightarrow \infty$.
\end{assumption}
Assumption \ref{ass:iters2} ensures that the optimization error from SEGA is asymptotically negligible relative to the $n^{-1/2}$ sampling error of the infeasible MLE. For example, using an iteration schedule like $k(n)\asymp n^{1/\beta}$ requires choosing $\alpha>\beta$ and satisfies Assumption \ref{ass:iters2}.
To derive our ultimate result in as direct a manner as possible, we also maintain the following condition on the initial value. This condition can be weakened at the cost of imposing additional assumptions that are less interpretable.
\begin{assumption}\label{ass:starting}
There exists a random variable $C_n$ such that $\|\theta_n^{(0)}-\hat{\theta}_n\| \leq C_n$ where $\mathbb{E}(C_n^2)<\infty$.
\end{assumption}
The following result shows that the smoothed SEGA estimator $\overline\theta_n$ is first-order asymptotically equivalent to the infeasible MLE.
\begin{theorem}\label{thm:limit}
Assumptions \ref{ass:mds}-\ref{ass:starting} are satisfied. If $\sqrt{n}(\hat{\theta}_n-\theta_0)\Rightarrow N(0,V)$, then $\sqrt{n}(\overline\theta_n-\theta_0)\Rightarrow N(0,V)$.
\end{theorem}
Theorem \ref{thm:limit} is the main inferential result for SEGA: for a sufficiently large number of iterations, SEGA inherits the limiting distribution of the infeasible MLE, including the usual sandwich form under misspecification as defined in Section~\ref{sec:inference}.
\section{Numerical experiments}\label{sec:experiments}
We use numerical experiments to evaluate the finite-sample performance of SEGA in the multinomial probit model. We compare SEGA with the method of simulated scores (MSS), simulated maximum likelihood (SML), and the exact MLE. We examine whether the theoretical concentration and asymptotic equivalence results for SEGA translate into accurate and computationally scalable finite-sample inference. Experiments with the random effects Tobit model are reported in Supplemental Appendix~\ref{A:tobit_simulation}.
\subsection{Design and implementation}
We consider the MNP model described in Section~\ref{sec:examplesMNP}. The data is generated according to \eqref{eq:Y_i} and \eqref{eq:Z_i}. The predictor matrix $x_i$ has rows $x_{ij} = (e_j^\top, z_{ij})$, where $e_j$ denotes the $j$-th standard basis vector in $\mathbb{R}^J$ and $z_{ij} = \log p_{ij} - \log p_{i0}\sim N(0,1)$ denotes the log price difference of alternative $j$ relative to the base. Hence, the model contains $J$ alternative-specific intercepts and one price coefficient. Experiments with $J$ equal to 2, 3, and 6, set $\beta = [0.5, -0.5, -0.8]^\top$, $\beta = [0.5, -0.5, -0.2, -0.8]^\top$, and $\beta = [-0.5, -0.30,-0.1,0.1,0.3,0.5,-0.8]^\top$, respectively. All experiments assume $\Sigma = 0.5 I_J + 0.5 \iota_J \iota_J^\top$.
SEGA is implemented with 20,000 iterations with the final 1,000 iterations used for smoothing. We use $S=20,000$ draws for the variance estimator. SML is implemented using an off-the-shelf optimization routine with the analytical score of the simulated likelihood derived in \citet{bolduc1999practical}. Following \citet{gates2006mata}, the number of simulation draws is set to $(J+1)\times 50$. MSS is implemented using an off-the-shelf root finding algorithm where the score vector equation is evaluated using the Gibbs resampling simulator proposed in \cite{hajivassiliou1998method}. We vary the number of simulation draws \(R\in\{10,100\}\) and the number of Gibbs-resampling steps
\(G\in\{1,5,10,15,20\}\). When $J=2$, we compute the exact MLE by numerical integration, described in Supplemental Appendix~\ref{A:mnp_MLE}.
\subsection{Simulation noise and comparison with MSS}\label{sec:MSS_Emp}
In this section, we compare the simulation noise inherent in SEGA against that in MSS. As discussed in Section \ref{sec:MSS_discuss}, MSS can deliver unstable parameter estimates if the number of simulations $S$ is not chosen large enough, or if the Gibbs sampler used to estimate the scores does not draw from the stationary distribution. This latter behavior can occur even with small choice problems when the Gibbs sampler does not take enough steps, and, in general, as the dimension of the choice set increases we must run the Gibbs sampler longer to draw from the correct stationary distribution.
To isolate the simulation noise in the SEGA and MSS algorithms, we consider one fixed simulated dataset with $n = 10{,}000$ observations and repeatedly estimate the model using SEGA and MSS. Since the data are held fixed, variation across repetitions reflects algorithmic noise rather than sampling variation. We compare SEGA and MSS to the estimate obtained from SML, which provides an estimate close to the MLE in this setting.
Figure~\ref{fig:HMcompare_densities} shows the price coefficient estimates across 100 repetitions. The results are qualitatively similar across all parameters. The horizontal line denotes the SML estimate. The gold boxplots correspond to SEGA, while the grey boxplots correspond to MSS. Panels~(a) and~(b) report results for MSS with \(R=10\) simulation draws in the three-choice and four-choice models, respectively. Panels~(c) and~(d) repeat the comparison for MSS with \(R=100\) simulation draws.
\begin{figure}[tb!]
\caption{Simulation noise in SEGA and MSS estimates}
\centering
\includegraphics*[width=\textwidth,trim=1.5cm 1cm 1.5cm .8cm,clip]{figures_simulation/IntroFig.eps}
\begin{flushleft}
\footnotesize
Parameter estimate variation under MSS and SEGA. The first column reports estimates from the three-choice MNP model, and the second column reports estimates from the four-choice MNP model. The horizontal line denotes the SML estimate. Gold boxplots correspond to SEGA, while grey boxplots correspond to MSS. The first and second rows show MSS results using $R=10$ and $R=100$ simulations, respectively. All boxplots are constructed from 100 repetitions.
\end{flushleft}
\label{fig:HMcompare_densities}
\end{figure}
The figure shows that MSS is highly-sensitive to its tuning parameters. When \(R=10\), the MSS estimates are substantially dispersed across repetitions, indicating substantial Monte Carlo variability in the estimator of the score function used by MSS. Increasing the number of Gibbs-resampling steps reduces bias relative to the benchmark, but does not remove the dispersion. This pattern is especially visible in the four-choice model in Panel~(b), where insufficient resampling produces estimates that are far from the SML estimate. Increasing the number of simulation draws to \(R=100\) reduces the variability of MSS, as shown in Panels~(c) and~(d). However, the remaining bias still depends strongly on the number of Gibbs-resampling steps, especially in the higher-dimensional model.
In contrast, the SEGA estimates are tightly concentrated around the SML estimate in all panels. This stability arises from the sequential smoothing in SGA, which averages simulation noise progressively over iterations to reduce their variability. In contrast, MSS attempts to reduce simulation noise within each iteration, a strategy that substantially limits applicability in large-scale settings.
The experiment illustrates that the instability of MSS noted in the literature can be directly attributed to its sensitivity of its tuning parameters, together with its dependence on the initialization of the latent variables and model parameters. The design of this experiments is deliberately favorable to MSS: the method is initialized at the true parameter value, whereas SEGA is initialized at default values. When initialized at default values, MSS frequently encounters numerical failures in the root-finding algorithm; for example, with \(G=10\) and \(R=100\), approximately 30\% of runs fail. Given the implementation challenges associated with MSS, we do not include it in the following experiments.
\subsection{Sampling behavior and comparison with SML}
The second set of experiments studies sampling variability. We generate $1{,}000$ independent data sets and compare the sampling distributions of SEGA and SML, using the exact MLE as a benchmark when it is computationally available. Figure~\ref{fig:MCdens28} compares the repeated-sampling distributions of SEGA and SML for the price coefficient in the MNP model. Panels~(a) and~(c) report the low-dimensional case with \(J=2\), for which the exact MLE can also be computed by numerical integration. Panels~(b) and~(d) report the higher-dimensional case with \(J=6\), for which exact likelihood evaluation is computationally infeasible. The top row uses \(n=1{,}000\), while the bottom row uses \(n=20{,}000\). The vertical line indicates the true parameter value.
\begin{figure}[tb!]
\caption{Monte Carlo distributions coefficient estimates}
\centering
\includegraphics*[width=\textwidth]{figures_simulation/MCdens28.eps}
\begin{flushleft}
\footnotesize
This figure shows the Monte Carlo distributions of the estimates for $\beta_J$ by SEGA (solid line) and SML (dotted line). The panels correspond to Monte Carlo experiments with $J=2,6$ and $n=1000,20000$. The vertical lines indicate the parameter value in the data generating process.
\end{flushleft}
\label{fig:MCdens28}
\end{figure}
In the low-dimensional case, the SEGA, SML, and exact MLE distributions are virtually indistinguishable. This is visible in Panel~(a) and (c), where all three distributions are tightly concentrated around the true value. These results show that, when the exact MLE is available as a benchmark, SEGA reproduces its sampling behavior closely. The higher-dimensional case is more challenging. In Panel~(b), where \(J=6\) and \(n=1{,}000\), the SEGA and SML distributions are both more dispersed and display skewness. The two distributions are not identical in this small-sample, higher-dimensional setting, reflecting the greater difficulty of estimating the MNP model when the latent utility dimension increases. However, Panel~(d) shows that this difference largely disappears when the sample size increases to \(n=20{,}000\). The SEGA and SML distributions then become tightly concentrated and nearly coincide. This pattern is consistent with the asymptotic equivalence result in Theorem~\ref{thm:limit}: as the sample size grows and the number of SEGA iterations is sufficiently large, the sampling behavior of SEGA matches that of MLE.
Table~\ref{tab:MCstat28} summarizes the same comparison in terms of absolute bias, root mean squared error, and computation time. The bias and RMSE results confirm the message from Figure~\ref{fig:MCdens28}. In the \(J=2\) designs, SEGA, SML, and the exact MLE have nearly identical accuracy. In the \(J=6\) designs, SEGA and SML also deliver similar RMSEs, with only small differences in absolute bias. Thus, the stochastic score approximation used by SEGA does not lead to a meaningful loss of statistical accuracy in these experiments.
\begin{table}[tb!]
\centering
\caption{Monte Carlo bias, RMSE, and computation time}
\begin{threeparttable}
\begin{tabular}{llrrrrrr}
\toprule \toprule
& & \multicolumn{3}{c}{$J=2$} & \multicolumn{3}{c}{$J=6$} \\
\cmidrule(lr){3-5} \cmidrule(lr){6-8}
\multicolumn{1}{l}{$n$} & & Absolute bias & RMSE & Time & Absolute bias & RMSE & \multicolumn{1}{l}{Time} \\
\midrule
1000 & SEGA & 0.006 & 0.062 & 0.085 & 0.046 & 0.191 & 0.424 \\
& SML & 0.003 & 0.062 & 0.028 & 0.008 & 0.205 & 2.175\\
& MLE & 0.002 & 0.062 & 0.275 & -& -& -\\
\cmidrule(lr){2-8}
10000 & SEGA & 0.001 & 0.014 & 0.552 & 0.012 & 0.047 & 3.0580 \\
& SML & 0.000 & 0.014 & 0.261 & 0.010 & 0.045 & 22.8210\\
& MLE & 0.000 & 0.014 & 7.485 & - & - & - \\
\bottomrule \bottomrule
\end{tabular}
\begin{tablenotes}
\footnotesize
\item This table shows the absolute bias and root mean squared error averaged across all parameters in the multinomial probit model. Furthermore, it reports the average computation time in minutes for each method across the replications in each experiment.
\end{tablenotes}
\end{threeparttable}
\label{tab:MCstat28}
\end{table}
The main difference between the methods is computational. When \(J=2\), SML is faster
than SEGA because the simulated likelihood is relatively cheap to evaluate. However, the
computational cost of SML grows rapidly with the number of alternatives. For \(J=6\) and
\(n=1{,}000\), SML is already more than five times slower than SEGA. For \(J=6\) and
\(n=20{,}000\), SML is more than seven times slower. This difference reflects the distinct way in which the methods handle integration over the latent utilities. SML approximates choice probabilities at each likelihood evaluation, and this becomes increasingly costly as the dimension of the latent utility vector grows. SEGA instead uses a single draw of the latent variables at each iteration and averages the resulting simulation noise over the stochastic-gradient path.
\subsection{Inferential accuracy}
The final experiment evaluates whether the standard errors implied by the asymptotic distribution in Theorem~\ref{thm:limit} deliver well-calibrated confidence intervals in finite samples. We compute coverage rates across 1,000 Monte Carlo replications for experiments with \(J=2\) and a correctly specified model, \(J=2\) and a misspecified model, and a correctly specified but higher dimensional setting with \(J=6\).
The misspecified \(J=2\) design generates from the same MNP model but the covariance matrix is made heteroskedastic: \(\Sigma_i = \left(0.1+0.5{z_{i1}^2}\right)^2 \left(0.5 I_J + 0.5 \iota_J\iota_J^\top\right)\). The estimated model nevertheless imposes the homoskedastic MNP covariance specification in \eqref{eq:Z_i}. Hence, the target parameter is the pseudo-true value \( \theta_0=\arg\max_{\theta\in\Theta}E\{\ell_i(\theta)\}\), rather than the parameter value used to generate the data. We approximate this pseudo-true parameter by computing the exact MLE on 100 independent samples of size 100,000 and averaging the resulting estimates.
Table~\ref{tab:coverage} reports the coverage rates for the experiments with $J=2$. For the correctly specified \(J=2\) model, coverage is close to the nominal 95\% level for all parameters. At \(n=1{,}000\), the coverage rates range from 0.91 to 0.95, and they improve further at \(n=10{,}000\), where all rates lie between 0.94 and 0.95.
The misspecified \(J=2\) design delivers similarly accurate coverage: the rates are close to 0.95 at both sample sizes, indicating that sandwich-based SEGA inference remains well calibrated when the target is interpreted as the pseudo-true parameter.
\begin{table}[tb!]
\centering
\caption{Empirical coverage of 95\% SEGA confidence intervals}
\label{tab:coverage}
\begin{threeparttable}
\begin{tabular}{lccccc}
\toprule\toprule
Sample size & \(\beta_1\) & \(\beta_2\) & \(\beta_p\) & \(\Sigma_{12}\) & \(\Sigma_{22}\) \\
\hline
&\multicolumn{5}{c}{Correct model specification} \\
\cline{2-6}
\(1{,}000\) & 0.94 & 0.93 & 0.95 & 0.93 & 0.91 \\
\(10{,}000\) & 0.95 & 0.94 & 0.95 & 0.95 & 0.94 \\
\hline
&\multicolumn{5}{c}{Model misspecification} \\
\cline{2-6}
\(1{,}000\) & 0.96 & 0.94 & 0.95 & 0.95 & 0.96 \\
\(10{,}000\) & 0.95 & 0.94 & 0.95 & 0.93 & 0.94 \\
\bottomrule \bottomrule
\end{tabular}
\begin{tablenotes}
\footnotesize
\item This table reports empirical coverage rates of nominal 95\% confidence intervals for the SEGA estimator across 1,000 Monte Carlo replications. The first reports results under correct specification and the second panel under misspecification.
\end{tablenotes}
\end{threeparttable}
\end{table}
In the high-dimensional setting, the Monte Carlo experiment indicates that larger sample sizes are required for accurate empirical coverage. For the \(J=6\) design, the model contains 27 parameters, and the individual coverages range between 0.83 and 0.99 with \(n=1{,}000\), and 0.89 and 0.96 with \(n=10{,}000\), while the median coverage is 0.94 with both sample sizes. As the number of alternatives increases, the likelihood involves higher-dimensional integrals, and accurate estimation of the sampling variability requires substantially more data. Our empirical application with the multinomial probit model contains more than one million observations, so the large-sample regime is the relevant one for that setting. However, for this highly-nonlinear model, Monte Carlo validation at that scale is computationally infeasible.
\section{Conclusion}\label{sec:conclusion}
Latent variable models are central to empirical work in economics, but likelihood-based inference in these models is often limited by the need to integrate over high-dimensional latent variables.
This paper proposes SEGA, a stochastic estimation framework that uses Fisher’s identity to construct an unbiased estimate of the likelihood score from a single draw of the latent variables. By embedding this score in a stochastic gradient ascent algorithm, SEGA averages simulation noise over the optimization path rather than requiring accurate likelihood or score approximations at each parameter value.
The theoretical results show that SEGA has the same first-order inferential properties as the infeasible MLE. In finite samples, the SEGA estimator concentrates around the infeasible MLE as the number of iterations increases. Asymptotically, when the number of iterations grows appropriately with the sample size, the SEGA estimator is asymptotically equivalent to the infeasible MLE and inherits its first-order limiting distribution, including the usual sandwich form under misspecification. The numerical experiments support these theoretical results.
The empirical applications illustrate the value of SEGA for applied researchers. In the multinomial probit application, SEGA allows us to estimate an unrestricted model of brand choice using more than one million purchase observations, construct confidence intervals for own- and cross-price elasticities, and test covariance restrictions that are commonly imposed for tractability. In the random effects Tobit application, SEGA allows us to estimate a household-level demand model with more than 200,000 households observed over 62 weeks, and to conduct inference on heterogeneity in baseline demand, price sensitivity, and state dependence.
These applications show that SEGA makes likelihood-based inference feasible in LDVMs at scales where existing methods are impractical. By removing this computational barrier, SEGA expands the set of empirical questions that applied researchers can study with flexible limited dependent variable models.
\newpage