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.
112,681 characters · 15 sections · 66 citation commands
Inference by Stochastic Optimization: A Free-Lunch Bootstrap
JEL Classification: C2, C3
Keywords: Stochastic gradient descent, Newton-Raphson, Simulation-Based Estimation.
\setcounter{page}{0} \thispagestyle{empty} \baselineskip=18.pt
Many questions of economic interest can be expressed as non-linear functions of unknown parameters $\theta$ that need to be estimated from a sample of data of size $n$. The typical econometric routine is to first obtain a consistent estimate $\hat\theta_n$ of the true value $\theta^0$ by minimizing an objective function $Q_n(\theta)$, after which its sampling uncertainty is assessed. Though gradient-free optimizers provide point estimates, its asymptotic variance is often analytically intractable. One remedy is to use bootstrap standard errors, but this requires solving the minimization problem each time the data is resampled, and for complex models, this is no simple task. There is a long-standing interest in finding `short-cuts' that can relieve the computation burden without sacrificing too much accuracy. Examples include Davidson1999, Andrews2002, Kline2012, Armstrong2014 and more recently Honore2017. These methods provide standard errors by taking a converged estimate $\hat\theta_n$ as given. As such, estimation always precedes inference.
This paper proposes a resampling scheme that will deliver both the point estimates of $\theta$ and its standard errors within the same optimization framework. Since the standard errors are obtained as a by-product of point estimation, we refer to the procedure as a `free-lunch bootstrap'.\footnote{In optimization, the no-free lunch theorem of wolpert-macready:97 states that, when averaged over all problems, the computation cost of finding a solution is the same across methods. We use the term to refer to the ability to compute the quantities for inference when the estimator is constructed.} The free-lunch is made possible by a specially designed stochastic optimization algorithm that resamples batches of data of size $m\le n$. Given an initial guess $\theta_0$, one updates $\theta_{b}$ for $b\ge 0$ to $\theta_{b+1}$ using the gradient, the inverse Hessian as conditioning matrix, and a suitably chosen learning rate. We first show that the average over $B$ draws of $\theta_{b}$ is equivalent to the mode $\hat\theta_n$ obtained by classical optimization up to order $\frac{1}{m}$. We then show that the distribution of $\sqrt{m}(\theta_b-\hat\theta_n)$ conditional on the original sample of data is first-order equivalent to that of $\sqrt{n}(\hat\theta_n-\theta^0)$ upon rescaling, making it a bootstrap distribution. Because the conditioning matrix is the inverse Hessian, the procedure is a resampled Newton-Raphson (rnr) algorithm. For other conditioning matrices, the draws from resampling still produce a consistent estimate but cannot be used for inference.
The main appeal of the proposed methodology is its simplicity. If the optimization problem can be solved by our stochastic optimizer, inference can be made immediately without further computations. Natural applications include two-step estimation when the parameters in the two steps are functionally dependent in a complicated way, as well as minimum distance estimation that compares the empirical moments with the model moments expressed as a function of the parameters. When such a mapping cannot be expressed in closed-form, simulation estimation makes progress by using Monte-Carlo methods to approximate the binding function, but computing standard errors of the simulation-based estimates remains a daunting task. Our algorithm provides an automated solution to compute standard errors and removes simulation noise, resulting in more accurate estimates. The algorithm also provides a convenient way to compute clustered standard errors and model diagnostics.
As compared to other stochastic optimization algorithms, we use a learning rate that is fixed rather than vanishing, and though a small $m$ is desirable from a pure computation perspective, valid inference necessitates that $m$ cannot be too small. As compared to conventional bootstrap methods, the simultaneous nature of estimation and inference means that a preliminary $\hat\theta_n$ is not needed for resampling. Though $\theta_{b}$ is a Markov chain, no prior distribution is required, nor are Bayesian computation tools employed. In simulated examples and applications, our bootstrap standard errors match up well with the asymptotic and bootstrap analogs, but at significantly lower computational costs.
The plan of the paper is as follows. Section (ref) begins with a review of classical and stochastic optimization. The proposed free-lunch algorithm is presented in Section (ref) and its relation to other resampling procedures is explained. The properties of the draws from the algorithm are derived in Section (ref). Simulated and empirical examples are presented in Section (ref). Section (ref) extends the main results to simulation-based estimation. Appendix (ref) provides derivations of the main results. An on-line supplement\footnote{The file is available for download at \textcolor{blue}{\url{www.columbia.edu/ sn2294/papers/freelunch-supp.pdf}}.} provides the r code to implement one of the applications considered, additional results with details for replications, as well as an analytical example for Section (ref).
Consider minimization of the objective function $Q_n(\theta)$ with respect to $\theta$ whose true value is $\theta^0$. The sample gradient and Hessian of $Q_n(\theta)$ are defined respectively by
The necessary conditions for a local minimum are $\|G_n(\hat\theta_n)\|=0$ and $H_n(\hat\theta_n)$ positive semi-definite. The sufficient conditions are $\|G_n(\hat\theta_n)\|=0$ and $H_n(\hat\theta_n)$ positive definite. To find the optimal solution, a generic rule for updating from the current estimate $\theta_k$ is \[ \theta_{k+1}=\theta_k -\gamma_k Z_n(\theta_k)\] where $\gamma_k$ is the step size and $Z_n=\frac{\partial \theta_{k+1}}{\partial \gamma_k}$ is the direction of change.
Gradient based methods specify $Z_n(\theta_k)= P_n(\theta_k) G_n(\theta_k)$ where $P_n(\theta_n)$ is a conditioning matrix. The updating rule then becomes
The method of gradient descent (gd) (also known as steepest descent) sets $ P_n=I_d$. Since gd does not involve the Hessian, it is a first order method and is less costly. Convergence of $\hat\theta_k$ to the minimizer $\hat\theta_n$ is linear under certain conditions,\footnote{In statistical computing, the convergence of $\theta_k$ to $\hat\theta_n$ is said to be linear if $\|\theta_{k+1}-\hat\theta_n\|/\|\theta_k-\hat\theta_n\|^q<r$ for some $r\in (0,1)$ if $q=1$ and quadratic if $q=2$. Convergence is superlinear if $\lim_{k\rightarrow\infty} \|\theta_{k+1}-\hat\theta_n\|/ \|\theta_k-\hat\theta_n\|=0$. See boyd-vandenberghe:04 Section 9.3.1 for linear convergence of gradient methods, and nocedal-wright:06 for quadratic convergence of Newton's method when $\gamma=1$ or $\gamma_k \to 1$ at an appropriate rate. `Damped Newton' updating with $\gamma_k \in (0,1)$ has a linear convergence rate, see boyd-vandenberghe:04 Section 9.5.3 and nesterov2018 Section 1.2.4.} but the rate depends on $I_d-\gamma H_n(\theta_k)$ being in the restricted region of $(-1,1)$, and can be slow when the ratio of the maximum to the minimum eigenvalue of $H_n$ is large. The Newton-Raphson algorithm puts $P_n=H_n(\theta_k)^{-1}$. It is a second-order method since it involves the Hessian matrix. When $\gamma_k=1$, the algorithm converges quadratically if $Q_n$ satisfies certain conditions. A drawback of Newton's algorithm is that it requires computation of the inverse of the Hessian. When strong convexity fails, the Hessian could be non-positive definite for $\theta$ away from the minimum. In such cases, it is not uncommon to replace the Hessian by $H_n(\theta_k)+c\cdot I_d$ for some $c>0$, or specify $P_n=(H_n(\theta_k)^\prime H_n(\theta_k))^{1/2}$ to restore positive definiteness around saddle-points, see nocedal-wright:06. Quasi-Newton methods bypass direct computation of the Hessian or its inverse, but analytical convergence results are more difficult to obtain. We focus our theoretical analysis on gradient descent and Newton-Raphson based algorithms but will consider quasi-Newton methods in some simulations.
Stochastic optimization finds the optima in noisy observations using carefully designed recursive algorithms. The idea can be traced to the theory of stochastic approximation when the goal is to minimize some function $Q(\theta)$ with gradient $G(\theta)$, which is equivalent to the root-finding problem $G(\theta)=0$ whose the true value is $\theta^0$. A classical optimizer would perform $ \theta_{k+1}=\theta_{k}-\gamma_{k} G(\theta_{k})$. Robbins1951 considers the situation when we only observe $G(\theta_{k})+e_k$ with $\mathbb E(e_k)=0$ and suggests to update according to \[ \theta_{k+1}=\theta_{k}- \gamma_{k} (G(\theta_{k})+e_k).\] Robbins1951 proved that $\theta_k\smash{\mathop{\longrightarrow}\limits^{as}} \theta^0$ for $G$ non-decreasing with step size sequence $\gamma_k \geq 0$ satisfying
The first condition ensures that all possible solutions will be reached with high probability regardless of the starting value, while the second ensures convergence to the true value. Building on the Robbins-Monro algorithm, the Kiefer-Wolfowitz algorithm uses a finite difference approximation $G(\theta_k)\approx G_n(\theta_k)=\frac{1}{2\epsilon_k}\bigg[ Q_n(\theta_k+\epsilon_k)-Q_n(\theta_k-\epsilon_k)\bigg]$. This is often recognized as the first implementation of stochastic gradient descent. Kiefer1952Stochastic proves convergence of $\theta_k$ to the maxiumum likelihood estimate $\hat\theta_n$ assuming that the likelihood $Q_n$ is convex, $\epsilon_k$ goes to zero, and that the two conditions stated in ((ref)) hold.
Modern stochastic gradient descent updates according to
where $G_m(\theta_k)=\frac{1}{m}\sum_{i=1}^m G(\theta_k;x_i)$ is an estimate of $G(\theta)$. It can be seen as Monte-Carlo based since the $m$ observations used to compute $G_m(\theta_k)$ are usually chosen from $\{1,\dots,n\}$ randomly. Though $m=1$ is computationally inexpensive and is a popular choice, a small $\gamma_k$ is often needed to compensate for the higher variation. A common rule is to choose $\gamma_k = \gamma k^{-\delta}$, where $\delta \in (1/2,1]$ and $\gamma>0$ are the choice parameters. Depending on $\delta$, convergence as measured by $\mathbb{E}(\|\theta_k-\theta^0\|^2)$ can occur at a $1/k$ rate or slower. To reduce sensitivity to the tuning parameters, ruppert1988efficient and polyak1992acceleration propose to accelerate convergence using what is now known as Polyak-Ruppert averaging: $\overline{\theta}_k = \frac{1}{k}\sum_{i=1}^k \theta_i$. Importantly, $\overline{\theta}_k$ converges at the fastest $1/k$ rate for all choices of $\delta \in (1/2,1]$. Bach2011 shows that the improvements hold even for a finite number of iterations $k$. We will return to Polyak-Ruppert averaging below.
Stochastic optimization presents an interesting alternative to classical optimization as it approximates the gradient on {\it minibatches} of the original data. This is particularly helpful in large scale learning problems such Lasso, support-vector machines and K-means clustering when non-linear optimization can be challenging. Improvements to sgd with $P_m=I_d$ include momentum polyak:64 and accelerated gradient nesterov:83 methods. Besides its computational appeal, stochastic optimization can improve upon its classical counterpart in non-convex settings.\footnote{See Goodfellow-et-al-2016 for an overview of sgd. ge2015escaping shows that noisy gradient descent can escape all saddle points in polynomial time under a strict saddle property whereas classical gradient methods converge at saddle points where the gradient is zero. jin2017escape shows that the dimension of $\theta$ has a negligible effect on the number of iterations needed to escape saddle points, making it an effective solution even in large optimization problems.}
A variation of sgd, known as Stochastic gradient Langevin dynamics (sgld) incorporates Langevin dynamics into a Bayesian sampler. As will be discussed further below, the update is based on the gradient of the posterior distribution. Of note now is that sgld has two types of noises: an injected noise, and the stochastic gradient noise based on $m \ll n$ observations. They play different roles in the algorithm. In the early phase of sgld, the stochastic gradient dominates and the algorithm performs optimization. In the second phase, the injected noise dominates and the algorithm behaves like a posterior sampler. The algorithm seamlessly switches from one phase to another for an appropriate choice of the learning rate.
Unlike classical optimization, stochastic Newton-Raphson with the inverse Hessian $H_m(\theta)$ as conditioning matrix is not popular because the Hessian is often noisy and near singular for $m$ small, rendering the algorithm unstable. We will show that using a variation of stochastic Newton-Raphson with larger batches of data can produce draws that not only provide an accurate estimate of $\theta^0$ but also yields frequentist assessment of sampling uncertainty. It thus integrates numerical optimization with statistical inference. In contrast, other conditioning matrices will yield consistent estimates but would not provide valid inference in our setup.
Consider extremum estimation of parameters $\theta\in\Theta \subset \mathbb{R}^d$ from data $x=(x_1,\ldots, x_n)$. Let $\theta^0$ be the minimizer of a twice differentiable population objective function $Q(\theta)$ whose sample analog is $ Q_n(\theta)\equiv Q_n(\theta;x)$. The sample extremum estimator is \[ \hat\theta_n=\text{argmin}_{\theta \in \Theta} Q_n(\theta).\] For likelihood estimation, $Q_n(\theta)=-\sum_{i=1}^n \ell_i(\theta)$ where $\ell_i$ is the likelihood of $\theta$ at observation $x_i$. For least squares estimation, $Q_n(\theta)$ is the sum of squared residuals $\sum_{i=1}^n e_i^2(\theta)$. For GMM estimate with positive weighting matrix $W_n$, $Q_n(\theta)=\overline{g}_n(\theta)^\prime W_n \overline{g}_n(\theta)$ where $\mathbb E [g_i(\theta^0)]=0$. Under regularity conditions stated in Theorem 2.1 of newey-mcfadden-handbook $\hat\theta_n$ is consistent for $\theta^0$. If, in addition, the assumptions in Theorem 3.1 of newey-mcfadden-handbook hold, then $\hat\theta_n$ is also $\sqrt{n}$-asymptotically normal: \[\sqrt{n}(\mathbb{V}^{0})^{-1/2}(\hat\theta_n-\theta^0)\smash{\mathop{\longrightarrow}\limits^d} N(0,I_d) \] where $\mathbb V^0=[H(\theta^0)]^{-1} \text{var}(\sqrt{n}G_n(\theta^0)) [H(\theta^0)]^{-1}$. Finite sample inference is typically based on an estimate of $\mathbb V^0$ which can be analytically intractable or costly to compute on the full sample. It is not uncommon to resort to bootstrap inference. We consider the $m$ out of $n$ bootstrap with $m \to \infty, m/n \to c\in [0,1]$ and samples $(x_1^{(b)},\dots,x_{m}^{(b)})$ with replacement from the data $(x_1,\dots,x_n)$ for $b=1,\dots,B$ and solve $B$ minimization problems: \[ \hat \theta_m^{(b)} =\text{argmin}_{\theta \in \Theta} Q_m^{(b)}(\theta), \] where the resampled objective $Q_m^{(b)}(\theta)=Q_m^{(b)}(\theta,x^{(b)})$ is computed over the sample $x^{(b)}$.
Let $\mathbb{E}^\star$ and $\text{var}^\star$ denote the bootstrap expectation and variance which are taken conditional on the sample data $(x_1,\dots,x_n)$, and $\overset{d^\star}{\to}$ denotes the convergence in distribution conditional on the data. Since we only consider correctly specified regular estimators, the desired convergence is to a Gaussian limit. The $m$ out of $n$ bootstrap can allow for different sampling schemes. Assuming that the resampling scheme is chosen to reflect the dependence structure of the data, it holds in a variety of settings that\footnote{For two-way clustering, a recommended procedure is to resample over one cluster dimension and reweigh along the other roodman2019fast. See also cameron2011robust on multiway clustering. For time-series data, block resampling is needed to preserve the dependence structure. For correctly specified GMM models the bootstrap described above is valid hahn1996. } \[ \sqrt{m} (\mathbb{V}_m)^{-1/2} \left( \hat\theta_m^{(b)} - \hat\theta_n \right) \overset{d^\star}{\to} \mathcal{N} \left( 0, I_d \right),\] where $\mathbb V_m=[H_n(\hat\theta_n)]^{-1} \text{var}^\star(\sqrt{m}G_m^{(b)}(\hat\theta_n)) [H_n(\hat\theta_n)]^{-1}$ depends on the inverse Hessian as well as the variance of the resampled score. The result implies that the resampled distribution of $\hat\theta_m^{(b)}$ can be used to approximate the sampling distribution of $\hat\theta_n$.
An $m$ out of $n$ bootstrap with $m<n$ uses smaller samples and performs as well as a $n$ out of $n$ bootstrap in simulations while requiring similar or weaker conditions, bickel2012. Nonetheless, it still makes multiple calls to a classical optimizer. As will be seen below, our proposed algorithm only requires one call to the optimizer.
We propose the following two algorithms and use `$b$' to index the iterates of the algorithm. In this notation, $G_m^{(b+1)}(\theta_{b})$ is the gradient computed on the $(b+1)$-th batch of resampled data of size $m$ and evaluated at $\theta_{b}$, the parameter value in the previous draw, and $H_m^{(b+1)}(\theta_{b})$ is similarly defined.
Algorithm (ref) produces an estimate $\overline\theta_{\textsc{re}}$ by resampling, hence the acronym re. It works for any conditioning matrix $P_b$ satisfying assumptions to be made precise in Theorem (ref). Algorithm (ref) produces draws using the inverse Hessian as $P_b$ as in Newton-Raphson, hence the acronym rnr. The free-lunch aspect relates to the fact that we get both an estimate $\overline{\theta}_{r\textsc{nr}}$ and its standard error in one run. The bootstrap aspect comes from the fact that under the assumptions of Theorem (ref), $\sqrt{m}(\theta_b-\hat\theta_n)$ has the same asymptotic distribution as $\sqrt{n}(\hat\theta_n-\theta^0)$ after an adjustment of $\frac{m}{\phi(\gamma)}$. The quantity $V_{r\textsc{nr}}$ is an estimate of the sandwich variance that is computationally costly for classical estimation. A Wald test for $H_0:\theta=\theta^\dag$ can be constructed as \[ \text{wald}= n(\overline\theta_{r\textsc{nr}}-\theta^\dag)^\prime V_{r\textsc{nr}}^ {-1} (\overline\theta_{r\textsc{nr}}-\theta^\dag)\] which has an asymptotic Chi-squared distribution under the null hypothesis. A 95% level bootstrap confidence interval can also be constructed after adjusting for $\frac{m}{n}$ and $\phi(\gamma)$ by taking the (0.025,0.975) quantiles of $\big\{ \overline{\theta}_{r\textsc{nr}}+\sqrt{\frac{m}{n\phi(\gamma)}}(\theta_{b} - \overline{\theta}_{r\textsc{nr}})\big\}_{b \geq 1}$.
Algorithms (ref) and (ref) have three features that distinguish them from existing gradient-based stochastic optimizers. First, $\gamma\in (0,1]$ does not change with $b$. Fixing $\gamma$ rather than letting $\gamma_b\rightarrow 0$ potentially permits faster convergence. Second, we sample $m$ out of $n$ observations with $m/n\rightarrow c\in [0,1]$ and $\sqrt{n}/m \to 0$. This precludes the popular choice in stochastic optimization of $m=1$, but admits $m=n$. We thus accept a higher computation cost to accommodate inference. Third, compared to sgd\, Algorithm (ref) uses the inverse Hessian as conditioning matrix.
This subsection uses the linear regression model to gain intuition of the Free-Lunch bootstrap. The model is $y_i=x_i^\prime\theta+e_i$. Let $\hat e_n = y_n-X_n \hat\theta_n$ be the $n\times 1$ vector of least squares residuals evaluated at the solution $\hat\theta_n$. $X_n$ denote the $n\times K$ matrix of regressors. The linear model is of interest because the objective function is quadratic and the quantities required for updating are analytically tractable. The gradient and Hessian of the full sample objective function $ Q_n(\theta) =(y_n-X_n\theta)^\prime(y_n-X_n\theta)/(2n)$ are $ G_n(\theta) = -X_n^\prime (y_n-X_n\theta)/n$ and $H_n(\theta) = X_n^\prime X_n/n$. The updates for this linear model evolve as
Convergence of $\theta_k$ for a given conditioning matrix $P_k$ can be studied by subtracting $\hat\theta_n$ from both sides of the updating equation and re-arranging terms (see Appendix (ref) for details). Table (ref) summarizes convergence of gd, nr, sgd, rgd\, and rnr. snr\, is not considered because $X_m^\prime X_m/m$ is singular for $m=1$ so $\theta_b$ is not well defined. The left panel of the table gives the updating rule in closed form and the right panel expresses the deviation of the draws from $\hat\theta_n$ as the sum of a deterministic and a stochastic component.
As seen from Table (ref), gd updates do not depend on the Hessian but convergence does, while for nr\, the opposite is true. Convergence of nr\; can be achieved after one iteration if $\gamma=1$. In sgd, rgd\, and rnr, batch resampling adds a stochastic component to the updates and convergence is no longer deterministic. The deviations for s\textsc{gd}\, and r\textsc{gd}\, $\theta_{b+1}-\hat\theta_n$ follow a VAR(1) process with varying and fixed coefficient matrices $I_d-\gamma_b H_b$ and $I_d-\gamma H_b$, respectively. In contrast, the r\textsc{nr}\ draws have an AR(1) representation with a fixed coefficient $(1-\gamma)$ that is dimension-free and independent of the Hessian. Note that r\textsc{gd}\, and r\textsc{nr}\, keep $\gamma$ fixed and rely on averaging over $b$ for convergence.
Our main result pertains to rnr\; so it is useful to have a deeper understanding of how it works. Unlike stochastic optimizers which require $\gamma_b$ vanishing, the learning rate $\gamma$ used to generate the rnr\, draws is constant. The draws evolve according to
where $\hat\theta_m^{(b+1)}=(X_m^{(b+1)\prime} X_m^{(b+1)})^{-1}X_m^{(b+1)\prime} y_m^{(b+1)}$ is obtained by classical optimization using the $(b+1)$-th bootstrap sample $(y_i^{(b+1)},x_{i}^{(b+1)})_{i=1,\dots,m}$. Being a bootstrap estimate, it holds under regularity conditions that the distribution of $\sqrt{m}(\hat\theta_m^{(b+1)}-\hat\theta_n)$ conditional on the data approximates the sampling distribution of $\sqrt{n}(\hat\theta_n-\theta^0) \smash{\mathop{\longrightarrow}\limits^d} N(0,\mathbb V^0)$.
Clearly when $\gamma=1$, ((ref)) implies $\theta_{b} = \hat\theta_m^{(b)}$, meaning that each rnr\; draw equals the bootstrap estimate $\hat\theta_m^{(b)}$. We want to show that the draws are still bootstrap estimates when $\gamma\in(0,1)$. For such $\gamma$, $\theta_{b+1}-\hat\theta_n$ is an AR(1) process where for each $b$, the innovations $ \gamma(\hat\theta_m^{(b+1)}-\hat\theta_n)$ are iid conditional on the original sample. Iterating the AR(1) formula backwards to the initial value $\theta_0$, we can decompose the draws $\theta_{b+1}$ into two terms:
where $\{\hat\theta_m^{(b+1-j)}\}_{j \geq 0}$ are the bootstrap estimates in the previous iterations. The constant learning rate is crucial in achieving this representation.
To show that our estimator $\overline\theta_{r\textsc{nr}}$ is $\sqrt{n}$-consistent for $\hat\theta_n$, i.e. $\overline \theta_{r\textsc{nr}}=\hat\theta_n+o_{p^\star}(\frac{1}{\sqrt{n}})$, we need to evaluate the average of the two terms in ((ref)) over $b$. The initialization bias in ((ref)) is due to taking an arbitrary starting value $\theta_0$ and is identical to the optimization error in classical Newton-Raphson. For $\gamma\in(0,1]$, $\frac{1}{B}\sum_{b=1}^B (1-\gamma)^{b+1}=O(\frac{1}{B})$ because $\{ (1-\gamma)^{b} \}_{b \ge 1}$ is a summable geometric series. Another bias term of order $O(\frac{1}{m})$ arises when $\mathbb{E}^\star(\hat\theta_m^{(b)}-\hat\theta_n)=O(\frac{1}{m})$. Since $\hat\theta_n$ is fixed as $b$ varies, we now have $\mathbb E^\star (\overline{\theta}_{r\textsc{nr}})=\hat\theta_n+O(\frac{1}{B}) +O(\frac{1}{m}).$ Thus $\mathbb E^\star (\overline{\theta}_{r\textsc{nr}})=\hat\theta_n +o(\frac{1}{\sqrt{n}})$ as required, assuming $\frac{\sqrt{n}}{\text{min}(m,B)}\rightarrow 0$. Turning to the variance, first note that by virtue of bootstrapping, $\{\hat\theta_m^{(b+1-j)}-\hat\theta_n\}_{j \geq 0}$ constitutes a sequence of conditionally iid errors each with variance that is $O(\frac{1}{m})$. Since $\{ (1-\gamma)^{b} \}_{b \geq 1}$ is summable, the variance in $\overline{\theta}_{r\textsc{nr}}$ due to resampling is $O(\frac{1}{mB})$. This becomes $o(\frac{1}{n})$ when $\frac{n}{mB}\rightarrow 0$, a sufficient condition being $\frac{\sqrt{n}}{\text{min}(m,B)} \to 0$, which is also required for the bias to be negligible. We have thus shown that $\overline{\theta}_{r\textsc{nr}} = \hat\theta_n + o_{p^\star}(\frac{1}{\sqrt{n}})$ for $\frac{\sqrt{n}}{\text{min}(m,B)} \to 0$, which is a simplified version of Theorem (ref) below. Though the result has the flavor of Polyak-Ruppert averaging in stochastic optimization, $\gamma$ is fixed here and $m$ increases with $n$.
To show bootstrap validity of rnr, we need to establish that, conditional on the sample of data, the distribution of $\sqrt{m}(\theta_{b+1}-\hat\theta_n)$ is asympotically equal, up to a constant scaling factor, to that of $\sqrt{m}(\hat\theta_m^{(b+1)}-\hat\theta_n)$. This requires that the initialization bias in each $\theta_{b}$ is $o(\frac{1}{\sqrt{m}})$, which holds when $\frac{\log(m)}{b} \to 0$. From ((ref)), $\sqrt{m}(\theta_{b+1}-\hat\theta_n)$ has variance $\gamma^2 \mathbb{V}_{m}$ conditional on $\theta_{b}$ and unconditional variance \[\text{var}\bigg(\sqrt{m}(\theta_{b+1}-\hat\theta_n)\bigg)=\frac{\gamma^2 + O([1-\gamma]^{b+2})}{1-[1-\gamma]^2}\mathbb{V}_{m}\approx \phi(\gamma)\mathbb V_m\] where $\phi(\gamma) = \frac{\gamma^2}{1-[1-\gamma]^2}$, and $\mathbb V_m=\text{var}(\sqrt{m}(\hat\theta^{(b+1)}_m-\hat\theta_n))$ is the bootstrap estimate of the sandwich variance $\mathbb{V}^0$ defined above. This establishes that the variance of $\theta_{b}$ is proportional to that of the bootstrap estimate. As shown in gonccalveswhite2005, $\mathbb V_m$ is consistent for $\mathbb V^0$ under certain moment conditions. This implies that, up to the scaling factor $\phi(\gamma)$, the co-variance of $\theta_b$ is also consistent for $\mathbb V^0$. Combined with asympotic normality of each $\sqrt{m}(\hat\theta_m^{(b+1)}-\hat\theta_n)$ for each $b$ and additional conditions to be made precise in Theorem (ref), we have \[ \bigg(\phi(\gamma)\mathbb{V}_{m}\bigg)^{-1/2}\sqrt{m} \left( \theta_{b+1} - \hat \theta_n \right) \overset{d^\star}{\to} \mathcal{N} \left( 0, I_d\right). \] But asymptotic theory gives the distribution of $\sqrt{n}(\hat\theta_n-\theta^0)$ with sample size $n$, not $m$. An adjustment for $\phi(\gamma)$ and $m$ is needed. Let $\mathbb V_{r\textsc{nr}}= \frac{m}{\phi(\gamma)}\text{var}^\star(\theta_{b}) = \mathbb{V}_m + o(1)$. For appropriate choice of $m$ and $\gamma$, $ \mathbb V_{r\textsc{nr}}^{-1/2} \sqrt{n}(\overline\theta_{r\textsc{nr}}-\theta^0) \overset{d^\star}{\to} \mathcal N(0,I_d) $ and Algorithm (ref) proposes a plug-in estimate of $\mathbb V_{r\textsc{nr}}$.
This section studies the properties of draws $\theta_{b}$ produced by Algorithms (ref) and (ref) for non-linear models. The proofs are more involved for two reasons. First, an arbitrary $\gamma \in (0,1]$ may not lead to convergence even for classical optimizers. Second, whereas in quadratic problems the draws $\theta_{b}$ have a tractable AR(1) representation, for non-quadratic objectives the draws $\theta_{b}$ follow a non-linear process which is more difficult to study. Hence we need to first show, under strong convexity conditions, that there exist fixed values of $\gamma \in(0,1]$ such that optimization of $Q_n(\theta)$ has a globally convergent solution. We then show, using the idea of coupling, for appropriate choices of $m$ and $B$ that $\theta_{b+1}$ can be made very close to a linear AR(1) sequence $\theta_{b+1}^\star$ that is constructed as if the objective were quadratic. This allow us to establish consistency of $\overline\theta_{\textsc{re}}$ for $\hat \theta_n$ in Theorem (ref) for a large class of $P_b$, and a distribution result in Theorem (ref) that validates inference for a particular choice of $P_b$.
Econometric theory typically studies the conditions under which $\hat\theta_n$ is consistent for $\theta^0$, taking as given that a numerical optimizer exists to produce a convergent solution $\hat\theta_n$. From newey-mcfadden-handbook, the regularity conditions for consistent estimation of $\theta$ are continuity of $Q(\theta)$ and uniform convergence of $Q_n(\theta)$ to $Q(\theta)$. Asymptotic normality further requires smoothness of $Q_n(\theta)$, $\theta^0$ being in the interior of the support, and non-singularity of $H(\theta^0)$. But classical Newton-type algorithms may only converge to a local minimum and a global convergent solution is guaranteed only when the objective function is strongly convex on the parameter space $\Theta$. For gradient-based optimizers to deliver such a solution, the following provides the required conditions.
Condition i. implies strong convexity of $Q_n$ on $\Theta$.\footnote{See boyd-vandenberghe:04, Chapter 9.1. A function $Q_n$ is strongly convex on $\Theta$ if for all $\theta\in\Theta$, there exists some $\underline{\lambda}>0$ such that $\nabla^2 Q_n(\theta) \ge \underline{\lambda} I_d$. For bounded $\Theta$, there also exists $\overline{\lambda}$ such that $\nabla^2 Q_n(\theta) \le \overline{\lambda} I_d$. Then $\overline{\lambda}/\underline{\lambda}$ is an upper bound on the condition number of $\nabla^2 Q_n(\theta)$.} Condition ii. imposes Lipschitz continuity of the Hessian. Assumption (ref) implies the following two inequalities which are known as the Polyak-\Lojasiewicz inequalities:
where $\tilde{\theta}_n$ is an intermediate value between $\theta$ and $\hat{\theta}_n$. Inequality ((ref)), due to lojasiewicz1963 and polyak1963, follows from the positive definiteness of $H_n(\tilde \theta_n)$. Together, ((ref)) and ((ref)) ensure that $\hat \theta_n$ is a unique (or global) minimizer of $Q_n$.
Assumption (ref) also implies that there exists $\gamma$ such that gradient based optimization is globally convergent. To see why, consider
where the last inequality is implied by Assumption (ref) i. and iii. Since a contraction occurs if $A(\gamma)\in [0,1)$, global convergence follows. Now at $\gamma=0$, $A(0)=1$ and $\partial_\gamma A(0)<0$, so by continuity and local monotonicity of $A(\cdot)$, there exists a nonempty subinterval of the form $(0,\tilde \gamma]$ with $\tilde \gamma \in (0,1]$ such that $A(\gamma) \in [0,1)$ for all $\gamma \in (0,\tilde \gamma]$. This establishes existence of an interval of values for $\gamma$ close to zero such that the gradient-based optimizer is globally convergent. But depending on $\underline{\lambda}_P \underline{\lambda}_H$ and $\overline{\lambda}_P\overline{\lambda}_H$, there may exist larger values of $\gamma \in (0,1]$ with $A(1) \geq 1$ that could frustrate convergence. The following Lemma shows that $\sqrt{A(\gamma)}$ is the global convergence rate of $\theta_k$ to $\hat\theta_n$.
{\bf Proof of Lemma (ref):} As discussed above, there exists $\gamma$ such that $A(\gamma) \in [0,1)$. For such $\gamma$, let $\overline{\gamma}(\underline \lambda_P,\underline \lambda_H,\overline\lambda_K,\overline\lambda_H) \in (0,1]$ independent of $k$ be such that: $A(\gamma) = (1-\overline{\gamma})^2 \in [0,1)$. It follows that
In general, a larger value of $\overline{\gamma}$ would result in faster convergence of $\|\theta_k-\hat\theta_n\|_2$ to zero. The choice of $\gamma$ and the implied $\overline{\gamma}$ in Lemma (ref) are typically data-dependent, but further insights can be gained in two special cases. For gd, the largest globally convergent $\gamma$ is $\underline{\lambda}_H/\overline{\lambda}_H^2$. In ill-conditionned problems when this ratio is small, convergence will be slow since $(1-\overline{\gamma})^2 = (1-[\underline{\lambda}_H/\overline{\lambda}_H]^2)$ will be large. For nr\; when $P(\theta)=H(\theta)^{-1}$, we can use $\overline\lambda_{PH}\leq \overline \lambda_{P}\overline \lambda_{H}$ and $\underline \lambda_{PH} \geq \underline\lambda_{P}\underline \lambda_H$ to obtain a tighter bound. The globally convergent $\gamma$ that minimizes $1-2\gamma \underline\lambda_{PH}+\gamma^2 \overline\lambda^2_{PH}$ is then $\gamma = \underline{\lambda}_{PH}/[\overline{\lambda}_{PH}]^2$ which is strictly less than $1$ for non-quadratic objectives. Since the $(1-\overline{\gamma})^2 = (1-[\underline{\lambda}_{PH}/\overline{\lambda}_{PH}]^2)$ associated with nr\; is typically smaller than for gd, nr\; will converge faster.
Resampling is usually used for inference, but Algorithm (ref) uses resampling for estimation. Unlike classical optimizers, the resampled gradient is noisy. As a consequence, the draws $(\theta_{b})_{b\geq 1}$ constructed by Algorithm (ref) no longer converge deterministically. The following conditions will be imposed on the resampled objective $Q_m^{(b)}$.
Assumption (ref) i. bounds the remainder term in the Taylor expansion of each resampled gradient around the sample minimizer $\hat \theta_n$. Assumption (ref) ii. implies that each resampled objective is also strongly convex. Conditions iii.-v. are tightness condition on the resampled gradient and Hessian empirical process. It implies uniform convergence over $\Theta$ at a $\sqrt{m}$-rate.\footnote{This is implied by a conditional uniform Central Limit Theorem. See VanderVaart1996 and kosorok2007 for iid data. Chen2003 provide high-level conditions for resampling two-step estimators when the first-step estimator can be nonparametric.} Condition iv. is satisfied with $C_3^\prime = 0$ for MLE and NLS estimators because $G_n$ is a sample mean. For over-identified GMM, $\overline{g}_n(\hat\theta_n) \neq 0$ and the gradient $G_n(\hat\theta_n) = 2 \partial_\theta \overline{g}_n(\hat\theta_n)^\prime W_n \overline{g}_n(\hat\theta_n)$ is not a sample mean. Condition iv. requires correct specification in GMM so that $\|\overline{g}_n(\hat\theta_n)\|_2$ goes to zero sufficiently fast as $n \to \infty$.
The following lemma shows that $\theta_b$ will converge in probability to and stays within a $\frac{1}{\sqrt{m}}$ neighborhood of $\hat \theta_n$ as $b$ increases.
{\bf Proof of Lemma (ref):} For any $\theta \in \Theta$, let $\mathbb G_b(\theta) = \sqrt{m} \left( G_m^{(b+1)}(\theta)-G_n(\theta) \right)$. By construction of $\theta_b$, we have $\theta_{b+1}-\hat\theta_n = \theta_{b}-\hat\theta_n-\gamma P_b G_b$. It follows that
Taking the $\|\cdot\|_2$ norm on both sides, applying the triangular inequality and using arguments analogous to Lemma (ref), we have for $\gamma \in (0,1]$ small enough such that the same $A(\gamma)\in [0,1)$:
Taking expectations on both sides:
The desired result is then obtained with $C_5 = \gamma\overline{\lambda}_P C_3$.\qed
Lemma (ref) shows stochastic convergence of $\theta_b$ to $\hat\theta_n$. To study the properties of our estimator $\overline \theta_{\textsc{re}}$, we will use a concept known as {\em coupling}. A coupling between two distributions $\mu$ and $\nu$ on an (unrestricted) common probability space is a pair of random variables $X$ and $Y$ such that $X\sim \nu$ and $Y \sim \mu$, and are equal, on average, up to Wasserstein distance of order $p\geq 1$. Precisely, the Wasserstein-Fr\'{e}chet-Kantorovich coupling distance between two distributions $\nu$ and $\mu$ is defined as: $ W_p(\nu,\mu)^p = \inf_{(X,Y), X \sim \nu, Y \sim \mu} \mathbb{E}(\|X-Y\|^p), p \geq 1. $
Of interest here is the coupling between $\theta_b$ and $\theta_b^\star$, where $\theta_b^\star$ is a linearized sequence of $\theta_b$ defined below. They have different marginal distributions because one is a linear and the other is a non-linear process. Nonetheless, they live on the same probability space because they rely on the same source of randomness originating from the resampled objective $Q_m^{(b)}$. Hence if we can show that $\|\theta_{b}-\theta_{b}^\star\|$ is small in probability, then we can work with the distribution of $\theta_{b}^\star$ which is more tractable.
Precisely, we are interested in a linearized sequence defined as
where $\Psi(\hat\theta_n) = I_d- \gamma \overline P_m H_n(\hat\theta_n)$ and $\overline P_m=I_d$ for rgd\ and $\overline P_m=[H_n]^{-1}$ for rnr. We saw earlier from ((ref)) in the linear regression model that $\Psi(\hat\theta_n)=(1-\gamma)I_d$ for rnr. We now provide conditions on $P_b$ for the draws produced in Algorithm (ref) to be close to those defined in ((ref)) in non-quadratic settings.
Assumption (ref) ii. is needed to ensure that the resampled conditioning matrix $P_b$ converges to $\overline{P}_m$ used in ((ref)) and Assumption (ref) i. ensures stability of the linearized process ((ref)). These assumptions allow us to study $\theta_{b+1}^\star-\hat\theta_n$ as a VAR process with parameters that depend on the Hessian, the conditioning matrix and the learning rate as in the OLS example.
For rgd\; with $P_b=\overline{P}_m=I_d$, Condition ii. holds automatically, while Condition i. requires $\gamma<2/\overline{\lambda}_H$. For rnr\; with $\overline{P}_m = [H_n(\hat\theta_n)]^{-1}$, it will be shown in Theorem (ref) that Conditions i.-ii. hold for any $\gamma \in (0,1]$ such that $A(\gamma) \in [0,1)$ under the assumptions of Lemmas (ref), (ref). This implies that $\theta_{b}^\star$ constructed in ((ref)) for rnr\; is an AR(1) process with autoregressive coefficient $1-\gamma$ as in the OLS case. Now define
The autoregressive structure of $\theta_b^\star-\hat\theta_n$ and $\theta_{b}-\hat\theta_n$ together with the assumed convergence of $P_b$ to $\overline{P}_m$ lead to the following result on the coupling distance between $\theta_{b}$ and $\theta_b^\star$.
The statement holds for any conditioning matrix $P_b$ evaluated on the subsamples satisfying Assumption (ref). Since $\sum_{b=1}^B \overline \rho^b\le \frac{1}{1-\overline{\rho}}$, Lemma (ref) implies
The result is useful because it implies that our estimator $\overline{\theta}_{\textsc{re}}$ equals $\overline{\theta}_{\textsc{re}}^\star = \frac{1}{B}\sum_{b=1}^B \theta_{b}^\star$ up to vanishing terms. By the triangular inequality:
The first term can be bounded by Lemma (ref) as discussed above, and $\mathbb{E}^\star(\theta_{b}^\star) = \hat\theta_n$ by construction of $\theta_{b}^\star$ in ((ref)). Furthermore, Assumption (ref) i. implies that the difference $\overline{\theta}_{\textsc{re}}^\star - \hat\theta_n$ is a $O_{p^\star}(\frac{1}{\sqrt{mB}})$ since $\theta_{b}^\star$ is asymptotically ergodic and its innovations have variance of order $\frac{1}{m}$.
Theorem (ref) says that the average of draws $\overline{\theta}_{\textsc{re}}$ is a consistent estimate of $\hat\theta_n$ for any choice of conditioning satisfying Assumption (ref). The inverse Hessian (rnr) and the identity matrix (rgd) are examples of such conditioning matrices $P_b$.
Theorem (ref) is valid for $P_b$ satisfying the assumptions of the analysis. This subsection specializes to rnr\; produced by Algorithm (ref) which uses the inverse Hessian as conditioning matrix. There are two reasons for this choice. First, it implies a faster decline in the initialization bias compared to e.g. $\overline{P}_m = I_d$ used in gradient descent. Second, such a conditioning matrix has a limit $\overline{P}_m = [H_n(\hat \theta_n)]^{-1}$. For $\overline{P}_m \neq [H_n(\hat \theta_n)]^{-1}$ the dynamics are approximated by a VAR(1) instead of a simple AR(1). While the variance of the AR(1) is proportional to the desired $\mathbb{V}_m$, up to a simple adjustment, this is generally not the case for the VAR(1).
Once Assumption (ref) is granted with $\overline{P}_m = [H_n(\hat \theta_n)]^{-1}$, the general idea of deriving the limiting distribution of the rnr\; draws is to ensure the increasing sum $\theta_{b}^\star = \hat\theta_n - \gamma \sum_{j=0}^{b-1} (1-\gamma)^j H_n(\hat\theta_n)^{-1}G_{m}^{(b-j)}(\hat\theta_n)$ preserves the convergence of each resampled $G_{m}^{(b-j)}(\hat\theta_n)$. The autoregressive nature of $\theta_{b}$ makes the argument somewhat different from the standard setting where each resampled minimizer $\hat\theta_m^{(b)}$ can usually be expressed as a function of a single $G_m^{(b)}(\hat\theta_n)$ plus negligible terms. In such cases, distributional statements about $\hat\theta_m^{(b)}$ follow from the convergence of each resampled $G_m^{(b)}$. Here, the increasing sum $\theta_{b}^\star$ depends on the entire history of the independently resampled $\{G_{m}^{(b-j)}(\hat\theta_n)\}_{j=0,\dots,b-1}$ for which we need to prove convergence.
Assumption (ref) requires non-degeneracy of the variance-covariance matrix which is required for Central Limit Theorems white2014. Assumption (ref) provides higher-order conditions to ensure that the bootstrap converges in distribution at a sufficiently fast rate. It can be understood as requiring the resampled data to have an Edgeworth expansion, the first term being the characteristic function of the standard normal distribution. This occurs with $\beta=1/2$, $\alpha=1$ for averages of iid data with finite third moment lahiri2006. By Assumption (ref), the error in the Gaussian approximation of $\sqrt{n}[H_n(\hat\theta_n)]^{-1}G_m^{(b)}(\hat\theta_n)$ depends on $\alpha$ through $r_m(\tau)$ and on $\beta$ through the inflation factor $1+\frac{r_m(\tau)}{m^\beta}$. These two parameters are of significance because the error in the asymptotic approximation for $\theta_{b}^\star$ inherits the error in $\sqrt{n}[H_n(\hat\theta_n)]^{-1} G_m^{(b)}$. The following theorem takes as given the validity of bootstrap standard errors, i.e. $\sqrt{m} (\mathbb{V}_m)^{-1/2} \Big( \hat\theta_n - \theta^0 \Big) \overset{d}{\to} \mathcal{N}\left( 0, I_d \right)$.
The thrust of the Theorem is that $[H_n(\hat\theta_n)]^{-1} G_m^{(b)}(\hat\theta_n)$ is approximately normal when properly standardized by $(\mathbb{V}_m)^{-1/2}$ and scaled by $\sqrt{m}$. The summation in the AR(1) representation ((ref)) preserves this property under the stated assumption but inflates the variance by a factor $\phi(\gamma)$ which needs to be adjusted. As pointed out above, the error in the Gaussian approximation of $\theta_{b}^\star$ is of the same order as $\sqrt{n}[H_n(\hat\theta_n)]^{-1}G_m^{(b)}(\hat\theta_n)$, which depends on $\alpha,\beta$ according to Assumption (ref). But $r_m(\tau)$ is inflated by a factor of $\frac{(2-\gamma)^\alpha}{1-[1-\gamma]^\alpha}$ which is $1$ when $\gamma=1$ and goes to infinity as $\gamma \to 0$. The Gaussian approximation is better for larger $\gamma$.
An implication of Theorems (ref) and (ref) is that \[ \mathbb V_{r\textsc{nr}}^{-1/2}\sqrt{n}\left( \overline{\theta}_{r\textsc{nr}} - \theta^0 \right) = \mathbb V_{r\textsc{nr}}^{-1/2}\sqrt{n}\left( \hat\theta_n - \theta^0 \right) + o_{p^\star}(1) \overset{d}{\to} \mathcal{N}\left( 0, I_d \right), \] where $\mathbb{V}_{r\textsc{nr}} = \frac{m}{\phi(\gamma)}\text{var}^\star(\theta_{b}-\hat\theta_n)$, and a plug-in estimator $V_{r\textsc{nr}}$ is defined in Algorithm (ref). This implies that standard errors and quantiles computed from the draws $\theta_{b}$, after adjusting for $m$ and $\phi(\gamma)$, can be used to make asymptotically valid inference. Confidence intervals can be constructed to test linear and non-linear hypotheses.
Theorem (ref) specializes to rnr\; where $P_b=[H_m^{(b+1)}(\theta_{b})]^{-1}$. Because the AR(1) representation in ((ref)) does not hold for rgd, simple adjustments cannot be designed that would allow rgd\; to provide valid inference. Furthermore, when $\gamma \in (0,1]$ is small enough such that rgd\, converges, it is not uncommon that $\lambda_{\text{max}}(\Psi(\hat\theta_n)\Psi(\hat\theta_n)^\prime) \simeq 1$ because of ill-conditioning, and when $\theta_{b}$ is very persistent, a much larger $B$ will be required.
Given the Markov chain nature of our $\theta_{b}$, convergence of the chain can be diagnosed using the standard tools from the MCMC literature such as convergence diagnostics considered in gelman1992,Gelman1998. As seen from ((ref)), the draws $\theta_{b}$ approximately follow $d$ univariate AR(1) processes with the same persistence parameter $(1-\gamma)$ which is user-chosen. This can be used to gauge the quality of our large sample approximation in the data for a given pair $(\gamma,m)$. We will illustrate this feature below.
It is noteworthy that while the appeal of stochastic optimization is the savings from using $m \ll n$, our rnr\; requires $m$ not to be too small. This can be seen as the cost of valid inference. Nonetheless, several additional shortcuts could improve the numerical performance of rnr. Our algorithm can be modified so that the Hessian is updated every few iterations rather than at each iteration. The draws would still be valid since the assumptions of Theorem (ref) would still hold. The Hessian could also be approximated using quasi-Newton methods which only requires computing gradients. However, as shown in dennis1977, nocedal-wright:06, the analytical properties of the Hessian approximated by bfgs can only be guaranteed under strong conditions for quadratic objectives. ren1983 show that the bfgs estimate $P_k$ may not converge to the Hessian even for quadratic objectives. Theoretical guarantees can be given for less popular but more tractable methods such as Broyden's method or the Symmetric Rank-1 (SR1) update Conn1991. These, unfortunately, tend to be less stable than bfgs even in classical optimization. Though an extension of Theorems (ref) and (ref) to quasi-Newton updating is left to future work, the results based on a resampled bfgs procedure are promising, as will be seen below.
Our algorithm is related to several other fast bootstrap methods. As discussed in the introduction, solving the minimization problem $B$ times can be computationally challenging or infeasible. Some shortcuts have been proposed to generate bootstrap draws for inference at a lower cost. Davidson1999 (hereafter dmk) proposes a $n$ out of $n$ {\em approximate} bootstrap that replaces non-linear estimation on each batch of re-sampled data by a small number of Newton steps using $\hat\theta_n$ as starting value. In our notation, they perform Newton-Raphson updating $\theta^{(b)}_{\textsc{dmk},j+1} = \theta^{(b)}_{\textsc{dmk},j} - [H_n^{(b)}(\theta^{(b)}_{\textsc{dmk},j})]^{-1} G_n^{(b)}(\theta^{(b)}_{\textsc{dmk},j})$ with $\theta^{(b)}_{\textsc{dmk},0} = \hat\theta_n$ and $j=0,\dots,k-1$ times for each $b=1,\ldots, B$ and report the draws $\theta^{(b)}_{\textsc{dmk},k}$. Armstrong2014 extends this approach for two-step estimation with a finite dimensional or nonparametric first-step estimator. Kline2012 (hereafter, ks) suggests a {\em score bootstrap} that uses random weights to perturb the score while holding the Hessian at the sample estimate. If the random weights are $\{ \omega^{(b)}_i\}$ with $\mathbb E[ \omega_i]=0, \mathbb E[\omega_i^2]=1$, then the distribution $[H_n(\hat\theta_n)]^{-1}\frac{1}{\sqrt{n}} \sum_{i=1}^n \omega^{(b)}_i G_i(\hat\theta_n; y_i,x_i)$ conditional on the data is used to approximate that of $\sqrt{n}(\hat\theta_n-\theta^0)$. The appeal is that the Hessian only needs to be computed once. Honore2017 proposes an approach where the resampled objective is minimized only in a scalar direction for a class of models.
The methods above all rely on a preliminary converged estimate, $\hat\theta_n$ and hence estimation precedes inference. We compute $\overline\theta_{r\textsc{nr}}$ and an estimate of its sampling uncertainty in the same loop, so no further computation is needed once $\overline\theta_{r\textsc{nr}}$ is available. Under our assumptions, the initialization bias will vanish. The practical implication is that for $B$ large enough, the initial values of rnr\ can be far away from the global minimum $\hat \theta_n$, and the algorithm will not be sensitive to the usual stopping criteria used in optimization to find $\hat \theta_n$.
liang-su:19 suggests a `moment-adjusted' algorithm (masgrad) that, in our notation, updates according to $\gamma_b \to 0$ with $P_b = \text{var}(\sqrt{n}G_n(\theta))^{-1/2}$, which is the asymptotic variance-covariance matrix of the sample gradient. In practice, they recommend to evaluate this quantity using the full sample. Under the information matrix equality, we have $\mathbb E[(G_n(\theta_{b}) G_n(\theta_{b})']= \mathbb E [H(\theta_{b})]$ so that the difference amounts to using $H(\theta_{b})^{-1/2}$ instead of $H(\theta_{b})^{-1}$. While such a conditioning matrix would result in consistent estimates, it would not provide asymptotically valid bootstrap draws, which requires $P_b = H(\theta_{b})^{-1}$ and $\gamma$ fixed as shown in our Theorem (ref).
The sgld algorithm proposed in welling-teh:11 updates according to
where $v_{b}\sim N(0,\gamma_b I_d)$ is an injected noise, $\gamma_b$ satisfies ((ref)) and $p(\theta_{b})$ is the prior distribution evaluated at $\theta_{b}$ while $\log p(x^{(b)}_i|\theta_{b})$ is the log-likelihood of a resampled observation $x^{(b)}_i$ evalutated a $\theta_{b}$. The update is thus based on the gradient of the log posterior distribution. Like sgld draws, the draws of our free-lunch bootstrap involve two phases: optimization and sampling. First, in the optimization phase, the shape of the objective function dominates the resampling noise until $\theta_{b}$ attains a neighborhood of $\hat\theta_n$. Then, in the sampling the resampling phase, the noise dominates and rnr\, draws have bootstrap properties. Compared to sgld, the noise is not injected exogenously and our $\gamma$ is fixed. welling-teh:11 shows that with carefully chosen step size $\gamma_b$ and noise variance $\sigma^2_v$, sgld draws can be used for Bayesian inference. Our free-lunch algorithm does not involve any prior and the goal is frequentist inference, as in the Laplace-type inference proposed in Chernozhukov2003 (hereafter CH).
Like CH, our goal is also to simplify the estimation of complex models. CH tackles non-smooth and non-convex objective functions by combining a prior with a transformation of the objective function. In principle, we can also handle non-convex objective functions through regularization, but smoothness is an assumption we need to maintain. CH relies on a Laplace approximation to validate the theory while we use the idea of coupling. By nature of the Metropolis-Hastings algorithm, not all CH draws are accepted and the Markov chain is better described as a threshold autoregressive process. All our draws are accepted and they constitute a nonlinear but smooth autoregressive process. Valid quasi-Bayes inference requires the optimal weighting matrix $W_n=\text{var}(\sqrt{n}(\overline g_n(\hat\theta_n))$ which needs to be estimated. Continuously updating $W_n(\theta)$ can result in local optima so that the MCMC chain can take significantly more time to converge. Whether or not convexification is required, our approach does not require a specific weighting matrix.
In terms of tuning parameters, CH requires as input the proposal distribution in the Metropolis-Hastings algorithm and the associated hyper-parameters. Our tuning parameters are confined to the fixed learning rate $\gamma$ and the resampling size $m$, which do not depend on the dimension of $\theta$. The complexity of the problem also affects the two algorithms in different ways. As seen from Lemma (ref), nr\; converges at a dimension-free linear rate of $(1-\overline{\gamma})$, whereas MCMC converges more slowly as the dimension of $\theta$ increases. For instance, the number of iterations needed for the random walk Metropolis-Hastings to converge increases quadratically with the condition number of the Hessian of the log-density and linearly in the dimension $d$ of $\theta$. To alleviate this issue, several samplers exploit gradient information roberts1996,girolami2011,neal2011,welling-teh:11. While these methods improve upon random walk Metropolis-Hastings, ill-conditioning can still render slow convergence. Scaling the proposal using Hessian information can reduce the effect of ill-conditioning but requires a preliminary estimate. See dwivedi2019log, Table 1, for an overview of mixing times in Metropolis-Hastings algorithms.
This section illustrates the properties of the rnr\; draws using simulated data and data used in published work. Throughout, we use a burn-in period of $\textsc{burn} = 1+\text{round}( \log(0.01)/\log(1-\gamma))$ so that the bias is approximately less than 1% of the initialization error $\|\theta_0-\hat\theta_n\|$. Additional implentation details are given in Appendix (ref). The set of $\gamma$ values satisfying the conditions for Lemma (ref) are data dependent, but in all simulated and empirical examples, $\gamma \in [0.1,0.3]$ performed well.
\paragraph{Example 1: OLS} We simulate data from the linear model with intercept $\beta_0=1$, slope $\beta_1=1$, $x_i \sim \mathcal E(2)$, $e_i \sim t(6)$, $n=200$. We set $B=1000$ plus burn-in draws. Homoskedastic standard errors with a degree of freedom adjustment are computed. Table (ref) reports estimates and standard errors for one simulated sample. We consider three values of batch size $m=200,50,10$ and for each batch size, three values of the learning rate $\gamma$. The results are denoted rnr$_\gamma$ for $\gamma=1, 0.1$, and $0.01$. The smaller the $\gamma$, the more persistent are the draws. Thus $\gamma=0.01$ is a case of extreme persistence, and as seen from the analysis of the linear model, the variance of the draws are larger the smaller $\gamma$ is.
The OLS estimator takes the value $\hat \beta_1 = 1.230$ for this simulated sample. We see from Table (ref) that the rnr\, estimate is very close to OLS when $m=200$ ($=n$) and the choice of $\gamma$ makes little difference. Theorem (ref) suggests that the estimation error should be of order $\frac{1}{\sqrt{m}}$. The large bias associated with a $m$ small is most visible at $m=10$, which is less than $\sqrt{n}$. The difference between the OLS and rnr\, estimates is nearly a third of a standard error for $\gamma=1$. The $m$ out of $n$ Bootstrap and rnr\; standard errors are also less accurate with $m=10$.
\paragraph{ Example 2: MA(1)} Consider the estimation of a MA(1) model by non-linear least squares (nlls). The data is generated as $y_t = \mu + e_t + \psi e_{t-1} $. We set $\mu = 0, \psi = 0.8$, $n = 500$ and $B=2,000$. In this example, $ Q_n(\theta) = \sum_{t=1}^n e_t(\theta)^2 $ where $e_t (\theta)$ are the nlls filtered residuals computed as described in Appendix (ref). In estimation, the gradient and Hessian are computed analytically. For the standard bootstrap, we implement a state-space resampling algorithm described in Appendix (ref). For rnr, we initialize at $\theta_0 = (0,0)$ with a learning rate set to $\gamma=0.6,0.1$ and $0.01$, noting that $\gamma=1$ was too large to get stable results in this example.
In this synthetic data, the nlls estimator is $\hat \psi = 0.816$. Table (ref) shows that when $m=n$, rnr\; produces a $\overline{\theta}_{r\textsc{nr}}$ that is very close to the full sample nlls estimate for all three values of $\gamma$. As in the OLS example above, the bias and standard errors are larger when $m$ is smaller, as suggested by Theorem (ref). The rnr\, standard errors are very similar to those obtained by the $m$ out of $n$ bootstrap for all values of $m$. The standard errors are quite poor for $\gamma=0.01$, most likely because of the strong persistence of the draws.
This subsection considers three examples, the first concerns probit estimation of labor force participation, the second is covariance structure estimation of earnings dynamics, and the third is structural estimation of a BLP model. \paragraph{Application 1: Labor Force Participation} The probit model is of interest because the objective function is strictly convex. To illustrate, we estimate the model for female labor force participation considered in mroz:87. The data consist of $n=753$ observations assumed iid. We set $B=1000$ and $\gamma=0.3$. Three values of $m$ are considered: $m=n$, $200$, $100$. Appendix (ref) provides r code for replicating rnr\; in this example. There are 8 parameters in this exercise and to conserve space, we only report 4 to get a flavor of the results. Table (ref) in the on-line Appendix reports all coefficients. As seen from Table (ref), the rnr\; estimates are close to the MLE ones. Furthermore, the rnr\; standard errors are close to the bootstrap standard errors. Table (ref) also shows results for resampled bfgs which is labeled rqn. Evidently, the r\textsc{qn} estimates are similar to r\textsc{nr}; but is much faster to compute because the Hessian is not computed directly.
Panel (a) of Figure (ref) illustrates the behavior of the draws produced by rnr. The dashed red line corresponds to the MLE estimate $\hat\theta_n$. The black line corresponds to rnr\; draws based on resampling the data with replacement. The blue line shows iterates from classical nr\; with the same $\gamma=0.3$. The top left panel shows the first $20$ draws in the convergence phase when the classical nr\; and the proposed rnr\; should behave similarly. While in this example, rnr\; converges after $5$ draws, \textsc{nr} requires $10$ to $15$ iterations to achieve convergence. The top right panel plots the next $200$ draws. Since convergence is achieved after 5 draws, these draws are in the re-sampling phase. Evidently, the transition between the convergence and the resampling phase of r\textsc{nr}\; is seamless. The AR(1) coefficient on $\theta_{b,educ}$ based on the converged draws (after discarding the first five) is estimated to be $0.673$ with a standard error $0.016$, which is not significantly different from $0.7=1-\gamma$ predicted by Lemma (ref).
Panel (b) of Figure (ref) uses the mroz:87 data to further illustrate Lemma (ref). We compare the rnr\; draws with two AR(1) sequences generated according to coupling theory in ((ref)), ie. $\theta_{b+1}^\star = \hat\theta_n + (1-\gamma)(\theta_{b}^\star-\hat\theta_n) - \gamma H_n^{-1}G_m^{(b+1)}(\hat\theta_n)$ with $\theta^\star_0 = \theta_0$. For $m=n$ shown in the left panel, the coupling result is very accurate after the short initial convergence phase as the two series are nearly indistinguishable. The right panel shows that coupling distance is noticeably greater when $m=100$.
Panel (c) of Figure (ref) illustrates Theorem (ref) by comparing the asymptotic Gaussian distribution with the bootstrap, rnr, dmk\, and ks\, distributions for the education coefficient with $m=n$ and $\gamma=0.3$. The distribution of the rnr\, draws is rescaled using the simple adjustment: $\overline{\theta}_{r\textsc{nr}} + \sqrt{\frac{m}{n\phi(\gamma)}}(\theta_{b}-\overline{\theta}_{r\textsc{nr}})$ after discarding a burn-in period of $10$ draws. The rnr\, distribution approximates the bootstrap distribution quite well.
\paragraph{Application 2: Earnings Dynamics}
moffitt-zhang:18 estimates earnings volatility using a subsample of $3508$ males in the Panel Study of Income Dynamics (PSID) dataset between 1970 and 2014 for a total of $36403$ observations. Let $y_{iat}$ denote an individual's earnings $i$ in age group $a$ (between 24 and 54) at time $t$. Earnings are assumed to be the sum of a permanent $\mu_{ia}$ and a transitory $\nu_{iat}$ component:
In moffitt-zhang:18, the variances are modeled via 11 parameters and estimated by a sequential quadratic programming algorithm (SQP). The Hessian in their example has both positive and negative eigenvalues, suggesting that the solution could be a saddle point. To abstract from identification issues, we estimate $\theta=(\nu_0,\delta_0,\gamma_0,\gamma_1)$ and fix the remaining 7 parameters.\footnote{Specifically, $ \text{var}(\mu_{i,0})$ by $\nu_0$, $ \text{var}(\omega_{ir})$ by $\delta_0, \delta_1$, $ \text{var}(\varepsilon_{ir})$ by $ \gamma_0,\gamma_1,k$, and $ \psi_{a,a-r}$ by $ \pi, \lambda_1, \eta_1,\eta_2,\eta_3$. We set $k=1,\lambda=5$, and all remaining parameters to zero.} We specify the conditioning matrix as $P_b = (H_b^\primeH_b)^{-1/2}$ to ensure positive definiteness. Algorithm (ref) converges using the starting values $\theta_0 = ( 0.054, -10.257, -4.355, 0.012 )$ as in the original paper, with $m=n$, $\gamma=0.2$ $B=2000$, and resampling at the age-cohort level. We also consider re-weighting instead of resampling which is denoted as rnr$_w$. Though our theory does not cover resampled quasi-Newton methods, we also use an implementation of bfgs that sets $P_b = (H_{b,\textsc{bfgs}}^\prime H_{b,\textsc{bfgs}})^{-1/2}$, where $H_{b,\textsc{bfgs}}$ is the bfgs approximation of the Hessian matrix, and report the results as rqn.
Table (ref) shows that the rnr, rnr$_w$ and rqn estimates are very close to $\hat\theta_n$ obtained by SQP. The rnr$_w$\; standard errors are larger than the rnr\; ones, which are in turn larger than the bootstrap ones, but the differences are not enough to change the conclusion that all four parameters are statistically different from zero. However, bootstrap inference of $\hat\theta_n$ is time-consuming, requiring 5h48m to produce 2000 draws, even after the original Matlab code was ported to r and C++ using Rcpp to get a greater than $10$ times speedup in computation time. In contrast, r\textsc{nr}\; produces estimates and standard errors in 1h4m, and the r\textsc{qn} in 38m. The time needed for \textsc{dmk}\ to produce standard errors is comparable to r\textsc{nr}, while r\textsc{qn} is more comparable to \textsc{ks}. However, both have an overhead of having to first obtain $\hat\theta_n$, which entails minimization of the objective function by SQP.
In addition to providing standard errors, an additional by-product of Algorithm (ref) is that the draws can be used for model diagnostics. Gentzkow2017 provides statistics to assess the sensitivity of the parameter estimates to assumptions of the model, taking the data as given. Our algorithm takes the model assumptions as given, but takes advantage of resampling to shed light on the sensitivity of the estimates to features of the data themselves. As pointed out in Chatterjee1986, influential observations could be outliers, or could be points of high leverage. If no such observations exist, removing them in the resampled data should not significantly affect the Markov chain. If their presence is influential, we should witness a `break' in the draws.
With this motivation in mind, we examine whether the parameter estimates of the earnings model are sensitive to data of a particular age group.
Figure (ref) presents results for estimation based on resampled data that exclude one age group at a time. The parameter that is least sensitive to age-groups appears to be $\gamma_1$. The parameter $\nu_0$ tends to be lower when the age group (29-33) is excluded, while $\gamma_0$ is higher when the age group (44-48) is excluded. The parameter most sensitive to age is $\gamma_0$, which is evidently smaller in absolute magnitude when the younger age groups are excluded. For example, it is -1.2 when the youngest age group is dropped, but is around -2.0 when the oldest age group is excluded.
\paragraph{Application 3: Demand for Cereal}
We evaluate the rnr\, algorithm on the BLP model of Berry1995 using the cereal data generated in nevo2000. The sample consists of market shares for $24$ cereal products across $94$ markets for a total of $2256$ market/product observations. This example is of interest because bootstrapping the BLP model is demanding. We use the BLPestimatoR R package which builds on C++ functions to evaluate the GMM objective and analytical gradient. The data consists of market shares for 24 products in 94 markets. In the BLP framework, parameters on terms that enter linearly can projected out by 2SLS. Of the remaining parameters that need to be estimated, we drop some interaction terms from the original paper that may be difficult to identify. This allows us to focus on coefficients that enter the moment conditions non-linearly since the BLP procedure requires a fixed-point inversion for these, making the moment conditions costly to evaluate. The parameter dimension including market fixed effects is $d=33$. To control for possible correlations in the unobservables at the market level, we compute cluster-robust standard errors. See Appendix (ref) for details
The data consists of market shares $s_{gj}$ in market $g \in \{1,\dots,94\}$ for product $j \in \{1,\dots,24\}$ with characteristic matrix $X_{gj}$. To resample at the market level, for each $b=1,\dots,B$ we draw markets $g^{(b)}_{1},\dots,g_{94}^{(b)}$ from $\{1,\dots,94\}$ with replacement and take the associated shares and characteristics $\{ s_{g^{(b)}j},X_{g^{(b)}j}\}_{j=1,\dots,24}$ as observations within each market. Since the number of clusers is relatively small, we only consider $m=n$. We set $\gamma = 0.2$ and a burn-in period of $10$ draws. Additional values of $\gamma$ are considered in Table (ref). Both rnr\, and rqn deliver estimates similar to those obtained from classical optimization and the standard errors are similar to the bootstrap ones. However, there is a significant difference. To generate $B=1000$ draws the standard bootstrap requires 4h45m while the rnr\, runs in 1h04m which is almost 5x faster. The rqn estimate only requires 13m, which is about 20x faster than the bootstrap. While dmk\, is similar to rnr, a preliminary estimate of $\theta$ is needed.\footnote{\textsc{ks} \; is not reported here because it needs significant rewriting of the \textsc{BLPestimatoR} package. Furthermore, \textsc{ks} only consider just-identified GMM models as indicated in their footnote 8.} Estimation of the AR(1) coefficient for the r\textsc{nr}\, draws associated with the parameters reported in Table (ref) finds that the estimates range from $0.78$ to $0.82$ with 95% confidence levels that always include the value of $(1-\gamma)=0.8$ predicted by theory.
{For this example we also consider the CH quasi-Bayesian estimator implemented using a random walk Metropolis-Hastings MCMC algorithm. The prior for each of the $33$ parameters is a $\mathcal{N}(0,100)$ distribution, and the inverse of the clustered variance-covariance matrix of the moments evaluated at $\hat\theta_n$ is used as the optimal weighting matrix $W_n$. The Markov chain is initialized at $\theta_0=\hat\theta_n$, the proposal in the random walk step is scaled by $0.2H_n(\hat\theta_n)^{-1/2}/\sqrt{n}$ which yields an acceptance rate of $0.335$. Though we generate a large number of draws ($B=50000$), the Markov chain is strongly persistent and the effective sample size as defined in gelman2013 is typically less than $100$. The CH estimates are further away from $\hat\theta_n$ than the rnr\; estimates, and the standard errors are also smaller than the $m$ out of $n$ bootstrap and the rnr\; or the dmk.
Simulation-based estimation is routinely used to analyze structural models associated with analytically intractable likelihoods. Estimators in this class include the simulated method of moments, indirect inference and efficient method of moments, and following jjng-15, we will generically refer to them as simulated minimum distance (smd) estimators. We now show that rnr\ will still provide valid inference. This is useful because computing standard errors for the smd estimates is not always straightforward.
The minimum-distance (md) estimator minimizes the distance between a sample auxiliary statistic $\hat\psi_n=\hat\psi_n(\theta^0)$ and it's expected value $\psi(\theta)$ and is defined as \[ \hat\theta_{n,\textsc{md}}=\text{argmin}_\theta\| \overline g_n(\theta)\|^2_{W_n}, \quad \overline g_n(\theta) = \hat\psi_n-\psi(\theta) \quad \tag{\textsc{md}} \] where $W_n$ is a weighting matrix. In cases when the binding function $\psi(\cdot)$ that maps $\theta$ to the auxiliary statistic is tractable, Algorithm (ref) provides a convenient way to compute standard errors for $\hat\theta_{n,\textsc{md}}$. When $\psi(\theta)$ is not tractable, smd simulates data $y_{i,s}(\theta)$ for given $\theta$ using iid errors $e_{i,s}$ and estimates $\psi(\theta)$ by $\hat \psi_{n,S}(\theta)=\frac{1}{S}\sum_{s=1}^S \psi_{n,s}(y_{n,s}(\theta))$. The estimator is \[ \hat\theta_{n,\textsc{smd}}=\text{argmin}_\theta \|\overline g_{n,S}(\theta)\|^2_{W_n} \quad \overline g_{n,S}(\theta)= \hat{\psi}_n - \hat\psi_{n,S}(\theta) \quad \tag{\textsc{smd}} \]
To motivate our simulation based rnr, consider the exactly identified case when it holds that $\hat\psi_n-\psi(\hat\theta_{n,\textsc{md}})=0$. Note that while the vector of auxiliary statistics $\hat\psi_m^{(b)}$ computed using resampled data satisfies $ \mathbb{E}^\star ( \hat\psi_m^{(b)} ) = \hat\psi_n$, the statistics $\hat\psi_{m,S}^{(b)}(\theta)$ computed by smd\; satisfies $ \mathbb{E}^\star [ \hat\psi_{m,S}^{(b)}(\theta)] = \psi(\theta)$, for all $\theta$. The two results together imply that $ \mathbb{E}^\star [ \hat\psi_m^{(b)} - \hat\psi_{m,S}^{(b)}(\theta) ]$ when evaluated at $\theta=\hat\theta_{n,\textsc{md}}$ is $\hat\psi_n-\psi(\hat\theta_{n,\textsc{md}})$ which takes the value of zero as in md\; estimation. This suggests a simulation based resampled objective function defined as: \[ Q_{m,S}^{(b)}(\theta)=\| \hat\psi_m^{(b)} -\hat \psi^{(b)}_{m,S}(\theta) \|^2_{W_n}, \quad \overline g_{m,S}(\theta)=\hat\psi_{m}^{(b)}-\hat\psi_{m,S}^{(b)}(\theta)\quad\tag{r\textsc{nr,s}}\] will have the same minimizer as the infeasible md, at least to a first-order. Let the draws be generated according $\theta_{b+1,S}=\theta_{b,S} - \gamma P_{b+1,S} G_{m,S}^{(b)}(\theta_{b})$ with gradient
By Theorem (ref), the mean $\overline \theta_{r\textsc{nr},S}$ is consistent for $\hat\theta_{n,\textsc{md}}$. By implication, $\theta_{\textsc{re},S}$ will also be more efficient than $\hat\theta_{n,\textsc{smd}}$, which we will verify in simulations below. To analyze $\overline\theta_{\textsc{re},S}$, we need the following:
Assumption (ref).iii' implies Assumption (ref).iii where $G_n(\theta)$ is the gradient of md\; by taking the difference $G_{m,S}^{(b)}(\theta)-G_n(\theta)$ and using the Cauchy-Schwarz inequality.
It remains to construct the variance of $\overline\theta_{r\textsc{nr},S}$. The foregoing analysis would suggest that valid inference would follow after the variance adjustment defined in Algorithm (ref). However, this is not the case. Intuitively, the estimator $\overline{\theta}_{r\textsc{nr},S}$ is consistent for $\hat\theta_{n,\textsc{md}}$ whose variance $\mathbb V^0$ does not involve simulation noise. But the quantity $ V_{r\textsc{nr}}$ defined in Algorithm (ref) presumes the presence of simulation noise in the estimate $\overline{\theta}_{r\textsc{nr},S}$ and will give standard errors that will, in general, be too large. The are many ways to overcome this problem, and most involve running a second chain of draws in parallel with the one used to compute the estimator. For example, taking the difference of two chains with the same simulated samples would work as the simulation noise will offset.
Our preferred approach is to use a second chain that directly estimates the variance of $\overline\theta_{r\textsc{nr},S}$. As in the first chain, this second chain is generated as $\theta^2_{b+1,S}=(1-\gamma)\theta^2_{b,S}-\gamma P_{b+1,S} \tilde{G}_{m,S}^{(b)}(\theta_{b,S}^2)$ as defined in ((ref)), but the gradient is
Compared to the first chain defined by ((ref)), the second chain replaces the simulated auxiliary statistics $\hat \psi_{m,S}^{(b)}(\hat\theta_n)$ by the sample estimates $\hat\psi_{n}$ which is already computed. As all other quantities involved in computing ((ref)) are taken from ((ref)), the computation overhead of generating $\theta^2_{b,S}$ is thus negligible.
jjng-15,jjng-14 shows that a weighted average of smd estimates with independent simulation draws constitutes a posterior mean which is asymptotically equivalent to the infeasible md estimator. This requires solving the optimization problem as many times (ie. $S>1$). Part i. of the proposition shows that this type of statistical efficiency can be achieved by rnr\; in a single run, ie($S=1$). Resampling by rnr\; involves taking draws from the joint distribution $F_n \times F_{\text{shocks}}$ to produce $\hat \psi^{(b)}_{m,S}(\theta)$, which is an estimate of population mapping $\psi(\theta)$. In practice, the simulation and resampling noise in $\hat \psi^{(b)}_{m,S}(\theta)$ and $\hat\psi_m^{(b)}$ are averaged out so that the variance of $\overline{\theta}_{r\textsc{nr},S}$ does not depend on $S$ asymptotically. This contrasts with the smd estimator $\hat \theta_{n,\textsc{smd}}$ which has vanishing simulation noise only when $S\rightarrow\infty$ as $n \to \infty$.
Part ii of the Proposition involves a second sequence $\theta^2_{b,S}$ which, as noted above, is used to compute the variance of the estimator. To understand its underpinnings, recall that the sandwich variance for $\hat\theta_{n,\textsc{md}}$ has a meat component that is the variance of the score $-2\partial_{\theta^\prime} \psi ( \hat \theta_{n,\textsc{md}} ) W_n ( \hat\psi_n-\psi (\hat\theta_{n,\textsc{md}}) )$. If $\psi$ were tractable, a bootstrap draw of this score would be $-2\partial_{\theta^\prime} \psi ( \hat \theta_{n,\textsc{md}} ) W_n ( \hat\psi_m^{(b)}-\psi (\hat\theta_{n,\textsc{md}}) )$. But this is approximately $-2\partial_{\theta^\prime} \hat\psi_{m,S}^{(b)}( \theta_{b,S}^1 ) W_n ( \hat\psi_m^{(b)}-\hat\psi_n )$ which is precisely the gradient ((ref)) used to generate $\theta_{b,S}^2$. Hence it provides a correct approximation of the variance of the scores. Though two chains are needed in the case of simulation estimation, it only needs $S=1$. These arguments are further illustrated using a simple example in Appendix (ref).
\paragraph{Example 3: Dynamic Panel} Consider the dynamic panel regression: \[ y_{it} = \alpha_i + \rho y_{it-1} + x_{it}^\prime\beta +\sigma_e e_{it}, \] with $\rho = 0.6,\beta=1,\sigma_e=1$, $x_{it} \sim \mathcal{N}(0,1), e \sim \mathcal{N}(0,1)$, $n=1000$ and $T=5$. Let $A = I_T - 1_T 1_T^\prime /T$, a matrix which computes the time de-meaned $y_{it} - \overline y_i$. The Least-Squares Dummy Variable (LSDV) estimator is obtained by regressing $A y_{2:T}$ on $A y_{1:T-1}$ and $A x_{2:T}$. The estimator is inconsistent for fixed $T$ as $n \to \infty$.
The LSDV estimator is inconsistent when $n\to\infty$ and $T$ is fixed. Gourieroux2010 shows that indirect inference, which has an automatic bias correction property, is consistent for fixed $T$. The idea is to match the sample LSDV estimator $\hat\psi_n = \hat\theta_{n,LSDV}$ with a simulated $\hat\psi_{n,S}(\theta) = \hat \theta^{sim}_{n,LSDV}(\theta)$ using $S\geq 1$ simulated samples.
To generate rnr\, draws, we resample $(y_{i1},\dots,y_{iT},x_{i1},\dots,x_{iT})_{i=1,\dots,m}$ with replacement over $i$ for given $m$ and compute $\hat\psi_m^{(b)}=\hat \theta_{m,LSDV}^{(b)}$, our resampled moments. Using the new simulation draws $e^{(b+1)}_{it}$ at each $b$, we simulate $S \geq 1$ panels: $y^{(b+1)}_{it,s} = \rho y_{it-1,s} + x_{it,s}^{ (b+1) \prime}\beta +\sigma_e e^{(b+1)}_{it,s}$, for $t=1,\dots,T$ and $i=1,\dots,m$ and compute the simulated moments $\hat\psi_{m,S}^{(b)}=\tilde \theta_{m,LSDV}^{(b)}(\theta_{b})$. An addional moment is needed to estimate $\sigma_e$; we use the standard deviation of the OLS residuals in the LSVD regression. The gradient and Hessian are computed using finite differences. We illustrate with $m=n,100,50$ for $n=1000$.
The LSDV estimate is 0.329 which is significantly downward biased. However, the indirect inference (ind) estimator corrects the bias as shown in Gourieroux2010. The estimate of $0.619$ in Table (ref) for $S=1$ bears this out. The rnr\, estimates are closer to $\theta^0$ than ind for $m=n$, $100$ and is similar for $m=50$. The ind estimates with $S=10$ are very close to the rnr\; estimates obtained over all $\gamma, m$ and $S$ including $S=1$. This implies that rnr\; achieves the efficiency of \textsc{ind} with large $S$ using just $S=1$. The standard errors are smaller than other methods except for $S=10$. Results for $S=2,5$ are reported in Table (ref).
We close the analysis with two remarks about the examples. As noted earlier, an ill-conditioned Hessian can render slow convergence of gradient-based optimizers. The values of $\frac{\lambda_\text{min}(H_n)}{\lambda_\text{max}(H_n)}$ evaluated at $\theta=\hat\theta_n$, are $ 10^{-7},8 \cdot 10^{-4}$ and $7 \cdot 10^{-6}$ for the probit, earnings dynamics, and BLP examples, respectively. Classical gd\ should be slow in converging in these cases, and the applications bear this out. Second, to reinforce the main result that Algorithm (ref) provides valid inference, we evaluate the coverage of rnr\; in all of the simulated examples considered. As seen from Table (ref), rnr\; delivers a 5% size in almost all cases. Details are given in Appendix (ref) of the online supplement.
In this paper, we design two algorithms to produce draws that, upon averaging, is asymptotically equivalent to the full-sample estimate produced by a classical optimizer. By using the inverse Hessian as conditioning matrix, the draws of Algorithm 2 immediately provide valid standard errors for inference, hence a free lunch. In problems that require $S$ simulations to approximate the binding function, our algorithm achieves the level of efficiency of smd with a large $S$, but at the cost of $S=1$. Numerical evaluations show that Algorithm 2 produces accurate estimates and standard errors but runs significantly faster than the conventional bootstrap and most of the `short-cut' methods.
\baselineskip=12.0pt