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.
88,341 characters · 17 sections · 80 citation commands
Estimation and Inference by Stochastic Optimization
JEL Classification: C2, C3
Keywords: Stochastic gradient descent, M-estimation, $m$ out of $n$ and multiplier Bootstrap.
\setcounter{page}{0} \thispagestyle{empty} \baselineskip=18.pt
Many empirical economic analyses involve the process of estimating parameters $\theta$ from a sample of $n$ observations and assessing sampling uncertainty. The typical routine is to first produce a consistent M-estimator $\hat\theta_n$ of the true $\theta^\dagger$ by minimizing the sample objective $Q_n$. Then, standard errors and confidence intervals are computed to perform inferences. Because the sandwich estimator is harder to implement for more complicated models, it is very common to use a bootstrap at that stage. For more computationally intensive nonlinear estimations, it is common to report standard errors based on few bootstrap draws, e.g. 20, 25, or 50.\footnote{71 of all papers published in the American Economic Review between 2012 and 2016 relied on bootstrap inference, 49 of them report bootstrap standard errors. More than one in five of those 49 papers use $100$ bootstrap replications or fewer to compute standard errors. Survey conducted by the author, details available upon request.} The main bottleneck is that each replication requires another estimation on resampled data, which is nontrivial for these complex models. There is a long-standing interest in finding shortcuts to relieve this computational burden without sacrificing too much accuracy. Examples include Davidson1999, Andrews2002, Kline2012, and more recently Honore2017. These methods compute standard errors taking a converged estimate $\hat\theta_n$ as given. As such, estimation always precedes inference.
This paper proposes to combine resampling with optimization to produce estimates, confidence intervals, and standard errors in a single step. This is very close in spirit to Markov-chain Monte Carlo (MCMC) methods that produce draws from a Bayesian posterior to be used for both estimation and inference. Unlike MCMC, the information matrix equality is not required for valid inference. The main loop of the Algorithm is a stochastic optimization routine which evaluates the objective on batches of resampled data with size $m\leq n$. Starting from an initial guess $\theta_0$, one updates $\theta_b$ to $\theta_{b+1}$, $b\geq 0$, using the resampled gradient and inverse Hessian (or an approximation) as a conditioning matrix, with a fixed learning rate $\gamma$. The Algorithm can also be used with reweighted samples of data. For a suitable choice of $\gamma$ and a strictly convex objective $Q_n$, it is shown that the average over $B$ draws $\theta_b$ is equal to $\hat\theta_n$ up to order $1/\text{min}(m,B)$. Even with $m$ and $B \ll n$, the estimates are consistent and asymptotically normal as long as $m$ and $B \gg \sqrt{n}$. The results also show that the distribution, conditional on the sample, of $\sqrt{m}(\theta_b-\hat\theta_n)$ is first-order equivalent to that of $\sqrt{n}(\hat\theta_n-\theta^\dagger)$, up to scale, and standard errors are consistent. Inference results are valid for a resampled Newton-Raphson (rnr) and a resampled quasi-Newton (rqn) algorithm based on a new quasi-Newton update. For other conditioning matrices, the draws can produce consistent estimates but they need not be valid for inference.
The main theoretical insight of the paper is that rnr\, and rqn\, draws are well approximated by a simple AR(1) process when the resampling size $m$ is sufficiently large. The process is centered around the full sample estimate $\hat\theta_n$, and its innovations have variance proportional to the sandwich formula. Hence, after burning-in the first couple of iterations, the distribution of the draws can be used as a bootstrap distribution -- up to a simple adjustment which accounts for the Markov-chain properties of the algorithms. The main results assume the sample and resampled objective functions are smooth and strictly convex. Additional results and illustrations are given in the Supplement for some non-strictly convex settings. In particular, the algorithms can perform well in situations where their non-resampled counterparts are typically challenged and return inconsistent estimates. This is another advantage of combining optimization with resampling.
After introducing the algorithms, a small empirical example illustrates the properties of the methods. Then one simulated and two empirical examples illustrate different ways in which they can be applied in empirical work. First, the Monte-Carlo example considers the estimation of a dynamic discrete choice model with heterogeneous agents. Estimates are significantly biased because of non-linearity, resulting in poor inferences due to size distortions. A bootstrap bias correction would be computationally intensive. rqn, combined with a simple split-panel bias-correction technique, produces accurate estimates and valid inferences in reasonable time. Next, using data from helpman2008, rnr\, and rqn\, are benchmarked against MCMC and stochastic gradient-descent (sgd) on a large-scale probit. rnr\, and rqn\, converge faster and have better mixing so that fewer draws are needed. Allowing for dependence to compute cluster robust standard errors is straightforward, which is not the case with MCMC. Finally, a more complex model is considered. A replication of donaldson2018 shows how a computationally intensive grid search used for both estimation and bootstrapping on a 100-core cluster environment can be replaced with r\textsc{nr}, running in less than six hours on a desktop computer. The time saved allows to explore other model specifications. The results also highlight identification issues which are easily diagnosed with r\textsc{nr}. Additional empirical and simulated examples are given in a companion paper jjng-aerpp:21. Finally, while the results in this paper are targeted at a class of smooth and convex M-estimation problems, the challenges of non-smooth and non-convex estimation constitute a serious impediment to empirical research. Another companion paper provides algorithms and theoretical results for finite-sample generalized method of moments estimation without assuming smoothness or convexity Forneron2022.
\paragraph{Outline of the paper.} Section (ref) begins with a review of classical and stochastic optimization. The proposed algorithms are presented in Section (ref), related methods are discussed. A simple analytical and empirical example highlight the theoretical aspects of the paper. Section (ref) provides the theoretical results. Section (ref) covers the simulated and empirical examples. Appendices (ref)-(ref) provide proofs, sample R codes, and primitive conditions for some of the assumptions. Issues related to non-convexity are covered in Appendix (ref).
Consider minimizing a twice-differentiable and convex sample objective function $Q_n(\theta)$ with respect to parameters $\theta \in \Theta \subset \mathbb{R}^{d_\theta}$: \[ \hat\theta_n = \text{argmin}_{\theta \in \Theta} Q_n(\theta),\quad Q_n(\theta) = \frac{1}{n} \sum_{i=1}^n q(z_i;\theta), \] using a sample of $n$ iid observations $z_i = (y_i,x_i)$. Examples include non-linear least-squares (NLS) for which $q(z_i;\theta) = [y_i - f(x_i;\theta)]^2$ where $f$ is the regression function, known up to the parameters of interest $\theta$. For (pseudo) maximum likelihood estimation (MLE), $q(z_i;\theta) = -\ell(y_i;\theta|x_i)$ is the negative of the log-likelihood. The true value, denoted by $\theta^\dagger$, minimizes the limiting objective $Q = \text{plim}_{n\to \infty} Q_n$. Under the conditions in Theorem 2.1 of newey-mcfadden-handbook, $\hat\theta_n$ is consistent for $\theta^\dagger$. If, in addition, the assumptions in Theorem 3.1 of newey-mcfadden-handbook hold, then $\hat\theta_n$ is also asymptotically normal: \[\sqrt{n}V_n^{-1/2}(\hat\theta_n-\theta^\dagger)\overset{d}{\to} N(0,I_d), \] where $V_n =[H_n(\hat\theta_n)]^{-1} \Sigma_n [H_n(\hat\theta_n)]^{-1}$ is a consistent estimator of the sandwich variance, $[H_n(\hat\theta_n)]^{-1} \overset{p}{\to} [H(\theta^\dagger)]^{-1}$ approximates the bread, and $\Sigma_n \overset{p}{\to} \lim_{n\to\infty}\text{var}[\sqrt{n}G_n(\theta^\dagger)]$ approximates the meat. Practitionners are often interested in conducting inference for a function of the parameters $h(\theta) \in \mathbb{R}$. Examples include inference on a single coefficient or a counterfactual based on the structural estimates $\hat\theta_n$. When $h$ is continuously differentiable, the delta-method gives the standard error formula $\text{se}[h(\hat\theta_n)] = \sqrt{ \nabla h(\hat\theta_n) V_n \nabla h(\hat\theta_n)^\prime /n}$ and the $95\%$ level confidence interval $h(\hat\theta_n) \pm 1.96 \times \text{se}[h(\hat\theta_n)]$. It is also common to compute confidence intervals from bootstrap draws using quantiles or the standard deviation of $h(\theta_n^{(b)})$, where $\hat\theta_n^{(b)}$ minimizes the resampled objectives $Q_n^{(b)}$.
For most problems, the estimator $\hat\theta_n$ does not have closed-form and is usually obtained by numerical optimization. The methods below rely on the sample gradient and Hessian of $Q_n$, defined respectively by:
When $Q_n$ is quadratic, as in OLS or IV regressions, the solution has closed form and $\hat\theta_n = \theta_0-[H_n(\theta_0)]^{-1} G_n(\theta_0)$, for any initial guess $\theta_0$. For general non-quadratic objectives, the solution does not have closed-form. A typical strategy is to consider a sequence of approximations:
where $P_k$ is a symmetric positive definite conditioning matrix, $\gamma_k \in (0,1]$ is called the learning rate, it penalizes for the quality of the approximation ((ref)). For $\gamma_k=1$ and $P_k^{-1} = H_n(\theta_k)$, the approximation is exact for $Q_n$ quadratic. Otherwise, it coincides with a second-order Taylor expansion around $\theta=\theta_k$, which is not exact. From a guess $\theta_k$, minimizing ((ref)) leads to:
Different choices of conditioning matrix $P_k$ lead to different algorithms. Newton-Raphson (nr) iterations use $P_k = [H_n(\theta_k)]^{-1}$ and $\gamma_k=1$; a damped update uses $\gamma_k<1$. For strongly convex problems, nr\, enjoys a fast quadratic convergence property, but it comes at the cost of computing both the gradient and the Hessian at each iteration $k$.\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.} To avoid the cost of computing the $(d_\theta+1)d_\theta/2$ second-order derivatives in $H_n$. Other updates rely on different choices of $P_k$. These include quasi-Newton (qn) and gradient-descent (gd).
The most popular quasi-Newton update is bfgs for which $P_k$ approximates $[H_n(\theta_k)]^{-1}$ iteratively using only gradients. Although not as fast as nr, convergence for quasi-Newton with an appropriate choice of $\gamma_k$ can be superlinear.\footnote{See dennis1977 and nocedal-wright:06 for an overview of quasi-Newton updates and their properties.} Quasi-Newton methods bypass Hessian computations, but analytical convergence results are more difficult to obtain. The gradient-descent algorithm (gd) sets $P_k=I_d$, it only requires computing a single gradient $G_n(\theta_k)$. Convergence, however, requires choosing $\gamma_k$ small enough that all eigenvalues of $I_d-\gamma_k H_n(\tilde \theta_k)$ lie in the restricted region $(-1,1)$, for some intermediate value $\tilde\theta_k$. This results in very slow convergence when the ratio of the largest to smallest eigenvalue of $H_n$ is large.
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. Note that for non-convex objectives, any $\theta_k$ with $G_n(\theta_k)=0$ is a fixed-point of the update rule ((ref)). This includes local and global minima and maxima as well as locally suboptimal solutions, called saddle points, for which $H_n$ has both positive and negative eigenvalues.
In their seminal paper on stochastic approximation, Robbins1951 consider the situation where only noisy gradients are available $G_n(\theta_{k})+e_k$ with $e_k$ iid, mean-zero. Using the notation above, they propose the gradient-descent update : $ \theta_{k+1}=\theta_{k}- \gamma_{k} [G(\theta_{k})+e_k]$, and show under regularity conditions that $\theta_k \overset{a.s.}{\to} \hat\theta_n$ when $\gamma_k > 0$ satisfy
The first condition ensures that all possible solutions will be reached with high probability regardless of the $\theta_0$ while the second ensures convergence to the true value. Kiefer1952Stochastic extend the results to non-linear least squares with gradients computed by finite-differences. Feasible choices of learning rate schedules include $\gamma_k = \gamma_0 k^{-\delta}$, with $\delta \in (1/2,1]$ and $\gamma_0>0$ as choice parameters. Depending on $\delta$, convergence as measured by $\mathbb{E}(\|\theta_k-\hat\theta_n\|^2)$ can occur at a $1/k$ rate or slower. ruppert1988efficient and polyak1992acceleration show that the averaged values $\overline{\theta}_k = \frac{1}{k}\sum_{i=1}^k \theta_i$, converge at the fastest $1/k$ rate for all choices of $\delta \in (1/2,1]$. This is commonly refered to as Polyak-Ruppert averaging. Building on these results, modern stochastic gradient-descent updates are given by
where $G_m(\theta_k)=\frac{1}{m}\sum_{i=1}^m \nabla q(z_i^{(k)},\theta_k)$ is an unbiased estimate of $G_n(\theta_k)$ based on $m<n$ randomly chosen observations $z_1^{(k)},\dots,z_m^{(k)}$ drawn at each $k$. Setting $e_k = G_m(\theta_k) - G_n(\theta_k)$ results in the same setting as stochastic approximation above. Though $m=1$ is computationally inexpensive and is a popular choice, the stochastic gradients $G_m(\theta_k)$ are very noisy and a small $\gamma_k$ satisfying ((ref)) is needed to compensate. This results in slow convergence making it well suited for problems where the gradients $G_m$ are very cheap to computed and running very many iterations is feasible. Unfortunately, this is typically not the case for more complicated economic models. A number of methods can accelerate convergence,\footnote{Improvements to sgd\, include momentum polyak:64 and accelerated gradient-descent nesterov:83. As for gd, using $P_k=I_d$ in sgd\, can lead to very slow convergence. Popular extensions of stochastic gradient descent include adagrad duchi-acm:11, and adam kingma2014.} but they tend to converge slower than nr, or \textsc{bfgs}. More closely related to this paper, some study stochastic \textsc{bfgs}, but compute the full sample $G_n$ many times which is more demanding.\footnote{See e.g. Algorithm 1, step 4, in pmlr-v51-moritz16. Note that stochastic Newton-Raphson (s\textsc{nr}), with a small $m =1$, is not used in practice because $H_m^{(k)}$ is often not invertible when $m<d_\theta$.}
Beyond reducing the computional cost, stochastic optimization can improve on classical methods in non-convex settings. For instance, suboptimal solutions aren't fixed points for sgd\, ge2015escaping,jin2017escape. Though sgd\, requires $\gamma_k \to 0$, it is common to use $\gamma_k = \gamma$ fixed in practice. dieuleveut-aos:20 consider $Q_n$ quadratic with $\gamma$ and $m$ fixed and find that it leads to estimation bias which needs to be corrected. Few results are available for inference, notably chen-aos:20 consider the computation of standard errors for sgd. This paper considers $\gamma$ fixed with $Q_n$ non-quadratic but requires $m \to \infty$. For inference, $P_k$ is either a Hessian or the new quasi-Newton update. The connection with bootstrap inference and the AR(1) representation are new.
Because the iterations are used both for estimation and bootstrap inference, they will be indexed by $b=0,1,2,\dots$, rather than $k=0,1,2,\dots$ used for estimation only. The main algorithm below relies on an resampled objective $Q_m^{(b)}$, and its derivatives $G_m^{(b)},H_m^{(b)}$, computed using either $m \leq n$ resampled observations, or using all $m=n$ observations but with random reweighting. Unless otherwise stated, $Q_m^{(b)}$ will refer to either a resampled or reweighted objective. Using these quantities, Algorithm (ref) describes a way to perform both estimation and inference for a function $h$ of $\theta$.
Similar to MCMC, after discarding an initial burn-in period, estimates are computed by averaging over draws. Also, standard errors and confidence intervals are computed from the standard deviation and quantiles of the draws but here with an adjustment $\sqrt{\frac{m}{n\phi(\gamma)}}$ which accounts for $m$ potentially smaller than $n$, and $\phi(\gamma)$ adjusts for the Markov-chain properties of $\theta_b$ as explained below. Compared to sgd, Algorithm (ref) requires $m \gg 1$ but allows for $m/n \to 0$ with the restriction that $\sqrt{n}/m \to 0$. This is the cost of using a fixed learning rate $\gamma$ rather than a decaying sequence $\gamma_k$ satisfying ((ref)), used in nearly all of the literature.
The learning rate $\gamma$ controls the optimization behaviour of the Algorithm, a larger value of $\gamma$ is associated with faster convergence. For quadratic $Q_n$, any $\gamma \in (0,1]$ is feasible, otherwise a smaller $\gamma < 1$ is needed, $\gamma \in [0.1,0.3]$ performs well in the Examples. The parameter $\gamma$ also affects the distributional properties of $\theta_b$ which can be more concentrated around $\hat\theta_n$ than conventional bootstrap draws. For instance, with $m=n$ and $\gamma=0.1$ we have $1/\phi(\gamma) = 19$: the variance of $\theta_b$ is $19$ times smaller than the variance of $\hat\theta_n$.\footnote{$\phi$ is decreasing in $\gamma$: the smaller $\gamma$, the more localized the Markov-chain will be. For comparison $\phi(1)=1$, $\phi(0.4)=4$, $\phi(0.2)=9$, and $\phi(0.01)=199$.} A number of different implementations are possible, a resampled Gradient-Descent (rgd) would use the same $P_b = I_d$ for each iteration $b$. Valid inference, however, requires implementations for which $P_b$ approximates the inverse Hessian sufficiently well. This is the case for the resampled Newton-Raphson and resampled quasi-Newton algorithms. The first is conceptually simple and provides good intuition for the results. For most applications, the second is more computationally efficient and will be preferred.
Depending on the choice of resampling or reweighting scheme, $G_m^{(b)},H_m^{(b)}$ will be computed differently. First, the resampled objective, gradient, and Hessian are:
where $z_1^{(b)},\dots,z_m^{(b)}$ are drawn independently with replacement from the empirical data $(z_1,\dots,z_n)$ at each iteration $b$. Unlike sgd, the resampling should match the dependence structure of the data. For instance, clustered data should be resampled at the cluster rather than individual level. Second, the re-weighted objective, gradient, and Hessian are:
where $w_i^{(b)}$ are iid with mean and variance equal to $1$. Examples include Gaussian weights $w_i^{(b)} \sim \mathcal{N}(1,1)$, exponential $w_i^{(b)} \sim \exp(1)$, or Poisson $w_i^{(b)} \sim \text{Pois}(1)$. This is known as a multiplier, or exchangeably weighted, bootstrap.\footnote{See VanderVaart1996 and kosorok2007. The Poisson bootstrap approximates, with computational advantages, the standard $n$ out of $n$ bootstrap chamandy2012.} Similar to resampling, the weights should vary only between clusters in clustered data, and more generally match the dependence of the data. Depending on the setting, reweighting can be preferred to resampling. For instance, if $\theta$ includes fixed effects, resampling can produce samples with no observation that is informative about one of the fixed effects.\footnote{This is also an issue for censored models where resampling can result in datasets with no censored observations. This also affects sgd; with $m=1$ a gradient is only informative about at one fixed effect at a time. This makes convergence very slow, see Appendix (ref) for an illustration.}
The resampled Newton-Raphson update described in Algorithm (ref) below can be used within Algorithm (ref) to produce consistent estimates and valid inferences.
The implementation is very simple, as it only requires computing a Hessian. In many settings, however, repeatedly computing this matrix of second-order derivatives can be very demanding so that a quasi-Newton update is often preferred. Nevertheless, the rnr\, implementation gives a good first intuition about why and how Algorithm (ref) produces estimates and valid inferences. This is illustrated with a simple analytical and a small empirical example below.
\paragraph{A pen and pencil example.} The following illustrates the main optimization and bootstrap properties of rnr\, using a tractable OLS regression of $y_i$ on $x_i$: $Q_n(\theta) = \frac{1}{2n} \sum_{i=1}^n (y_i - x_i^\prime \theta)^2$. Using $m$ out of $n$ resampling, we have $G_m^{(b+1)}(\theta) = -\frac{1}{m} \sum_{i=1}^m x_i^{(b+1)} (y_i^{(b+1)} - x_i^{(b+1)\prime} \theta)$, and $H_m^{(b)}(\theta) = \frac{1}{m} \sum_{i=1}^m x_i^{(b+1)} x_i^{(b+1)\prime}$. In more compact notation: $G_m^{(b+1)}(\theta) = - \frac{1}{m} X_m^{(b+1)\prime} ( y_m^{(b+1)} - X_m^{(b+1)}\theta)$, $H_m^{(b+1)}(\theta) = \frac{1}{m} X_m^{(b+1)\prime} X_m^{(b+1)}$. For $b \geq 0$, one rnr\, update takes the form:\footnote{Derivations for rnr\, and other methods are given in Appendix (ref).} \[ \theta_{b+1} = \theta_b + \gamma \left( X_m^{(b+1)\prime} X_m^{(b+1)} \right)^{-1} X_m^{(b+1)\prime} ( y_m^{(b+1)} - X_m^{(b+1)}\theta_b). \] Notice that $\left( X_m^{(b+1)\prime} X_m^{(b+1)} \right)^{-1} X_m^{(b+1)\prime} y_m^{(b+1)} = \hat\theta_m^{(b+1)}$ is the bootstrap OLS estimate. Substract $\hat\theta_n$ on both sides and re-arrange terms to find that: \[ \theta_{b+1} -\hat\theta_n = (1-\gamma)(\theta_b-\hat\theta_n) + \gamma ( \hat\theta_m^{(b+1)}-\hat\theta_n) \] is an AR(1) process with conditionally iid innovations based on the centered bootstrap estimates: $\hat\theta_m^{(b)}-\hat\theta_n$. Notice that neither $\hat\theta_n$ nor $\hat\theta_m^{(b)}$ are computed explicitly in the rnr\, update. Using a backwards recursion, we have: \[ \theta_{b+1} -\hat\theta_n = \underbrace{(1-\gamma)^b(\theta_0-\hat\theta_n) \vphantom{\gamma \sum_{j=0}^{b} (1-\gamma)^j ( \hat\theta_m^{(b-j)}-\hat\theta_n)}}_{\text{Starting value bias}} + \underbrace{\gamma \sum_{j=0}^{b} (1-\gamma)^j ( \hat\theta_m^{(b+1-j)}-\hat\theta_n)}_{\text{Resampling noise}} \] Suppose, for simplicity, that under the bootstrap expectation $\mathbb{E}^\star(\hat\theta_m^{(b+1)})=\hat\theta_n$. Then the starting value bias due to $\theta_0 \neq \hat\theta_n$ is \[ \mathbb{E}^\star(\theta_{b+1}-\hat\theta_n) = (1-\gamma)^b(\theta_0-\hat\theta_n),\] which converges exponentially fast to zero. Very quickly, the bias becomes negligible and the average draw consistently estimates $\hat\theta_n$. A seamless transition occurs from optimization to resampling. Using the formula for the variance of an AR(1) process, we have \[ \text{var}^\star ( \theta_{b+1} ) = \gamma^2\frac{1- (1-\gamma)^{b+1}}{1-(1-\gamma)^2} \text{var}^\star(\hat\theta_m^{(b)}) \simeq \phi(\gamma) \text{var}^\star ( \hat\theta_{m}^{(b)} ),\] which also converges exponentially fast. The factor $\phi(\gamma) = \gamma^2/[1-(1-\gamma)^2]$ is used in Algorithm (ref). Assuming $\text{var}^\star ( \hat\theta_{m}^{(b)} )$ consistently estimates standard errors, then so does rnr\, with the suggested adjustment. The AR(1) representation only holds if $P_b$ approximates the inverse Hessian sufficiently well, this is not the case for gd\, and several quasi-Newton updates.
\paragraph{Illustration with real data.} The main results extend the properties described for OLS above to more general non-quadratic strictly convex objectives. The main insight is that even though the AR(1) representation above does not hold exactly, there exists a sequence $\theta_b^\star$, called coupling, which follows an AR(1) process and is uniformly close to the Markov-chain $\theta_b$ up to an error of order $1/m$. A good approximation requires $m \to \infty$ but still allows for $m/n \to 0$.
To illustrate this feature, Figure (ref) shows the draws for the coefficient on education for a probit model fitted on data from mroz:87. This example uses $\gamma=0.3$ and resamples $m=n/2$ out of $n=753$ observations with replacement. The top left panel shows the first $20$ iterations, which corresponds to the initial convergence phase. The resampled rnr\, (black) and the deterministic nr\, (blue) have very similar paths, moving quickly to the true value (dashed black). Then, the top right panel shows the next $200$ iterations -- the resampling phase. Even with $m<n$, the rnr\, draws (black) are very close to the AR(1) process $\theta_b^\star$ (blue) used to prove the results. The bottom panel illustrates the bootstrap properties of the draws. The left panel compares a $\mathcal{N}(\hat\theta_n,V_n/n)$ with the standard $n$ out of $n$ bootstrap, and rnr\, draws adjusted for $m$ and $\phi(\gamma)$. The distributions are very close. The right panel further confirms this with a more detailed q-q plot view. Replication code for this example is provided in Appendix (ref) with a detailed comparison of rnr, rqn, for different choices of $m$, with standard errors computed using the sandwich formula, the standard bootstrap, and methods by Davidson1999, Kline2012. Additional comparisons with these methods can be found in jjng-aerpp:21.
In practice, quasi-Newton methods are often preferred to Newton-Raphson because the latter requires computing the full Hessian matrix many times, which can be costly.\footnote{In some settings, it is possible to compute the Hessian `for free' once the gradient is computed; see Example 3, Section (ref) and Appendix (ref).} The main idea behind the quasi-Newton update in Algorithm (ref) below is that the $d_\theta(d_\theta+1)/2$ unique elements in the Hessian can be recovered using only $d_\theta$ scalar derivatives. Consider the Hessian-vector product $y_b = H_m^{(b+1)}(\theta_b)s_b$, for a direction $s_b \neq 0$. Notice that $y_b = \lim_{\varepsilon \to 0} [G_m^{(b+1)}(\theta_b + \varepsilon s_b)-G_m^{(b+1)}(\theta_b - \varepsilon s_b)]/2\varepsilon$ is a scalar derivative which can be computed without the full matrix $H_m^{(b+1)}(\theta_b)$. The full Hessian matrix can be recovered from Hessian-vector products computed in $d_\theta$ linearly independent directions $s_b,s_{b-1},\dots$ by running a simple linear regression of $y_b$ on $s_b$. Using an overdetermined system with $L > d_\theta$ directions ensures the linear regression is well conditionned. In practice, the different $y_b$ are computed at different values of $\theta_b$ so that the Hessian recovered by least-squares need not be symmetric nor positive definite. A simple transformation is used to enforce these two features, while retaining its approximation properties. The resulting quasi-Newton update is new. The while loop ensures the regression is well-conditionned. It never ran in the examples of Section (ref) but should included as a safeguard. A more detailed discussion of difference with existing updates and the choice of tuning parameters is deferred to Section (ref).
To generate $B$ draws, the standard bootstrap with resampling or reweighting requires solving the minimization problem $\hat\theta_m^{(b)} = \text{argmin}_{\theta \in \Theta} Q_m^{(b)}(\theta)$ $B$ times which can be computationally intensive. Faster alternatives have been proposed in the literature. Davidson1999 consider a $k$-step nr\, update: $\hat\theta_{n,k+1}^{(b)} = \hat\theta_{n,k}^{(b)} - [H_n^{(b)}(\hat\theta_{n,k}^{(b)})]^{-1}G_n^{(b)}(\hat\theta_{n,k}^{(b)})$ starting at the same full sample estimate $\hat\theta_{n,0}^{(b)} = \hat\theta_n$ for all $b$. Andrews2002 shows that $k=1$ is enough for asymptotically valid inference, and $k>1$ provides asymptotic refinements. dmk\, also propose a $k$-step quasi-Newton update which requires $k>1$ steps.\footnote{Following earlier work by robinson1988, the results are based on quadratic and super-linear convergence properties of nr\, and bfgs, respectively, when $\gamma=1$ which holds only if $\|\hat\theta_{n}^{(b)}-\hat\theta_n\|$ is sufficiently small and if (and only if) the quasi-Newton matrix $P_k$ satisfies certain properties dennis1977.} In comparison, using the same $m=n$, a $k=1$-step implementation of dmk\, has the same cost as r\textsc{nr}, but is higher than r\textsc{nr}\, for $k>1$ steps. The $k$-step quasi-Newton update requires $k>1$ making it $k$ times more costly than rq\textsc{n}\, for $m=n$. Using $m<n$ further reduces the cost of both estimation and inference for r\textsc{nr} , rq\textsc{n}. The original idea of \textsc{dmk}\, has been extended to fast subsampling hong2006, two-step estimation Armstrong2014, and $\ell_1$-penalized estimation li2021.
Building on hu1995 and hu2000, Kline2012 propose a wild score bootstrap which can be written as a one-step nr\, update: $\hat\theta_{n}^{(b)} = \hat\theta_{n} - [H_n(\hat\theta_{n})]^{-1}G_n^{(b)}(\hat\theta_{n})$ where $G_n^{(b)}$ is reweighted using de-meaned weights. If both the estimation and using $m=n$ are not too demanding, this approach is computationally attractive. fastwild:19 point out some caveats with this approach. By only evaluating $G_m^{(b)}$ at $\hat\theta_n$, it does not capture sampling variation in $H_n(\hat\theta_{n})$ or non-linearities in the objective. The latter also applies to dmk\, with $k=1$. By design, rnr\, and rqn\, reflect variation in $H_m^{(b)},G_m^{(b)}$ around $\hat\theta_n$. Notice that if the optimizer returns a suboptimal solution $\tilde\theta_n$, then dmk\, with $k>1$ may not be centered around $\tilde\theta_n$ which is a red flag. By de-meaning the weights, \textsc{ks}\, produces draws centered around $\tilde\theta_n$ so that the issue goes unnoticed.\footnote{Appendix (ref), Table (ref), illustrates failed optimizations for a dynamic discrete choice model, in comparison rq\textsc{n}\, produces accurate estimates. It would be problematic to apply \textsc{dmk}\, or \textsc{ks}\, to the inconsistent estimates.} Honore2017 propose to approximate the sandwich $V_n$ using $ d_\theta(d_\theta+1)/2 \times B$ scalar minimizations of $Q_n^{(b)}$. When $d_\theta$ increases, this quickly becomes taxing.\footnote{In Section (ref), Example 2, $55945 \times B$ scalar optimizations would be needed which is quite large.} Using a preliminary estimate $\hat\theta_n$ to resample is the defining characteristic of these methods compared to r\textsc{nr}\, and rq\textsc{n}\, which never compute $\hat\theta_n$ explicitly.
For likelihood problems, it is common to use Bayesian inference and deploy the MCMC toolkit. Unlike frequentist estimation, Bayesian analyses rely on sampling from a posterior distribution rather than optimization. Well known samplers include Gibbs and Metropolis-Hastings. Because they can be slow to converge, gradient-based samplers are increasingly popular with Metropolis-adjusted Langevin dynamics (MALA), Hamiltonian Monte-Carlo (HMC), or stochastic gradient Langevin dynamics (SGLD) algorithms. A comparison between rnr, rqn, and MALA is given in Section (ref). Unlike Bayesian inference, this paper does not require the information matrix equality for valid inference. Convergence diagnostics used for MCMC algorithms such as trace plots or formal tests gelman1992,Gelman1998,gelman2013 can also be used to monitor convergence of rnr\, and rqn. Lastly, while the detailed balance condition guarantees the posterior is a stationary solution of a Metropolis-Hastings algorithm, convergence rate results typically assume a (nearly) concave log-posterior -- similar to the convexity assumption in this paper.\footnote{Recall that the detailed balance condition ensures that the random walk Metropolis-Hastings algorithm is ergodic but deriving the rate of convergence is more difficult. See e.g. mengersen1996, brooks1998, belloni2009.} jjng-16 show how optimization can be used to sample from a Bayesian posterior with an intractable likelihood. Here, optimization is used for resampling.
This section has three parts. The first subsection presents the main assumptions on $Q_n$, $Q_m^{(b)}$, $P_b$, and derives convergence results for classical optimizers and Algorithm (ref) as well as the pivotal coupling result with the AR(1) process $\theta_b^\star$. These results are non-asymptotic. The second subsection focuses on large sample estimation and inference. The theorems can be applied to any implementation of Algorithm (ref) for which $P_b$ has the required properties. The third subsection specializes them to rnr\, and rqn\, described in Algorithms (ref) and (ref).
The econometric theory for extremum estimation in e.g. newey-mcfadden-handbook takes as a given that an optimizer finds the unique solution. But the econometric assumptions are not enough to ensure a numerical optimizer will provide a globally convergent solution to the sample problem $Q_n(\theta)$. The following is standard in gradient-based optimization theory:
Condition i. implies that $Q_n$ is strongly convex with probability $1$ while condition ii. requires Lipschitz continuity of the Hessian.\footnote{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)$.} Strong convexity of the sample objective function can be less restrictive than it appears as some non-convex problems can be convexified by penalization.\footnote{In some cases, the strict convexity assumption restricts the parameter space to be bounded in the absence of penalization. See e.g. the discussion of Assumption (H3) in Bach2011.} Assumption (ref) implies the so-called Polyak-\Lojasiewicz inequalities:
where $\tilde{\theta}_n$ is an intermediate value, between $\theta$ and $\hat{\theta}_n$. The first inequality follows from the positive definiteness of $H_n(\hat\theta_n)$. To see the usefulness of Assumption (ref), note that for any symmetric positive definite conditioning matrix $P_k$:
where the last inequality follows from Assumption (ref) and $0<\underline{\lambda}_P \leq \overline{\lambda}_P < \infty$ bound the eigenvalues of $P_k$ (see Assumption (ref) iii). Now a contraction occurs if the choice of $\gamma$ is such that $A(\gamma)\in [0,1)$. To see that this is feasible, note that $A(0)=1$ and $\partial _\gamma A(0)<0$. Thus by continuity and local monotonicity of $A(\cdot)$, there exists a $\gamma^\star \in(0,1]$ such that $A(\gamma)\in[0,1)$ for all $\gamma\in(0,\gamma^\star]$ as desired. For such choice of $\gamma$, define $\overline{\gamma}(\underline \lambda_P,\underline \lambda_H,\overline\lambda_P,\overline\lambda_H) \in (0,1]$ independent of $k$ be such that: $A(\gamma) = (1-\overline{\gamma})^2 \in [0,1)$. It follows that:
Global convergence follows from iterating on the contraction at each $k$. The derivations above rely on the positive definiteness of both $H_n$ and $P_k$. A non-definite $P_k$ can be problematic for convergence. While Assumption (ref) i-ii require a smooth and strictly convex objective function, Assumption (ref) iii restricts the choice of conditioning matrix.
The finite upper bound in Assumption (ref) iii is particularly important to ensure the optimizer is well behaved. It is automatically satisfied for gd\, since $P_b = I_d$ but requires the lower bound in Assumption (ref) i for nr. For quasi-Newton, it needs to be enforced as in Algorithm (ref). Note that for nr, with an adjustement similar to the $\tau_b$ in Algorithm (ref), it possible to get convergence even when $H_n$ is singular away from $\hat\theta_n$.\footnote{This is illustrated in Appendix (ref). Note that modifications to the resampled hessian is also present in gonccalveswhite2005 when deriving consistency of bootstrap standard errors in OLS regressions.}
The first Lemma derives the deterministic global convergence of a classical gradient-based optimizer for an appropriate choice of $\gamma$ at a dimension-free linear rate of $(1-\overline{\gamma})$. Convergence of Newton-Raphson with $P_k=H_k^{-1}$ is implied by the Lemma under Assumption (ref) only. This first result is useful for deriving the stochastic convergence of the resampled Algorithm (ref). In particular, notice that the resampled update is: \[ \theta_{b+1} = \underbrace{ \vphantom{\theta_b - \gamma P_b G_n(\theta_b)} \vphantom{\gamma P_b [ G_n(\theta_b) - G_m^{(b+1)}(\theta_b) ]}\theta_b - \gamma P_b G_m^{(b+1)}(\theta_b)}_{\text{Resampled update}} = \underbrace{ \vphantom{\gamma P_b [ G_n(\theta_b) - G_m^{(b+1)}(\theta_b) ]}\vphantom{\theta_b - \gamma P_b G_m^{(b+1)}(\theta_b)} \theta_b - \gamma P_b G_n(\theta_b)}_{\text{Deterministic update}} + \underbrace{ \vphantom{\theta_b - \gamma P_b G_n(\theta_b)} \vphantom{\theta_b - \gamma P_b G_m^{(b+1)}(\theta_b)} \gamma P_b [ G_n(\theta_b) - G_m^{(b+1)}(\theta_b) ]}_{\text{Resampling noise}}, \] which is the sum of the deterministic update studied in Lemma (ref) with resampling noise. The convergence has a deterministic and a stochastic component. The deterministic part is covered by Assumption (ref) and the condition on $\gamma$ in Lemma (ref). The following assumption will be used to study the stochastic part of the update.
These high-level conditions will be used to prove convergence of $\theta_{b}$ in $L_p$-norm, $p \geq 2$. Appendix (ref) provides sufficient conditions for i-ii. using either $m$ out of $n$ resampling with replacement or reweighting with Gaussian multiplier weights.\footnote{Dependence of $C_2,C_3$ on $n$ is omitted to simplify notation. The conditions under Gaussian reweighting are the same for $p=2$ and $p>2$ but not with resampling. The results can be extended to other weighting distributions with the same conditions as resampling by adapting the proof of Theorem 3.6.13 in VanderVaart1996 similar to cheng2015.} Consistent estimation and valid inference using quantiles in Algorithm (ref) (steps 4 and 5), will require $p=2$ in Assumption (ref). Consistency of standard errors further requires $p=4$. Consistency of moments for the standard bootstrap is also associated with particular moment conditions, see gonccalveswhite2005, Kato:11, cheng2015. Since the estimator in Algorithm (ref) is computed by averaging draws, the derivations are fundamentally moment-based.
Assumption (ref) iv reiterates the importance of the bounds for $P_b$ for convergence but in the stochastic case. Under Assumptions (ref), (ref), and the same choice of $\gamma$ used in Lemma (ref): \[ \|\theta_{b+1}-\hat\theta_n\|_2 \leq (1-\overline{\gamma})\|\theta_{b}-\hat\theta_n\|_2 + \gamma \overline{\lambda}_P \sup_{\theta \in \Theta} \| G_m^{(b+1)}(\theta)-G_n(\theta) \|_2, \] with probability 1. Iterating the recursion reveals the general behaviour of $\theta_b$: \[ \|\theta_{b+1}-\hat\theta_n\|_2 \leq (1-\overline{\gamma})^b\|\theta_{0}-\hat\theta_n\|_2 + \gamma \overline{\lambda}_P \sum_{j=0}^{b-1} (1-\overline{\gamma})^j\sup_{\theta \in \Theta} \| G_m^{(b+1-j)}(\theta)-G_n(\theta) \|_2, \] with probability 1. The leading term converges deterministically and exponentially fast to zero. It corresponds to a starting values bias. When Assumption (ref) i holds for some $p \geq 2$ and $(1-\overline{\gamma})\in [0,1)$, the second term has bounded $L_p$-norm, of order $1/\sqrt{m}$. This characterizes the stochastic convergence of $\theta_b$ as show in the following Lemma.
Lemma (ref) shows that $\theta_b$ converges in $L_p$-norm to a $1/\sqrt{m}$-neighborhood of $\hat\theta_n$. Unlike classical optimizers or sgd\, with $\gamma_k \to 0$ where $\theta_k \to \hat\theta_n$ either deterministically or in probability, here $\theta_{b} \not\to \hat\theta_n$ as $b$ increases -- this is similar to MCMC which produces draws around $\hat\theta_n$. Except for quadratic objectives, the Markov-process followed by $\theta_b$ around $\hat\theta_n$ is intractable. Nevertheless, the following will show that
gives a good approximation of $\theta_b$ and its distribution. The autoregressive coefficient $\Psi_n$ and $\overline{P}_m$ depend on the choice of update $P_b$. Table (ref) makes them explicit for rgd, rnr, and rqn.
In the literature, $\theta_{b}^\star$ is known as a coupling of $\theta_{b}$. The two have different marginal distributions because one follows a linear and the other a non-linear process. Nonetheless, they live on the same probability space because they are generated from the same resampled objective $Q_m^{(b)}$. Hence if we can show that $\|\overline{\theta}_{r\textsc{nr}}-\overline{\theta}_{r\textsc{nr}}^\star\|$ is negligible, we can work with the distribution of $\theta_{b}^\star$ which is more tractable.
The constant $\overline{\gamma} \in (0,1]$ is the global rate of convergence in Lemma (ref). Putting $p=2$, the lemma states that that the expected deviation $\mathbb{E}^\star(\|\theta_{b}-\theta_{b}^\star\|)$ is of order $1/m$ plus a term that vanishes exponentially fast in $b$. This in turn implies that for $p=2$, upon averaging, $\mathbb{E}^\star(\|\overline{\theta}_{\textsc{re}}-\overline{\theta}_{\textsc{re}}^\star\|_2)$ is of order $1/m$ plus a $1/B$ term due to $d_{0,n}$. Using these non-asymptotic bounds, the large sample estimation and inference properties of Algorithm (ref) can be derived.
Proposition (ref) is pivotal for what follows. Suppose for simplicity that there is no approximation error, i.e. $\theta_b=\theta_b^\star$, then: \[ \mathbb{E}^\star(\theta_b) = \hat\theta_n + [\Psi_n]^b (\theta_0-\hat\theta_n), \] converges to $\hat\theta_n$ exponentially fast under Assumption (ref) i. This implies that rgd, rnr, and rqn\, deliver consistent estimates. This should be of interest for practitionners using sgd\, with $\gamma_k = \gamma$ fixed and Polyak-Ruppert averaging, which is common in practice, since the theory gives restrictions on $m$ for the bias to be negligible. Furthermore, $\theta_b^\star$ being an ergodic first-order vector autoregressive process: \[ m \times \text{var}^\star(\theta_b) \to \gamma^2 \Psi_n \text{var}^\star[\sqrt{m}G_m^{b}(\hat\theta_n)] \Psi_n^\prime , \] converging exponentially fast. It consistently estimates the asymptotic variance, up to scale, when $\text{var}^\star[\sqrt{m}G_m^{b}(\hat\theta_n)]$ consistently estimates the meat and $\Psi_n$ approximates the bread $[H_n(\hat\theta_n)]^{-1}$ (up to scale). This is the case for rnr\, and rqn, but not r\textsc{gd}. Suppose $\sqrt{m}G_m^{b}(\hat\theta_n) \sim \mathcal{N}(0,\Sigma_n)$ then, for large $b$ we approximately have $\theta_b \sim \mathcal{N}(\hat\theta_n,\gamma^2 \Psi_n \Sigma_n \Psi_n^\prime/n)$ which matches the asymptotic distribution, after re-scaling the variance, if $\Psi_n$ approximates the bread (up to scale). The following Theorem gives a non-asympotic bound for the approximation error $\|\overline{\theta}_{\textsc{re}}-\hat\theta_n\|_2$ and the large-sample consistency result.
First-order asymptotic equivalence of $\overline{\theta}_{\textsc{re}}$ and $\hat\theta_n$ requires both $B$ and $m$ to be large relative to $\sqrt{n}$. Using fixed $m$ and $\gamma$ does not produce consistent estimates. This is relevant for sgd\, with $m\ll n$ and $\gamma_k=\gamma$. For inference, consider the standard bootstrap as a benchmark. Let $\overset{d^\star}{\to}$ denote convergence in distribution conditional on $(z_1,\dots,z_n)$. Under regularity conditions, an $n$ out of $n$ bootstrap estimate $\hat \theta_n^{(b)} =\text{argmin}_{\theta \in \Theta} Q_n^{(b)}(\theta)$ is such that
If $\sqrt{n}G_n^{(b)}(\hat\theta_n)$ satisfies a conditional CLT then the bootstrap distribution targets the correct asymptotic distribution. The bootstrap is often used to compute standard errors. Consistency of moments further requires uniform integrability.\footnote{See billingsley2013. Specifically, to consistently estimate the variance $V_n$, the condition requires $\sup_{n} \mathbb{E}^\star(\|\sqrt{n}(\hat \theta_n^{(b)} - \hat\theta_n)\|^{2+\varepsilon})$ finite for some $\varepsilon >0$. See gonccalveswhite2005, Kato:11, cheng2015, for applications to moment convergence of the bootstrap in OLS regressions and M-estimation.} Here, convergence in distribution of $\theta_b$ comes from Proposition (ref) and a conditional CLT for $\theta_b^\star$, i.e. the growing sum $\sum_{j=0}^b (1-\gamma)^j [H_n(\hat\theta_n)]^{-1}\sqrt{m}G_m^{(b)}(\hat\theta_n)$. Assumption (ref) is used to ensure that distributional convergence for $\sqrt{m}G_m^{(b)}(\hat\theta_n)$ extends to $\theta_b^\star$. Primitive conditions are given in Appendix (ref) for $m$ out of $n$ resampling, and reweighting with Gaussian weights.
Assumption (ref) is standard. It requires asymptotic normality of $\hat\theta_n$, consistency of the sandwich estimator, and validity of the delta method for inference on a function $h$ of $\theta$. Theorem (ref) presents the large sample inference. First, the asympotic validity quantile-based confidence intervals. Second, consistency of standard errors computed in Algorithm (ref). Both results require a choice of $P_b$ satisfying $\overline{P}_m = [H_n(\hat\theta_n)]^{-1}$.
Having the general estimation and inference results for Algorithm (ref), the following verifies the implementation-specific Assumptions (ref) iv, (ref) to show that Theorems (ref)-(ref) apply to rnr\, and rqn. For rqn, it is assumed, without loss of generality, that $\underline{\lambda} = \underline{\lambda}_H/2$.\footnote{The proof only relies on $\underline{\lambda} < \lambda{\text{min}}(H_n(\theta))$ around $\hat\theta_n$.}
Consistency of standard errors implies that they can be used to compute t-statistics, confidence intervals of the form $h(\overline{\theta}_{\textsc{re}}) \pm 1.96 \sqrt{m/[n\phi(\gamma)] \widehat{\text{var}}[h(\theta_b)]}$, where $\widehat{\text{var}}[h(\theta_b)]$ is the sample variance of $h(\theta_b)$, after discarging the first $\textsc{burn} \gg \log(m)$ iterations. A simple rule-of-thumb to gauge the order of magnitude for $\textsc{burn}$ is that $(1-\overline{\gamma}) \geq (1-\gamma)$ and the initial value bias is of order $(1-\overline{\gamma})^b d_{0,n}$. Suppose, we target a bias of order $d_{0,n}/100$, then we need $b > \log(1/100)/\log(1-\overline{\gamma}) \geq \log(1/100)/\log(1-\gamma) \simeq 43$ for $\gamma=0.1$. This implies that the burn-in should include at least $40-50$ draws. A burn-in of $225$ draws corresponds to setting $\overline{\gamma}=0.02$ in the same calculations.
The quasi-Newton update proposed in Algorithm (ref) differs from bfgs and other methods in the literature to ensure the following three properties are satisfied: i) $P_b$ is symmetric, and ii) positive definite with strictly positive eigenvalues, iii) Assumption (ref) holds when $\|\theta_b - \hat\theta_n\|$ is small. Conditions i-ii) are required to ensure optimization is stable and apply Lemma (ref). Condition iii) is needed for valid inference. Other methods do not satisfy i), ii), and iii) simultaneously. Amongst the most widely used methods, bfgs enforces i) and ii) plus the so-called secant equation: $P_b y_{k} = s_k$ so that $P_b H_n(\theta_k)s_{k} = s_k$ always holds but is only guaranteed in this direction so that iii) need not hold for non-quadratic $Q_n$.\footnote{See ren1983 for counter-examples. The BFGS matrix $P_k$ can have negative eigenvalues when $Q_n$ is non-convex. This is problematic for convergence as illustrated in Section (ref).} The sr1 update, used with trust-region algorithms, satisfies i) and iii) but not ii) which makes it unstable in estimation using ((ref)).\footnote{See fiacco1990, Conn1991 for derivations of iii) with SR1.} Both methods update $P_k$ from $P_{k-1}$ using a rank-one update based on the secant equation above. Here $P_b$ is fully updated using OLS, this is related to the multi-secant update in schnabel:83 which can satisfy iii) but typically not i-ii), and requires $L \leq d_{\theta}$.
In the proof of the Corollary, the OLS estimate $\hat H_b$ is shown to satisfy iii) but not necessarily i)-ii). The transformation $\hat H_b^\prime \hat H_b$ enforces i), and the regularization $\tau_b$ enforces ii).\footnote{These two transformations can also be useful for nr\, and rnr\, as illustrated in Section (ref).} Algorithm (ref) requires three inputs. At least $L \geq d_\theta$ secant updates are needed to compute $\hat H_b$ by OLS; $L \geq \text{max}(25,1.5 \times d_\theta)$ works well in the examples. The cutoff $\underline{\lambda}_S$ ensures the least-squares problem is well-conditionned. With $\underline{\lambda}_S=10^{-6}$, the while loop never runs in the applications, but it should be included as a safeguard. The requirement for $\underline{\lambda}$ is to be small enough to satisfy $0<\underline{\lambda} < \lambda_{\text{min}}(H_n)$ around $\hat\theta_n$.
The following provides one Monte-Carlo and two empirical examples which illustrate the properties of rnr\, and rqn\, for estimation and inference.
The first example considers likelihood estimation and inference in a single-agent dynamic discrete choice model with unobserved heterogeneity.\footnote{See e.g. aguirregabiria2010 and arcidiacono-ellickson:11 for reviews on the estimation of dynamic discrete choice models.} It is common in this class of models to use a parametric mixture distribution to model heterogeneity since it is parsimonious relative to a more computationally demanding non-linear fixed effect estimation. There are two challenges in this example. First, the likelihood is non-convex because of the unobserved mixture component. Second, the non-linear transformation used to integrate out the unobserved heterogeneity results in finite sample bias and sizable distortion for inference. Analytical bias correction can be challenging to implement, because the terms are not easily tractable. Standard bootstrap bias correction would be computationally intensive. Using a large-T linearization of the bias, it is shown that the split panel jackknife of dhaene2015 can reduce this bias, and straightforward to use with rqn.\footnote{See Appendix (ref) for the derivations under simplifying assumptions.} A comparison of rnr\, and rqn\, with MLE and the standard bootstrap for a simpler homogeneous agents model are given in Appendix (ref), Table (ref).
The setting is similar to the seminar model of rust:87, each agent $i$ solves a dynamic programming problem at each period $t$ with a flow utility for the set of actions $a_{it}=0$ or $1$: \[ U_{i0}(x_{it},a_{it-1}) = \varepsilon_{it}(0), \quad U_{i1}(x_{it},a_{it-1}) = \beta_{0i} + \beta_{1i}x_{it} -(1-a_{it-1})\delta_1 + \varepsilon_{it}(1), \] where $\delta_1$ is the entry cost, $\varepsilon_{it}(0)$ and $\varepsilon_{it}(1)$ follow an extreme-value type I distribution, $x_{it}$ is an exogenous regressor which has finite support and follows a first-order Markov process with transition probability $\Pi$. Both $a_{it}$ and $x_{it}$ are observed for $i=1,\dots,n$ and $t=1,\dots,T$. Here the heterogeneity is modeled on the intercept and slope parameters of the utility function. The intercept $\beta_{0i}$ equals $\mu_{0}^1$ with probability $\omega \in [0,1]$, and $\mu_{0}^2$ with probability $1-\omega$. The slope is continuously distributed $\beta_{1i} \sim \omega \mathcal{N}(\mu_1,\sigma_1^2) + (1-\omega)\mathcal{N}(\mu_2,\sigma_2^2)$. The parameters of interest are $\theta=(\mu_0^1,\mu_1^1,\mu_0^2,\mu_1^2,\sigma_1,\sigma_2,\omega,\delta_1)$, the discount factor $\rho=0.98$ is fixed. The matrix $\Pi$ is estimated by sample average and taken as input in the log-likelihood below: \[ \ell_{nT}(\theta,\Pi) = \frac{1}{nT}\sum_{i=1}^n \log \left( \int \exp[\sum_{t=1}^T \ell_{it}(\theta,\Pi,\beta_i)]f(\beta_i|\theta)d\beta_i \right), \] where $\ell_{it}(\theta,\Pi,\beta_i) = \log[\mathbb{P}(a_{it}|x_{it},a_{it-1},\theta,\Pi,\beta_i)]$ is the conditional probability of choice $a_{it} \in \{0,1\}$ computed by solving the dynamic programming problem by fixed point iterations, and $f(\cdot|\theta)$ is the mixture distribution for $\beta_i = (\beta_{0i},\beta_{1i})$.
The split panel jacknife is implemented using a full sample implementation of rqn\, $\theta_{b,nT}$ and two half panel implementations $\theta_{b,nT/2}^1,\theta_{b,nT/2}^2$ based on the first and last $T/2$ time-observations each. All three use the same $m=n/2$ out of $n$ resampled observations at each $b$. Bias-corrected draws are given by $\tilde{\theta}_{b,nT}=2\theta_{b,nT}-[\theta_{b,nT/2}^1+\theta_{b,nT/2}^2]/2$. Table (ref) summarizes the results from $400$ Monte Carlo replications. Full panel rqn\, estimates (denoted as rqn) are close to $\theta^\dagger$. Some coefficients have a bias comparable in magnitude to sampling uncertainty. Inference for these coefficient is far from the nominal $5\%$ size using either quantiles or standard errors. The split panel correction (rqn-bc) reduces bias by an order of magnitude and bias-corrected inference is much closer to nominal size. Table (ref) provides additional results with $T=50$ and a comparison with R's bfgs optimizer. In that configuration, bfgs estimates are very inaccurate, and far less reliable than rq\textsc{n}.
It is common to use MCMC methods for large-scale estimation problems. The second example compares the properties of rnr\, and rqn\, with a Metropolis-Hastings algorithm called MALA (Metropolis-adjusted Langevin Algorithm). The idea is to compare the convergence and mixing properties of the draws when there are hundreds of parameters to estimate.\footnote{A description of and implementation details for MALA are given in the Appendix.}
The model is a simple Probit using the data from helpman2008. There are 10 regressors plus 324 exporter and importer fixed effects for a total of 334 parameters to be estimated using 248,060 country/year observations. rnr\, and rqn\, are implemented with $\gamma = 0.1$ by re-weighting with Gaussian multiplier weights. Implementation details for the different algorithms are given in the Appendix. While rqn\, accepts all draws by design, MALA only accepts 57% of them which makes the Markov-chain more persistent as illustrated in Figure (ref) below. The top panel shows the faster convergence rnr\, and rqn. After discarding the initial burn-in, the autocorrelation coefficients are $0.889$ for rnr\, and rq\textsc{n}; $0.995$ for MALA, which indicates better mixing.\footnote{Note that $0.995^{25} \simeq 0.889$ this is why $25$ times more draws are generated for MCMC than rq\textsc{n}.} The bottom panel shows that r\textsc{nr}\, and rq\textsc{n}\, draws are more localized around $\hat\theta_{n}$ than MCMC. Indeed, $\gamma=0.1$ implies $\phi(\gamma)=19$: the variance of $\theta^{rq\textsc{n}}_b$ is 19 times smaller than $\theta^{\text{\textsc{mcmc}}}_b$.
Table (ref) in the Appendix further reports estimates and standard errors for the 10 regressors and the intercept. MLE and its standard errors are computed using R's glm function and vcovHC for robust standard errors. Results for rqn\, are close to MLE for all 11 coefficients. Bayesian posterior standard deviations tend to underestimate. Clustered re-weighting with exponential weights are also reported for comparison, it only requires modifying the line of code generating the weights. In terms of computation, the likelihood, gradient, and hessian are computed in C++ using Rcpp for performance. R's built-in glm takes about 2min for estimation using iteratively reweighted least-squares, rqn\, 4min, rnr\, 10min, and MCMC 1h30min. As for most applications of gradient-based optimization on large datasets with many parameters, a good implementation of the code computing the gradient and Hessian is crutial for fast computation: unoptimized code runs is about 2 hours for rq\textsc{n}, and 1.5 days for MCMC. Results for s\textsc{gd}\, and an infeasible implementation of s\textsc{nr}, which relies on the full sample $H_n(\hat\theta_n)$, are reported in the Appendix. Even after 5 million iterations ($>2$ hours), s\textsc{gd}\, is still far from the solution $\hat\theta_n$.
The third and final example replicates an estimation in donaldson2018 of the transportation costs of salt in India between 1861-1930. The model is a simple NLS regression using the log-price of salt of type $o$ at destination $d$ in year $t$: \[ \log(\text{price}_{dt}) = \beta_0 + \delta \log( \text{LCRED}(R_t,\alpha)_{odt} ) + \text{controls} + \varepsilon_{odt}. \] The main regressor $\text{LCRED}(R_t,\alpha)_{odt}$, measures the lowest-cost route effective distance between origin $o$ and destination $d$ in year $t$ given the transportation network $R_t$ and the transportation costs $\alpha = (\alpha_{\text{rail}},\alpha_{\text{road}},\alpha_{\text{river}},\alpha_{\text{coast}})$ for each mode of transportation for the salt products.\footnote{The normalization $\alpha_{\text{rail}}=1$ is used so that all costs are relative to rail.} The regressor $\text{LCRED}(R_t,\alpha)_{odt}$ is not observed directly. For a given value $\alpha$, it is constructed using Dijkstra’s shortest-path algorithm applied to each triplet $o,d,t$ to find the most cost-effective route along the transportation network $R_t$. The parameters of interest are transportation costs $\alpha$, inferred from the regression and the price elasticity $\delta$.
This particular example is interesting for two reasons. First, the model is fairly difficult and time-consuming to estimate. donaldson2018 uses a grid search for both estimation and bootstrap inference in a 100 core cluster environment.\footnote{The estimation relied on a coarse grid $(\alpha_{\text{road}},\alpha_{\text{river}},\alpha_{\text{coast}}) \in [1,10]^3$ with step size $0.125$ for a total of 373248 grid points, and bootstrap inference was conducted using 200 replications on a coarser grid with step size 0.5 totalling in 5832 grid points.} rnr\, runs in less than 6 hours on a desktop computer. Second, the model is potentially unidentified for several values. This is the case if $\delta=0$, or for $\delta \neq 0$ when a transportation mode systematically dominates the others, e.g. if $\alpha_{\text{road}}$ is small and $\alpha_{\text{river}},\alpha_{\text{coast}}$ are large. The specific values for which this occurs depend on the structure of the transportation network $R_t$. A quadratic penalty is added to the objective: $\text{pen}(\theta) = \lambda ( \|\alpha-\overline{\alpha}\|_2^2 + \delta^2)$, with $\lambda = 0.1/n$ and $\overline{\alpha} = (\overline{\alpha}_{\text{road}},\overline{\alpha}_{\text{river}},\overline{\alpha}_{\text{coast}}) = (4.5,3.0,2.25)$ are observed historical relative freight rate estimates reported in donaldson2018. Notice that if e.g. $\alpha_{\text{river}}$ is unidentified, then $\partial_{\alpha_{\text{river}}}Q_m^{(b)}(\theta)=0$ around $\hat\theta_n$ so that only the penalty determines the path of $\alpha_{b,\text{river}}$, and, after convergence, the draws are degenerate $\alpha_{b,\text{river}}=\overline{\alpha}_{\text{river}}$. With a penalty, degeneracy is indicative that the coefficient is unindentified. This is not the case for a Bayesian posterior.
Figure (ref) shows the draws and the adjusted rnr\, distribution for $(\alpha,\delta)$. Baseline OLS results evaluated for $\alpha= \overline{\alpha}$ are reported with bootstrap standard errors clustered at the district level -- indexed by $d$. There are also $142$ origin/destination fixed effects, a linear time trend and the intercept, totalling in $147$ parameters to be estimated using 7345 observations.\footnote{donaldson2018 includes origin/time fixed effects instead of a linear time trend but the design matrix is singular with these controls in the baseline OLS regression. Hence, they are replaced with the linear trend.} rnr\, converges quickly in about 75 iterations. 1000 draws are produced in 5h40m on a desktop computer with 7000 calls to Dijkstra’s algorithm, the main computational bottleneck.\footnote{For reference, a grid search using a small set of 19 equally spaced points for each $\alpha_j$ requires $19^3=6859$ calls to the algorithm.} The bottom panel shows the distribution of the draws after discarding the first 75 iterations, and adjusting for $\phi(\gamma)$. The distributions for coast and road are degenerate, suggesting they are not identified. The results are not sensitive to the choice of $\theta_0$ and $\lambda$. Table (ref) reports estimates and standard errors with and without fixed effects. The price elasticity is very similar to the one estimated by OLS when fixing $\alpha=\overline{\alpha}$. The estimated coefficient $\alpha_{\text{river}}$ is less than half the historical average $\overline{\alpha}_{\text{river}}$. The null that river transportation is more cost effective than rail, $\alpha_{\text{river}}<1$, cannot be rejected at the 5% significance level.
This paper proposes Algorithms that produce estimates, bootstrap confidence intervals and standard errors in a single run. While the theory relies on convexity of the sample and resample objectives, Appendix (ref) discusses implementation and convergence in some non-convex settings. In Examples 1 and 3, rnr\, and rqn\, perform well without convexity. Extending the Algorithms to handle general non-convex problems would be a useful avenue for future research. Forneron2022 shows how to modify a local gradient-based algorithm to design a fast globally convergent one without using multiple starts. This could also be very useful for sampling, as considered here.
\baselineskip=12.0pt