Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
145,224 characters · 19 sections · 140 citation commands
Approximate Minimax Estimation of a Bounded Normal Mean via Stochastic Mirror Ascent
\doublespacing
We propose a numerical algorithm to find an approximately minimax estimator (under quadratic loss) of the mean of a standard normal distribution, when such mean is known to belong to the interval $[-m,m]$. casella_strawderman_1981 refer to this classical problem in statistical decision theory as the Bounded Normal Mean estimation problem. As noted by johnstone1992minimax, the minimax risk in this problem plays a central role in comparing the behavior of linear to nonlinear estimators in an important class of nonparametric estimation problems; see ibragimov1985nonparametric and donoho_liu_macgibbon_1990. The minimax risk of the Bounded Normal Mean problem has also been used recently by armstrong_kline_sun_2025 as an input to construct estimators that adapt to misspecification.
Our suggested algorithm finds an approximately minimax estimator by first obtaining an approximately least-favorable distribution. This strategy is common in the literature; see kempthorne1987numerical and chamberlain2000econometric. One point of departure relative to previous work is that we use a stochastic mirror ascent routine (with negative entropy as a mirror map) to approximately solve the following concave optimization problem: we maximize Bayes risk over the space of probability distributions supported on a fixed (equally-spaced) grid of $I$ points over the interval $[-m,m]$. In what follows, we let $\Delta^{I-1}$ denote the space of such distributions, without making explicit reference to the support points.
Before providing more specific details about our contribution, it is worth noting that our numerical strategy has two sources of approximation error, both of which are explicitly taken into account in our analysis. The first source of error is an obvious discretization error, which comes from the fact that we optimize Bayes risk over the simplex $\Delta^{I-1}$ and not over all probability distributions supported on $[-m,m]$. We show that we can control this discretization error thanks to a special (Lipschitz) continuity property of the risk function (see Lemma (ref) and Lemma (ref)). The second source of error is what we term optimization error, which comes from the fact that the mirror ascent routine for concave maximization is designed to find an approximate solution, and not an exact one. An explicit analysis of this type of error is common in theoretical computer science, operations research, and optimization where it is typical “to relax the requirement of finding an optimal solution, and instead settle for a solution that is good enough” shmoys2011design. We can control this optimization error thanks to the rich literature that has analyzed the theoretical properties of mirror descent (or ascent).
The methods of mirror descent nemirovski_yudin_1983 are a family of iterative procedures recommended in the optimization literature for approximately solving convex problems of high dimensions. These methods are iterative, first-order optimization algorithms, in that they typically require repeated evaluations of the objective function and an estimator of its (sub or super) gradient, but do not exploit any further smoothness information about the objective function.\footnote{The algorithm we recommend only requires repeated evaluations of an estimator of the supergradient of the objective function, but not of the objective function itself. See Algorithm (ref).} The optimization error associated with these algorithms comes from the fact that the procedure is designed to stop after finitely many iterations. Applications of the methods of mirror descent to other problems in statistics are given in juditsky2008 and fernandez2025approximate.
There are at least four reasons that explain why mirror ascent was the algorithm we chose for finding an approximately least-favorable distribution. First, even in the simple Bounded Normal Mean problem, it is difficult to get an explicit characterization of higher-order derivatives of the objective function defining the least-favorable distribution, which means that first-order methods (which only use information concerning the gradient and do not require further smoothness assumptions) are particularly attractive. Second, there is also a sense in which our algorithm is essentially the best first-order iterative algorithm for approximately minimizing (maximizing) a convex (concave), Lipschitz function over the simplex (as in the problem defining the approximately least-favorable distribution); see Proposition 4.2 in ben2001ordered.\footnote{It is also known that the rate of convergence of mirror descent (or ascent) improves over regular subgradient descent for convex minimization (maximization) problems in the probability simplex bubeck_2015.} Third, as our results show, the supergradient of the objective function at any distribution $\pi \in \Delta^{I-1}$ depends on the risk function $R(d_{\pi},\theta)$, where $d_\pi$ is the Bayes estimator associated with $\pi$. This means that even though evaluating the supergradient is potentially costly (as it requires numerical integration), we can rely on the well-known stochastic version of mirror ascent that replaces the supergradient by an unbiased estimator (and in our case, this unbiased estimator can be constructed using only one Monte Carlo draw per grid point per iteration, further reducing the computational burden of the algorithm). Fourth, since there is a rich theory regarding the convergence analysis of stochastic mirror descent (ascent) bubeck_2015, we can provide explicit theoretical guarantees about the performance of our suggested algorithm.\footnote{While it is entirely possible that off-the-shelf methods implemented in numeric computing platforms (for example, second-order optimization methods like sequential quadratic programming via fmincon in Matlab) perform as satisfactorily as our recommended algorithm, it is not clear to us how to provide theoretical guarantees without being more explicit about the tuning parameters used by the algorithm to generate approximate minimax solutions.}
The first theoretical result in this paper (Theorem (ref)) shows that our stochastic mirror ascent algorithm provably finds an approximately least-favorable distribution for the Bounded Normal Mean problem after finitely many iterations. The theoretical guarantees we provide are expressed in terms of a scalar parameter $\epsilon>0$. This parameter controls the overall approximation error of our numerical algorithm. Let $\lceil \cdot \rceil$ denote the ceiling function.\footnote{The ceiling function returns the smallest integer that is greater than or equal to a given number.} The theorem states that if the number of iterations, $J(\epsilon)$, and the number of points in the grid, $I(\epsilon)$, are chosen to be \[ J(\epsilon) := \left \lceil \frac{2 M^2 \ln (I(\epsilon))}{(\epsilon/2)^2} \right \rceil, \quad I(\epsilon):= \left \lceil 1 + \frac{3m^3 + 4m^2}{\epsilon/2} \right \rceil, \quad \textrm{ where }M:= 4m^2, \] then, for any $\alpha \in (0,1)$, with probability at least $1-\alpha$, the Bayes risk associated with the prior distribution $\pi^\epsilon$ found by the algorithm is close to the maximin value of the Bounded Normal Mean problem. And, in particular, it differs at most by an additive factor \[ \epsilon \left(1 + \sqrt{\frac{\ln(1/\alpha)} {\ln(I(\epsilon))}} \right).\] With the exception of $\epsilon$, our theorem makes explicit recommendations about all of the parameters needed to run the algorithm (such as the initial condition and the learning rate). It is worthwhile to make three remarks about our result. First, the probabilistic statement that appears in the theorem comes from the fact that our theoretical analysis acknowledges the fact that the gradient of the objective function needs to be estimated by Monte Carlo integration. This estimation induces---just as any other stochastic approximation method---an additional source of variation that will make our approximately least-favorable distribution change for different runs of the algorithm if the random seed is not fixed. Second, the total number of iterations (which scales logarithmically in the number of grid points) should be viewed as \emph{sufficient} to guarantee a good performance of the algorithm, but it is by no means \emph{necessary}. In fact, our simulation results show that the approximation error is likely to be of size $\epsilon$ for a considerably smaller number of iterations (which is important, given that $J(\epsilon)$ is already in the order of millions for values of $m$ as small as $m=3$). \emph{Third}, the convergence guarantee we provide in Theorem (ref) is expressed in terms of how well we can approximate the \emph{maximin} value of the Bounded Normal Mean problem, and not in terms of how close the approximate least-favorable distribution is relative to the true least-favorable distribution.\footnote{The least-favorable distribution is known to exist and be unique; see Proposition 4.19 in johnstone_2019. See also delatte26. }
The second theoretical result of this paper (Theorem (ref)) complements Theorem (ref) by showing that $\pi^{\epsilon}$ will indeed be close (in terms of the 1-Wasserstein distance) to the least-favorable distribution of the Bounded Normal Mean problem, provided $\epsilon$ is small enough. To the best of our knowledge, none of the other existing papers in the literature that have suggested algorithms for the Bounded Normal Mean problem---such as feldman1989minimax, donoho_liu_macgibbon_1990, and gourdin1994global---have derived results like ours. This is also true for other existing papers that have suggested algorithms for finding approximate solutions for more general minimax/maximin problems---such as kempthorne1987numerical, johnstone1992minimax, chamberlain2000econometric, montielolea_aradillasfernandez_blanchet_qiu_stoye_tan_2024, fernandez2025approximate, elliott2015nearly, and guggenberger2025numerical.
The third theoretical result of this paper (Theorem (ref)) shows that there is a sense in which the Bayes estimator associated with the prior distribution $\pi^\epsilon$ found by the stochastic mirror ascent algorithm is approximately minimax. The notion of approximate minimaxity that we use in this paper extends the classical definition of minimax sequence given in the seminal work of ghosh1964uniform, p. 1037. Broadly speaking, ghosh1964uniform defines a deterministic sequence of decision rules to be a minimax sequence if their worst-case risk converges to the minimax risk of the problem of interest. Since, as we have explained before, our stochastic mirror ascent algorithm introduces noise beyond the statistical uncertainty contained in the estimation problem, we need to allow for the possibility of a stochastic sequence of decision rules that approximate the minimax value of interest. In light of this observation, we define a sequence of decision rules to be a stochastic minimax sequence if their worst-case risk converges in probability to the minimax value of interest. This means that the tail of a stochastic minimax sequence contains decision rules that are $\epsilon$-minimax---in the sense of ferguson_1967, Chapter 1, Definition 4, p. 33---with probability at least $1-\delta$, for any $\epsilon,\delta>0$.
The last theoretical result of the paper (Theorem (ref)) shows that the Bayes estimator associated with the prior distribution $\pi^\epsilon$ approximates the exact minimax estimator of a bounded normal mean in a stronger sense: for any compact set $K$ of data realizations, there is a small enough $\epsilon$ (that could depend on the set) such that the approximate and the exact minimax estimators are uniformly close to each other over $K$, with probability approaching one. This means that our approximate minimax estimator converges---uniformly over compact sets---to the minimax estimator of the Bounded Normal Mean problem.
{ Simulations:} It is known that in the Bounded Normal Mean problem the worst-case risk of the best minimax linear estimator is at most 25% above the best, unconstrained, minimax estimator. This is the so-called ibragimov1985nonparametric constant, see Theorem 4.7 in johnstone_2019. The minimax linear estimator thus provides a computationally feasible (and trivial to compute) benchmark for the Bounded Normal Mean problem. The value of the ibragimov1985nonparametric constant implies that the largest reduction in risk that we could expect to see from our approximately minimax estimator (relative to the minimax linear estimator) is 20%. Figures (ref) and (ref) show that our stochastic mirror ascent routine generates risk improvements over the minimax linear estimators that range from about 6% to 18%, as $m$ ranges from 1 to 4 (with the largest improvement occurring at $m=1.6$). In our simulations we choose $\epsilon$ to be $\epsilon(m) := (1/5)m^2/(1+m^2)$. We provide a rationale for this choice in Section (ref).
In terms of the computational effort required to achieve the 6% to 18% improvement reported above, we note that the largest reduction in risk that we see takes about 45 seconds on a personal computer. The computational time increases with $m$, and implementing our algorithm for $m=3$ takes about 2 hours (see Figure (ref)). The main bottleneck we face in the implementation of our algorithm is the number of iterations recommended in Theorem (ref). We report simulations where we bound the wall-clock time of our procedure to 5 and 15 minutes. With the 15-minute limit, we still see improvements between 6% and 18%, but the improvement over linear estimators when $m=3$ drops from about 14% without constraints, to almost 6%.
{ Application:} A common problem in macroeconometrics is the estimation of the dynamic causal effects---or impulse response functions---of aggregate shocks to macroeconomic policies or economic fundamentals. The two most common empirical methods currently used by applied macroeconomists for doing so are structural vector autoregressions (henceforth interchangeably referred to as VARs or SVARs), going back to the seminal work of sims_1980, and the so-called local projections (LPs), introduced by jorda_2005. It is now understood that there is a sense in which VARs and LPs share the same estimand plagborgmoller_wolf_2021,xu2026local, but that in finite samples these procedures lie on opposite ends of a bias-variance trade-off li2024local: LP estimators have lower bias than VAR estimators, but they also have substantially higher variance at intermediate and long horizons.
The question of interest is how to best aggregate the information in LP and VAR estimators to estimate a particular impulse response coefficient. We follow closely the work of armstrong_kline_sun_2025 and consider the problem of constructing a minimax estimator based on quadratic loss, and a bivariate Gaussian model for the VAR and LP estimators.\footnote{In the bivariate Gaussian model, the LP estimator is unbiased but has potentially high variance, and the VAR estimator is biased but with smaller variance than the LP estimator; see Equation (ref).} We impose an a priori bound on the bias of the VAR. As we explain in Section (ref), both the bivariate normality and the bound on the bias can be motivated using the locally misspecified vector autoregression model used in montielolea_plagborgmoller_qian_wolf_forthcoming. In particular, the bound on the bias can be obtained from restricting the fraction of the residual variance that is explained by the lagged shocks in the misspecified part of the VAR model; see montielolea_qian_wolf_plagborgmoller_2025, p. 16.
Under a convenient invariance assumption, the minimax invariant estimator is a debiased version of the VAR estimator. The optimal debiasing function depends on the data only through the standardized difference between the VAR and LP estimators, which can be interpreted as a hausman1978specification statistic for misspecification and which we denote as $\Delta$. As shown in armstrong_kline_sun_2025, the optimal debiasing function turns out to correspond to the minimax estimator of a Bounded Normal Mean problem, where we use $\Delta$ to try to estimate the magnitude of the bias of the VAR estimator.
We report our estimator for four influential macroeconomic papers, where we allow the fraction of the residual variance that is explained by the lagged shocks in the misspecified part of the VAR model to be at most 1%. The four papers we consider are gertler_karadi_2015, focusing on the impulse response of log industrial production to a monetary policy shock; romer_romer_2010, focusing on the impulse response of log of real GDP per capita to a tax shock; ramey_2011, focusing on the impulse response of log of real GDP per capita to a government spending shock (based on military news); and francis_owyang_roush_dicecio_2014, focusing on the impulse response of log of labor productivity to a technology shock. Even though the amount of misspecification that we allow is relatively small (which should favor the VAR estimator), in all applications the approximately minimax estimator is very close to the LP estimator. This is also consistent with the minimax linear estimator reported by montielolea_plagborgmoller_qian_wolf_forthcoming (Corollary 4.2), which under our assumptions, places about 70% of the weight on the LP estimator. In our applications, for medium and long horizons the shrinkage towards the LP estimator is a bit more aggressive than the one implied by the linear estimator (and the gains in excess risk relative to an oracle that knows the bias are about 15%).
The remainder of the paper is organized as follows. Section (ref) introduces the Bounded Normal Mean problem and its minimax formulation. Section (ref) presents the algorithm and the main theoretical results. Section (ref) reports simulation results. Section (ref) applies the method to estimation of impulse responses.
Fix $m > 0$ and consider the statistical model \[ Y \sim \mathcal N(\theta, 1), \quad \theta \in [-m,m].\] The problem of interest is to estimate $\theta$ under the loss function: \[\mathcal{L} (a, \theta) = (a- \theta)^2.\] casella_strawderman_1981 refer to this problem as the bounded normal mean estimation problem. An excellent and detailed treatment of the problem is given in Chapter 4.6 of johnstone_2019.
A large literature in statistics has studied the theoretical properties of different estimators for the Bounded Normal Mean problem; see, for example, bickel1981minimax, levit1981asymptotic, donoho_liu_macgibbon_1990, and the references therein. As we will explain below, an important difference of our paper relative to existing literature is that we make it our goal to suggest a concrete numerical strategy---based on the methods of stochastic mirror descent (ascent)---to obtain an approximate minimax estimator (in a sense we make precise). The main theoretical results in our paper will be focused on showing that our numerical strategy indeed achieves the desired objective.
We note that we have deliberately set the variance of the outcome $Y$ to equal one. As we discuss later, this is without loss of generality as the Bounded Normal Mean problem treats the variance as known. See also the scale-invariance relation in Equation 4.35, Chapter 4.6 of johnstone_2019. It is also without loss of generality to restrict the action space to be $[-m,m]$. The reason is that, under squared-error loss, truncating any real-valued estimator to $[-m,m]$ weakly reduces risk pointwise in $\theta$.\footnote{One benchmark is the projection of $Y$ onto $[-m,m]$, $d_{\textrm{trunc}}(y) = \min \{\max \{y,-m\}, m\},$ which is measurable, valued in $[-m,m]$, and, by truncation, weakly dominates the unconstrained estimator $d(y) = y$ under squared loss as $(d_{\textrm{trunc}}(y) - \theta)^2 \le (y-\theta)^2$ for all $\theta \in \Theta$ pointwise.}
Let $\mathcal{D}$ be the class of all measurable decision rules $d: \mathbb{R} \to [-m,m]$. As usual, the risk of a rule $d$ is \[R (d, \theta) := \mathbb{E}_{\theta} [(d(Y)- \theta)^2],\] where the expectation is taken assuming $Y \sim \mathcal{N}(\theta,1)$. The first result of this paper establishes an important property of the risk function that will be exploited by our numerical strategy; namely, that $R(d,\cdot)$ is Lipschitz continuous in $\theta$ for any $d \in \mathcal{D}$.
This first lemma shows that the Lipschitz constant is independent of $d$, but dependent on the bound $m$. To the best of our knowledge this result has not appeared previously in the literature. The key argument we use to establish this result is an upper bound on the total variation distance between two normal distributions with unit variance; see Theorem 1.3 in devroye2018total.\footnote{For normal distributions with common covariance, Theorem 1 in barsov1987proximity implies the sharper bound: \[ \frac{1}{2} \int_{-\infty}^{\infty} \left| \phi(y;\theta) - \phi(y;\theta') \right| dy \le \frac{ |\theta-\theta'| }{\sqrt{2\pi}}.\] Using this bound in Appendix (ref) would give the sharper Lipschitz constant equal to $3 \sqrt{\frac{2}{\pi}} m^2 + 4 m$.}
Given $m>0$, we say that an estimator $d^*$ is minimax---or $m$-minimax, adopting the terminology of armstrong_kline_sun_2025---if it satisfies
We refer to $\bar{v}^*(m)$ as the minimax value of the bounded normal mean estimation problem associated with the bound $m>0$.
It is known that when the bound $m$ is sufficiently small, there is a closed-form solution for the problem in (ref); see Theorems 3.1 and 4.2 in casella_strawderman_1981. More generally, the problem of finding a minimax estimator can be viewed as an infinite dimensional convex programming problem; see pp. 1419 and 1420 in donoho_liu_macgibbon_1990. This means that finding a minimax estimator must necessarily involve some type of numerical approximation.
To be explicit about the fact that any given numerical strategy will only approximately solve the minimax problem in (ref), we build on the seminal work of ghosh1964uniform and define a stochastic minimax sequence as follows. Let $\{P_n\}_{n=1}^{\infty}$ be a sequence of probability distributions over the space $\mathcal{D}$ and let $\{d^*_{n}\}_{n=1}^{\infty}$ be a stochastic sequence of decision rules in which each decision rule $d^*_n$ is drawn from $P_n$.
Definition (ref) extends the classical definition of minimax sequence given in ghosh1964uniform, p. 1037. ghosh1964uniform defines a deterministic sequence of decision rules $\{d_n\}_{n=1}^{\infty} \subseteq \mathcal{D}$ to be a minimax sequence if
Since the numerical procedure that we will recommend (which is akin to stochastic gradient descent for convex minimization problems) introduces noise beyond the statistical uncertainty contained in the estimation problem, we need to allow for the possibility that $d_n$ is stochastic. To this end, Definition (ref) replaces the $``\rightarrow''$ relation in (ref) by convergence in probability, where convergence in probability is assessed relative to the distributions $P_n$.
Definition (ref) is also related to the notion of $\epsilon$-minimax decision rules, which is a common way of formalizing the notion of an approximate minimax solution in statistical decision problems; see ferguson_1967, Chapter 1, Definition 4, p. 33.\footnote{See also montielolea_aradillasfernandez_blanchet_qiu_stoye_tan_2024 for recent applications and an algorithmic approach for finding $\epsilon$-minimax decision rules in a class of statistical decision problems.} Note that if $\{d^*_n\}_{n=1}^{\infty}$ is a stochastic minimax sequence, then for every pair $\epsilon,\delta > 0$ there exists $n(\epsilon,\delta)$ such that whenever $n \geq n(\epsilon,\delta)$ \[ P_n \left( \bar{v}^*(m) \leq \sup_{\theta \in [-m,m]} R(d^*_n,\theta) \leq \bar{v}^*(m) + \epsilon \right) \geq 1-\delta. \] This means that the tail of a stochastic minimax sequence contains decision rules that are $\epsilon$-minimax with probability at least $1-\delta$.
Our strategy for constructing a stochastic minimax sequence for the problem in (ref) consists of finding a sequence of priors that approximate the problem's least-favorable distribution. We do this by discretizing the parameter space. The strategy of discretizing the parameter space to find an approximate least-favorable distribution has been recommended before in the econometrics literature by chamberlain2000econometric and it has been used recently by armstrong_kline_sun_2025. To the best of our knowledge, neither paper provides theoretical results that guarantee that such a strategy indeed succeeds in finding an approximate minimax rule (or even an approximate least-favorable distribution). The closest result that we are aware of is Lemma 1 of guggenberger2025numerical, which provides a set of sufficient conditions that guarantee that the value of a discretized version of a minimax problem converges to the minimax value of the original problem.
For the sake of exposition, this section presents a quick overview of the typical relation between a minimax decision rule and a least-favorable distribution. We also present the discretized Bounded Normal Mean problem and the definition of an approximate least-favorable distribution for this problem.
Let $\tau$ denote the standard (subspace) topology over $[-m,m]$. It is known that the topological space $([-m,m],\tau)$ is a Polish space.\footnote{This follows from the fact that $\tau$ is metrizable by the absolute value metric, and the space $([-m,m],|\cdot|)$ is a complete and separable metric space.} In order to define probability distributions over $[-m,m]$, we endow this space with the $\sigma$-algebra generated by $\tau$; i.e., its Borel $\sigma$-algebra. Let $\Delta([-m,m])$ denote the set of all (Borel) probability measures over $[-m,m]$.
Let $\pi$ be an element of $\Delta([-m,m])$. The Bayes risk of $d$ with respect to $\pi$ is defined as \[r(d, \pi) := \mathbb{E}_{\pi}[\mathbb{E}_{\theta} (d(Y)- \theta)^2] = \int_{[-m,m]} R(d,\theta) d\pi.\] The rule $d_{\pi}$ is said to be a Bayes rule with respect to $\pi$ if it satisfies
Define a least-favorable distribution as a distribution $\pi^*$ that satisfies
In Appendix (ref) we show that the maximin problem that defines the least-favorable distribution can be viewed as a concave maximization problem over the infinite-dimensional space $\Delta([-m,m])$. The objective function in this problem---which corresponds to the function $f(\pi) = \inf_{d \in \mathcal{D}} r(d,\pi)$---is Lipschitz continuous in $\pi$, with respect to the 1-Wasserstein metric over $\Delta([-m,m])$.
It is known that
See Proposition 4.19 in johnstone_2019 and also their Theorem A.5.\footnote{There are two additional interesting properties that we will reference later. First, the least-favorable distribution is unique. Second, the least-favorable distribution is supported on finitely many points (and all support points have the same risk). See, again, Proposition 4.19 in johnstone_2019.} It is also known that if (ref) holds, any minimax rule $d^*$ must be a Bayes rule with respect to the least-favorable distribution $\pi^*$.
Discretized Maximin Problem: As mentioned before, our numerical strategy to construct a stochastic minimax sequence starts by discretizing the parameter space $\Theta := [-m, m]$ using an equally-spaced grid of $I>1$ points $\{\theta_1,\dots,\theta_{I}\}$ in the interval $[-m,m]$. This means that the $i$-th point in the grid equals: \[ \theta_i: = -m + \frac{2m(i-1)}{I-1}. \] The gap between any consecutive points $\theta_{i+1}, \theta_i$ is
Define the $(I-1)$-dimensional simplex $\Delta^{I-1} := \left \{\pi \in [0,1]^I: \sum_i \pi_i = 1 \right\}$ and then consider the discretized version of the maximin problem in (ref)
There is an evident gap between the value of the original maximin problem in (ref)---which we denoted as $\underline{v}^*(m)$---and the value of the discretized maximin problem in (ref)---which we are denoting as $\underline{v}(m;I)$. Obviously, $\underline{v}(m;I) \leq \underline{v}^*(m)$, since the feasible set of the outer maximization in the discretized maximin problem is a subset of $\Delta([-m,m])$. The second result of this paper establishes a concrete quantitative relation between the maximin values of the Bounded Normal Mean problem and its discretized version.
The key argument to establish Lemma (ref) is the Lipschitz continuity of the risk function (which was established in Lemma (ref)) and the minimax theorem applied to both the discretized and the original problem. As expected, the gap between the discretized maximin problem and the original problem decreases as a function of $I$, and larger values of $m$ will require a finer grid to attain a good approximation. We note that Lemma (ref) readily yields a formula to select $I$ in order to guarantee that the gap between the two maximin problems is small. For example, to guarantee a gap of at most $\epsilon/2$ it suffices to set the number of points in the uniform grid equal to
The function $\lceil \cdot \rceil$ is the ceiling function: the function that returns the smallest integer that is greater than or equal to a given number. In light of Lemma (ref) our numerical strategy will focus on obtaining an approximate solution to the discretized maximin problem. The following definition is entirely analogous to the notion of $\epsilon$-minimax decision rules---see ferguson_1967, Chapter 1, Definition 4, p. 33---but applied to the maximin problem.
Let $I(\epsilon)$ be defined as in (ref). Note that an important implication of Lemma (ref) is that if $\pi^{(\epsilon/2,I(\epsilon))}$ is an $(\epsilon/2)$-least-favorable distribution for the discretized problem with $I(\epsilon)$ equally-spaced points, then $\pi^{(\epsilon/2,I(\epsilon))}$ is an $\epsilon$-least-favorable distribution of the original Bounded Normal Mean problem. Our strategy to construct a stochastic minimax sequence for the Bounded Normal Mean problem relies on this observation. More concretely, we will consider a sequence $\{(\epsilon_n,I_n)\}_{n=1}^{\infty}$, where $\epsilon_n \rightarrow0$ and $I_n \rightarrow \infty$. We will show that we can build an associated stochastic sequence of priors $\pi_n \equiv \pi^{\epsilon_n}$ that are guaranteed to be---with high probability---approximately the least-favorable distribution for the original problem. In the next section we will show how to obtain such a stochastic sequence of priors. Our main result will be that the sequence of Bayes decision rules corresponding to the sequence of priors $\{\pi_n\}_{n=1}^{\infty}$ is a stochastic minimax sequence.
Note also that without an upper bound on $\epsilon$, any distribution $\pi \in \Delta^{I-1}$ is approximately least-favorable. The risk function in the Bounded Normal Mean problem is bounded by $4m^2$, so we are only interested in considering approximation errors below this upper bound.
This section presents our main results. First, we show that a stochastic mirror ascent algorithm for concave maximization over the $(I(\epsilon)-1)$-dimensional simplex can be used to provably obtain---with high probability---an approximate least-favorable distribution for the maximin problem in (ref). Second, we show that the Bayes rules associated with the sequence of approximately least-favorable priors obtained via stochastic mirror ascent constitute a stochastic minimax sequence.
For any $\pi \in \Delta([-m,m])$, define the function $f: \Delta([-m,m]) \rightarrow \mathbb{R}_{+}$ given by
The first key observation here is that the function $f(\cdot)$ is concave in $\pi$ over $\Delta([-m,m])$, and also over $\Delta^{I-1}$ for any $I > 1$. This result has been established more generally before; see, for example, chamberlain2000econometric, Equation 5, p. 630, and the discussion therein. For the sake of completeness, we formalize this observation in Lemma (ref) in Appendix (ref). Lemma (ref) also shows that for any prior $\pi \in \Delta^{I-1}$, the vector $g(\pi):=(R(d_{\pi}, \theta_1), \ldots, R(d_{\pi},\theta_{I}))$ is a supergradient of $f(\cdot)$ at $\pi$ (where $d_{\pi}$ is the Bayes rule corresponding to $\pi$). Lastly, since $|R(d, \theta)| \leq 4m^2$ for all $(d,\theta) \in \mathcal{D} \times \Theta$, the function $f:\Delta^{I-1} \rightarrow \mathbb{R}$ is Lipschitz (with respect to the $L_1$-norm) with Lipschitz constant $M:=4m^2$.
A popular algorithmic approach for approximately solving concave maximization problems---in particular, those of high dimensions---is to use the methods of (stochastic) mirror descent of nemirovski_yudin_1983. The methods of mirror descent are a family of iterative procedures designed for finding an approximate solution of convex minimization (or concave maximization) problems in high dimensions. These methods require repeated (but finitely many) evaluations of the subgradient of the objective function or a noisy but unbiased estimator of the subgradient. These methods exploit the geometry of the optimization domain, but they do not exploit any higher-order smoothness information about the optimization problem. Mirror descent is known to outperform gradient descent in some optimization domains, such as the probability simplex; see Section 4.3 of bubeck_2015. See also Section 5.1.1 in nemirovski2009robust.
We consider the iterative algorithm described below for maximizing $f(\cdot)$ over the simplex $\Delta^{I-1}$. The stochastic component of this algorithm comes from the evaluation of the supergradient. The pseudo-code that describes Algorithm (ref) below takes the grid $\{\theta_i\}_{i=1}^{I}$ as given. Another input is the step size (denoted as $\eta$), which controls the adjustment of the prior at each epoch. We also take the number of epochs ($J$) as given.
The next theorem shows that---with high probability---the stochastic mirror ascent routine outputs an approximate least-favorable distribution for the Bounded Normal Mean problem.
Theorem (ref) shows that our suggested stochastic mirror ascent algorithm---with the choice of tuning parameters in (ref)---succeeds in finding an approximate least-favorable distribution for the maximin problem (ref). Since Algorithm (ref) relies on an estimator of the supergradient (see Step 7), it is important to quantify this uncertainty in our theoretical guarantees. The theorem provides an explicit description of the effects of the stochastic gradient on the optimization procedure. In particular, the theorem says that with probability at least $1-\alpha$, the Bayes risk of $\pi^\epsilon$ is at most within the additive factor \[ \epsilon \left(1 + \sqrt{\frac{\ln(1/\alpha)} {\ln(I(\epsilon))}} \right)\] of the maximin value of interest. Note that this term is increasing in $\epsilon$. However, a smaller value of $\epsilon$ increases the computational complexity of the procedure because it increases the dimension of the concave optimization problem via its effect on $I(\epsilon)$.
To complement Theorem (ref) we also show that there is a sense in which $\pi^{\epsilon}$ is close to the (unique) least-favorable distribution of the original maximin problem, which we have denoted as $\pi^*$. To show this, take an arbitrary sequence $\{\epsilon_n\}_{n=1}^{\infty} \subseteq \mathbb{R}_{++}$ such that $\epsilon_n \rightarrow 0$. For each $\epsilon_n$ in the sequence, we obtain the corresponding approximate least-favorable prior $\pi_n : = \pi^{\epsilon_n}$ as described in Theorem (ref). We show that $\pi_n$ converges (in 1-Wasserstein distance, as defined below) to the least-favorable distribution $\pi^*$.
The 1-Wasserstein Distance. For arbitrary Borel probability distributions $\pi, \pi' \in \Delta([-m,m])$ define $\Gamma(\pi,\pi')$ to be the set of all couplings of $(\pi,\pi')$. Define the 1-Wasserstein distance (with cost function $|\cdot|$) as \[ W(\pi,\pi') := \inf_{\gamma \in \Gamma(\pi,\pi')} \int_{[-m,m]^2} |\theta-\theta'|d\gamma. \] It is known that $W(\pi,\pi')$ defines a metric over all Borel probability measures on $\Delta([-m,m])$; see Theorem 7.3 in villani2021topics. Moreover, it is common to refer to $W(\cdot,\cdot)$ as the Kantorovich-Rubinstein distance. It is also known that the Kantorovich-Rubinstein duality formula implies that
where \[ \| \varphi \|_{\textrm{Lip}} := \sup_{\theta \neq \theta'} \frac{|\varphi(\theta)-\varphi(\theta')|}{|\theta-\theta'|}. \] See Remark 7.5, p. 207 in villani2021topics.
Theorems (ref) and (ref) have shown that the stochastic mirror ascent algorithm can be used to generate an approximate least-favorable prior. In this section, we use these results to construct a stochastic minimax sequence.
Theorems (ref) and (ref) translate the convergence of the approximate least-favorable distribution into convergence of the approximately minimax estimator. Theorem (ref) shows that the sequence of Bayes rules corresponding to the approximately least-favorable distribution forms a bona fide stochastic minimax sequence. The stochastic convergence reflects the randomness from the Monte Carlo estimation of supergradients. Theorem (ref) complements Theorem (ref) by showing that the sequence of approximately minimax rules converges (uniformly over compacts) with high probability to the problem's minimax rule.
This section explores the numerical performance of the stochastic mirror ascent algorithm in the Bounded Normal Mean problem. We compare the risk of the approximately minimax estimator to the linear minimax benchmark. We also assess the computational cost of the suggested procedure.
As discussed in the introduction---with the exception of $\epsilon$---our Theorem (ref) makes explicit recommendations about all of the parameters that need to be specified in order to run the stochastic mirror ascent algorithm. We now argue that in the Bounded Normal Mean problem, there is a simple way to choose $\epsilon$ (as a function of $m$).
Note first that the minimax risk over linear estimators in the Bounded Normal Mean problem is
see Equation 4.28 in johnstone_2019. It is also known that
where the upper bound is due to donoho_liu_macgibbon_1990; see also Theorem 4.17 in johnstone_2019.\footnote{The factor 1.25 is typically referred to as the ibragimov1985nonparametric constant.}
If $\epsilon$ were chosen to satisfy $\epsilon > 0.25 \bar{v}^*(m)$, then the minimax linear estimator would be $\epsilon$-minimax: $\bar{v}_{\mathrm{linear}}(m) \leq 1.25 \bar{v}^*(m) < \bar{v}^*(m) + \epsilon$. We are interested in values of $\epsilon$ associated with smaller approximation errors (to avoid having provable $\epsilon$-minimax guarantees for the minimax linear estimator). Thus, using (ref) we recommend choosing $\epsilon$ as
Figure (ref) displays both the required number of iterations of the stochastic mirror ascent algorithm and the implied grid size when we choose $\epsilon(m)$ as in (ref) and $m \in [0.1,9]$. The figure illustrates that both the number of iterations and the grid size grow (polynomially) with $m$. The dominant cost is the number of iterations, since $J(m):=J(\epsilon(m))$ grows as $m^4\ln(m)$ while the number of grid points $I(m):=I(\epsilon(m))$ grows as $m^3$. For example, at $m=1$, $I(m) = 141$ and $J(m) \approx 63,345$ (yielding a total time of about 3 seconds for the algorithm to run); while at $m=3$, $I(m) = 1,301$ and $J(m) \approx 2.3$ million (which makes the runtime of the algorithm about 2 hours).
Figure (ref) illustrates the total time required to run the stochastic mirror ascent routine and evaluate the risk with 1,000 MC draws and the values of $\epsilon(m)$ for $m \in \{1, \ldots, 3\}$, implied by the formula in (ref).
Figure (ref) compares the risk of the minimax estimator obtained via stochastic mirror ascent and the linear minimax estimator for $m \in \{1,\ldots,4\}$. The upper bound (gray dashed line) is $\bar{v}_{\mathrm{linear}}(m)=m^2/(1+m^2)$, the worst-case risk of the minimax linear estimator. The lower bound (orange dashed line) is $0.8 \bar{v}_{\mathrm{linear}}(m)$, obtained from the ibragimov1985nonparametric constant. The difference between these upper and lower bounds equals $(1/5)m^2/(1+m^2)$, which is our choice of $\epsilon(m)$.
We obtain the worst-case risk of the estimator generated by the stochastic mirror ascent algorithm (red line, with circle markers) in two steps. In the first step, we evaluate the risk function of the approximately minimax estimator over a finite grid of 200,000 equally-spaced points in the interval $[-m,m]$. The risk function is evaluated by Monte Carlo simulation using $1,000$ draws.\footnote{We simulate $N=1,000$ independent draws $Z_{j} \sim N(0,1)$ and compute, for each grid point $\theta_i$,
} In the second step, we acknowledge the fact that the worst-case risk evaluated on the finite grid understates the estimator's worst-case risk over the parameter space (since the worst-case risk may be attained between grid points). In Appendix (ref) we show how the Lipschitz continuity of the risk function in Lemma (ref) can be used to translate the evaluation of the risk over a finite grid into a certifiable upper bound for the worst-case risk of the approximately minimax estimator---and, consequently, for the problem's minimax value. The red line (with circle markers) in Figure (ref) is such an upper bound.
The blue line (with circle markers) is the average risk---computed according to the approximate least-favorable distribution---of the estimator generated by the stochastic mirror ascent algorithm. Since, by construction, the approximately minimax estimator is a Bayes decision rule with respect to the approximately least-favorable prior, the blue line (with circle markers) is a legitimate lower bound for both the problem's minimax value and for the estimator's worst-case risk.
Since the worst-case risk of $d_{\pi^\epsilon}$---the estimator generated by the stochastic mirror ascent algorithm---lies in between the red and blue lines (both with circle markers), we conclude that $d_{\pi^\epsilon}$ is $\epsilon$-minimax optimal for a value of $\epsilon$ considerably smaller (the difference between the lines with circle markers) than the target $\epsilon(m)$ (the difference between the dashed lines). Figure (ref) is thus graphical evidence of the success of our stochastic mirror ascent algorithm.
To further illustrate the success of our stochastic mirror ascent algorithm, Figure (ref) reports a lower bound on the percentage improvement in the worst-case risk of $d_{\pi^\epsilon}$ relative to the worst-case risk of the minimax linear estimator. That is, we would like to report the risk improvement computed as \[100 \times \left(1 - \frac{\sup_{\theta \in [-m,m]} R(d_{\pi^\epsilon}, \theta)}{\bar{v}_{\mathrm{linear}}(m)} \right),\] but we replace $\sup_{\theta \in [-m,m]} R(d_{\pi^\epsilon}, \theta)$ by the certified upper bound reported in Figure (ref).
For comparison, in Appendix (ref) we also report the worst-case risk of two alternative estimators: the truncated maximum likelihood estimator $d_{\textrm{trunc}} (y) = \max\{\min\{y,m\}, -m\}$ and the clipped linear estimator (purple line) $d_{\textrm{linear, clipped}} = \max\{\min\{c \cdot y,m\}, -m\}$, where $c=m^2/(1+m^2)$.
Figure (ref) plots the approximately minimax estimator and the minimax linear estimator for $m \in \{1,2,4\}$. Note that for $m=2$ and $m=4$, the most visible differences between $d_{\pi^\epsilon}$ and the minimax linear estimator concentrate near the boundaries.
Suppose that we are only willing to run the algorithm for at most $x$ minutes. That is, we will run the algorithm either until the required number of iterations is reached or the time budget is exhausted, whichever happens first. These types of simulations are typically referred to as fixed wall-clock time simulations.\footnote{See Section 7 of JMLR:v24:21-0377 for an example of numerical simulations that also fix the wall-clock time. }
Figure (ref) displays the number of iterations completed within a five-minute limit for each $m \in \{1, \dots, 4\}$, versus the theoretical target implied by $\epsilon(m) = \frac{1}{5} \frac{m^2}{1+m^2}$. As illustrated by this figure, the five-minute time limit starts binding at $m=2.2$.
Figure (ref) displays different risk values: analytical bounds ($\bar{v}$ and $\frac{1}{1.25}\bar{v}$), and the grid-evaluated worst-case risk of the stochastic mirror ascent algorithm.
Figure (ref) illustrates the percentage improvements in risk, analogous to Figure (ref), but under time restrictions: a five-minute limit (solid red line) and a fifteen-minute limit (dashed red line). Even with the five-minute limit binding for $m > 2.2$, the proposed algorithm still delivers a meaningful improvement in the worst-case risk up to $m \approx 3$. Beyond this point, the five-minute constraint is no longer sufficient to certify any improvement, and the SMA-based decision rule is dominated by the linear minimax rule. Increasing the limit to fifteen minutes pushes the break-even point out to $m \approx 3.3$ and recovers additional risk improvement throughout the upper range of $m$. However, at large $m$ the problem becomes hard enough that this time is not sufficient to beat the linear rule.
We conclude this section by making a brief comment about adaptive stochastic gradient methods---see leijordan2020 for a recent discussion. Our interest is centered on discussing potential applications of adaptive algorithms to the Bounded Normal Mean problem. As mentioned earlier, the algorithm we developed in the paper only exploits the fact that the objective function defining the least-favorable distribution is concave and Lipschitz. More concretely, our algorithm used a simple characterization of the supergradient and a bound on the Lipschitz constant ($M:=4m^2$). However, it is entirely possible that the objective function has additional properties that could be exploited by an iterative optimization routine; for instance, the function could be strongly concave, or it could possess well-behaved higher-order derivatives, or it could have a tighter modulus of Lipschitz continuity.
For strong concavity, it would be interesting to study the extent to which the algorithms in leijordan2020---which adapt to an unknown modulus of strong convexity---could be used in our problem.
For adaptation to the problem's Lipschitz constant, one could explore the possibility of using the algorithms in duchi2011, in which the step sizes are computed from the observed gradients. Interestingly, these algorithms could yield theoretical guarantees that are comparable to the ones obtained with the optimal parameters chosen a priori. This type of adaptation seems particularly promising in our problem, as fewer iterations could suffice to achieve the same accuracy. In Appendix (ref) we present some numerical evidence that the improvement is non-negligible. To show this, we run Algorithm (ref) for $m \in \{1, 1.6, 2\}$ and record the magnitude of the stochastic supergradient at every iteration. This exercise shows that replacing the worst-case Lipschitz bound with its empirical counterpart would reduce the number of iterations by factors of approximately 2, 3.5, and 5 (see Figure (ref) in Appendix (ref)), with the improvement increasing in $m$. We note, however, that extending the high-probability guarantees in Theorem (ref) to data-dependent step sizes would require a different type of concentration argument, which is beyond the scope of this paper.
A common problem in macroeconometrics is the estimation of the dynamic causal effects---or impulse response functions---of aggregate shocks to macroeconomic policies or economic fundamentals. The two most common empirical methods currently used by applied macroeconomists for doing so are structural vector autoregressions (henceforth interchangeably referred to as VARs or SVARs), going back to the seminal work of sims_1980, and the so-called local projections (LPs), introduced by jorda_2005.
Understanding the properties of these two estimation procedures has been an active area of recent research; see the references and discussion in montielolea_plagborgmoller_2021, montielolea_qian_wolf_plagborgmoller_2025, and jorda2025local. It is now understood that there is a sense in which VARs and LPs share the same estimand plagborgmoller_wolf_2021,xu2026local, but that in finite samples these procedures lie on opposite ends of a bias-variance trade-off li2024local: LP estimators have lower bias than VAR estimators, but they also have substantially higher variance at intermediate and long horizons.
Let $\theta_h$ denote the impulse response coefficient of interest at horizon $h$. For the sake of notation, we omit any explicit reference to the horizon of interest in the estimators and the covariance matrix (although both of these clearly depend on $h$). We consider a set-up where an applied macroeconomist has access to both LP and VAR estimators of $\theta_h$. We assume that these estimators follow the simple parametric model
We treat the matrix $\Sigma$ as known, and assume that the variance of the LP estimator $(\sigma^2_{\textrm{LP}})$ is strictly larger than the variance of the VAR estimator $(\sigma^2_{\textrm{VAR}})$. The unknown parameter $b$ is used to represent the potential bias of the VAR estimator.
While the exact bivariate normality of the estimators is unlikely to hold in finite samples, the asymptotic normality of the estimators can be established by using the locally misspecified vector autoregression model used in montielolea_plagborgmoller_qian_wolf_forthcoming.\footnote{montielolea_plagborgmoller_qian_wolf_forthcoming consider stationary data generating processes that are well approximated by a finite-order SVAR model, but subject to local misspecification in the form of an asymptotically vanishing moving average process of potentially infinite order.} In particular, their Proposition 3.1 provides an asymptotic expansion that shows that the large-sample distribution of the LP estimator is unaffected by local misspecification (and so, there is a sense in which LPs are asymptotically unbiased). In contrast, their Proposition 3.2 shows that local misspecification typically leads to asymptotic bias in the VAR estimator. This bias is captured by the scalar parameter $b$ in (ref), which will depend on the magnitude of the local misspecification. When the local misspecification is of order $O_{P}(1/\sqrt{T})$, a martingale central limit theorem can then be applied to show that the right-hand side of (ref) provides a valid asymptotic approximation to the finite-sample distribution of the VAR and LP estimators under the locally misspecified VAR model.
We make two additional remarks about the structure of the bivariate normal distribution in (ref). First, when the bias parameter $b$ equals zero, the estimator $\alpha \hat{\theta}_{\textrm{LP}} + (1-\alpha) \hat{\theta}_{\textrm{VAR}}$ is also unbiased for the parameter of interest $\theta_h$ (this holds for any positive or negative $\alpha \in \mathbb{R}$). Since we want to treat the VAR estimator as the efficient estimator in the absence of potential bias, this means that the linear estimator with the smallest possible variance must coincide with the VAR estimator. This is enforced by setting $\rho = \sigma_{\textrm{VAR}}/\sigma_{\textrm{LP}}$ in (ref), which implies that the covariance between the two estimators is $\sigma^2_{\textrm{VAR}}$. We note that this result can also be established more generally under the locally misspecified vector autoregression model used in montielolea_plagborgmoller_qian_wolf_forthcoming; see their Corollary A.2.
Second, throughout the rest of this section, we will assume that there is a known bound on the bias parameter $b$ that enters the model (ref). The motivation for this bound is the following. Let $T$ denote the sample size used to produce the VAR estimator. montielolea_plagborgmoller_qian_wolf_forthcoming showed that if the researcher is willing to place an upper bound on the fraction of the variance of the moving average residual that is explained by the moving average component of their locally misspecified VAR (denote this parameter by $\mathcal{M}$), then the worst-case bias (in absolute value) for the VAR estimator takes the following form:
See Proposition 4.1 of montielolea_plagborgmoller_qian_wolf_forthcoming and also Equation 3.3 in montielolea_qian_wolf_plagborgmoller_2025.
The question of interest is how to best aggregate the information in LP and VAR estimators to estimate the impulse response coefficient $\theta_h$. This question has been considered recently by nemtyrev2026targeted and chen_pesavento_vonnak_2026. More generally, the question of how to best aggregate the information in two estimators---one of which is potentially biased but more efficient---has been analyzed by armstrong_kline_sun_2025 and lin2026introducing. A related problem is also analyzed in wuetal2026.
We follow closely the work of armstrong_kline_sun_2025 and consider the problem of constructing a minimax estimator (or decision rule) based on quadratic loss, and the bivariate Gaussian model in (ref). We want to impose the bound (ref) on the bias parameter $|b| \leq m \cdot \sigma_{\Delta}$.
It is without loss of generality to focus on estimators $d_h(\hat{\theta}_{\textrm{LP}}, \Delta)$ that depend on the VAR and LP estimators only through the statistics
where the absolute value of $\Delta$ can be viewed as a hausman1978specification test of correct specification of the VAR model that compares the VAR and LP impulse response estimates. A test of this kind was proposed by stock2018identification for testing the invertibility of a vector autoregressive model. It will be convenient to introduce an auxiliary random variable $V \sim \mathcal{N}(0,\sigma^2_{\textrm{LP}}-\sigma^2_{\Delta})$, assumed to be independent of $\Delta$. We can use this random variable to construct a simple representation of the estimation error in the LP estimator. In particular, if we define $U := V - \sigma_{\Delta}(\Delta - b/\sigma_{\Delta}) \sim \mathcal{N}(0, \sigma^2_{\textrm{LP}} )$, then $(U,\Delta)$ has the same joint distribution as $(\hat{\theta}_{\textrm{LP}}-\theta_h,\Delta)$.
We follow armstrong_kline_sun_2025 and consider the following invariance restriction on the estimators $d: \mathbb{R}^{2} \rightarrow \mathbb{R}$:
This restriction means that if we shift the unbiased LP estimator by $s'$ units (but we keep the Hausman statistic constant), then the estimator $d$ should also be adjusted by $s'$ units. Note that, by definition of $\Delta$, any estimator of $\theta_h$ that tries to additively “bias-adjust” the VAR estimator by means of a function $\hat{b}(\cdot)$ that only depends on $\Delta$ satisfies the invariance restriction in (ref). To see this, note that
Note also that any invariant estimator $d_h(\hat{\theta}_{\textrm{LP}},\Delta)$ can be written as a “bias-adjusted” VAR estimator: \[d_h(\hat{\theta}_{\textrm{LP}},\Delta)=d_h(0,\Delta)+\hat{\theta}_{\textrm{LP}}=\hat{\theta}_{\textrm{VAR}}+\underbrace{(-\sigma_{\Delta}\Delta+d_h(0,\Delta))}_{:=-\hat{b}(\Delta)}.\] The risk function of estimators satisfying the invariance restriction in (ref) admits a convenient simplification:
Consequently, the worst-case risk of any invariant estimator does not depend on $\theta_h$, and can be expressed in terms of the minimax risk of the Bounded Normal Mean problem we discussed in Section (ref). In particular, let $\mathcal{D}_{I}$ be the set of all invariant estimators for $\theta_h$. Then,
Moreover, if $d^*_{\textrm{BNM-m}}$ denotes the minimax estimator for the Bounded Normal Mean problem with bound equal to $m$, the minimax invariant estimator of the parameter $\theta_h$ becomes
That is, the minimax invariant estimator for the impulse response coefficient $\theta_h$ is simply a debiased VAR estimator, where the debiasing function is based on a minimax estimator of $b/\sigma_{\Delta}$ (imposing the bound $m$). The estimator of $b/\sigma_{\Delta}$ relies only on the marginal distribution of $\Delta$.
By the same logic as in Section (ref), the stochastic mirror ascent algorithm applied to the Bounded Normal Mean problem with bound $m$ induces, via (ref), a stochastic minimax sequence of invariant estimators for $\theta_h$, with worst-case risk converging to $\bar{v}^*_h(m;\Sigma)$. A formal statement is given in Appendix (ref).
To illustrate our results, we consider four examples of impulse response estimation from four influential macroeconomic papers. In each application, impulse response estimates are obtained using both local projection and vector autoregression methods, including the same set of observable variables and controls. These estimates are based on the replication files provided by montielolea_plagborgmoller_qian_wolf_forthcoming.\footnote{Source: \href{https://github.com/ckwolf92/lp_var_inference}{GitHub repository with replication files for montielolea_plagborgmoller_qian_wolf_forthcoming}. The original data are from ramey_2016.}
Figure (ref) shows the estimated impulse responses along with bootstrapped standard errors for both the LP and VAR estimators.
We now apply the estimator in (ref) to the four applications described above. For each application and horizon $h$, we compute: (1) the hausman1978specification specification test statistic $|\Delta_h|$, (2) the bias correction $\sigma_{\Delta_h} d_{\textrm{BNM-m}}^*(\Delta_h)$, (3) the resulting minimax estimator as described by Equation (ref), and (4) the percentage risk improvement relative to the minimax-optimal linear combination of LP and VAR. For reference, Table (ref) reports the implied values of $m$ and $\epsilon(m)$ across applications. We use those values to run the stochastic mirror ascent algorithm as described in Section (ref) and construct an approximate minimax decision rule corresponding to $\Delta_h$.
Figure (ref) displays the value of the Hausman specification test statistic $|\Delta_h|$ across horizons in each application. We set the Hausman test statistic to be equal to zero at horizon $h=0$ (although it theoretically equals 0 divided by 0, since the LP and VAR estimators coincide on impact).
For each application, we compute the approximate rule $d_{\textrm{BNM-m}}^*$ in the corresponding minimax Bounded Normal Mean problem using stochastic mirror ascent. The bias correction of the VAR estimator is nonlinear in $\Delta_h$ and is given by $\sigma_{\Delta_h} d_{\textrm{BNM-m}}^*(\Delta_h)$.
Figure (ref) plots the implied nonlinear bias correction across horizons for the four different applications.
We then construct the bias-corrected VAR estimator for impulse responses in each of these applications. Figure (ref) illustrates the comparison of this estimator to the underlying VAR and LP estimators.
To benchmark the nonlinear bias-corrected VAR, we compare it with the minimax-optimal linear combination of LP and VAR. That is, a linear combination of VAR and LP estimators, with weights chosen to minimize the worst-case risk. As shown in montielolea_plagborgmoller_qian_wolf_forthcoming, the worst-case risk for such a linear combination is given by \[(1-\omega)^2m^2(\sigma_{\textrm{LP}}^2 - \sigma_{\textrm{VAR}}^2) + \sigma_{\textrm{VAR}}^2 + \omega^2(\sigma_{\textrm{LP}}^2 - \sigma_{\textrm{VAR}}^2).\] The optimal weight that minimizes this expression is \[\omega^* = \frac{m^2}{1+m^2},\] leading to the estimator \[d^*_{\textrm{linear}}=\hat{\theta}_{\textrm{VAR}} - \sigma_\Delta \omega^* \Delta.\]
In the following exercise, we compare these values with $\mathcal{R}(d_h,\theta_h,b)$, which uses the (grid-evaluated) worst-case risk of the optimal decision rule in the Bounded Normal Mean problem (1,000 Monte Carlo draws).
Figure (ref) compares the estimated worst-case risk of the nonlinear minimax correction with the worst-case risk of this optimal linear combination. Positive values indicate that the nonlinear rule improves on the best linear rule.
We presented a computational approach to find an approximately minimax estimator (under quadratic loss) for the Bounded Normal Mean problem of casella_strawderman_1981. Our suggested algorithm considers a fixed (equally-spaced) grid of $I$ points that discretize the parameter space $[-m,m]$. We then approximately maximize Bayes risk over the space of probability distributions supported on the $I$ points of the grid using a stochastic mirror ascent routine for concave maximization. We use this approximately least-favorable distribution to find the corresponding Bayes estimator. Our main results (Theorems (ref)-(ref)) showed that the approximately least-favorable distribution and the corresponding Bayes estimator are indeed close (in a sense we made precise) to their desired targets.
As we explained in the paper, our algorithm has two sources of approximation error: discretization and optimization error. The Lipschitz continuity property of the risk function that we established in Lemma (ref) controls the first component. The well-known theoretical properties of stochastic mirror descent allow us to control the second component. Taking into account these two sources of approximation error, we showed that the approximately least-favorable distribution obtained from our algorithm converges in $1-$Wasserstein distance to the exact least-favorable distribution of the Bounded Normal Mean problem. In addition, the approximately minimax estimator converges (uniformly over compacts) to the exact minimax estimator.
We think that several questions remain open. First, as we have explained in the paper, our algorithm relies on an ex ante discretization of the parameter space $[-m,m]$. This makes the problem tractable and allows the use of mirror ascent on a finite-dimensional probability simplex, but it also introduces discretization error that needs to be accounted for. A natural direction for future work would be to explore the possibility of approximating the least-favorable distribution by looking directly at the concave problem of maximizing Bayes risk over all probability distributions on the interval $[-m,m]$. An approach for this problem could be based on the algorithm presented in the work of kent2021modified for optimizing infinite-dimensional functionals of probability measures.\footnote{In Appendix (ref) we derived different properties of the general maximin problem that could be helpful for this purpose.} Another possibility would be to use recent advances in online learning with continuous action spaces; see krichene2015hedge and negrea2021minimax. Second, some of the theoretical results we derived in order to control optimization error are based on constants/bounds that are conservative. For example, the Lipschitz constant for the risk function we derived in Lemma (ref) could perhaps be tighter (which could decrease the number of iterations in Theorem (ref)). Similarly, the number of iterations for the mirror-ascent algorithm is designed to work for any (Lipschitz) concave optimization problem over the simplex. Since we are only interested in solving the Bounded Normal Mean problem, it is perhaps possible to sharpen the analysis and reduce the number of iterations. Third, the least-favorable distribution in the Bounded Normal Mean problem is known to have a very particular structure: the support of the distribution is finite. The stochastic mirror ascent algorithm that we use is initialized at the uniform distribution in $\Delta^{I-1}$, and by design, the output of the algorithm places positive mass on every point (although as our simulations show, the mass associated with some of these points is numerically very close to zero). It could be useful to consider other algorithms for finding the least-favorable distribution that increase the number of support points in every iteration, as in the recent work of guggenberger2025numerical.
We view our paper as part of a larger research program that emphasizes the importance of providing a careful analysis of the computational aspects surrounding problems in econometrics and statistics; see, for example, the recent work of cattaneo2024onimplicitbias,cattaneo2025howmemory,cattaneo2025modifiedlossmomentum. As mentioned by chandrasekaran2013computational, “Although classical statistics gave little consideration to computational complexity, computational issues have come increasingly to the fore in modern high-dimensional statistics.” There are many examples now within statistics and machine learning that illustrate the benefits of bringing computational awareness to the forefront of these fields.\footnote{See, for illustration, the lectures of the “Computational Complexity of Statistical Inference Boot Camp” of the Simons Institute for the Theory of Computing: \href{https://www.youtube.com/playlist?list=PLgKuh-lKre13obZNGYdY-25ihfD0PIyJ-}{video lectures}.} Advances in econometrics in this regard have not followed the same pace, but as mentioned recently by Patrick Kline in one of his interventions at the Online Chamberlain Seminar reflecting on the future of econometrics, “There is probably going to be a more important role for computation moving forward. I think we have not rewarded computational advances enough in the profession, and we ought to be more serious about that ... And it is important for procedures that are proposed to actually be implementable. And in some cases ... Having some guarantees on computability of different procedures is important and I think we ought to learn more computer science and we will.”\footnote{See minute 19:58 of \href{https://www.youtube.com/watch?v=ifp_BUKl3ig}{this talk}}
Finally, we hope that the methods outlined in this paper could make the application of statistical decision theory in econometrics more frequent. As noted by manski2021econometrics, the “primary challenge to use of statistical decision theory is computational.”
\doublespacing
{20pt}