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.
84,428 characters
SAUSS: Stochastic Approximation with Unbiased Simulated Scores for Limited Dependent Variable Models
\maketitle
\begin{abstract}
Multinomial choice models allow flexible substitution patterns but become computationally demanding with many alternatives or observations. With a fixed per-observation simulation budget, simulated maximum likelihood introduces simulation bias, while each optimization step requires a full-sample likelihood evaluation. We propose Stochastic Approximation with Unbiased Simulated Scores (SAUSS), an averaged stochastic approximation based on conditionally unbiased mini-batch score estimates. Each iteration uses a fixed mini-batch regardless of sample size. For multinomial probit, accept--reject sampling provides exact conditional draws and unbiased score estimates for any fixed number of accepted draws. Under local conditions, asymptotic theory for the averaged estimator and the partial-sum process of the SAUSS iterates incorporates mini-batch and simulation variability and supports random-scaling and plug-in inference. In simulations and an application, SAUSS gives comparable results in less than 1\% of the computation time of simulated maximum likelihood. SAUSS extends to limited dependent variable models with conditional-expectation score representations and exact conditional sampling.
\end{abstract}
\noindent\textbf{Keywords:}
Accept--reject simulation; limited dependent variable models; multinomial
probit; simulated scores; stochastic gradient descent.
\onehalfspacing
\section{Introduction}\label{sec:introduction}
{Statistical modeling of choices among multiple alternatives across a variety
of fields requires balancing flexible dependence structures against
computational tractability.} The multinomial logit
model remains a workhorse in part because its choice probabilities are available
in closed form \citep{mcfadden1974conditional,Train2009}. Its
independence-of-irrelevant-alternatives structure, however, restricts how
substitution across alternatives can respond to unobserved heterogeneity.
Multinomial probit (MNP) allows the latent utility shocks to have a flexible
covariance structure and thereby accommodates richer substitution patterns. The
price of this flexibility is computation: each likelihood contribution is a
multivariate normal region probability, with dimension increasing in the number
of alternatives. This integration problem has long limited routine likelihood
estimation of richly parameterized MNP models
\citep[e.g.,][]{gewekeKeaneRunkle1994alternative,Train2009}.
The computational burden has two distinct dimensions. As the number of
alternatives \(J\) grows, {each observation requires a
higher-dimensional region probability and, after scale normalization, a
flexible covariance specification contains a growing number of parameters.} As
the sample size \(n\) grows, conventional likelihood optimization repeatedly
evaluates these costly contributions for the full sample. {Conventional
full-sample simulated maximum likelihood mitigates the first problem by
replacing each probability with a smooth simulator, most commonly one based on the}
Geweke--Hajivassiliou--Keane (GHK) construction
\citep[e.g.,][]{borschSupanHajivassiliou1993,hajivassiliou1996simulation}.
{It does not by itself address the second problem created by a
large sample.} Moreover, placing an estimated probability inside a
logarithm creates finite-simulation bias even when the probability simulator
itself is unbiased. {Under conventional simulated-likelihood
asymptotics, the number of draws per observation must therefore increase with
\(n\), while every optimization step continues to process all \(n\)
observations.} Simulation
effort per contribution to the log likelihood and the number of contributions per step can thus grow
together.
{Classical simulation-based work addresses this integration
problem in several ways. \citet{mcfadden1989method} introduces the method of
simulated moments for discrete-response models, avoiding direct numerical
integration of choice probabilities. \citet{keane1994computational} develops a
computationally practical extension for panel limited dependent variable (LDV)
models by factorizing moment conditions in terms of simulated transition
probabilities. \citet{gewekeKeaneRunkle1994alternative} compare simulated
maximum likelihood, simulated moments, and Bayesian data augmentation directly
for MNP. Together with the development and evaluation of GHK and related
multivariate-normal simulators
\citep{borschSupanHajivassiliou1993,hajivassiliou1996simulation}, these
contributions make simulation-based estimation of flexible probit models
practical at moderate sample sizes but leave a basic large-sample tension:
accurate simulation and repeated full-sample computation remain costly as
\(n\) grows.}
{The closest antecedent to our approach is the method of simulated
scores (MSS) of \citet{hajivassiliou1998method}. Rather than inserting an
estimated probability into a criterion, MSS represents the likelihood score as
a conditional expectation of a complete-data score. Exact conditional draws
obtained by accept--reject simulation yield an unbiased simulated score for any
fixed number of accepted draws. MSS also develops continuous score simulators
based on recursive conditioning and Gibbs resampling. The unbiased
accept--reject score simulator is generally
discontinuous in the parameter, whereas the continuous constructions may
require increasing simulation effort to eliminate approximation bias.
{A finite Gibbs chain provides an exact conditional draw if
initialized in stationarity; from a nonstationary initialization, its
distribution converges to the target as the chain length tends to infinity.}
{Because MSS is formulated as a full-sample simulated-score
estimating equation, conventional numerical solution methods face a tradeoff
between unbiased score simulation and smoothness.}}
{This paper develops SAUSS---Stochastic Approximation with Unbiased
Simulated Scores---to sidestep that tradeoff. {SAUSS uses exact
conditional accept--reject draws to construct an unbiased, potentially
discontinuous score estimate inside a stochastic approximation recursion.}
This supplies a
conditionally unbiased update direction with suitable moment properties
without requiring each realization of the simulator to be smooth in the
parameter. At each iteration, SAUSS draws a mini-batch of \(m\) observations,
where \(m\) can be much smaller than \(n\) and may equal one, generates fresh
conditional score draws for those observations, and updates the parameter using
the resulting average negative score. The proposed estimator averages the
iterates in the manner of \citet{polyak1992acceleration}. With \(m\) fixed, the
number of observations processed per update does not grow with \(n\). SAUSS
thus retains the exact likelihood score as its target while avoiding both the
logarithm of a simulated probability and repeated full-sample updates.}
{Relative to existing methods, SAUSS makes three contributions.
First, it embeds the exact conditional simulation and unbiased score
construction of MSS in mini-batch stochastic approximation with iterate averaging.
{Second, it gives an explicit MNP implementation with a
scale-normalized, otherwise unrestricted covariance parameterization and
evaluates its computational performance as the number of alternatives grows.}
Third, it develops inference that jointly accounts for mini-batch and
score-simulation variation. The framework extends to other LDV models admitting
conditional-expectation score representations and exact conditional
simulation.}
{Under the reported configurations, SAUSS delivers results broadly
comparable to GHK-based simulated maximum likelihood, while GHK requires more
than 100 times as much computation time in both the Monte Carlo benchmark and
the empirical application. With the same tuning constants, SAUSS remains computationally feasible
as the choice set expands far beyond the benchmark design.
{Although the reported runtime ratios depend on the simulation
budgets and implementations, their consistent pattern illustrates the
practical value of the SAUSS architecture.}}
{The statistical theory treats the variation created by that
architecture as part of the inferential problem. Under local regularity and
initialization conditions, we establish a linear representation and functional
central limit theorem for the partial-sum path in the following
\(\epsilon\)-sense: for any \(\epsilon>0\), the learning-rate constant can be
chosen so that the recursion remains in the local neighborhood with probability
at least \(1-\epsilon-o(1)\), and the conclusions hold on that event. The
endpoint yields, in the same \(\epsilon\)-sense, the corresponding normal
approximation for the averaged estimator. The limit separates variation
generated by mini-batch sampling from variation generated by score simulation
and supports random-scaling and plug-in procedures that account for both
sources of algorithmic variation.}
{The closest contemporaneous work is
\citet{frazierLoaizaMayaNibbering2026scalable}, who propose Stochastically
Estimated Gradient Ascent (SEGA). Both SEGA and SAUSS combine conditional-score
simulation and iterate averaging. Drawing on two
simulation routes already developed in MSS, SEGA uses Gibbs-based score
simulation in \emph{full-sample} updates, whereas SAUSS uses score simulation
based on exact accept--reject conditional draws in \emph{mini-batch} updates.
SEGA's theory assumes conditional unbiasedness, which is satisfied by exact
conditional draws or by a Gibbs chain initialized in stationarity. The SAUSS
inferential limit retains variation from mini-batch sampling and score
simulation. Under SEGA's
regularity conditions and an iteration schedule that makes optimization error
negligible relative to \(n^{-1/2}\) sampling error, its averaged estimator is
first-order asymptotically equivalent to the infeasible maximum likelihood estimator (MLE) based on the exact
likelihood and therefore inherits its limiting distribution. The two approaches
thus use the same conditional-score lineage from MSS but adopt different computational
regimes and asymptotic approximations.}
{A parallel Bayesian literature avoids direct likelihood evaluation
through latent-data augmentation. \citet{albertChib1993bayesian} provide the
foundational framework for binary and polychotomous probit models. For MNP,
\citet{mccullochRossi1994exact} use latent Gaussian utilities in a Gibbs
posterior simulator, and \citet{imaiVanDyk2005bayesian} develop marginal data
augmentation. Recent work develops scalable Bayesian estimation using
structured covariance models \citep{loaizaMayaNibbering2022scalable} and fast
variational approximations \citep{loaizaMayaNibbering2023fast}. These methods
target posterior inference, whereas SAUSS targets the exact likelihood score
and delivers frequentist inference for the averaged stochastic approximation.}
{Other recent approaches modify the likelihood approximation or
the computation of choice probabilities, including approximate composite
marginal likelihood \citep{bhat2011macml}, expectation propagation within an EM
algorithm \citep{dingImbensQuYe2024probit}, and a neural-network emulator for
correlated discrete-choice probabilities \citep{huchKeane2026amortized}.
SAUSS instead retains the exact likelihood score as its target and changes the
computational and inferential architecture in which unbiased simulated-score
contributions are subsampled across observations, accumulated over iterations,
and used for inference.}
{Finally, SAUSS connects the classical stochastic approximation
recursion of \citet{robbinsMonro1951stochastic} and iterate averaging of
\citet{polyak1992acceleration} to simulation-based likelihood estimation.
Recent stochastic-approximation work provides scalable estimation and inference
for smooth-loss models, with applications to linear mean and logistic regression
\citep{leeLiaoSeoShin2022randomScaling}; for large-scale linear quantile
regression \citep{leeLiaoSeoShin2025quantile}; for semiparametric models
\citep{chenTamerYao2026fast}; and for overidentified moment condition models through
SGMM \citep{chenLeeLiaoSeoShinSong2025sgmm} and SLIM
\citep{chenKimLeeSeoSong2025slim}. In SAUSS, however, each update contains
additional conditional-simulation variation because the exact likelihood score
is not directly available. The functional limit and variance procedures
explicitly accommodate this additional source of variation without smoothing
the accept--reject mechanism.}
The rest of the paper is organized as follows.
Section~\ref{sec:sauss-framework} presents the SAUSS framework, and
Sections~\ref{sec:ldv-identity} and~\ref{sec:mnp-sauss} develop the LDV score
identity and its MNP implementation. Sections~\ref{sec:asymptotic-theory}
and~\ref{sec:inference} establish the asymptotic theory and inference procedures.
Sections~\ref{sec:monte-carlo} and~\ref{sec:application} report the Monte Carlo
evidence and empirical application, and
Section~\ref{sec:conclusion} concludes.
Proofs of the asymptotic results and the supporting lemmas are collected in
Appendix~\ref{app:proof-roadmap}.
\section{The SAUSS Framework}\label{sec:sauss-framework}
{We propose estimation and inference using stochastic gradient descent, which recursively updates the estimator as follows:}
\[
\theta_t= \theta_{t-1} - \gamma_t \hat G_t(\theta_{t-1})
\]
with a suitably chosen learning rate $\gamma_t$. Here, {$\hat G_t(\theta_{t-1})$ is the negative of a simulated likelihood score.} Unlike the usual stochastic gradient descent inference framework, directly evaluating the score from the log-likelihood is a challenging task in limited dependent variable models.
This section presents the conceptual framework for SAUSS. It first
reviews why smooth simulated likelihoods are attractive but biased at a
fixed simulation budget, then explains how MSS supplies unbiased
simulated score contributions, and finally sets out the resulting
stochastic approximation recursion in descent form.
\subsection{Smooth Simulated Likelihood and Bias}
Let \(P_i(\theta)\) denote the probability of the observed outcome for
observation \(i\), and let
\[
\mathcal{L}_n(\theta)
=
\sum_{i=1}^n \log P_i(\theta)
\]
be the exact log-likelihood. Simulated maximum likelihood replaces
\(P_i(\theta)\) with a simulated probability estimator
\(\hat P_i(\theta)\) and maximizes
\[
\hat{\mathcal{L}}_n(\theta)
=
\sum_{i=1}^n \log \hat P_i(\theta).
\]
{GHK-based simulated maximum likelihood is a prominent
implementation of this strategy for MNP \citep[e.g.,][]{bolduc1999practical}.}
Much of the classical simulation literature therefore places heavy
weight on simulators that are continuous, differentiable, bounded away
from zero, and suitable for full-sample numerical optimization
\citep[e.g.,][]{Train2009}. Smoothness is valuable because the simulated
object is treated as an objective function to be maximized.
The difficulty is that the smooth simulated objective is generally not
an unbiased version of the exact log-likelihood. Even when
\(\hat P_i(\theta)\) is unbiased for \(P_i(\theta)\), the nonlinear
transformation \(\log \hat P_i(\theta)\) is generally biased. A related
problem arises for the simulated score. If
\(\hat P_i(\theta;\mathbf{u})\) denotes a simulated probability computed
from random draws \(\mathbf{u}\), then, for a fixed finite simulation
budget,
\begin{equation}\label{eq:sml-score-bias}
\mathbb{E}_{\mathbf{u}}
\left[
\nabla_\theta \log \hat P_i(\theta;\mathbf{u})
\right]
\neq
\nabla_\theta \log P_i(\theta)
\qquad
\text{in general},
\end{equation}
whenever the derivative and expectation are well defined. Here
\(\mathbb{E}_{\mathbf{u}}\) denotes expectation over the simulation
draws conditional on the observation.
In full-sample simulated maximum likelihood, this approximation error
is controlled by increasing the simulation budget as the sample size
grows. That strategy is less attractive for stochastic approximation
because a growing simulation budget also raises the cost of every
mini-batch evaluation. The alternative pursued here is to simulate an
unbiased likelihood-score contribution directly.
The number of simulated score draws used in each evaluation can then
remain fixed. In accept--reject implementations, the number of proposals
required to obtain those draws may nevertheless be random.
\subsection{Unbiased Simulated Scores}
The method of simulated scores (MSS) proposed by
\citet{hajivassiliou1998method} changes the object being simulated.
Instead of simulating \(P_i(\theta)\) and then taking a logarithm, it
simulates the score
\[
s_i(\theta)
=
\nabla_\theta \log P_i(\theta)
\]
directly.
When the simulator supplies exact conditional draws, or otherwise
provides an unbiased score construction, the resulting estimate can
satisfy
\[
\mathbb{E}_{\mathbf{u}}
\{\hat s_i(\theta)\mid D_i\}
=
s_i(\theta).
\]
In \citet{hajivassiliou1998method}, the MSS estimator is formulated as a
full-sample Z-estimator based on simulated scores, that is, as a solution
to
\begin{align}\label{def:mss-z}
n^{-1} \sum_{i=1}^n \hat s_i(\theta) = 0.
\end{align}
This root-finding problem can be difficult when
\(\theta\mapsto\hat s_i(\theta)\) is discontinuous, as often occurs with
accept--reject constructions.
SAUSS retains the unbiased simulated score but changes how it is used.
At each iteration, the negative simulated score provides a conditionally
unbiased stochastic direction for the gradient of the negative
log-likelihood. This direction need not itself be the derivative of a
smooth simulated objective.
\subsection{SAUSS as Stochastic Approximation}\label{sec:sgd}
We write the recursion in descent form for the negative log-likelihood.
Given the observed sample \(D_{1:n}\), define the individual loss
\(\ell_i(\theta)=-\log P_i(\theta)\) and the empirical criterion
\(Q_n(\theta)=n^{-1}\sum_{i=1}^n\ell_i(\theta)\). Let
\(\hat\theta_n\) denote a solution of the sample score equation
\(n^{-1}\sum_{i=1}^n s_i(\hat\theta_n)=0\).
Let \(\psi_i(\theta,\mathbf{u})\) be a one-draw simulated score satisfying
\(\mathbb{E}_{\mathbf{u}}\{\psi_i(\theta,\mathbf{u})\mid D_i\}
=s_i(\theta)\), and write \(g_i(\theta,\mathbf{u})
=-\psi_i(\theta,\mathbf{u})\). Thus \(g_i\) is a conditionally unbiased
stochastic direction for the individual loss gradient
\(\nabla_\theta\ell_i(\theta)=-s_i(\theta)\).
Set \(\theta_1=\theta_{\mathrm{init}}\). For iterations
\(t=2,\ldots,T\), draw
ordered mini-batch indices \(I_{t1},\ldots,I_{tm}\) independently and
uniformly from \(\{1,\ldots,n\}\), and generate the random
inputs \(\mathbf{u}_{tbr}\), \(r=1,\ldots,R\), independently for
every batch position, score draw, and iteration. The mini-batch
direction and SAUSS update are
\begin{equation}\label{eq:sa-update}
\begin{aligned}
\hat G_t(\theta_{t-1})
&=
\frac{1}{mR}
\sum_{b=1}^{m}\sum_{r=1}^{R}
g_{I_{tb}}(\theta_{t-1},\mathbf{u}_{tbr}), \\
\theta_t
&=
\theta_{t-1}
-\gamma_t\hat G_t(\theta_{t-1}).
\end{aligned}
\end{equation}
Here \(\gamma_t\) is a decreasing step size. If \(\mathcal F_{t-1}\)
denotes the information available
immediately before the draws at iteration \(t\), then
\(\mathbb{E}\{\hat G_t(\theta_{t-1})\mid
D_{1:n},\mathcal F_{t-1}\}
=\nabla_\theta Q_n(\theta_{t-1})\). Conditional on the sample, the ideal
recursion therefore targets \(\hat\theta_n\).
The reported estimator is the Polyak--Ruppert average
\begin{equation}\label{eq:pr-average}
\bar\theta_T
=
\frac{1}{T-T_0}
\sum_{t=T_0+1}^{T}
\theta_t,
\end{equation}
where \(T_0\) is the averaging start index. Thus \(T\) is the final
iterate index, the run contains \(T\) stored iterates including the
initial estimator, and the recursion performs \(T-1\) stochastic
updates.
The key statistical implication is that \(R\) need not diverge in
order to eliminate simulation bias from the ideal recursion: for every
fixed \(R\geq1\), the stochastic direction remains conditionally
unbiased. Holding \(R\) fixed does not, however, eliminate simulation
variance. The magnitude of the resulting algorithmic uncertainty
depends jointly on \(T\), \(m\), and \(R\).
\begin{remark}[Accepted draws and computational cost]
Fixing \(R\) fixes the number of accepted score draws, not the number
of proposals. For an accept--reject implementation with acceptance
probability \(a_i(\theta)\), the expected number of proposals required
to obtain \(R\) accepted draws is \(R/a_i(\theta)\). Thus the
exact recursion has random computational cost even when \(R\) is
fixed.
\end{remark}
The next section states the LDV score identity that delivers the
unbiased score contributions.
\section{The LDV Score Identity}\label{sec:ldv-identity}
This section states the identity used in
\citet{hajivassiliou1998method}. The key point is to represent the
observed outcome as a fixed latent-region event, so that differentiating
the log probability produces a conditional expectation of a complete-data
score.
Write the observed data as \(D_i=(Y_i,X_i)\), where \(Y_i\) is the
limited dependent variable outcome and \(X_i\) collects the observed
conditioning variables. Conditional on \(X_i\), suppose that \(Y_i\) is
generated by an underlying latent vector \(Y_i^*\). For the realized
outcome \(Y_i=y_i\), choose latent coordinates in which
\(\{Y_i=y_i\}\) is equivalent to the fixed-region event
\(Y_i^*\in\mathcal{D}_i\), where
\(\mathcal{D}_i=\mathcal{D}(y_i,X_i)\). The region may depend on the
observed outcome and conditioning variables, but it does not vary with
\(\theta\). All probabilities and densities below are conditional on
\(X_i\), with this conditioning suppressed to simplify notation. The
parameter enters through the conditional latent density
\(f_i(y;\theta)\), and the likelihood contribution is
\begin{equation}\label{eq:ldv-prob}
P_i(\theta)
=
\int_{\mathcal{D}_i}
f_i(y;\theta)
\,dy.
\end{equation}
\begin{lemma}[Hajivassiliou--McFadden identity]\label{lem:hm-identity}
Fix \(i\), condition on \(X_i\), and fix \(\theta\) in the interior of
the parameter space.
Assume that the observed outcome \(Y_i=y_i\) is equivalent to
\(Y_i^*\in\mathcal{D}_i\) for a measurable set \(\mathcal{D}_i\) that
does not depend on \(\theta\). Assume also that \(P_i(\theta)>0\) and
that \(f_i(y;\theta)>0\) for almost every \(y\in\mathcal{D}_i\).
Finally, suppose that there is a neighborhood \(\mathcal{N}_\theta\) of
\(\theta\) such that, for almost every \(y\in\mathcal{D}_i\), the map
\(\vartheta\mapsto f_i(y;\vartheta)\) is differentiable on
\(\mathcal{N}_\theta\), and there exists an integrable envelope \(M_i\)
such that
\[
\sup_{\vartheta\in\mathcal{N}_\theta}
\left\|
\nabla_\vartheta f_i(y;\vartheta)
\right\|
\leq
M_i(y),
\qquad
\int_{\mathcal{D}_i} M_i(y)\,dy<\infty .
\]
Then
\begin{equation}\label{eq:ldv-score-identity}
\nabla_\theta \log P_i(\theta)
=
\mathbb{E}_\theta
\left[
\nabla_\theta \log f_i(Y_i^*;\theta)
\mid
Y_i^* \in \mathcal{D}_i
\right].
\end{equation}
\end{lemma}
\begin{proof}
By the stated assumptions,
\begin{align}
\nabla_\theta P_i(\theta)
&=
\int_{\mathcal{D}_i}
\nabla_\theta f_i(y;\theta)
\,dy \notag \\
&=
\int_{\mathcal{D}_i}
\left[
\nabla_\theta \log f_i(y;\theta)
\right]
f_i(y;\theta)
\,dy.
\label{eq:ldv-score-derivation}
\end{align}
Dividing \eqref{eq:ldv-score-derivation} by \(P_i(\theta)\) gives
\[
\nabla_\theta \log P_i(\theta)
=
\int_{\mathcal{D}_i}
\nabla_\theta \log f_i(y;\theta)
\frac{f_i(y;\theta)}{P_i(\theta)}
\,dy,
\]
which is \eqref{eq:ldv-score-identity} because \(f_i(y;\theta)/P_i(\theta)\) is the conditional density of \(Y_i^*\) given \(Y_i^*\in\mathcal{D}_i\).
\end{proof}
\begin{remark}[Binary probit and fixed regions]
For binary probit, let
\(Y_i^*=X_i^\top\beta+\varepsilon_i\), with
\(\varepsilon_i\sim\mathcal{N}(0,1)\), and let
\(Y_i=\mathbb{I}\{Y_i^*>0\}\). In latent-utility coordinates, the
observed outcome corresponds to the fixed region \((0,\infty)\) when
\(Y_i=1\) and to \((-\infty,0]\) when \(Y_i=0\). In error coordinates,
the same event has a boundary that depends on \(\beta\). The identity is
therefore applied in latent-utility coordinates, where the region is
fixed and the parameter enters through the density.
\end{remark}
Lemma~\ref{lem:hm-identity} identifies the simulation target. If we can
draw independent exact realizations \(Y_{ir}^*(\theta)\) from the
conditional distribution of \(Y_i^*\) given
\(Y_i^* \in \mathcal{D}_i\), then
\begin{equation}\label{eq:ldv-unbiased-score}
\hat s_i(\theta)
=
\frac{1}{R}
\sum_{r=1}^{R}
\left.
\nabla_\theta \log f_i(y;\theta)
\right|_{y=Y_{ir}^*(\theta)}
\end{equation}
is unbiased for \(s_i(\theta)\) for any fixed \(R \geq 1\).
For an accept--reject simulator, when the same underlying proposal
sequence is used across parameter values, the simulated-score map may be
discontinuous in \(\theta\) because the identity of the accepted
proposal can change discretely as the parameter moves.
For full-sample numerical root finding, this discontinuity can create a
difficult computational problem. For stochastic approximation,
conditional unbiasedness is the central property, together with the
moment and stability conditions stated later.
The identity applies directly to the region-probability contributions in
canonical LDV settings emphasized by
\citet{hajivassiliou1998method}, including panel probit and multinomial
choice models with correlated latent errors. In Tobit-type models, it
applies to censored contributions, while uncensored density
contributions can be scored directly.
We build on the Hajivassiliou--McFadden identity and focus on LDV models
in which the score admits a conditional-expectation representation and
for which either an exact unbiased score simulator is available or a
practical simulator is available whose approximation error can be
diagnosed. Section~\ref{sec:mnp-sauss} verifies this construction
explicitly for MNP under the covariance normalization used by the
algorithm.
\section{SAUSS for Multinomial Probit}\label{sec:mnp-sauss}
This section specializes the LDV score identity to MNP and
describes how SAUSS can be implemented for MNP.
\subsection{Model and Choice Regions}
Consider \(n\) agents choosing among \(J\) alternatives. The latent
utility for agent \(i\) and alternative \(j\) is
\begin{equation}\label{eq:mnp-utility}
U_{ij} = x_{ij}^\top \beta + \varepsilon_{ij},
\quad j=1,\ldots,J,
\end{equation}
where \(x_{ij} \in \mathbb{R}^K\), \(\beta \in \mathbb{R}^K\), and
\(\varepsilon_i:=(\varepsilon_{i1},\ldots,\varepsilon_{iJ})^\top\)
satisfies \(\varepsilon_i\sim\mathcal{N}(0,\Sigma)\). Agent \(i\)
chooses \(y_i=j\) if \(U_{ij}>U_{ik}\) for all \(k \neq j\).
Utility levels are not identified, so we work in differences relative
to alternative 1. For \(j=2,\ldots,J\), define
\begin{equation}\label{eq:mnp-diff}
\Delta U_{ij}
:=
U_{ij}-U_{i1}
=
(x_{ij}-x_{i1})^\top\beta+\nu_{ij},
\qquad
\nu_{ij}:=\varepsilon_{ij}-\varepsilon_{i1}.
\end{equation}
Let \(d:=J-1\),
\(\Delta U_i:=(\Delta U_{i2},\ldots,\Delta U_{iJ})^\top\),
\(\nu_i:=(\nu_{i2},\ldots,\nu_{iJ})^\top\),
and let \(\Delta X_i\) be the \(d \times K\) matrix whose
\((j-1)\)th row, corresponding to alternative \(j\),
is \((x_{ij}-x_{i1})^\top\). Then
\[
\Delta U_i = \Delta X_i\beta+\nu_i,
\qquad
\nu_i \sim \mathcal{N}(0,\Omega).
\]
The covariance matrix of differenced errors satisfies
\[
\Omega_{j-1,k-1}
=
\Sigma_{jk}
+
\Sigma_{11}
-
\Sigma_{j1}
-
\Sigma_{1k},
\qquad
j,k=2,\ldots,J.
\]
Thus the first coordinate of \(\nu_i\) is \(\nu_{i2}\), and the first
row and column of \(\Omega\) correspond to alternative 2 relative to
the base alternative.
The base alternative is a labeling convention. It fixes the coordinate
system for utility differences but does not restrict substitution
patterns.
Let \(W_i:=\Delta U_i\), let
\(w:=(w_2,\ldots,w_J)^\top\in\mathbb{R}^d\) denote a realization of
\(W_i\), and set \(w_1:=0\) for the base alternative.
The observed choice corresponds to a fixed region of the differenced
latent-utility space. If the base alternative is chosen,
\[
\mathcal{D}_1
:=
\left\{
w \in \mathbb{R}^d
\mid
w_j < 0 \ \text{for all } j=2,\ldots,J
\right\}.
\]
For a non-base alternative \(j=2,\ldots,J\),
\[
\mathcal{D}_j
:=
\left\{
w \in \mathbb{R}^d
\mid
w_j > 0
\ \text{and}\
w_j > w_k \ \text{for all } k=2,\ldots,J,\ k\neq j
\right\}.
\]
Thus the probability of choice \(j\) is the integral of the
\(\mathcal{N}(\Delta X_i\beta,\Omega)\) density of
\(\Delta U_i\) over the fixed region \(\mathcal{D}_j\):
\begin{equation}\label{eq:mnp-prob}
P_{ij}(\theta)
:=
\int_{\mathcal{D}_j}
\phi_d(w;\Delta X_i\beta,\Omega)
\,dw,
\end{equation}
where \(\theta\) collects \(\beta\) and the free covariance parameters
specified below, and \(\phi_d(\cdot;\mu,\Omega)\) denotes the
\(d\)-variate normal density.
\subsection{Covariance Normalization}
The model is identified only up to scale. For any \(c>0\), the
transformation
\[
(\beta,\Omega)
\mapsto
(c\beta,c^2\Omega)
\]
leaves all choice probabilities unchanged. The main SAUSS
implementation uses the covariance Cholesky normalization
\begin{equation}\label{eq:covariance-cholesky}
\Omega:=LL^\top,
\qquad
L :=
\begin{pmatrix}
1 & & \\
a_{21} & e^{\lambda_2} & \\
\vdots & \ddots & \ddots \\
a_{d1} & \cdots & a_{d,d-1} & e^{\lambda_d}
\end{pmatrix}.
\end{equation}
Thus \(L_{11}=1\) fixes the standard deviation of \(\nu_{i2}\) at one
and thereby fixes the scale; the remaining diagonal elements are
positive by construction, and the strictly lower-triangular elements
are unrestricted. When \(d=1\), this normalization reduces to
\(\Omega=1\), as in binary probit.
Let
\[
\theta:=(\beta^\top,\theta_L^\top)^\top,
\qquad
\theta_L:=\operatorname{free}(L).
\]
Here \(\operatorname{free}(L)\) stacks the strictly lower-triangular
entries of \(L\) and the log-diagonal coordinates
\(\lambda_k:=\log L_{kk}\), \(k=2,\ldots,d\), excluding the fixed
element \(L_{11}\).
When \(d=1\), \(\theta_L\) is empty and the MNP normalization reduces
to the usual binary-probit scale normalization.
\subsection{Likelihood and Descent Target}
Let \(d_{ij}:=\mathbb{I}\{y_i=j\}\), where \(\mathbb{I}\{\cdot\}\)
denotes the indicator function. The sample log-likelihood is
\begin{equation}\label{eq:mnp-loglik}
\mathcal{L}_n(\theta)
:=
\sum_{i=1}^n \sum_{j=1}^J
d_{ij}\log P_{ij}(\theta).
\end{equation}
Because \(P_{ij}(\theta)\) lacks a closed form, simulated maximum
likelihood replaces it by a simulated probability estimate. If
\(\tilde P_{ij}^{(r)}(\theta;u_i^{(r)})\) is the contribution from
the \(r\)-th sequence of random draws, the standard simulated
objective is
\begin{equation}\label{eq:mnp-sim-loglik}
\hat{\mathcal{L}}_n(\theta)
:=
\sum_{i=1}^n \sum_{j=1}^J
d_{ij}
\log
\left(
\frac{1}{R}
\sum_{r=1}^{R}
\tilde P_{ij}^{(r)}(\theta;u_i^{(r)})
\right).
\end{equation}
SAUSS does not maximize \eqref{eq:mnp-sim-loglik}. It targets the
score of \eqref{eq:mnp-loglik} directly through simulation and uses
the negative simulated score as a descent direction for the negative
log-likelihood.
\subsection{Unbiased Simulated Scores}
The score for one observed choice is
\[
s_i(\theta)
:=
\nabla_\theta \log P_{iy_i}(\theta).
\]
By Lemma~\ref{lem:hm-identity},
\[
s_i(\theta)
=
\mathbb{E}_\theta
\left[
\nabla_\theta
\log \phi_d(W_i;\Delta X_i\beta,\Omega)
\mid
W_i\in\mathcal{D}_{y_i}
\right].
\]
The accept--reject simulator draws
\(\nu_i\sim\mathcal{N}(0,\Omega)\), forms
\(W_i:=\Delta X_i\beta+\nu_i\), and keeps the draw if
\(W_i\in\mathcal{D}_{y_i}\). Thus an accepted draw has distribution
\[
\nu_i
\mid
\{\Delta X_i\beta+\nu_i\in\mathcal{D}_{y_i}\}.
\]
Conditional on acceptance, \(W_i\) has the latent normal distribution
restricted to the observed choice region, so the average latent-score
contribution over accepted draws is an unbiased estimate of
\(s_i(\theta)\). This statement is exact for any fixed number
\(R \geq 1\) of accepted draws, provided the simulator runs until
those accepted draws are obtained.
In the one-draw score formulas below, \(\nu_i\) denotes a generic
accepted draw from this observation-specific conditional distribution,
equivalently \(\nu_i=W_i-\Delta X_i\beta\). Let \(Q:=\Omega^{-1}\) for
notational convenience. We suppress the dependence of this accepted
draw on \(\beta\) and \(\Omega\) when no confusion arises.
The one-draw complete-data score contribution for \(\beta\) is
\begin{equation}\label{eq:mnp-beta-score}
s_{\beta i}^{\mathrm{c}}
:=
\Delta X_i^\top Q \nu_i.
\end{equation}
The matrix score with respect to \(\Omega\) is
\begin{equation}\label{eq:mnp-omega-score}
S_{\Omega i}^{\mathrm{c}}
:=
\frac{1}{2}
Q(\nu_i\nu_i^\top-\Omega)Q.
\end{equation}
Because \(\Omega=LL^\top\), the corresponding score for the
Cholesky factor is
\begin{equation}\label{eq:mnp-L-score}
S_{L i}^{\mathrm{c}}
:=
2S_{\Omega i}^{\mathrm{c}}L
=
Q(\nu_i\nu_i^\top-\Omega)QL.
\end{equation}
The descent implementation uses the negative of these score
contributions. For any \(d\times d\) matrix \(A\), write
\(\mathcal{C}_{\theta_L}(A;L)\) for the vector that keeps the free
lower-triangular entries of \(A\), replaces each free diagonal entry
\(A_{kk}\), \(k=2,\ldots,d\), by \(A_{kk}L_{kk}\) to account for the
log-diagonal parameterization
\(\lambda_k:=\log L_{kk}\), \(k=2,\ldots,d\), and drops the fixed coordinate
\((1,1)\).
Let \(\nu_{ir}\), \(r=1,\ldots,R\), denote independent accepted
draws from the conditional distribution above for observation \(i\).
Averaging the negative score contributions over these draws gives
\begin{equation}\label{eq:mnp-one-observation-gradient}
\widehat g_{\beta i}(\theta)
:=
-\frac{1}{R}
\sum_{r=1}^{R}
\Delta X_i^\top Q\nu_{ir},
\qquad
\widehat g_{L i}(\theta)
:=
-\frac{1}{R}
\sum_{r=1}^{R}
\mathcal{C}_{\theta_L}
\{Q(\nu_{ir}\nu_{ir}^\top-\Omega)QL;L\}.
\end{equation}
Write
\[
\widehat g_i(\theta)
:=
\bigl(
\widehat g_{\beta i}(\theta)^\top,
\widehat g_{L i}(\theta)^\top
\bigr)^\top
\]
for the loss-gradient contribution in the full parameter vector.
\subsection{SAUSS Algorithm}
The main recursion uses a scalar decreasing learning rate,
\begin{equation}\label{eq:mnp-learning-rates}
\gamma_t
:=
{\gamma_0}(t-1)^{-a},
\qquad
a\in(1/2,1),
\qquad t\geq2.
\end{equation}
This is the standard decreasing-rate form used by the theory.
\begin{algorithm}[!htbp]
\DontPrintSemicolon
\caption{SAUSS for Multinomial Probit}\label{alg:sauss-mnp}
\KwIn{Data \(\{(\Delta X_i,y_i)\}_{i=1}^n\); \(\theta_{\mathrm{init}}\); tuning parameters \(\{T,T_0,m,R,{\gamma_0},a\}\), with \(T\geq2\) and \(T_0\in\{0,\ldots,T-1\}\)}
\KwOut{Averaged parameter \(\bar\theta_T\), reported as \(\hat\beta\) and \(\hat\Omega=\hat L\hat L^\top\)}
Set \(\theta_1\leftarrow\theta_{\mathrm{init}}\), \(N_{\mathrm{avg}}\leftarrow 0\), and \(\bar\theta\leftarrow 0\)\;
\If{\(T_0=0\)}{
\(N_{\mathrm{avg}}\leftarrow 1\) and \(\bar\theta\leftarrow\theta_1\)\;
}
\For{\(t=2,\ldots,T\)}{
Draw \(I_{t1},\ldots,I_{tm}\) independently and uniformly from \(\{1,\ldots,n\}\), and set \(\mathcal{B}_t\leftarrow(I_{t1},\ldots,I_{tm})\)\;
Construct \(\beta\) and \(L\) from \(\theta_{t-1}\), with \(L_{11}=1\)\;
Set \(\Omega\leftarrow LL^\top\) and \(Q\leftarrow\Omega^{-1}\)\;
Set \(\gamma_t\leftarrow{\gamma_0}(t-1)^{-a}\)\;
\For{\(b=1,\ldots,m\)}{
Set \(i\leftarrow I_{tb}\)\;
Run the accept--reject simulator until \(R\) accepted draws are obtained for observation \(i\)\;
Compute \(\widehat g_{tb}(\theta_{t-1})\) from \eqref{eq:mnp-one-observation-gradient} for observation \(i\), using fresh accepted draws\;
}
\(\bar g_t\leftarrow m^{-1}\sum_{b=1}^m\widehat g_{tb}(\theta_{t-1})\)\;
\(\theta_t\leftarrow \theta_{t-1}-\gamma_t\bar g_t\)\;
\If{\(t>T_0\)}{
\(N_{\mathrm{avg}}\leftarrow N_{\mathrm{avg}}+1\)\;
\(\bar\theta\leftarrow\bar\theta+(\theta_t-\bar\theta)/N_{\mathrm{avg}}\)\;
}
}
Set \(\bar\theta_T\leftarrow\bar\theta\)\;
Construct \(\hat\beta\) and \(\hat L\) from \(\bar\theta_T\), and set \(\hat\Omega\leftarrow\hat L\hat L^\top\)\;
\end{algorithm}
SAUSS resamples the simulation draws whenever an observation enters a
score evaluation. This differs from simulated likelihood implementations
that use common random numbers to stabilize a fixed simulated likelihood
objective. Here the algorithm does not optimize such a fixed simulated
objective. It uses simulation to produce an unbiased
stochastic direction for the gradient of the exact average negative
log-likelihood, so the simulation randomness is renewed along the
stochastic approximation path.
Algorithm~\ref{alg:sauss-mnp} takes as data input the
differenced covariates and observed choices
\(\{(\Delta X_i,y_i)\}_{i=1}^n\). The algorithmic parameter vector is
\(\theta=(\beta^\top,\theta_L^\top)^\top\). Although the algorithm is
written in terms of \(\theta\), each score evaluation reconstructs
\(\beta\) and the Cholesky factor \(L\) from \(\theta\), imposes
\(L_{11}=1\), and sets \(\Omega=LL^\top\).
The initial value \(\theta_{\mathrm{init}}\) is therefore
an input to the procedure.
The computational budget is summarized by \(T\), \(m\), and \(R\).
Here \(T\) is the final iterate index, so the algorithm performs
\(T-1\) stochastic approximation updates; \(m\) is
the number of observations in each mini-batch, and \(R\) is the
number of accepted simulation draws used to form each
observation-level loss-gradient. With iid mini-batches sampled from
the empirical distribution, \(m(T-1)/n\) is the expected number of
passes through the sample.
The averaging input \(T_0\) determines when Polyak--Ruppert averaging
begins. Algorithm~\ref{alg:sauss-mnp} reports the average of the
iterates \(\theta_{T_0+1},\ldots,\theta_T\), consistent with the
framework in Section~\ref{sec:sauss-framework}; when \(T_0=0\), this
average includes the initial estimator \(\theta_1\).
{The asymptotic results in
Section~\ref{sec:asymptotic-theory} apply when \(T_0=T_0(T)\) satisfies
\(T_0/T\to0\).}
The learning-rate inputs in Algorithm~\ref{alg:sauss-mnp} are
\({\gamma_0}\) and \(a\). The exponent \(a\in(1/2,1)\) gives the
decreasing-rate form used in the theory, while \({\gamma_0}\) controls the
overall scale of the stochastic approximation step. For simplicity,
the main algorithm uses a single scalar learning rate.
{In computation, fixed block-specific scale constants for the
coefficient and Cholesky components can be used as a diagonal
preconditioner, as described in Remark~\ref{rem:minibatch}.}
Algorithm~\ref{alg:sauss-mnp} describes the exact recursion: each
observation-level simulator runs until \(R\) accepted draws are
obtained, so the simulated score contribution is unbiased for
\(s_i(\theta)\) for every fixed \(R\geq1\); equivalently,
\(\widehat g_i(\theta)\) is unbiased for \(-s_i(\theta)\). Thus the
ideal recursion uses a fixed number \(R\) of accepted draws, but the
number of proposal draws, and hence runtime, is random.
\begin{remark}[Capped computation]
In practice, we impose a large cap \(K_{\max}\) on the number of proposal
draws. If the cap is reached after at least one acceptance but before
\(R\) acceptances, averaging the available accepted draws remains
centered on the exact loss-gradient, but uses a random, smaller
simulation size and has a different conditional variance. If no draw
is accepted, returning a zero contribution may bias the capped score
estimator. We therefore monitor trials per accepted draw and the
frequencies of partial-cap and no-acceptance events.
\end{remark}
\section{Asymptotic Theory}\label{sec:asymptotic-theory}
This section states the asymptotic theory for the SAUSS recursion with the exact simulator, which runs until \(R\) accepted draws are obtained for each observation. Inference based on these results is developed in Section~\ref{sec:inference}. Throughout, \(n\) denotes the sample size, \(T\) is the final iterate index, \(m\) is the mini-batch size, and \(R\) is the number of accepted simulation draws per observation. The initial estimator is \(\theta_1\), and updates produce \(\theta_t\) for \(t=2,\ldots,T\), so the run contains \(T-1\) stochastic approximation updates.
{All statements indexed by the run length are taken along a
sequence in which \(T\to\infty\) and the pilot size
\(n_0=n_0(T)\to\infty\); no relative rate between \(n_0\) and \(T\) is
imposed. The mini-batch size \(m\), accepted-draw count \(R\), and
learning-rate constants \((\gamma_0,a)\), as well as the parameter
dimension \(p\) (and hence \(J\) in the MNP specialization), are held
fixed.}
The recursion is \eqref{eq:sa-update} of Section~\ref{sec:sgd}, with mini-batch direction \(\hat G_t(\theta)\) built from the one-draw loss-gradient simulator \(g_i(\theta,U)=-\psi_i(\theta,U)\). Let \(\mathbb{E}_t^*\) denote expectation over the simulation draws in iteration \(t\), conditional on the past and the selected mini-batch, and let \(\mathcal F_t\) be the sigma-field generated by \(\theta_1\) and all mini-batch and simulation draws through iteration \(t\). For \(t\geq2\), define
\[
G_t(\theta)
=
\mathbb{E}_t^*\{\hat G_t(\theta)\},
\qquad
\mathcal R(\theta)
=
\mathbb{E}\{G_t(\theta)\mid\mathcal F_{t-1}\},
\qquad
H=\nabla_\theta\mathcal R(\theta_0).
\]
Under population sampling, \(\mathcal R(\theta)=\mathbb{E}\{g_i(\theta)\}=-\mathbb{E}\{s_i(\theta)\}\), so the target \(\theta_0\) satisfies \(\mathcal R(\theta_0)=0\), and \(H\) is the negative Jacobian of the population score, that is, the population information matrix under correct specification.
\subsection{Assumptions}
The assumptions below are stated for the recursion with the exact simulator; the proof appendix uses them with \(m=1\) (see Remark~\ref{rem:minibatch}). Throughout, \(C<\infty\) and \(c>0\) denote generic constants.
\begin{assumption}[Local identification and information]\label{ass:population-info}
The normalized parameter space \(\Theta\subset\mathbb{R}^p\) contains \(\theta_0\) as an interior point, and \(\theta_0\) is the only zero of \(\mathcal R(\theta)\) in a neighborhood \(\mathcal N_0\) of \(\theta_0\). The Jacobian \(H(\theta)=\nabla_\theta\mathcal R(\theta)\) exists, is Lipschitz continuous, and satisfies \(\|H(\theta)\|\leq C\) on \(\mathcal N_0\); \(H=H(\theta_0)\) is symmetric with \(\lambda_{\min}(H)>c\).
\end{assumption}
Assumption~\ref{ass:population-info} gives the local drift condition used for consistency: by the integral form of the mean value theorem, there is a ball \(B_\rho=\{\theta:\|\theta-\theta_0\|\leq\rho\}\subset\mathcal N_0\) on which
\[
(\theta-\theta_0)'\mathcal R(\theta)
\geq
c_0\|\theta-\theta_0\|^2
\quad\text{and}\quad
\|\mathcal R(\theta)-H(\theta-\theta_0)\|\leq C\|\theta-\theta_0\|^2
\]
for some \(c_0>0\). This is the descent-form analogue of local concavity of the population log-likelihood; symmetry of \(H\) holds for likelihood problems, where \(H\) is the population information matrix.
\begin{assumption}[Local initialization]\label{ass:algorithm-stability}
The recursion is initialized by a preliminary estimator \(\theta_1=\tilde\theta_{n_0}\) satisfying \(\tilde\theta_{n_0}\xrightarrow{p}\theta_0\) as \(n_0\to\infty\).
\end{assumption}
\begin{assumption}[Smoothness of the exact loss-gradient]\label{ass:smooth-G}
For each observation, the exact loss-gradient \(G_t(\theta)=\mathbb{E}_t^*\{\hat G_t(\theta)\}\) is continuously differentiable in \(\theta\), and
\[
\sup_{\theta\in\mathcal N_0}\mathbb{E}\big(\|\nabla_\theta G_t(\theta)\|^2\mid\mathcal F_{t-1}\big)\leq C .
\]
\end{assumption}
The simulated loss-gradient \(g_i(\theta,U)\) may be discontinuous in \(\theta\); Assumption~\ref{ass:smooth-G} concerns only its conditional expectation given the data, which for LDV models is the exact score contribution and is smooth.
\begin{assumption}[Unbiased simulated loss-gradients and moments]\label{ass:unbiased-scores}
For every \(\theta\in\mathcal N_0\), the one-draw simulator satisfies
\[
\mathbb{E}\{g_i(\theta,U)\mid D_i\}
=
g_i(\theta)
=
-s_i(\theta),
\]
and the simulation draws are fresh across observations, accepted draws, and stochastic approximation iterations. Let \(a\in(1/2,1)\) be the learning-rate exponent of Assumption~\ref{ass:sampling-steps}. For some \(\delta>0\) and some \(p_0>(1-a)^{-1}\), with \(q=\max\{2+\delta,2p_0\}\),
\[
\sup_{\theta\in\mathcal N_0}
\mathbb{E}\big(\|g_i(\theta,U)\|^{q}\mid\mathcal F_{t-1}\big)
\leq C,
\qquad
\sup_{\theta\in\mathcal N_0}
\mathbb{E}\big(\|G_t(\theta)\|^{q}\mid\mathcal F_{t-1}\big)
\leq C .
\]
\end{assumption}
The \(2+\delta\) moments are used for the Lindeberg condition; the \(2p_0\) moments are used to control the Polyak--Ruppert weighted remainder, and \(p_0\) grows as \(a\to1\).
\begin{assumption}[Continuity of the simulation variance]\label{ass:sim-variance}
Let \(A_t(\theta)=\operatorname{Var}\{g_i(\theta,U)\mid D_i\}\) be the conditional covariance of the one-draw simulated loss-gradient given the observation. Then \(\theta\mapsto\mathbb{E}\{A_t(\theta)\mid\mathcal F_{t-1}\}\) is Lipschitz continuous on \(\mathcal N_0\) with a bounded Lipschitz constant.
\end{assumption}
Assumption~\ref{ass:sim-variance} replaces a mean-square continuity condition on the simulator itself: it restricts only the conditional variance of the simulated loss-gradient, not a coupling of draws at different parameter values, and it is what the variance-consistency step of the martingale functional central limit theorem uses.
\begin{assumption}[Mini-batches and learning rates]\label{ass:sampling-steps}
At each update \(t\geq2\), the mini-batch indices \(I_{t1},\ldots,I_{tm}\) are sampled independently with replacement from the population distribution, or from the empirical distribution in the fixed-sample version, independently of the past conditional on the current iterate. The learning rate is
\[
\gamma_t=\gamma_0(t-1)^{-a},
\qquad
a\in(1/2,1),
\qquad t\geq2.
\]
Hence \(\sum_{t\geq2}\gamma_t=\infty\), \(\sum_{t\geq2}\gamma_t^2<\infty\), and \(\sum_{t\geq2}\gamma_t/\sqrt{t}<\infty\).
\end{assumption}
\begin{assumption}[Conditional variance limit]\label{ass:variance-limit}
For \(t\geq2\), let
\[
\zeta_t(\theta)=\hat G_t(\theta)-G_t(\theta),
\qquad
\xi_t(\theta)=G_t(\theta)-\mathcal R(\theta),
\]
and write \(\varepsilon_t=\zeta_t(\theta_{t-1})+\xi_t(\theta_{t-1})\) for the noise evaluated at the current iterate, which is a martingale difference with respect to \(\mathcal F_t\). The limit
\[
V_{m,R}
=
\operatorname{plim}_{t\to\infty}
\Big[
\operatorname{Var}\{G_t(\theta_0)\mid\mathcal F_{t-1}\}
+
\frac{1}{mR}\mathbb{E}\{A_t(\theta_0)\mid\mathcal F_{t-1}\}
\Big]
\]
exists and is positive definite.
\end{assumption}
For iid mini-batches with replacement and iid simulation draws, the variance in Assumption~\ref{ass:variance-limit} has the explicit decomposition
\begin{equation}\label{eq:vmr}
V_{m,R}
=
\frac{1}{m}\operatorname{Var}\{g_i(\theta_0)\}
+
\frac{1}{mR}
\mathbb{E}
\left[
\operatorname{Var}\{g_i(\theta_0,U)\mid D_i\}
\right].
\end{equation}
The first term is mini-batch sampling variation in the exact loss-gradient. The second term is simulation variation conditional on the observation. Since \(g_i(\theta_0)=-s_i(\theta_0)\), the first term equals \(\operatorname{Var}\{s_i(\theta_0)\}/m\).
\begin{assumption}[Smooth transformations]\label{ass:transformation}
{When inference concerns a scalar parameter \(\kappa=h(\theta_0)\), the map \(h:\Theta\to\mathbb R\) is twice continuously differentiable on \(\mathcal N_0\) and satisfies
\[
\sup_{\theta\in\mathcal N_0}\|\nabla^2h(\theta)\|
+\|\nabla h(\theta_0)\|
\leq C.
\]}
\end{assumption}
\subsection{Main stochastic approximation result}
\begin{theorem}[Consistency, linear representation, and FCLT]\label{thm:sa-linear}\label{thcflc}
Let \(m=1\), so that each update uses one observation and \(R\) accepted draws, and suppose Assumptions~\ref{ass:population-info}--\ref{ass:variance-limit} hold. Define the average and, for \(r\in[0,1]\), the partial-sum and limiting processes by
\[
\begin{aligned}
\bar\theta_T
&=
\frac{1}{T}\sum_{t=1}^{T}\theta_t, \\
X_T(r)
&=
\frac{1}{\sqrt{T}}
\sum_{t=1}^{\lfloor Tr\rfloor}(\theta_t-\theta_0), \\
X(r)
&=
\big(H^{-1}V_{1,R}H^{-1}\big)^{1/2}W(r),
\end{aligned}
\]
where \(W\) is a \(p\)-dimensional standard Wiener process.
{Let \(\tau\) be the first exit time from \(B_\rho\) and set
\(\mathcal E:=\{\tau=\infty\}\). Then, for every \(\epsilon>0\), there
exists \(\bar\gamma(\epsilon)>0\) such that, if
\(\gamma_0\leq\bar\gamma(\epsilon)\),
\(\mathbb P(\mathcal E^c)\leq\epsilon+o(1)\), and the following
conclusions hold:}
\begin{enumerate}
\item[(i)] (consistency) {\(\|\theta_t-\theta_0\|\mathbb I\{\mathcal E\}\xrightarrow{p}0\)} as \(t\to\infty\);
\item[(ii)] (linear representation)
\begin{equation}\label{eq:pr-linear-rep}
\sqrt{T}(\bar\theta_T-\theta_0)
=
-H^{-1}
\frac{1}{\sqrt{T}}
\sum_{t=2}^{T}
\{\zeta_t(\theta_{t-1})+\xi_t(\theta_{t-1})\}
+
Z_T,
\qquad {Z_T\mathbb I\{\mathcal E\}\xrightarrow{p}0},
\end{equation}
and the same representation holds uniformly in \(r\) for the partial sums \(X_T(r)\);
\item[(iii)] (FCLT) for every bounded continuous \(f\) on \(D[0,1]^p\),
\[
\limsup_{T\to\infty}\Big|\mathbb{E}f(X_T)-\mathbb{E}f(X)\Big|\leq 2\epsilon\sup|f| ;
\]
in particular \(\sqrt{T}(\bar\theta_T-\theta_0)=X_T(1)\) is asymptotically \(\mathcal N(0,H^{-1}V_{1,R}H^{-1})\) up to the same \(\epsilon\)-error.
\end{enumerate}
\end{theorem}
\begin{remark}[On the \(\epsilon\)-formulation]\label{rem:epsilon}
{The tolerance \(\epsilon\) originates in a single step of the proof. For any \(b\in(0,\rho)\), the stopped-process argument and a supermartingale maximal inequality give
\[
\mathbb P(\tau<\infty)
\leq
\mathbb P(\|\theta_1-\theta_0\|>b)
+
\frac{b^2+C\gamma_0^2\sum_{s\geq1}s^{-2a}}{\rho^2}.
\]
For a given \(\epsilon\in(0,1)\), set \(b_\epsilon=\rho\sqrt{\epsilon/2}\) and choose \(\gamma_0\) sufficiently small that
\[
C\gamma_0^2\sum_{s\geq1}s^{-2a}
\leq
\frac{\rho^2\epsilon}{2}.
\]
Then
\[
\mathbb P(\tau<\infty)
\leq
\mathbb P(\|\theta_1-\theta_0\|>b_\epsilon)
+\epsilon
=
\epsilon+o(1),
\]
where the final equality follows from pilot consistency as \(n_0=n_0(T)\to\infty\). All conclusions of the theorem are established on the event \(\{\tau=\infty\}\), and \(\epsilon\) bounds the probability of its complement.
Ordinary weak convergence without the \(\epsilon\)-error is not available under the stated local assumptions because the unbounded simulation noise gives the recursion a positive probability of leaving the neighborhood on which those assumptions hold. Establishing such convergence would require additional stability conditions controlling the recursion outside that neighborhood. These conditions are generally unsuitable for MNP because the population objective can flatten as latent variances increase. We therefore retain the local \(\epsilon\)-formulation, which applies to the diminishing learning-rate schedule used in implementation, with \(\gamma_0\) and \(a\) held fixed as \(T\) increases.}
\end{remark}
\begin{remark}[Mini-batches, step-size shift, and burn-in]\label{rem:minibatch}
{Theorem~\ref{thm:sa-linear} is stated and proved for \(m=1\). For any fixed \(m\), suppose that the mini-batch observations are sampled independently with replacement and that their simulation draws are conditionally independent across contributions. The mini-batch direction \(\hat G_t(\theta)\) in \eqref{eq:sa-update} is then the average of \(m\) conditionally independent copies of the one-observation direction. It has the same conditional mean, and its conditional variance is \(1/m\) times the one-observation variance, yielding \(V_{m,R}=V_{1,R}/m\). Its moment and Lipschitz bounds follow from those of the individual contributions with adjusted constants. The same argument therefore applies with \(V_{1,R}\) replaced by \(V_{m,R}\).
{The same result covers a fixed symmetric positive-definite
preconditioner \(P\) in the update
\[
\theta_t
=
\theta_{t-1}
-
\gamma_tP\hat G_t(\theta_{t-1}),
\qquad
\gamma_t=\gamma_0(t-1)^{-a}.
\]
In the centered coordinate
\(\phi=P^{-1/2}(\theta-\theta_0)\), the drift and innovation covariance
are \(P^{1/2}HP^{1/2}\) and \(P^{1/2}V_{m,R}P^{1/2}\),
respectively. Applying the theorem in this coordinate and mapping back
preserves the Polyak--Ruppert covariance
\(H^{-1}V_{m,R}H^{-1}\), with an adjusted small-step threshold.}
A fixed shift \(\gamma_t=\gamma_0(t-1+t_0)^{-a}\) also leaves the limits unchanged. For the averaged estimator, let \(T_0=T_0(T)\) be deterministic and satisfy \(T_0/T\to0\). The functional limit theorem then gives \(\{\sum_{t=1}^{T_0}(\theta_t-\theta_0)\}\mathbb I\{{\mathcal E}\}=o_{\mathbb P}(\sqrt T)\), and \(T/(T-T_0)\to1\), so the linear representation and functional limit remain unchanged in the same \(\epsilon\)-sense as Theorem~\ref{thm:sa-linear}.}
\end{remark}
{Many scalar parameters of interest are nonlinear transformations of \(\theta_0\), for example a particular substitution elasticity or choice probability in MNP. Let \(h\) be such a transformation and define \(\kappa_t=h(\theta_t)\) and \(\kappa_0=h(\theta_0)\).}
\begin{theorem}[Smooth transformations]\label{thm:transformation}\label{thnonlinear}
Under the conditions of Theorem~\ref{thm:sa-linear}, suppose in addition that Assumption~\ref{ass:transformation} holds. Then, in the same sense as in Theorem~\ref{thm:sa-linear}(iii),
{
\[
\frac{1}{\sqrt{T}}\sum_{t=1}^{\lfloor Tr\rfloor}(\kappa_t-\kappa_0)
\Rightarrow
\nabla h(\theta_0)'\big(H^{-1}V_{1,R}H^{-1}\big)^{1/2}W(r),
\qquad r\in[0,1].
\]
}
\end{theorem}
Theorem~\ref{thm:transformation} says that the scalar transformed path \(\kappa_t\) obeys the same functional limit as a linear functional of \(\theta_t\), so the path-based inference of Section~\ref{sec:inference} applies to \(\kappa_t\) directly, without a delta-method variance estimator.
\section{Inference}\label{sec:inference}
This section collects the inference consequences of
Theorems~\ref{thm:sa-linear} and \ref{thm:transformation}. Three features
of SAUSS shape it. {First, the recursion has two sources of
algorithmic variation, mini-batch selection and score simulation. Their
magnitude relative to sampling uncertainty from the data depends on \(n\),
\(T\), \(m\), and \(R\).} Second, random-scaling inference uses only
the solution path, so it does not require estimates of \(H\) or of the
simulation variance. {Third, many scalar objects of interest, such as a particular substitution elasticity or choice probability, are nonlinear transformations of \(\theta_0\), and the transformed path can be used directly.}
\subsection{Sampling and algorithmic uncertainty}\label{sec:inference-regimes}
{For a sample of size \(n\), the same recursion can be viewed as a
randomized algorithm targeting the exact sample MLE.} Let
\[
S_n(\theta)=\frac{1}{n}\sum_{i=1}^{n}s_i(\theta),
\]
and let \(\hat\theta_n\) solve \(S_n(\hat\theta_n)=0\), suppressing standard existence qualifications. Under the usual likelihood expansion,
\[
\sqrt{n}(\hat\theta_n-\theta_0)
=
H^{-1}
\frac{1}{\sqrt{n}}
\sum_{i=1}^{n}s_i(\theta_0)
+
o_p(1).
\]
{For a fixed sample, the exact sample estimator \(\hat\theta_n\) is the algorithmic target. The sampling and algorithmic errors satisfy the exact decomposition}
\[
{
\sqrt{n}(\bar\theta_T-\theta_0)
=
\sqrt{n}(\hat\theta_n-\theta_0)
+
\sqrt{n}(\bar\theta_T-\hat\theta_n).
}
\]
{
\begin{remark}[Finite-run variance approximation]\label{rem:finite-run-variance}
For a finite sample and a finite number of stochastic approximation iterations, a useful first-order approximation is
\[
\operatorname{Var}(\bar\theta_T-\theta_0)
\approx
\frac{1}{n}H^{-1}\Sigma_sH^{-1}
+
\frac{1}{T}H^{-1}V_{m,R}H^{-1}.
\]
{Equivalently, this approximation suggests}
\begin{equation}\label{eq:sample-algorithm-combined}
\sqrt{n}(\bar\theta_T-\theta_0)
\overset{\mathrm{approx.}}{\sim}
\mathcal{N}
\left(
0,
H^{-1}\Sigma_sH^{-1}
+
{\frac{n}{T}}H^{-1}V_{m,R}H^{-1}
\right).
\end{equation}
Here \(\Sigma_s=\operatorname{Var}\{s_i(\theta_0)\}\). The first term represents sampling uncertainty in the exact sample estimator, while the second represents the additional uncertainty from terminating the randomized algorithm after \(T\) iterations. Thus, the sum can provide a more informative finite-run variance approximation than the conventional likelihood variance alone when \(T\) is not large relative to \(n\). By Remark~\ref{rem:minibatch}, \(V_{m,R}=V_{1,R}/m\), so the approximation also provides the reductions from larger mini-batches and additional accepted simulation draws.
A corresponding plug-in approximation is
\[
\widehat{\operatorname{Var}}(\bar\theta_T)
:=
\frac{1}{n}\hat H^{-1}\hat\Sigma_s\hat H^{-1}
+
\frac{1}{T}\hat H^{-1}\hat V_{m,R}\hat H^{-1}.
\]
The fixed-sample implementation underlying this approximation uses iid sampling with replacement from the empirical distribution and fresh simulation draws at every score evaluation. When \(m(T-1)>n\), observations are sampled multiple times on average, and \(m(T-1)/n\) is the expected number of passes.
\end{remark}
}
\subsection{Random scaling and sandwich inference}\label{sec:inference-rs}
The functional central limit theorem in Theorem~\ref{thm:sa-linear} gives two routes to inference on the algorithmic target. Let \(A\) be an \(\ell\times p\) matrix with full row rank. For random-scaling inference, consider the null restriction \(A\theta_0=a_0\).
\begin{corollary}[Inference consequences of the FCLT]\label{cor:inference}
Under the FCLT in Theorem~\ref{thm:sa-linear}, with the fixed-\(m\)
extension in Remark~\ref{rem:minibatch} when \(m>1\), let
\(\bar\theta_T=T^{-1}\sum_{t=1}^{T}\theta_t\) and define
\[
\widehat V_{\mathrm{RS}}(A)
=
\frac{1}{T^2}
\sum_{t=1}^{T}
\left\{
\sum_{s=1}^{t} A(\theta_s-\bar\theta_T)
\right\}
\left\{
\sum_{s=1}^{t} A(\theta_s-\bar\theta_T)
\right\}'.
\]
{Under the null restriction \(A\theta_0=a_0\), when
\(\widehat V_{\mathrm{RS}}(A)\) is nonsingular, the random-scaling statistic}
\[
T(A\bar\theta_T-a_0)'
\widehat V_{\mathrm{RS}}(A)^{-1}
(A\bar\theta_T-a_0)
\]
converges (in the sense of Theorem~\ref{thm:sa-linear}(iii)) to the standard self-normalized Brownian functional
\[
W_\ell(1)'
\left\{
\int_0^1
\bar W_\ell(r)\bar W_\ell(r)'dr
\right\}^{-1}
W_\ell(1),
\qquad
\bar W_\ell(r)=W_\ell(r)-rW_\ell(1).
\]
Alternatively, if \(\hat H\) consistently estimates \(H\) and \(\hat V_{m,R}\) consistently estimates \(V_{m,R}\), then
\[
\hat H^{-1}\hat V_{m,R}\hat H^{-1}
\xrightarrow{p}
H^{-1}V_{m,R}H^{-1}.
\]
Consequently,
\[
\frac{1}{T}
\hat H^{-1}\hat V_{m,R}\hat H^{-1}
\]
is a sandwich estimator for the first-order algorithmic covariance of
\(\bar\theta_T\).
\end{corollary}
{Under population sampling, the target in Corollary~\ref{cor:inference} is \(\theta_0\). Conditional on a fixed sample, path variation instead describes the remaining randomized-algorithm variation around the sample target \(\hat\theta_n\) and does not by itself account for sampling uncertainty. Moreover, \(\widehat V_{\mathrm{RS}}(A)\) is a self-normalizer with a nondegenerate Brownian-functional limit, not a consistent estimator of \(AH^{-1}V_{m,R}H^{-1}A'\).}
Random scaling uses the solution path and follows the online-inference logic developed in \citet{leeLiaoSeoShin2022randomScaling} and \citet{leeLiaoSeoShin2025quantile}. Sandwich inference is closer to conventional MLE inference and is useful when the algorithm is run long enough that computational uncertainty is intended to be small; the MNP form of \(V_{m,R}\) and its plug-in estimator are given next.
\paragraph{The MNP form of \(V_{m,R}\) and a plug-in estimator.}
For the MNP specialization, \(V_{m,R}\) can be written directly in
terms of accepted-draw loss-gradient contributions. {For an
accepted draw \(\nu\) of the differenced error vector defined relative to
alternative 1, define}
\[
{
\Gamma_i(\theta,\nu)
:=
\begin{pmatrix}
-\Delta X_i^\top Q\nu \\
-\mathcal{C}_{\theta_L}
\{Q(\nu\nu^\top-\Omega)QL;L\}
\end{pmatrix},
\qquad
\Omega:=LL^\top,
\qquad
Q:=\Omega^{-1}.
}
\]
Let \(\nu_i^\star(\theta_0)\) denote one accepted draw from the conditional distribution of the canonical error vector given \((\Delta X_i,y_i)\) at \(\theta_0\). Then
\[
g_i(\theta_0)
=
\mathbb{E}
\{
\Gamma_i(\theta_0,\nu_i^\star(\theta_0))
\mid
\Delta X_i,y_i
\},
\]
and
\begin{equation}\label{eq:mnp-vmr}
V_{m,R}^{\mathrm{MNP}}
=
\frac{1}{m}
\operatorname{Var}
\{g_i(\theta_0)\}
+
\frac{1}{mR}
\mathbb{E}
\left[
\operatorname{Var}
\{
\Gamma_i(\theta_0,\nu_i^\star(\theta_0))
\mid
\Delta X_i,y_i
\}
\right].
\end{equation}
Equivalently, if
\[
\widehat g_i^{R}(\theta_0)
=
\frac{1}{R}
\sum_{r=1}^{R}
\Gamma_i(\theta_0,\nu_{ir}^\star(\theta_0)),
\]
then
\begin{equation}\label{eq:mnp-vmr-total}
V_{m,R}^{\mathrm{MNP}}
=
\frac{1}{m}
\operatorname{Var}
\{
\widehat g_i^{R}(\theta_0)
\}.
\end{equation}
A direct plug-in estimator can be built from independent accepted-draw loss-gradient calls at {a consistent estimate \(\hat\theta\), where \(\hat\theta\xrightarrow{p}\theta_0\)}. Generate \(B\) independent \(R\)-accepted calls
\[
\widehat g_{ib}^{R}(\hat\theta)
=
\frac{1}{R}
\sum_{r=1}^{R}
\Gamma_i(\hat\theta,\hat\nu_{ibr}^\star),
\qquad
b=1,\ldots,B,
\]
set
\[
\bar g
=
\frac{1}{nB}
\sum_{i=1}^{n}
\sum_{b=1}^{B}
\widehat g_{ib}^{R}(\hat\theta),
\]
and use
\begin{equation}\label{eq:mnp-vmr-plugin}
\widehat V_{m,R}^{\mathrm{MNP}}
=
\frac{1}{m}
\frac{1}{nB-1}
\sum_{i=1}^{n}
\sum_{b=1}^{B}
\left\{
\widehat g_{ib}^{R}(\hat\theta)-\bar g
\right\}
\left\{
\widehat g_{ib}^{R}(\hat\theta)-\bar g
\right\}'.
\end{equation}
{For fixed \(m\), \(R\), and \(B\), under iid sampling,
fresh conditionally independent evaluation calls, and standard continuity
and local uniform-law conditions at
\(\hat\theta\xrightarrow{p}\theta_0\), the estimator in
\eqref{eq:mnp-vmr-plugin} is consistent for
\(V_{m,R}^{\mathrm{MNP}}\) as \(n\to\infty\). Calls sharing the same
observation \(i\) are conditionally independent given \(D_i\), rather than
unconditionally independent.}
With \(B=1\), this estimator uses one independent loss-gradient call per observation and estimates the total one-step variance in \eqref{eq:mnp-vmr-total}. {Larger \(B\) reduces Monte Carlo noise in estimating the sandwich input, and \(B\geq2\) is required to report the within- and between-observation decomposition in \eqref{eq:mnp-vmr}.}
\subsection{Nonlinear transformations}\label{sec:inference-nonlinear}
{Let the scalar parameter of interest be \(\kappa_0=h(\theta_0)\) for a smooth map \(h\) satisfying Assumption~\ref{ass:transformation}, and suppose \(\nabla h(\theta_0)\neq0\) so that its first-order limit is nondegenerate. Rather than estimating \(\nabla h(\hat\theta)\) and applying the delta method, we track the transformed path \(\kappa_t=h(\theta_t)\) inside the recursion and apply random scaling to \(\kappa_t\) directly. Theorem~\ref{thm:transformation} shows that the partial-sum process of \(\kappa_t-\kappa_0\) has the limit \(\nabla h(\theta_0)'(H^{-1}V_{1,R}H^{-1})^{1/2}W(r)\). Define
\[
\bar\kappa_T=\frac{1}{T}\sum_{t=1}^T\kappa_t,
\qquad
\widehat V_{\mathrm{RS}}^\kappa
=
\frac{1}{T^2}\sum_{t=1}^T
\left\{\sum_{s=1}^t(\kappa_s-\bar\kappa_T)\right\}^2.
\]
The scalar version of Corollary~\ref{cor:inference} therefore gives the self-normalized statistic
\[
\frac{T(\bar\kappa_T-\kappa_0)^2}
{\widehat V_{\mathrm{RS}}^\kappa}.
\]
}
\section{Monte Carlo Experiments}\label{sec:monte-carlo}
In this section we examine the finite-sample performance of the SAUSS estimator for the multinomial probit
model. We use the accept--reject simulator of Section~\ref{sec:mnp-sauss} with the covariance parameterization
\(\Omega=LL'\), identified by the normalization \(L_{11}=1\). A sequential accept--reject variant produced
similar results and we do not report it separately.
\subsection{Design and Implementation}\label{sec:mc-design}
The baseline design, labeled D1, is a multinomial probit model written in utility differences relative to
alternative~1. For individual \(i=1,\ldots,n\) and alternatives \(j=2,\ldots,J\),
\begin{equation}\label{eq:mc-dgp}
U_{ij}-U_{i1} = (x_{ij}-x_{i1})'\beta_0+\nu_{ij},
\qquad
\nu_i=(\nu_{i2},\ldots,\nu_{iJ})'\sim \mathcal{N}(0,\Omega_0),
\end{equation}
and the observed choice \(y_i\) is the alternative with the largest latent utility. Each alternative has \(K=3\)
covariates drawn independently from the standard normal distribution, and the true coefficient vector is
\(\beta_0=(-0.5,1,1)'\). The differenced errors are equicorrelated: \(\Omega_0\) has unit diagonal and all
off-diagonal elements equal to \(0.3\). The number of alternatives varies over \(J\in\{4,8,16,32,64\}\), so that
the dimension of \(\Omega_0\) ranges from \(3\times 3\) to \(63\times 63\). The sample size is fixed at
\(n=2000\) throughout, and each design cell is replicated \(1000\) times. We therefore focus on computational scaling with the number of alternatives and on finite-sample performance across different Monte Carlo designs.
The implementation follows Section~\ref{sec:mnp-sauss}. At each iteration a mini-batch of \(m=10\) observations
is drawn from the sample independently with replacement, and the accept--reject simulator is run until
\(R=5\) accepted draws are obtained for each observation in the mini-batch, subject to a trial cap of
\(K_{\max}=20{,}000\). If the simulator obtains no accepted draw for an observation within the cap, that
observation's contribution to the mini-batch gradient is omitted from the update.
{Let \(s=1,\ldots,80{,}000\) index the stochastic updates. The block-specific learning rates are \(\gamma_{\beta,s}=0.5s^{-0.501}\) and
\(\gamma_{L,s}=0.2s^{-0.501}\); the exponent lies in the interval \((1/2,1)\) required for the asymptotic
optimality of the averaged iterate \citep{polyak1992acceleration}. We use Polyak--Ruppert averaging to
construct the reported estimates.} The tuning constants were fixed in a small
pilot study at \(J=4\) and then held constant across all designs and all values of \(J\); a tuning sweep at
\(J=64\), not reported here, produced stable results in a neighborhood of these values.
Table~\ref{tab:mc-design} summarizes the design.
\begin{table}[t]
\centering
\caption{Monte Carlo design}
\label{tab:mc-design}
\small
\begin{tabular}{@{}ll@{}}
\toprule
Model & Multinomial probit in utility differences, equation~\eqref{eq:mc-dgp} \\
Alternatives & \(J\in\{4,8,16,32,64\}\) (D1); \(J=4\) (D2, D3) \\
Covariates & \(K=3\), \(x_{ijk}\stackrel{iid}{\sim}\mathcal{N}(0,1)\) \\
Coefficients & D1: \(\beta_0=(-0.5,1,1)'\); D2, D3: \(\beta_0=(-1.5,2,2)'\) \\
Error covariance & D1: unit diagonal, equicorrelation \(0.3\); D2, D3: equation~\eqref{eq:mc-omega-d2} \\
Identification & \(\Omega=LL'\), \(L_{11}=1\) \\
Sample size & \(n=2000\) \\
Replications & \(1000\) per design cell \\
Mini-batch size & \(m=10\), sampled i.i.d.\ with replacement \\
Accepted draws & \(R=5\) per observation, trial cap \(K_{\max}=20{,}000\) \\
Step size & \(0.5s^{-0.501}\) (\(\beta\)), \(0.2s^{-0.501}\) (\(L\)) \\
Updates & \(80{,}000\), with Polyak--Ruppert averaging \\
\bottomrule
\end{tabular}
\end{table}
We evaluate the estimates with three criteria. {Although the normalization \(L_{11}=1\) fixes the scale of the
parameters, we report criteria that are invariant to the scale normalization so that comparisons with
implementations adopting other normalizations remain transparent.} First, the
coefficient criterion measures the recovery of the direction of \(\beta_0\). Writing
\(\tilde\beta=\beta/\|\beta\|\) and letting \(\hat\beta_b\) denote the estimate in replication
\(b=1,\ldots,B\),
\[
\mathrm{RMSE}_{\beta}
=
\Bigl\{ B^{-1}\textstyle\sum_{b=1}^B \|\tilde{\hat\beta}_b-\tilde\beta_0\|^2 \Bigr\}^{1/2}.
\]
Second, the correlation criterion measures the recovery of the error dependence structure. Let \(C(\Omega)\)
denote the correlation matrix implied by \(\Omega\) and \(d=J-1\); then
\[
\mathrm{RMSE}_{C}
=
\Bigl\{ \tfrac{2}{Bd(d-1)}\textstyle\sum_{b=1}^B\sum_{k<\ell}
\bigl(C_{k\ell}(\hat\Omega_b)-C_{k\ell}(\Omega_0)\bigr)^2 \Bigr\}^{1/2}.
\]
Third, the probability criterion measures out-of-sample fit of the implied choice probabilities. For each \(J\)
we fix a test design of \(N_{\mathrm{test}}=5000\) covariate draws from the distribution above, independent of
all estimation samples, and approximate the choice probabilities \(P_{qj}(\theta)\) at any parameter value by
frequency simulation with \(S_{\mathrm{eval}}=5000\) common random draws. The criterion is
\[
\mathrm{RMSE}_{P}
=
\Bigl\{ \tfrac{1}{BN_{\mathrm{test}}J}\textstyle\sum_{b=1}^B\sum_{q=1}^{N_{\mathrm{test}}}\sum_{j=1}^J
\bigl(\hat P_{qj}(\hat\theta_b)-\hat P_{qj}(\theta_0)\bigr)^2 \Bigr\}^{1/2},
\]
where the same test design and the same evaluation draws are used for all replications and all estimators within
a given \(J\), so that simulation noise in the evaluation step is common across the comparison. Because average
choice probabilities scale as \(1/J\), levels of \(\mathrm{RMSE}_P\) are not comparable across different values
of \(J\); we therefore use \(\mathrm{RMSE}_P\) only for comparisons within a given \(J\).
\subsection{Main Results}\label{sec:mc-crossj}
{We benchmark SAUSS against simulated maximum likelihood based on the Geweke--Hajivassiliou--Keane (GHK) simulator at \(J=4\), where the conventional estimator is computationally reliable \citep{borschSupanHajivassiliou1993,hajivassiliou1996simulation,Train2009}. The simulated log-likelihood is maximized numerically. We use the implementation in the R package \texttt{mlogit} \citep{croissant2020mlogit}, with \(100\) GHK draws per observation and the package's default starting values, optimizer, and convergence tolerances. Both estimators are applied to the same \(1{,}000\) simulated data sets and use the same normalization \(L_{11}=1\). For the larger choice sets, we report SAUSS alone to examine its computational and statistical scaling as \(J\) increases.}
\begin{table}[t]
\centering
\caption{SAUSS performance across \(J\) and comparison with GHK-based simulated maximum likelihood at \(J=4\) (design D1)}
\label{tab:mc-ghk}
\small
\setlength{\tabcolsep}{4pt}
\begin{tabular}{rlcccr}
\toprule
\(J\) & Estimator & \(\mathrm{RMSE}_{\beta}\) & \(\mathrm{RMSE}_{C}\) & \(\mathrm{RMSE}_{P}\) & Median time (s) \\
\midrule
4 & SAUSS & 0.0205 & 0.0997 & 0.0120 & 0.5 \\
& GHK-SML & 0.0202 & 0.0939 & 0.0122 & 71.8 \\
\midrule
8 & SAUSS & 0.0171 & 0.2010 & 0.0121 & 2.6 \\
\midrule
16 & SAUSS & 0.0159 & 0.2864 & 0.0101 & 9.8 \\
\midrule
32 & SAUSS & 0.0148 & 0.3077 & 0.0075 & 43.4 \\
\midrule
64 & SAUSS & 0.0140 & 0.3091 & 0.0056 & 80.5 \\
\bottomrule
\end{tabular}
\par
\begin{minipage}{0.95\textwidth}
\vspace{1ex}
\footnotesize
\emph{Notes:} GHK-SML is \texttt{mlogit} with \texttt{probit=TRUE} and \(100\) GHK draws at default settings. At \(J=4\), both estimators are applied to the same \(1{,}000\) simulated data sets.
\end{minipage}
\end{table}
Consider first the behavior of SAUSS across \(J\). The coefficient criterion improves as \(J\) grows with \(n\)
fixed, from \(0.0205\) at \(J=4\) to \(0.0140\) at \(J=64\): each observation contributes a comparison among
\(J\) alternatives, so the information about the direction of \(\beta_0\) per observation increases with \(J\).
{The correlation criterion moves in the opposite direction, consistent with the growing number of covariance
parameters relative to the fixed sample size: the number of free correlation parameters, \(d(d-1)/2\), grows
quadratically in \(J\).} The acceptance rate of the accept--reject simulator declines steeply with \(J\),
though somewhat more slowly than the rate \(1/J\) at the larger values, from \(0.281\) at \(J=4\) to \(0.033\)
at \(J=64\), as the observed-choice region occupies a shrinking share of the error space. {The
declining acceptance rate increases simulation effort, but no-acceptance events remain rare. Even at \(J=64\),
the simulator failed to produce an accepted draw in only \(143\) of the \(800{,}000\) observation-level
evaluations per replication, less than \(0.02\) percent.} {Computation time increases with \(J\), but the median is only \(80.5\) seconds per replication at \(J=64\).}
{At \(J=4\), the two estimators have similar accuracy and the
absolute differences across the three criteria are at most \(0.0058\).
Under the reported configurations, SAUSS has a median runtime of
approximately \(0.5\) seconds, compared with \(71.8\) seconds for GHK-SML,
a factor of about \(140\). In this small-\(J\) setting, SAUSS delivers
similar aggregate accuracy to the conventional estimator. The larger
choice sets provide computational
stress tests with an unrestricted covariance parameterization. In
particular, the \(J=64\) design demonstrates computational feasibility
rather than precise recovery of the full covariance matrix, with the
tuning constants held fixed across \(J\).}
\subsection{Heterogeneous Covariance and a Rare Alternative}\label{sec:mc-robust}
{The equicorrelated design is favorable to the accept--reject simulator because no choice region is unusually
small. We therefore consider two additional designs at \(J=4\) that probe robustness beyond this baseline.} Design D2 increases the
coefficients to \(\beta_0=(-1.5,2,2)'\), making choices more deterministic conditional on covariates, and
replaces \(\Omega_0\) with the matrix
\begin{equation}\label{eq:mc-omega-d2}
\Omega_0^{\mathrm{D2}}
=
\begin{pmatrix}
1.00 & 0.85 & 0.75 \\
0.85 & 1.40 & 1.05 \\
0.75 & 1.05 & 2.00
\end{pmatrix},
\end{equation}
which has unequal variances and large, uneven correlations. Design D3 augments D2 with alternative-specific
intercepts \(\delta_0=(1.5,0.5,-2.0)'\) in the utility differences, so that alternative~4 is chosen with low
probability: its population choice share is \(9.2\) percent, against \(41\) percent for the most popular
alternative. The intercepts are estimated by absorbing them into the covariate matrix as alternative dummies,
with no change to the estimator, and we report the direction criterion separately for the slope and intercept
blocks.
\begin{table}[t]
\centering
\caption{SAUSS performance under the baseline and stress designs at \(J=4\)}
\label{tab:mc-robust}
\small
\setlength{\tabcolsep}{4pt}
\begin{tabular}{lcccccc}
\toprule
Design & \(\mathrm{RMSE}_{\beta}\) & \(\mathrm{RMSE}_{\delta}\) & \(\mathrm{RMSE}_{C}\) & \(\mathrm{RMSE}_{P}\)
& Acceptance rate & Median time (s) \\
\midrule
D1 & 0.0205 & -- & 0.0997 & 0.0120 & 0.281 & 0.5 \\
D2 & 0.0138 & -- & 0.1341 & 0.0157 & 0.414 & 0.6 \\
D3 & 0.0144 & 0.0514 & 0.1586 & 0.0205 & 0.453 & 0.5 \\
\bottomrule
\end{tabular}
\par
\begin{minipage}{0.95\textwidth}
\vspace{1ex}
\footnotesize
\emph{Notes:} \(J=4\), \(n=2000\), \(1000\) replications, tuning as in Table~\ref{tab:mc-design}.
D1 is the equicorrelated baseline. \(\mathrm{RMSE}_{\delta}\) is the direction criterion for the intercept block
(D3 only). The rare alternative in D3 has a \(9.2\) percent choice share.
\end{minipage}
\end{table}
{Table~\ref{tab:mc-robust} compares the two stress designs with the \(J=4\) baseline. The overall acceptance
rates in D2 and D3, \(0.414\) and \(0.453\), exceed the baseline rate of \(0.281\). D3 nevertheless creates
substantial heterogeneity in computational difficulty: the acceptance rate for observations choosing the rare
alternative is \(0.26\), roughly two-fifths of the rate for the most frequently chosen alternative~2
(\(0.66\)). The simulator remains effective: on average about two observation evaluations per replication
produced no accepted draw, out of \(800{,}000\), and the trial cap was reached fewer than five times per
replication. The accuracy measures remain in the same broad range across the three designs. The coefficient
criterion is smaller in D2 and D3 than at the baseline, whereas the correlation and probability criteria are
somewhat larger. Because the data-generating processes differ, these comparisons should be interpreted as
robustness diagnostics rather than as a ranking of the designs. Overall, SAUSS continues to perform well under
the heterogeneous covariance structure and the rare-alternative design.}
\section{Application: Maternal Labor-Supply Intentions}
\label{sec:application}
We revisit the maternal labor-supply model of
\citet{boneva2026}, which asks whether respondents' beliefs about the
consequences of maternal employment predict the labor-supply choice that
they or their partner would make. Column 2 of its Table~4 fits a
multinomial probit to the stated choices of $N = 2{,}873$ respondents among
$J = 3$ alternatives, which are `not working', `part-time', and `full-time', respectively. In this scenario,
full-time daycare is available. Utility is linear in five
alternative-specific belief and norm measures and six case-specific
demographic controls, so the utility-differenced model carries $K = 19$
parameters. The published estimates use GHK-based simulated maximum
likelihood with 2{,}000 draws.
{We re-estimate that column with SAUSS by setting 1{,}600 epochs, mini-batches
of $m = 10$, five accepted simulator draws per observation,
Polyak--Ruppert averaging, and a $t^{-0.501}$ step-size decay. We
compare it with GHK simulated maximum likelihood.
{Because the published \texttt{cmmprobit} estimates in \texttt{Stata} use the alternative normalization \(L_{11}=\sqrt{2}\), we divide their coefficients and standard errors by \(\sqrt{2}\) to express them on the \(L_{11}=1\) scale used by \texttt{mlogit} and SAUSS (both in \texttt{R}). On this common scale, our 2{,}000-draw \texttt{mlogit} GHK coefficients differ from the published \texttt{cmmprobit} coefficients by at most \(0.02\) published standard errors.}}
\begin{table}[t]
\centering
\caption{Multinomial probit estimates, daycare scenario of
\citet{boneva2026}, Table~4, column~2; $N = 2{,}873$. Standard errors in
parentheses. Times are one fit on a single core.}
\label{tab:bgr}
\begin{tabular}{lcc}
\toprule
& GHK-SML & SAUSS \\
\midrule
Child skills & $\underset{(0.1684)}{0.4553}$ & $\underset{(0.1809)}{0.4641}$ \\
Family outcomes & $\underset{(0.1876)}{1.1614}$ & $\underset{(0.1993)}{1.2722}$ \\
Maternal earnings (1{,}000 EUR) & $\underset{(0.0014)}{0.0033}$ & $\underset{(0.0016)}{0.0039}$ \\
Family's opinion & $\underset{(0.0365)}{0.2689}$ & $\underset{(0.0378)}{0.2872}$ \\
Friends' opinion & $\underset{(0.0388)}{0.2616}$ & $\underset{(0.0399)}{0.2843}$ \\
\addlinespace[3pt]
$\Omega_{21}$ & $\underset{(0.0946)}{-0.0022}$ & $\underset{(0.1000)}{-0.0418}$ \\
$\Omega_{22}$ & $\underset{(0.1379)}{0.4494}$ & $\underset{(0.1764)}{0.5840}$ \\
\midrule
Computation time (s) & 4180 & 33 \\
\bottomrule
\end{tabular}
\end{table}
{Table~\ref{tab:bgr} reports the five
alternative-specific coefficients emphasized in the original paper,
together with the error covariance and computation time.}
{{The five reported coefficients are broadly
similar across the two methods at the stated computational budgets. Each
SAUSS estimate lies within one GHK standard error of the corresponding GHK
estimate, and the signs and 5\% significance classifications agree. The SAUSS
coefficients are somewhat larger, accompanied by a fitted $\Omega_{22}$ of
$0.58$, compared with $0.45$ for GHK. This covariance difference is less than
one standard error under either fit, although the continued downward movement
of the SAUSS estimate near the end of the run indicates some finite-run
sensitivity at the reported budget.}}
{The SAUSS standard errors use the finite-run plug-in variance
approximation in Remark~\ref{rem:finite-run-variance} and are computed from
the reported run in $0.7$ seconds. After $4.6$ million observation-level
gradient evaluations, the estimated algorithmic component accounts for
$0.16\%$ of the total variance in this calculation, so the estimated
variance is dominated by sampling uncertainty. Across ten random seeds,
the coefficient estimates vary by only three to five percent of one
reported standard error, indicating limited seed sensitivity for $\beta$
at this budget.} {One caution concerns estimating $H$ using
the outer product of simulated scores. For fixed $R$, this outer product
estimates $\Sigma_s$ plus a simulation-variance component of order
$1/R$. Using it in place of $H$ therefore understates the standard
errors: in this application, the resulting standard errors are, in median,
$0.56$ times the simulation-adjusted plug-in standard errors at $R=5$
and $0.96$ times them at $R=200$. After accounting for simulation
variance, the plug-in standard errors are essentially stable across these
values of $R$.}
{The methods differ sharply in cost under these reported configurations:
$33$ seconds for SAUSS compared with $4{,}180$ seconds for the GHK likelihood.}
Some of that gap is
implementation rather than method, and a practitioner content with $200$
draws would face a much smaller one, since the GHK fit then takes $163$
seconds and moves by at most $0.16$ standard errors. {Because
$J=3$, this application provides a comparison in a small-choice-set setting
where both methods are computationally feasible.}
\section{Conclusion}\label{sec:conclusion}
This paper develops SAUSS, a stochastic approximation approach to
likelihood-based estimation when unbiased simulated scores are available
but simulated probabilities are costly or biased inside logarithms. The
main argument is that the method-of-simulated-scores construction, which
was difficult to combine with deterministic full-sample optimization
because of discontinuity, is well suited to mini-batch stochastic
approximation. {The MNP implementation shows how the general
LDV score identity can be turned into an explicit descent algorithm by
differencing relative to alternative 1 and imposing the \(L_{11}=1\) scale
normalization, with
accept--reject diagnostics separating the ideal unbiased recursion from
capped finite computation.}
{The numerical results show that SAUSS achieves accuracy comparable to
GHK-based simulated maximum likelihood in the small-choice-set benchmark
and remains computationally feasible as the number of alternatives
increases. In the empirical application, SAUSS produces broadly similar
estimates with substantially lower computation time under the reported
implementations.}