$\displaystyle=\frac{\sqrt{m}}{\sqrt{\phi(\gamma)}}V_{n}^{-1/2}(1-\gamma)^{b}(\theta_{0}-\hat{\theta}_{n})$
$\displaystyle\quad-\gamma\sum_{j=0}^{b-1}(1-\gamma)^{j}\frac{\sqrt{m}}{\sqrt{\phi(\gamma)}}V_{n}^{-1/2}[H_{n}(\hat{\theta}_{n})]^{-1}G_{m}^{(b-j)}(\hat{\theta}_{n}).$

The first term is a $o_{p}^{\star}(1)$ as long as $\log(m)/b\to 0$. For the
second term, since the
$[H_{n}(\hat{\theta}_{n})]^{-1}G_{m}^{(b-j)}(\hat{\theta}_{n})$ are
conditionally independent and identically distributed, we have by a
convolution argument:

$\displaystyle\mathbb{E}^{\star}\left(\exp(\mathbf{i}\tau^{\prime}\sqrt{m}\gamma\sum_{j=0}^{b-1}(1-\gamma)^{j}\frac{\sqrt{m}}{\sqrt{\phi(\gamma)}}V_{n}^{-1/2}[H_{n}(\hat{\theta}_{n})]^{-1}G_{m}^{(b-j)}(\hat{\theta}_{n}))\right)$
$\displaystyle=\prod_{j=0}^{b-1}\mathbb{E}^{\star}\left(\exp(\mathbf{i}\tau^{\prime}\sqrt{m}\gamma(1-\gamma)^{j}\frac{\sqrt{m}}{\sqrt{\phi(\gamma)}}V_{n}^{-1/2}[H_{n}(\hat{\theta}_{n})]^{-1}G_{m}^{(b-j)}(\hat{\theta}_{n}))\right)$
$\displaystyle=\prod_{j=0}^{b-1}\left[\exp\left(-\frac{\|\tau\|_{2}^{2}}{2}\frac{\gamma^{2}(1-\gamma)^{2j}}{\phi(\gamma)}\right)\left(1+\frac{r_{m}(\gamma(1-\gamma)^{j}\tau/\phi(\gamma))}{m^{\beta}}\right)\right]$
$\displaystyle=\underbrace{\vphantom{\prod_{j=0}^{b-1}\left[\left(1+\frac{r_{m}(\gamma(1-\gamma)^{j}\tau/\phi(\gamma))}{m^{\beta}}\right)\right]}\exp\left(-\frac{\|\tau\|_{2}^{2}}{2}\frac{\gamma^{2}[1-(1-\gamma)^{2b}]}{[1-(1-\gamma)^{2}]\phi(\gamma)}\right)}_{=\exp(-\|\tau\|_{2}^{2}/2\times[1+o(1)])}\underbrace{\prod_{j=0}^{b-1}\left[\left(1+\frac{r_{m}(\gamma(1-\gamma)^{j}\tau/\phi(\gamma))}{m^{\beta}}\right)\right]}_{(\text{error})}.$

To show that the last product is convergent under the stated assumptions, take
logs and use the inequality $\frac{x}{1+x}\leq\log(1+x)\leq x$ for $x>-1$.
Then

$\displaystyle\log\left(\|(\text{error})\|\right)$
$\displaystyle=\sum_{j=0}^{b-1}\log\left(1+\frac{|r_{m}(\gamma(1-\gamma)^{j}\tau/\phi(\gamma))|}{m^{\beta}}\right)\leq\sum_{j=0}^{b-1}\frac{|r_{m}(\gamma(1-\gamma)^{j}\tau/\phi(\gamma))|}{m^{\beta}}$
$\displaystyle\leq\sum_{j=0}^{b-1}\frac{\|\gamma\tau/\phi(\gamma)\|^{\kappa}(1-\gamma)^{\kappa
j}}{m^{\beta}}\leq\frac{\|\gamma\tau/\phi(\gamma)\|^{\kappa}}{[1-(1-\gamma)^{\kappa}]m^{\beta}}\to
0.$

Note that $\frac{\gamma}{\phi(\gamma)}=2-\gamma\geq 1$ for $\gamma\in(0,1]$.
Putting everything together we have:

$\displaystyle\mathbb{E}^{\star}\left(\exp(\mathbf{i}\tau^{\prime}\frac{\sqrt{m}}{\sqrt{\phi(\gamma)}}V_{n}^{-1/2}(\theta_{b}-\hat{\theta}_{n})\right)=\exp\left(-\frac{\|\tau\|_{2}^{2}}{2}\right)\left(1+O\left(\frac{\|\tau\|^{\kappa}}{m^{\beta}}\frac{(2-\gamma)^{\kappa}}{[1-(1-\gamma)^{\kappa}]}\right)\right).$

Now since convergence of the characteristic function implies weak convergence
(Billingsley, 2013, Th26.3, p349), which implies the desired convergence in
distribution.∎

##  Appendix E Implementing rnr in R

To illustrate how the rnr is implemented in a real data setting, the following
provides detailed commented R code for a probit model using the Mroz (1987)
data.

#### Illustration using Labor Force Participation

The table below presents the estimates and standard errors for all methods and
coefficients in the Mroz (1987) application.

Table E4: Labor Force Participation: Estimates and Standard Errors | Estimates
---|---
| mle |  |  |  | rnrn | rnr200 | rnr100 | rqnn | rqn200 | rqn100
nwifeinc | -0.012 | - | - | - | -0.012 | -0.013 | -0.014 | -0.012 | -0.011 | -0.012
educ | 0.131 | - | - | - | 0.132 | 0.138 | 0.143 | 0.131 | 0.129 | 0.129
exper | 0.123 | - | - | - | 0.123 | 0.124 | 0.123 | 0.123 | 0.124 | 0.125
exper2 | -0.002 | - | - | - | -0.002 | -0.002 | -0.002 | -0.002 | -0.002 | -0.002
age | -0.053 | - | - | - | -0.053 | -0.053 | -0.055 | -0.052 | -0.052 | -0.052
kidslt6 | -0.868 | - | - | - | -0.874 | -0.892 | -0.902 | -0.864 | -0.855 | -0.844
kidsge6 | 0.036 | - | - | - | 0.037 | 0.038 | 0.041 | 0.036 | 0.035 | 0.032
const. | 0.270 | - | - | - | 0.271 | 0.216 | 0.234 | 0.248 | 0.256 | 0.249
| Standard Errors
| ase | boot | dmk | ks | rnrn | rnr200 | rnr100 | rqnn | rqn200 | rqn100
nwifeinc | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.006 | 0.005 | 0.005 | 0.005 | 0.005
educ | 0.025 | 0.026 | 0.026 | 0.025 | 0.025 | 0.027 | 0.028 | 0.027 | 0.025 | 0.025
exper | 0.019 | 0.020 | 0.019 | 0.019 | 0.019 | 0.020 | 0.021 | 0.019 | 0.018 | 0.017
exper2 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001
age | 0.008 | 0.009 | 0.008 | 0.008 | 0.009 | 0.008 | 0.009 | 0.009 | 0.008 | 0.008
kidslt6 | 0.119 | 0.120 | 0.118 | 0.118 | 0.120 | 0.119 | 0.129 | 0.117 | 0.113 | 0.117
kidsge6 | 0.043 | 0.046 | 0.045 | 0.045 | 0.045 | 0.048 | 0.047 | 0.044 | 0.042 | 0.045
const. | 0.509 | 0.512 | 0.507 | 0.505 | 0.494 | 0.535 | 0.544 | 0.544 | 0.494 | 0.506

#### Sample R code to implement the Mroz Example.

⬇

set.seed(123) # set the seed

library(numDeriv) # compute numerical derivaties using finite differences,
alternative: library(pracma) is usually faster

library(foreign) # to load the data set in Stata dta format

data = read.dta(’mroz.dta’) # read the mroz data

y = data$inlf # outcome variable

X = cbind(data$nwifeinc,data$educ,data$exper, # regressors

data$exper^2,data$age,data$kidslt6,data$kidsge6,1)

colnames(X) = c(’nwifeinc’,’educ’,’exper’,’exper2’, # labels

’age’,’kidslt6’,’kidsge6’,’constant’)

n = 753 # sample size

index0 = 1:n # indices for the sample data

loglik <- function(coef,index=index0) {

# compute the log-likelihood for the Probit model on the observations indexed
by index (default 1:n, the original sample) at theta = coef

score = X[index,]%*%coef # compute the z-scores

ll = y[index]*log( pnorm(score) ) +

(1-y[index])*log( 1-pnorm(score) )

return( sum( ll ) )

}

d_loglik <- function(coef,index=index0) {

# compute the gradient of the log-likelihood for the Probit model on the
observations indexed by index (default 1:n, the original sample) at theta =
coef

# In this example, the gradient is analytically tractable, it could be
evaluated by finite differences by using the following:

# d_loglik <- function(coef,index=index0) {
return(jacobian(loglik,coef,index=index)) }

yy = y[index] # keep observations indexed by index

XX = X[index,] # keep observations indexed by index

score = XX%*%coef # compute the z-score

dll = 0 # initialize the gradient

for (i in 1:length(index)) {

dll = dll +

(yy[i]*XX[i,]*dnorm(score[i])/pnorm(score[i]) -

(1-yy[i])*XX[i,]*dnorm(score[i])/(1-pnorm(score[i])))

}

return(dll)

}

rNR <- function(coef0, learn = 0.1, iter = 500, m = n) {

# generate ’B = iter’ rNR draws with learning rate ’gamma = learn’ with m out
of n resampling

coefs = matrix(NA,iter,length(coef0)) # matrix where draws will be stored

coefs[1,] = coef0 # initialize the first-draw

for (i in 2:iter) {

index = sample(1:n,m,replace=TRUE) # sample m out of n observations with
replacement

G = d_loglik(coefs[i-1,],index=index) # compute the resampled gradient G using
analytical derivatives. Alternative using finite differences:

# G = jacobian(loglik,coefs[i-1,],index=index)

H = hessian(loglik,coefs[i-1,],index=index) # compute the resampled hessian H
using finite differences; we could also compute the jacobian of the gradient
d_loglik

coefs[i,] = coefs[i-1,] - learn*solve(H,G) # update

}

colnames(coefs) = colnames(X) # label the coefficients

return( list(coefs = coefs) ) # return draws

}

# estimates and standard errors (source: Introductory Econometrics, A Modern
Approach 2nd Edition, Wooldridge)

coef = c(-0.012,0.131,0.123,-0.0019,-0.053,-0.868,0.036,0.270)

ses = c( 0.005,0.025,0.019, 0.0006, 0.008, 0.119,0.043,0.509)

iter_rNR = 2e3 # number of rNR draws

learn = 0.3 # learning rate

coef0 = coef*3.25 # starting value

m1 = 753 # m = n

m2 = 200 # m = 200

m3 = 100 # m = 100

# adjustments to get valid standard errors

adj_rnr1 = sqrt(m1/n)*sqrt( (1-(1-learn)^2)/learn^2 )

adj_rnr2 = sqrt(m2/n)*sqrt( (1-(1-learn)^2)/learn^2 )

adj_rnr3 = sqrt(m3/n)*sqrt( (1-(1-learn)^2)/learn^2 )

b1 = 1 + round(log(0.01)/log(1-learn)) # burn-in sample size

# generate rNR draws

out_rNR1 = rNR(coef0,learn,b1 + iter_rNR, m1)

out_rNR2 = rNR(coef0,learn,b1 + iter_rNR, m2)

out_rNR3 = rNR(coef0,learn,b1 + iter_rNR, m3)

# format output

estimates =

rbind( coef,

apply(out_rNR1$coef[b1:(iter_rNR+b1),],2,mean),

apply(out_rNR2$coef[b1:(iter_rNR+b1),],2,mean),

apply(out_rNR3$coef[b1:(iter_rNR+b1),],2,mean))

std_errs =

rbind( ses,

apply(out_rNR1$coef[b1:(iter_rNR+b1),],2,sd)*adj_rnr1,

apply(out_rNR2$coef[b1:(iter_rNR+b1),],2,sd)*adj_rnr2,

apply(out_rNR3$coef[b1:(iter_rNR+b1),],2,sd)*adj_rnr3)

estimates = as.data.frame(estimates)

colnames(estimates) = colnames(X)

rownames(estimates) = c(’MLE’,’rNRn’,’rNR200’,’rNR100’)

std_errs = as.data.frame(std_errs)

colnames(std_errs) = colnames(X)

rownames(std_errs) = c(’ase’,’rNRn’,’rNR200’,’rNR100’)

# print results

print(round( cbind( t(estimates), t(std_errs) ), digits = 3 ))

# output printed below:

# MLE rNRn rNR200 rNR100 ase rNRn rNR200 rNR100

#nwifeinc -0.012 -0.012 -0.013 -0.014 0.005 0.005 0.005 0.005

#educ 0.131 0.132 0.136 0.140 0.025 0.026 0.026 0.028

#exper 0.123 0.123 0.123 0.125 0.019 0.019 0.020 0.021

#exper2 -0.002 -0.002 -0.002 -0.002 0.001 0.001 0.001 0.001

#age -0.053 -0.053 -0.054 -0.055 0.008 0.008 0.009 0.009

#kidslt6 -0.868 -0.872 -0.895 -0.917 0.119 0.121 0.121 0.126

#kidsge6 0.036 0.038 0.040 0.038 0.043 0.045 0.047 0.049

#constant 0.270 0.272 0.282 0.276 0.509 0.506 0.505 0.535

##  Appendix F Primitive Conditions for Assumption 2

The following provides primitive conditions for Assumption 2 i-ii. in the main
text. Although the dimension of $\theta$ is fixed in the main results, the
derivation make explicit the dependence of $C_{2},C_{3}$ on the dimension
$d_{\theta}$.

#### Setup and Notation.

Let $z_{i}=(y_{i},x_{i})$ be iid, the parameter space $\Theta$ is a compact,
convex subset of $\mathbb{R}^{d_{\theta}}$, the sample objective, gradient and
hessian are given by:

$Q_{n}(\theta)=\frac{1}{n}\sum_{i=1}^{n}q(z_{i},\theta),\quad
G_{n}(\theta)=\frac{1}{n}\sum_{i=1}^{n}\nabla q(z_{i},\theta),\quad
H_{n}(\theta)=\frac{1}{n}\sum_{i=1}^{n}\partial^{2}_{\theta,\theta^{\prime}}q(z_{i},\theta).$

The re-sampled objective relies on $m$ draws $z_{i}^{(b)}$ taken with equal
probability and with replacement from $(z_{1},\dots,z_{n})$, the re-sampled
objective, gradient and hessian are given by:

$Q_{m}^{(b)}(\theta)=\frac{1}{m}\sum_{i=1}^{m}q(z_{i}^{(b)},\theta),\quad
G_{m}^{(b)}(\theta)=\frac{1}{m}\sum_{i=1}^{m}\nabla
q(z_{i}^{(b)},\theta),\quad
H_{m}^{(b)}(\theta)=\frac{1}{m}\sum_{i=1}^{m}\nabla^{2}q(z_{i}^{(b)},\theta).$

Alternatively, the re-weighted objective relies on $m=n$ iid random weights
$w_{i}^{(b)}$ with mean and variance equal to one, the re-sampled objective,
gradient and hessian are given by:

$Q_{m}^{(b)}(\theta)=\frac{1}{m}\sum_{i=1}^{m}w_{i}^{(b)}q(z_{i},\theta),\quad
G_{m}^{(b)}(\theta)=\frac{1}{m}\sum_{i=1}^{m}w_{i}^{(b)}\nabla
q(z_{i},\theta),\quad
H_{m}^{(b)}(\theta)=\frac{1}{m}\sum_{i=1}^{m}w_{i}^{(b)}\nabla^{2}q(z_{i},\theta),$

for Gaussian multiplier weights we have $w_{i}^{(b)}\sim\mathcal{N}(1,1)$.

######  Assumption F6 (Lipschitz Derivatives).

There exists a measurable function $C_{q}$ such that for any two
$\theta_{1},\theta_{2}\in\Theta$ and any $z_{i}$:

$\displaystyle\|\nabla q(z_{i},\theta_{1})-\nabla q(z_{i},\theta_{2})\|$
$\displaystyle\leq C_{q}(z_{i})\|\theta_{1}-\theta_{2}\|,$ (F.5)
$\displaystyle\|\nabla^{2}q(z_{i},\theta_{1})-\nabla^{2}q(z_{i},\theta_{2})\|$
$\displaystyle\leq C_{q}(z_{i})\|\theta_{1}-\theta_{2}\|,$ (F.6)

Let $p\geq 2$, $C_{q}$ and the derivatives have finite $p$-th moment:

$\displaystyle\mathbb{E}(|C_{q}(z_{i})|^{p})<+\infty,\quad\mathbb{E}(\|\nabla
q(z_{i},\theta^{\dagger})\|^{p})<+\infty,\quad\mathbb{E}(\|\nabla^{2}q(z_{i},\theta^{\dagger})\|^{p})<+\infty$
(F.7)

The following subsections will prove the Lemma below under re-sampling and
then under re-weighting with Gaussian multiplier weights.

######  Lemma F3.

Suppose Assumption F6 holds, then we have:

  * i.

$\mathbb{E}^{\star}[G_{m}^{(b)}(\theta)]=G_{n}(\theta)$,

  * ii.

$(\mathbb{E}^{\star}[\sup_{\theta\in\Theta}\|G_{m}^{(b)}(\theta)-G_{n}(\theta)\|_{2}^{p}])^{1/p}\leq
C_{2n}m^{-1/2}$, where $C_{2n}\overset{p}{\to}C_{2}$ finite,

  * ii.

$(\mathbb{E}^{\star}[\sup_{\theta\in\Theta}\|H_{m}^{(b)}(\theta)-H_{n}(\theta)\|_{2}^{p}])^{1/p}\leq
C_{3n}m^{-1/2}$, where $C_{3n}\overset{p}{\to}C_{3}$ finite.

###  F.1 Re-sampled objective

###### Proof of Lemma F3.

Condition i. is immediate from iid resampling and the additivity over $i$ of
the sample and re-sampled gradients. Lemma F3 ii. and iii. are proved the same
way so the below will only focus on ii. for brievety.

Let
$F_{q}(z_{i},\theta)=C_{q}(z_{i})\|\theta-\theta^{\dagger}\|+\|\partial_{\theta}q(z_{i},\theta^{\dagger})\|$
be the enveloppe function which satisfies:
$\|\partial_{\theta}q(z_{i},\theta)\|\leq F_{q}(z_{i},\theta)$ for all
$\theta\in\Theta$. We have $F_{q}(z_{i},\theta)\leq
F_{q}(z_{i})=C_{q}(z_{i})\text{diam}(\Theta)+\|\partial_{\theta}q(z_{i},\theta^{\dagger})\|$,
which does not depend on $\theta$. Here
$\text{diam}(\Theta)=\sup_{\theta_{1},\theta_{2}\in\Theta}\|\theta_{1}-\theta_{2}\|$,
by compactness it is finite. The main idea is to apply Theorem 2.14.5 in van
der Vaart and Wellner (1996), under the distribution $\mathbb{P}^{\star}$, for
a $p\geq 2$ such that Assumption F6 holds:

$\displaystyle\left[\mathbb{E}^{\star}\left(\sup_{\theta\in\Theta}\sqrt{m}\|G_{m}^{(b)}(\theta)-G_{n}(\theta)\|^{p}\right)\right]^{1/p}$
$\displaystyle\lesssim\mathbb{E}^{\star}\left(\sup_{\theta\in\Theta}\sqrt{m}\|G_{m}^{(b)}(\theta)-G_{n}(\theta)\|\right)+m^{1/p-1/2}\left[\frac{1}{n}\sum_{i=1}^{n}|F_{q}(z_{i})|^{p}\right]^{1/p},$

where $\lesssim$ stands for less or equal than, up to universal constants. By
the strong law of large numbers and the finite moment Assumption F6 (F.7), we
have, as $n\to\infty$.
$\left[\frac{1}{n}\sum_{i=1}^{n}|F_{q}(z_{i})|^{p}\right]^{1/p}\overset{p}{\to}\left[\mathbb{E}(|F_{q}(z_{i})|^{p})\right]^{1/p}$.
Also for $p\geq 2$ and $m\geq 1$, $0\leq m^{1/p-1/2}\leq 1$. Now using Theorem
2.14.2 in van der Vaart and Wellner (1996), again under the distribution
$\mathbb{P}^{\star}$, we have:

$\displaystyle\mathbb{E}^{\star}\left(\sup_{\theta\in\Theta}\sqrt{m}\|G_{m}^{(b)}(\theta)-G_{n}(\theta)\|\right)\lesssim
J_{[\,]}(1,\mathcal{F},L_{2}(\mathbb{P}_{n}))\left[\frac{1}{n}\sum_{i=1}^{n}|F_{q}(z_{i})|^{2}\right]^{1/2},$

where $J_{[\,]}(1,\mathcal{F},L_{2}(\mathbb{P}_{n}))=\int_{0}^{1}\sqrt{1+\log
N_{[\,]}(\varepsilon\|F_{q}\|,\mathcal{F},L_{2}(\mathbb{P}_{n}))}d\varepsilon$
is the bracketing integral of the functions class
$\mathcal{F}=\\{\theta\to\partial_{\theta}q(z_{i},\theta)\\}$ and
$N_{[\,]}(\varepsilon\|F_{q}\|,\mathcal{F},L_{2}(\mathbb{P}_{n}))$ its
bracketing number. By Theorem 9.23 in Kosorok (2007), we have
$N_{[\,]}(\varepsilon\|F_{q}\|,\mathcal{F},L_{2}(\mathbb{P}_{n}))\leq
N(\varepsilon/2,\Theta,\|\cdot\|)\leq(6/\varepsilon)^{d_{\theta}}\text{vol}(\Theta)/\text{vol}(B)$
where $\text{vol}(\Theta),\text{vol}(B)$ are the volumes of $\Theta$ and the
unit sphere $B$ in $\mathbb{R}^{d_{\theta}}$, respectively. As a result, we
have:

$\displaystyle\mathbb{E}^{\star}\left(\sup_{\theta\in\Theta}\sqrt{m}\|G_{m}^{(b)}(\theta)-G_{n}(\theta)\|\right)\lesssim\sqrt{d_{\theta}}\left[\frac{1}{n}\sum_{i=1}^{n}|F_{q}(z_{i})|^{2}\right]^{1/2},$

and
$\left[\frac{1}{n}\sum_{i=1}^{n}|F_{q}(z_{i})|^{2}\right]^{1/2}\overset{p}{\to}\left[\mathbb{E}(|F_{q}(z_{i})|^{2})\right]^{1/2}$
finite. Using the inequality:
$\left[\frac{1}{n}\sum_{i=1}^{n}|F_{q}(z_{i})|^{2}\right]^{1/2}\leq\left[\frac{1}{n}\sum_{i=1}^{n}|F_{q}(z_{i})|^{p}\right]^{1/p}$,
we can conclude that:

$\left[\mathbb{E}^{\star}\left(\sup_{\theta\in\Theta}\sqrt{m}\|G_{m}^{(b)}(\theta)-G_{n}(\theta)\|^{p}\right)\right]^{1/p}\lesssim(1+\sqrt{d_{\theta}})\left[\frac{1}{n}\sum_{i=1}^{n}|F_{q}(z_{i})|^{p}\right]^{1/p},$

where the last term converges in probability to
$[\mathbb{E}(|F_{q}(z_{i})|^{p})]^{1/p}$, finite. ∎

###  F.2 Re-weighted objective with Gaussian multiplier weights.

The following specializes to the case where $w_{i}^{(b)}\sim\mathcal{N}(1,1)$,
iid. Notice that:

$\displaystyle G_{m}^{(b)}(\theta)$
$\displaystyle=G_{n}(\theta)+\frac{1}{m}\sum_{i=1}^{m}(w_{i}^{(b)}-1)\nabla
q(z_{i},\theta)\Rightarrow
G_{m}^{(b)}(\theta)\sim\mathcal{N}\left(G_{n}(\theta),\frac{1}{m}V_{n}(\theta)\right),$

where $V_{n}(\theta)=1/n\sum_{i=1}^{n}\nabla q(z_{i},\theta)\nabla
q(z_{i},\theta)^{\prime}$ and $w_{i}^{(b)}-1\sim\mathcal{N}(0,1)$. Note that
the proof below only requires Assumption F6 to hold with $p=2$, even if when
the desired result is stated for $p>2$.

###### Proof of Lemma F3.

Condition i. is immediate from the above. For any $\theta\in\Theta$,

$\sqrt{m}\left(G_{m}^{(b)}(\theta)-G_{n}(\theta)\right)\sim\mathcal{N}(0,V_{n}(\theta)),$

is a vector-valued Gaussian process, conditional on the sample of data
$z_{1},\dots,z_{n}$. For each $i=1,\dots,n$ we have:

$\|(w_{i}^{(b)}-1)\nabla
q(z_{i},\theta)\|\leq|w_{i}^{(b)}-1|[C_{q}(z_{i})\text{diam}(\Theta)+\|\nabla
q(z_{i},\theta^{\dagger})\|]=|w_{i}^{(b)}-1|F_{q}(z_{i}),$

which defines the enveloppe function for the re-weighted objective. For any
$a\in S^{d_{\theta}}=\\{a\in\mathbb{R}^{d_{\theta}},\|a\|=1\\}$, the surface
of the unit sphere, define
$\sqrt{m}a^{\prime}(G_{m}^{(b)}(\theta)-G_{n}(\theta))$, a scalar-value
Gaussian process defined on $S^{d_{\theta}}\times\Theta$ a compact subset of
$\mathbb{R}^{d_{\theta}}\times\mathbb{R}^{d_{\theta}}$. Let
$\sigma_{n}^{2}=\sup_{(a,\theta)\in
S^{d_{\theta}}\times\Theta}\mathbb{E}^{\star}(m|a^{\prime}(G_{m}^{(b)}(\theta)-G_{n}(\theta)|^{2})=\sup_{(a,\theta)\in
S^{d_{\theta}}\times\Theta}a^{\prime}V_{n}(\theta)a$ which is finite and
converges in probability to
$\sup_{\theta\in\Theta}[\lambda_{\text{max}}(\mathbb{E}[\nabla
q(z_{i},\theta)\nabla q(z_{i},\theta)^{\prime}])]$, also finite. We can now
apply results for scalar valued Gaussian processes. In particular, using
Theorem 5.8 in Boucheron et al. (2013), see also Proposition 3.19 in Massart
(2007), we have for any $u\geq 0$:

$\displaystyle\mathbb{P}^{\star}\left(\sup_{(a,\theta)\in
S^{d_{\theta}}\times\Theta}\sqrt{m}|a^{\prime}(G_{m}^{(b)}(\theta)-G_{n}(\theta)|\geq
M_{n}+u\right)\leq\exp\left(-\frac{u^{2}}{2\sigma_{n}^{2}}\right),$

where $M_{n}=\mathbb{E}^{\star}\left[\sup_{(a,\theta)\in
S^{d_{\theta}}\times\Theta}\sqrt{m}|a^{\prime}(G_{m}^{(b)}(\theta)-G_{n}(\theta)|\right]$.

Let $Z_{m}^{(b)}=\sup_{(a,\theta)\in
S^{d_{\theta}}\times\Theta}\sqrt{m}|a^{\prime}(G_{m}^{(b)}(\theta)-G_{n}(\theta)|$.
The main idea is to notice that $\sup_{(a,\theta)\in
S^{d_{\theta}}}\sqrt{m}|a^{\prime}(G_{m}^{(b)}(\theta)-G_{n}(\theta)|$ is
atained at
$a=[G_{m}^{(b)}(\theta)-G_{n}(\theta)]/\|G_{m}^{(b)}(\theta)-G_{n}(\theta)\|$
so that
$|Z_{m}^{(b)}|^{p}=\sup_{\theta\in\Theta}\|G_{m}^{(b)}(\theta)-G_{n}(\theta)\|^{p}$
which is the quantity we want to bound in expectations. Given that
$Z_{m}^{(b)}\geq 0$ we can write its $p$-th moment as:

$\displaystyle\mathbb{E}^{\star}(|Z_{m}^{(b)}|^{p})$
$\displaystyle=p\int_{0}^{\infty}z^{p-1}\mathbb{P}^{\star}(Z_{m}^{(b)}>z)dz$
$\displaystyle=p\int_{0}^{M_{n}}z^{p-1}\mathbb{P}^{\star}(Z_{m}^{(b)}>z)dz+p\int_{0}^{\infty}(M_{n}+z)^{p-1}\mathbb{P}^{\star}(Z_{m}^{(b)}>M_{n}+z)dz$
$\displaystyle\leq
M_{n}^{p}+p\int_{0}^{\infty}(M_{n}+z)^{p-1}\exp\left(-\frac{z^{2}}{2\sigma_{n}^{2}}\right)dz$
$\displaystyle=M_{n}^{p}+p\sigma_{n}\int_{0}^{\infty}(M_{n}+\sigma_{n}z)^{p-1}\exp\left(-\frac{z^{2}}{2}\right)dz.$

We already known that $\sigma_{n}$ is bounded and converges in probability to
a finite limit. It remains to bound the moment $M_{n}$ before we can conclude.

Let $X_{m}^{(b)}(\theta)=\sqrt{m}(G_{m}^{(b)}(\theta)-G_{n}(\theta)$, by
Gaussian re-weighting it is a Gaussian process. We will use an inequality for
separable sub-Gaussian processes, specifically Corollary 2.2.8 in van der
Vaart and Wellner (1996), to bound $M_{n}$. The following verifies the
assumptions required to apply the result. For any
$\theta_{1},\theta_{2}\in\Theta$, we have:

$\displaystyle
X_{m}^{(b)}(\theta_{1})-X_{m}^{(b)}(\theta_{2})=\frac{1}{\sqrt{m}}\sum_{i=1}^{m}(w_{i}^{b}-1)[\partial_{\theta}q(z_{i},\theta_{1})-\partial_{\theta}q(z_{i},\theta_{2})]\sim\mathcal{N}\left(0,V_{m}(\theta_{1},\theta_{2})\right),$

where
$V_{m}(\theta_{1},\theta_{2})=\frac{1}{m}\sum_{i=1}^{m}[\partial_{\theta}q(z_{i},\theta_{1})-\partial_{\theta}q(z_{i},\theta_{2})][\partial_{\theta}q(z_{i},\theta_{1})-\partial_{\theta}q(z_{i},\theta_{2})]^{\prime}$.
For any $\theta$, let $X_{m,j}^{(b)}(\theta)$ be the $j$-th row of
$X_{m}^{(b)}(\theta)$. For any $j\in\\{1,\dots,d_{\theta}\\}$, we have
$X_{m,j}^{(b)}(\theta_{1})-X_{m,j}^{(b)}(\theta_{2})\sim\mathcal{N}(0,V_{m,jj}(\theta_{1},\theta_{2}))$
where $V_{m,jj}(\theta_{1},\theta_{2})$ is the $j$-th diagonal element of
$V_{m}(\theta_{1},\theta_{2})$ which satisfies
$V_{m,jj}(\theta_{1},\theta_{2})\leq\text{trace}[V_{m}(\theta_{1},\theta_{2})]=\frac{1}{m}\sum_{i=1}^{m}\|\partial_{\theta}q(z_{i},\theta_{1})-\partial_{\theta}q(z_{i},\theta_{2})\|^{2}\leq\frac{1}{m}\sum_{i=1}^{m}F_{q}(z_{i})^{2}\|\theta_{1}-\theta_{2}\|^{2},$
by Lipschitz continuity. This implies the tail inequality:

$\mathbb{P}^{\star}\left(|X_{m,j}^{(b)}(\theta_{1})-X_{m,j}^{(b)}(\theta_{2})|>x\right)\leq
2\exp\left(\frac{-x^{2}}{2\frac{1}{m}\sum_{i=1}^{m}F_{q}(z_{i})^{2}\|\theta_{1}-\theta_{2}\|^{2}}\right),$

which implies that each $X_{m,j}^{(b)}$ is a separable sub-Gaussian process
under the semi-metric
$d(\theta_{1},\theta_{2})=[\frac{1}{m}\sum_{i=1}^{m}F_{q}(z_{i})^{2}]^{1/2}\|\theta_{1}-\theta_{2}\|$.
Using Corollary 2.2.8 in van der Vaart and Wellner (1996), we have for some
universal constant $K$:

$\mathbb{E}^{\star}\left(\sup_{\theta\in\Theta}|X_{m,j}^{(b)}(\theta)|\right)\leq\mathbb{E}^{\star}\left(|X_{m,j}^{(b)}(\theta^{\dagger})|\right)+K\int_{0}^{\infty}\sqrt{\log
D(\varepsilon,d)}d\varepsilon,$

where by Gaussianity of $X_{m,j}^{(b)}$, we have
$\mathbb{E}^{\star}\left(|X_{m,j}^{(b)}(\theta^{\dagger})|\right)=2/\pi
V_{n,jj}^{1/2}$, $V_{n,jj}$ is the $j$-th diagonal element of $V_{n}$ defined
earlier and $D(\varepsilon,d)\leq
3[\frac{1}{n}\sum_{i=1}^{n}F_{q}(z_{i})^{2}]^{1/2}\text{diam}(\Theta)/\varepsilon$.
Now notice that:

$\displaystyle
M_{n}=\mathbb{E}^{\star}\left(\sup_{\theta\in\Theta}\|X_{m}^{(b)}(\theta)\|\right)$
$\displaystyle\leq\sum_{j=1}^{d_{\theta}}\mathbb{E}^{\star}\left(\sup_{\theta\in\Theta}|X_{m,j}^{(b)}(\theta)|\right)$
$\displaystyle\leq
d_{\theta}\left[2/\pi(\text{max}_{j}V_{n,jj}^{1/2})+K\int_{0}^{\infty}\sqrt{\log
D(\varepsilon,d)}d\varepsilon\right].$

Without loss of generality, assume that $M_{n}\geq 1$ then we have:

$\displaystyle\mathbb{E}^{\star}(|Z_{m}^{(b)}|^{p})$ $\displaystyle\leq
M_{n}^{p}\left[1+p\sigma_{n}\int_{0}^{\infty}(1+\sigma_{n}z)^{p-1}\exp\left(-\frac{z^{2}}{2}\right)dz\right].$

Putting everything together, we have the following inequality for any $p\geq
2$:

$\displaystyle\left[\mathbb{E}^{\star}\left(\sup_{\theta}\|\sqrt{m}(G_{m}^{(b)}(\theta)-G_{n}(\theta))\|^{p}\right)\right]^{1/p}$
$\displaystyle\leq
d_{\theta}\left[2/\pi(\text{max}_{j}V_{n,jj}^{1/2})+K\int_{0}^{\infty}\sqrt{\log
D(\varepsilon,d)}d\varepsilon\right]\left[1+p\sigma_{n}\int_{0}^{\infty}(1+\sigma_{n}z)^{p-1}\exp\left(-\frac{z^{2}}{2}\right)dz\right]^{1/p}.$

where $\int\sqrt{\log D(\varepsilon,d)}d\varepsilon$ and $\sigma_{n}$ converge
in probability to a finite limit if Assumption F6 holds for $p=2$ since the
quantities involved only depend on second moments. ∎

##  Appendix G Primitive Conditions for Assumption 4

The derivations below use the same setup and notation as Appendix F above.

#### Re-sampled objective.

Suppose the gradient has finite third moment, i.e. $\mathbb{E}(\|\nabla
q(z_{i},\theta)\|^{3})<\infty$. The calculations below are based on the
derivations in Lahiri (2013, Ch6.2, pp147-148). Note that
$V_{n}^{-1/2}[H_{n}(\hat{\theta}_{n})]^{-1}G_{m}^{(b)}(\hat{\theta}_{n})\overset{d^{\star}}{=}\Sigma_{n}^{-1/2}G_{m}^{(b)}(\hat{\theta}_{n})$
so that we can focus on the latter. Also, note that
$\mathbb{E}^{\star}[\Sigma_{n}^{-1/2}G_{m}^{(b)}(\hat{\theta}_{n})]=0$,
$\text{var}^{\star}[\Sigma_{n}^{-1/2}G_{m}^{(b)}(\hat{\theta}_{n})]=I_{d}$.
Now compute the characteristic function under $\mathbb{E}^{\star}$:

$\varphi^{\star}_{m}(\tau)=\mathbb{E}^{\star}\left(\exp[\mathbf{i}\tau^{\prime}V_{n}^{-1/2}\sqrt{m}G_{m}^{(b)}(\hat{\theta}_{n})]\right)=\left[\mathbb{E}^{\star}\left(\exp\left[\mathbf{i}\tau^{\prime}V_{n}^{-1/2}\frac{\nabla
q(z_{i}^{(b)},\hat{\theta}_{n})}{\sqrt{m}}\right]\right)\right]^{m}.$

Using a Taylor expansion with remainder, we have for each $\tau$:

$\mathbb{E}^{\star}\left(\exp\left[\mathbf{i}\tau^{\prime}V_{n}^{-1/2}\frac{\nabla
q(z_{i}^{(b)},\hat{\theta}_{n})}{\sqrt{m}}\right]\right)=1-\frac{\|\tau\|_{2}^{2}}{2m}+R_{m}(\tau),$

where $|R_{m}(\tau)|\leq\mathbb{E}^{\star}(\|\nabla
q(z_{i}^{(b)},\hat{\theta}_{n})\|_{2}^{3})\|\tau\|_{2}^{3}m^{-3/2}$ because
all derivatives of the exponential term have modulus less than one. Apply a
Taylor expansion with integral remainder to the logarithm:

$\log[\varphi_{m}^{\star}(\tau)]=-\frac{\|\tau\|_{2}^{2}}{2}+mR_{m}(\tau)-m\int_{0}^{-\frac{\|\tau\|_{2}^{2}}{2m}+R_{m}(\tau)}\frac{(-\|\tau\|_{2}^{2}/[2m]+R_{m}(\tau)-t)}{(1+t)^{2}}dt,$

where $mR_{m}(\tau)\leq\mathbb{E}^{\star}(\|\nabla
q(z_{i}^{(b)},\hat{\theta}_{n})\|_{2}^{3})\|\tau\|_{2}^{3}m^{-1/2}$;
$\mathbb{E}^{\star}(\|\nabla
q(z_{i}^{(b)},\hat{\theta}_{n})\|_{2}^{3})\overset{p}{\to}\mathbb{E}(\|\nabla
q(z_{i},\theta^{\dagger})\|_{2}^{3})<\infty$ using a uniform law of large
numbers. For $m$ sufficiently large,
$|-\frac{\|\tau\|_{2}^{2}}{2m}+R_{m}(\tau)|<1/2$ with probability approaching
$1$, so we can use the following bound:

$m\Big{|}\int_{0}^{-\frac{\|\tau\|_{2}^{2}}{2m}+R_{m}(\tau)}\frac{(-\|\tau\|_{2}^{2}/[2m]+R_{m}(\tau)-t)}{(1+t)^{2}}dt\Big{|}\leq
6m\left[\frac{\|\tau\|_{2}^{4}}{4m^{2}}+|R_{m}(\tau)|^{2}\right],$

which satisfies the conditions required for Assumption 4.

#### Re-weighted objective with Gaussian multiplier weights.

Assumption 4 automatically holds because as noted in Appendix F, conditionally
on the sample: $G_{m}^{(b)}(\hat{\theta}_{n})\sim\mathcal{N}(0,\Sigma_{n})$
which implies
$V_{n}^{-1/2}[H_{n}(\hat{\theta}_{n})]^{-1}G_{m}^{(b)}(\hat{\theta}_{n})\sim\mathcal{N}(0,I_{d})$
for which the characteristic function is the same as in Assumption 4 with
$r_{m}(\tau)=0$.

##  Appendix H Additional Results for Section 5

###  H.1 Example 1: Dynamic Discrete Choice with Unobserved Heterogeneity

#### Baseline results with homogeneity.

Table H5 provides baseline results for homogeneous dynamic discrete choice
model where estimation is much faster so that a comparion with rnr, and the
standard bootstrap is feasible. The specification is taken from Abbring and
Klein (2020). The first three columns correspond to MLE estimates and
rejection rates using the standard $m$ out of $n$ bootstrap. The other 6
columns correspond to rnr and rqn.

Table H5: Baseline estimation and inference results with homogeneous agents

|  mle/bootstrap | rnr | rqn
---|---|---|---
m | $\beta_{0}$ | $\beta_{1}$ | $\delta_{1}$ | $\beta_{0}$ | $\beta_{1}$ | $\delta_{1}$ | $\beta_{0}$ | $\beta_{1}$ | $\delta_{1}$
| Average Estimate
1000 | -0.500 | 0.200 | 1.000 | -0.500 | 0.200 | 1.000 | -0.500 | 0.200 | 1.000
500 | - | - | - | -0.500 | 0.200 | 1.000 | -0.500 | 0.200 | 1.000
100 | - | - | - | -0.500 | 0.200 | 1.000 | -0.500 | 0.200 | 1.000
| Standard Deviation
1000 | 0.014 | 0.004 | 0.013 | 0.014 | 0.004 | 0.013 | 0.014 | 0.004 | 0.013
500 | - | - | - | 0.014 | 0.004 | 0.014 | 0.014 | 0.004 | 0.014
100 | - | - | - | 0.014 | 0.004 | 0.013 | 0.014 | 0.004 | 0.013
| Rejection Rates
1000 | 0.053 | 0.051 | 0.055 | 0.054 | 0.050 | 0.051 | 0.051 | 0.049 | 0.052
500 | 0.052 | 0.053 | 0.062 | 0.053 | 0.052 | 0.066 | 0.051 | 0.051 | 0.061
100 | 0.064 | 0.060 | 0.061 | 0.057 | 0.052 | 0.059 | 0.061 | 0.058 | 0.055

Legend: $n=1000$, $T=100$, $\gamma=0.1$, nominal size $=5\%$.

#### Results with $T=50$, comparison with BFGS

Table H6: Dynamic Discrete Choice Model with Heterogeneity ($T=50$)

| $\mu_{0}^{1}$ | $\mu_{1}^{1}$ | $\mu_{0}^{2}$ | $\mu_{1}^{2}$ | $100\sigma_{1}$ | $100\sigma_{2}$ | $\omega$ | $\delta_{1}$
---|---|---|---|---|---|---|---|---
$\theta^{\dagger}$ | -2.000 | 0.300 | -1.000 | 0.900 | 0.010 | 0.010 | 0.300 | 1.000
| Average Estimates
rqn | -2.032 | 0.302 | -0.983 | 0.906 | 1.132 | 0.907 | 0.294 | 0.954
rqn-bc | -2.014 | 0.303 | -1.000 | 0.894 | 1.140 | 0.793 | 0.294 | 1.000
bfgs | -1.823 | 0.222 | -0.950 | 0.945 | 4.029 | 3.620 | 0.329 | 0.860
bfgs-bc | -1.996 | 0.269 | -1.014 | 0.923 | 1.167 | 1.938 | 0.320 | 0.955
| Standard Deviation
rqn | 0.044 | 0.012 | 0.025 | 0.015 | 0.165 | 0.365 | 0.013 | 0.023
rqn-bc | 0.047 | 0.013 | 0.025 | 0.017 | 0.242 | 0.429 | 0.013 | 0.026
bfgs | 0.049 | 0.018 | 0.031 | 0.015 | 0.306 | 0.180 | 0.011 | 0.028
bfgs-bc | 0.056 | 0.021 | 0.032 | 0.016 | 0.368 | 0.258 | 0.013 | 0.031
| Rejection Rate
rqn | 0.035 | 0.035 | 0.068 | 0.035 | 0.033 | 0.393 | 0.018 | 0.255
rqn-bc | 0.022 | 0.033 | 0.013 | 0.015 | 0.030 | 0.370 | 0.018 | 0.043
rqn-bcse | 0.022 | 0.030 | 0.015 | 0.015 | 0.033 | 0.367 | 0.018 | 0.045

Legend: rqn-bc = split panel bias-corrected rqn; rqn, rqn-bc = quantile-based
CIs, rqn-bcse = std-error based CIs.
$\gamma=0.1,n=1000,T=50,m=n/2,B=2000,\textsc{burn}=250,$ nominal level $=5\%.$
bfgs: R’s optim with bfgs optimizer and tight convergence criteria, bfgs-bc:
split panel biased corrected.

#### Implementation details.

A grid of $60$ sobol points is used to integrate out the random coefficients
and solve the fixed point problem. During the burn-in of 250 draws, a smaller
grid of $30$ points is used. The split panel
$\theta_{b,nT/2}^{1},\theta_{b,nT/2}^{2}$ are only computed after
$\frac{2}{3}\textsc{burn}$ draws and initialized at $\theta_{b,nT}$ to save
time. The quasi-Newton matrix $P_{b}$ is pooled between split panel draws but
computed separately for the full panel ones. The regularization schedule is
$\lambda=20$ for the first $\textsc{burn}/2$ draws and then an exponentially
decaying schedule is used: $\lambda_{b+1}=0.9\lambda_{b}$.

#### Derivations for the second-order bias.

The following derives the second-order asymptotic bias of the random
coefficient dynamic discrete choice model in Section 5.1 under large-$T$
asymptotics and simplifying assumptions. In particular, the transition
probability $\Pi$ is assumed to be known and a quadratic assumption is used to
shorten integral derivations. Let $Q_{nT}(\theta)$ be the negative sample log-
likelihood:

$Q_{nT}(\theta,\Pi)=-\frac{1}{nT}\sum_{i=1}^{n}\log[\int\exp(T\ell_{iT}(\theta,\Pi,\beta))f(\beta|\theta)d\beta],$

where
$\ell_{iT}(\theta,\Pi,\beta)=\sum_{t=1}^{T}\ell_{it}(\theta,\Pi,\beta)/T$ is
the log-likelihood for individual $i$ using observations $t=1,\dots,T$ for a
given set of parameters $(\theta,\Pi,\beta)$, $f(\beta|\theta)$ is the mixture
distribution. The gradient $G_{nT}$ of the objective function $Q_{nT}$ is
given by:

$G_{nT}(\theta,\Pi)=-\frac{1}{nT}\sum_{i=1}^{n}\frac{\int(T\partial_{\theta}\ell_{iT}(\theta,\Pi,\beta)+\partial_{\theta}\log[f(\beta|\theta)])\exp(T\ell_{iT}(\theta,\Pi,\beta)+\log[f(\beta|\theta)])d\beta}{\int\exp(T\ell_{iT}(\theta,\Pi,\beta)+\log[f(\beta|\theta)])d\beta}.$

Under regularity conditions, the estimator $\hat{\theta}_{nT}$ is
$\sqrt{n}$-consistent and asymptotically normal, with $T$ fixed with second
order representation:111 The estimator is $\sqrt{n}$ not $\sqrt{nT}$
asympotically normal because of the dependence over $t$ for each $i$ due to
$\beta_{i}$. See Rilstone et al. (1996) for general results on second-order
bias in nonlinear models.

$\sqrt{n}\left(\hat{\theta}_{nT}-\theta^{\dagger}\right)=-[H_{nT}(\theta^{\dagger},\Pi^{\dagger})]^{-1}G_{nT}(\theta^{\dagger},\Pi^{\dagger})+\frac{1}{\sqrt{n}}B_{nT}(\theta^{\dagger},\Pi^{\dagger})+o_{p}(n^{-1/2}),$

and $B_{nT}(\theta^{\dagger},\Pi^{\dagger})$ converges in probability to some
limit which depends on $T$. The main idea for the following is to capture some
of the terms in $B_{nT}(\theta^{\dagger},\Pi^{\dagger})$ using a large-T
asymptotic framework where $T$ and $n$ grow at the same rate. The main idea
here is to approximate $\sqrt{n}G_{nT}(\theta^{\dagger},\Pi^{\dagger})$ as a
asymptotically normal term which does not involve the log/integral/exponential
transformation plus a $1/\sqrt{T}\asymp 1/\sqrt{n}$ term which accounts for
the transformation. The formula is very similar to the incidental parameter
bias but is of second rather than first order. To further simplify the
derivations below, the following assumption will be used.

######  Assumption H7.

The functions $\ell_{iT}(\theta,\Pi,\beta)$,
$\partial_{\theta}\ell_{iT}(\theta,\Pi,\beta)$, $\log[f(\beta|\theta)]$ are
quadratic in $\beta$ for each $\theta,\Pi$.

Let $\beta_{i}$ be the true value of $\beta$ for each individual $i$, using
the change of variable $\beta=\beta_{i}+hT^{-1/2}$, re-write:

$\displaystyle G_{nT}(\theta,\Pi)=$
$\displaystyle-\frac{1}{n}\sum_{i=1}^{n}\partial_{\theta}\ell_{iT}(\theta,\Pi,\beta_{i})$
(H.8)
$\displaystyle-\frac{1}{n\sqrt{T}}\sum_{i=1}^{n}\partial^{2}_{\theta,\beta}\ell_{iT}(\theta,\Pi,\beta_{i})\frac{\int
h\exp(T\ell_{iT}(\theta,\Pi,\beta_{i}+hT^{-1/2})+\log[f(\beta_{i}+hT^{-1/2})])dh}{\int\exp(T\ell_{iT}(\theta,\Pi,\beta_{i}+hT^{-1/2})+\log[f(\beta_{i}+hT^{-1/2})])dh}$
(H.9)
$\displaystyle-\frac{1}{2nT}\sum_{i=1}^{n}\sum_{j=1}^{d_{\theta}}\partial^{3}_{\theta,\beta,\beta_{j}}\ell_{iT}(\theta,\Pi,\beta_{i})\frac{\int
hh_{j}\exp(T\ell_{iT}(\theta,\Pi,\beta_{i}+hT^{-1/2})+\log[f(\beta_{i}+hT^{-1/2})])dh}{\int\exp(T\ell_{iT}(\theta,\Pi,\beta_{i}+hT^{-1/2})+\log[f(\beta_{i}+hT^{-1/2})])dh}.$
(H.10)

The first term, (H.8), is $\sqrt{n}$-asymptotically normal; it corresponds to
a likelihood where all the heterogeneity is fully observed. The second term,
(H.9), is the linear expansion term of the transformation and is also
$\sqrt{n}$-asymptotically normal; i.e. it contributes to the asymptotic
variance by taking into account the estimation of the unobserved heterogeneous
distribution. The third term, (H.9), is the quadratic expansion term of the
transformation which accounts for the non-linear effect of the transformation.
It will contribute to the second-order asymptotic bias. Notice that, with the
simplifying assumptions, the terms inside the exponential function are
quadratic so that the corresponding distribution is Gaussian. To derive the
mean and variance-covariance matrix of the Gaussian let
$Q_{iT}(\theta,\Pi,\beta_{i}+hT^{-1/2})=-\ell_{iT}(\theta,\Pi,\beta_{i}+hT^{-1/2})-\frac{1}{T}\log[f(\beta_{i}+hT^{-1/2})]$
and $G_{iT}(\theta,\Pi,\beta_{i}),H_{iT}(\theta,\Pi,\beta_{i})$ be its
gradient and hessian. Given the quadratic assumption, we have:

$\displaystyle\exp(-TQ_{iT}(\theta,\Pi,\beta_{i}+hT^{-1/2}))=$
$\displaystyle\exp\left(-TQ_{iT}(\theta,\Pi,\beta_{i})+\frac{T}{2}G_{iT}(\theta,\Pi,\beta_{i})^{\prime}H_{iT}(\theta,\Pi,\beta_{i})G_{iT}(\theta,\Pi,\beta_{i})\right)$
$\displaystyle\times\exp\left(-\frac{1}{2}\left[h+\sqrt{T}H_{iT}(\theta,\Pi,\beta_{i})^{-1}G_{iT}(\theta,\Pi,\beta_{i})\right]^{\prime}H_{iT}(\theta,\Pi,\beta_{i})\left[h+\sqrt{T}H_{iT}(\theta,\Pi,\beta_{i})^{-1}G_{iT}(\theta,\Pi,\beta_{i})\right]\right),$

which is, up to normalizing constants, equal to the density of a multivariate
Gaussian with mean
$-\sqrt{T}H_{iT}(\theta,\Pi,\beta_{i})^{-1}G_{iT}(\theta,\Pi,\beta_{i})$ and
variance $H_{iT}(\theta,\Pi,\beta_{i})^{-1}$. Notice that this is the
asymptotic distribution of $\sqrt{T}(\hat{\beta}_{i}-\beta_{i})$, the
individual fixed effect estimator of $\beta$ when the information matrix
equality holds. This implies that (H.9) simplifies to:

$(\ref{eq:expr2})=\frac{1}{n}\sum_{i=1}^{n}\partial^{2}_{\theta,\beta}\ell_{iT}(\theta,\Pi,\beta_{i})H_{iT}(\theta,\Pi,\beta_{i})^{-1}G_{iT}(\theta,\Pi,\beta_{i}).$

Likewise, (H.10) simplifies to:

$(\ref{eq:expr3})=-\frac{1}{2nT}\sum_{i=1}^{n}\sum_{j=1}^{d_{\theta}}\partial^{3}_{\theta,\beta,\beta_{j}}\ell_{iT}(\theta,\Pi,\beta_{i})[H_{iT}(\theta,Pi,\beta_{i})^{-1}]_{j},$

where $[H_{iT}(\theta,Pi,\beta_{i})^{-1}]_{j}$ is the $j$-th column of
$H_{iT}(\theta,Pi,\beta_{i})^{-1}$.

######  Assumption H8.

Suppose that as both $n,T\to\infty$:

$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\sum_{t=1}^{T}[\partial_{\theta}\ell_{it}(\theta^{\dagger},\Pi^{\dagger})+\partial^{2}_{\theta,\beta}\ell_{iT}(\theta^{\dagger},\Pi^{\dagger},\beta_{i})H_{iT}(\theta^{\dagger},\Pi^{\dagger},\beta_{i})^{-1}\sqrt{T}H_{iT}(\theta^{\dagger},\Pi^{\dagger},\beta_{i})]\overset{d}{\to}\mathcal{N}(0,\Sigma),$

and that:

$\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{d_{\theta}}\partial^{3}_{\theta^{\dagger},\beta,\beta_{j}}\ell_{iT}(\theta^{\dagger},\Pi^{\dagger},\beta_{i})[H_{iT}(\theta^{\dagger},Pi^{\dagger},\beta_{i})^{-1}]_{j}\overset{p}{\to}2\mathcal{B}(\theta^{\dagger},\Pi^{\dagger}),$

where
$\mathcal{B}(\theta^{\dagger},\Pi^{\dagger})=\mathbb{E}\left(\sum_{j=1}^{d_{\theta}}\partial^{3}_{\theta^{\dagger},\beta,\beta_{j}}\ell_{iT}(\theta^{\dagger},\Pi^{\dagger},\beta_{i})[H_{iT}(\theta^{\dagger},Pi^{\dagger},\beta_{i})^{-1}]_{j}\right)/2.$

######  Proposition H2 (Asymptotic Bias).

Suppose Assumptions H7-H8 hold and $n/T\to\kappa^{2}\in(0,+\infty)$, then:

$\displaystyle\sqrt{n}G_{nT}(\theta^{\dagger},\Pi^{\dagger})\overset{d}{\to}\mathcal{N}(0,\Sigma),$
(H.11)
$\displaystyle\sqrt{nT}\left(G_{nT}(\theta^{\dagger},\Pi^{\dagger})-(\ref{eq:expr1})-(\ref{eq:expr2})\right)\overset{p}{\to}-\kappa\mathcal{B}(\theta^{\dagger},\Pi^{\dagger}).$
(H.12)

Suppose $\hat{\theta}_{nT}$ satisfies the asymptotic expansion:

$\displaystyle\sqrt{n}(\hat{\theta}_{nT}-\theta^{\dagger})=-\sqrt{n}[H_{\infty}(\theta^{\dagger},\Pi^{\dagger})]^{-1}G_{nT}(\theta^{\dagger},\Pi^{\dagger})+\frac{1}{\sqrt{n}}B_{nT}(\theta^{\dagger},\Pi^{\dagger})+o_{p}(n^{-1/2}),$
(H.13)

where
$H_{nT}(\theta^{\dagger},\Pi^{\dagger})\overset{p}{\to}H_{\infty}(\theta^{\dagger},\Pi^{\dagger})$
positive definite. The estimator is asymptotically normal:

$\displaystyle\sqrt{n}\left(\hat{\theta}_{nT}-\theta^{\dagger}\right)\overset{d}{\to}\mathcal{N}(0,V),$
(H.14)

where
$V=[H_{\infty}(\theta^{\dagger},\Pi^{\dagger})]^{-1}\Sigma[H_{\infty}(\theta^{\dagger},\Pi^{\dagger})]^{-1}$.
The second-order bias can be decomposed in two:

$\displaystyle\sqrt{nT}\left(\hat{\theta}_{nT}-\theta^{\dagger}-[H_{nT}(\theta^{\dagger},\Pi^{\dagger})]^{-1}[(\ref{eq:expr1})+(\ref{eq:expr2})]\right)\overset{p}{\to}$
$\displaystyle\kappa[H_{\infty}(\theta^{\dagger},\Pi^{\dagger})]^{-1}\mathcal{B}(\theta^{\dagger},\Pi^{\dagger})$
$\displaystyle+\kappa^{-1}\text{plim}_{n,T\to\infty}B_{nT}(\theta^{\dagger},\Pi^{\dagger})$

The asymptotic expansion condition (H.13) can be verified using conditions in
Rilstone et al. (1996), the convergence in probability of the Hessian can be
derived from a law of large numbers. The positive definiteness is a local
identification condition.

###### Proof.

Note that
$\sqrt{n}G_{nT}(\theta^{\dagger},\Pi^{\dagger})=\sqrt{n}[(\ref{eq:expr1})+(\ref{eq:expr2})]+\sqrt{n}(\ref{eq:expr3})$.
Assumption H8 implies that the first term converges in distribution to a
Gaussian distribution with zero mean and asymptotic variance-covariance matrix
$\Sigma$. Using the same Assumption, the second term converges in probability
to zero. Using the same identity,
$\sqrt{nT}[G_{nT}(\theta^{\dagger},\Pi^{\dagger})-(\ref{eq:expr1})-(\ref{eq:expr2})]=\sqrt{nT}(\ref{eq:expr3})\overset{p}{\to}-\kappa\mathcal{B}(\theta^{\dagger},\Pi^{\dagger})$
which is the first result. Then (H.12) comes from Assumption H8. Pre-
multiplying (H.12) by
$H_{nT}(\theta^{\dagger},\Pi^{\dagger})^{-1}\overset{p}{\to}H_{\infty}(\theta^{\dagger},\Pi^{\dagger})^{-1}$
and using the asymptotic expansion (H.13) then implies the last result. ∎

######  Proposition H3 (Spit Panel Bias Reduction).

Let $\hat{\theta}_{nT/2}^{1},\hat{\theta}_{nT/2}^{2}$ be the estimators
computed using the first $\lfloor T/2\rfloor$ and the last $T-\lfloor
T/2\rfloor$ time-observations, both using all individuals $i=1,\dots,n$.
Suppose the associated objectives $Q_{nT/2}^{1},Q_{nT/2}^{2}$ satisfy
Assumptions (H7)-(H8) using $\lfloor T/2\rfloor$ and $T-\lfloor T/2\rfloor$
instead of $T$. Let $\hat{\theta}^{1,2}_{nT/2}$ be either estimator, suppose
they both satisfy the asymptotic expansions:

$\displaystyle\sqrt{n}(\hat{\theta}^{1,2}_{nT/2}-\theta^{\dagger})=-\sqrt{nT}[H^{1,2}_{\infty}(\theta^{\dagger},\Pi^{\dagger})]^{-1}G^{1,2}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})+\frac{1}{\sqrt{n}}B^{1,2}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})+o_{p}(n^{-1/2}),$
(H.15)

where
$H^{1,2}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})\overset{p}{\to}H_{\infty}(\theta^{\dagger},\Pi^{\dagger})$
positive definite. Then the bias corrected estimator
$\tilde{\theta}_{nT}=2\hat{\theta}_{nT}-(\hat{\theta}^{1}_{nT/2}+\hat{\theta}^{2}_{nT/2})/2$
satisfies:

$\displaystyle\sqrt{n}\left(\tilde{\theta}_{nT}-\theta^{\dagger}\right)\overset{d}{\to}\mathcal{N}\left(0,V\right),$
(H.16)
$\displaystyle\sqrt{nT}\left(\tilde{\theta}_{nT}-\theta^{\dagger}-[H_{\infty}(\theta^{\dagger},\Pi^{\dagger})]^{-1}[(\ref{eq:expr1})+(\ref{eq:expr2})]\right)\overset{p}{\to}\kappa^{-1}\text{plim}_{n,T\to\infty}B_{nT}(\theta^{\dagger},\Pi^{\dagger}).$
(H.17)

where $V$ is the same as the asymptotic variance of $\hat{\theta}_{nT}$.

The main takeaway of the Proposition is (H.17) which indicates that the part
of the second-order bias which is associated with the integral transformed can
be removed using the split panel jackknife.

###### Proof.

Combining asymptotic expansions (H.13), (H.15), we have:

$\displaystyle\sqrt{n}\left(\tilde{\theta}_{nT}-\theta^{\dagger}\right)$
$\displaystyle=-\sqrt{n}H_{\infty}(\theta^{\dagger},\Pi^{\dagger})^{-1}\left(2G_{nT}(\theta^{\dagger},\Pi^{\dagger})-[G^{1}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})+G^{2}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})]/2\right)$
$\displaystyle+\frac{1}{\sqrt{n}}\left(2B_{nT}-[B_{nT/2}^{1}+B_{nT/2}^{2}]/2\right).$

Notice that:

$\displaystyle
2G_{nT}(\theta^{\dagger},\Pi^{\dagger})-[G^{1}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})+G^{2}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})]/2=(\ref{eq:expr1})+(\ref{eq:expr2})+(\ref{eq:extra}),$

where (H.18) is defined below. The identity comes from additivity over $t$ of
$(\ref{eq:expr1}),(\ref{eq:expr2})$. The (H.18) term comes from combining the
(H.10) terms in each one of $G_{nT}$, $G^{1}_{nT/2}$, and $G^{2}_{nT/2}$:

$\displaystyle=-2\frac{1}{2nT}$
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{d_{\theta}}\partial^{3}_{\theta,\beta,\beta_{j}}\ell_{iT}(\theta,\Pi,\beta_{i})[H_{iT}(\theta,Pi,\beta_{i})^{-1}]_{j}$
(H.18)
$\displaystyle-\frac{1}{2}\frac{1}{2nT/2}\sum_{i=1}^{n}\sum_{j=1}^{d_{\theta}}\partial^{3}_{\theta,\beta,\beta_{j}}\ell^{1}_{iT/2}(\theta,\Pi,\beta_{i})[H^{2}_{iT/2}(\theta,Pi,\beta_{i})^{-1}]_{j}$
$\displaystyle-\frac{1}{2}\frac{1}{2nT/2}\sum_{i=1}^{n}\sum_{j=1}^{d_{\theta}}\partial^{3}_{\theta,\beta,\beta_{j}}\ell^{2}_{iT/2}(\theta,\Pi,\beta_{i})[H^{2}_{iT/2}(\theta,Pi,\beta_{i})^{-1}]_{j}.$

Using the assumptions and $\sqrt{n/[T/2]}\to\sqrt{2}\kappa$,
$\sqrt{nT}(\ref{eq:extra})$ has the following limit:

$\displaystyle\sqrt{nT}(\ref{eq:extra})$
$\displaystyle=-2\sqrt{n/T}\frac{1}{2n}\sum_{i=1}^{n}\sum_{j=1}^{d_{\theta}}\partial^{3}_{\theta,\beta,\beta_{j}}\ell_{iT}(\theta,\Pi,\beta_{i})[H_{iT}(\theta,Pi,\beta_{i})^{-1}]_{j}$
$\displaystyle-\frac{1}{2}\sqrt{2}\sqrt{n/[T/2]}\frac{1}{2nT/2}\sum_{i=1}^{n}\sum_{j=1}^{d_{\theta}}\partial^{3}_{\theta,\beta,\beta_{j}}\ell^{1}_{iT/2}(\theta,\Pi,\beta_{i})[H^{2}_{iT/2}(\theta,Pi,\beta_{i})^{-1}]_{j}$
$\displaystyle-\frac{1}{2}\sqrt{2}\sqrt{n/[T/2]}\frac{1}{2nT/2}\sum_{i=1}^{n}\sum_{j=1}^{d_{\theta}}\partial^{3}_{\theta,\beta,\beta_{j}}\ell^{2}_{iT/2}(\theta,\Pi,\beta_{i})[H^{2}_{iT/2}(\theta,Pi,\beta_{i})^{-1}]_{j}.$
$\displaystyle\overset{p}{\to}2[-2\kappa\mathcal{B}(\theta^{\dagger},\Pi^{\dagger})]-\sqrt{2}[-2\sqrt{2}\kappa\mathcal{B}(\theta^{\dagger},\Pi^{\dagger})]=0.$

Putting all the derivations together, we get the two sets of results:

$\displaystyle\sqrt{n}\left(2G_{nT}(\theta^{\dagger},\Pi^{\dagger})-[G^{1}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})+G^{2}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})]/2\right)=\sqrt{n}[H_{\infty}(\theta^{\dagger},\Pi^{\dagger})]^{-1}[(\ref{eq:expr1})+(\ref{eq:expr2})]+o_{p}(1),$

from which the asymptotic normality result follows. Also,

$\displaystyle\sqrt{nT}\left(2G_{nT}(\theta^{\dagger},\Pi^{\dagger})-[G^{1}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})+G^{2}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})]/2-[H_{\infty}(\theta^{\dagger},\Pi^{\dagger})]^{-1}[(\ref{eq:expr1})+(\ref{eq:expr2})]\right)$
$\displaystyle=o_{p}(1)+\kappa^{-1}\text{plim}_{n,T\to\infty}B_{nT}(\theta^{\dagger},\Pi^{\dagger}),$

where the $o_{p}(1)$ term includes the bias reduction in (H.18). ∎

#### A split panel implementation of rqn.

The split panel bias corrected estimator is
$\tilde{\theta}_{nT}=2\hat{\theta}_{nT}-[\hat{\theta}^{1}_{nT/2}+\hat{\theta}^{2}_{nT/2}]/2$.
The goal is to target $\tilde{\theta}_{nT}$ using rqn. Let
$\theta_{b,nT},\theta^{1}_{b,nT/2},\theta^{2}_{b,nT/2}$ be rqn draws computed
using the full and split panels, respectively but the same resampled
individuals $i$. Using Proposition 1, we have the three couplings:

$\displaystyle\theta^{\star}_{b+1,nT}$
$\displaystyle=\hat{\theta}_{nT}+(1-\gamma)(\theta^{\star}_{b,nT}-\hat{\theta}_{nT})-\gamma[H_{nT}(\theta^{\dagger},\Pi^{\dagger})]^{-1}G_{nT}(\theta^{\dagger},\Pi^{\dagger})$
$\displaystyle\theta^{\star,1}_{b+1,nT/2}$
$\displaystyle=\hat{\theta}^{1}_{nT/2}+(1-\gamma)(\theta^{\star,1}_{b,nT/2}-\hat{\theta}^{1}_{nT/2})-\gamma[H^{1}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})]^{-1}G^{1}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})$
$\displaystyle\theta^{\star,2}_{b+1,nT/2}$
$\displaystyle=\hat{\theta}^{2}_{nT/2}+(1-\gamma)(\theta^{\star,2}_{b,nT/2}-\hat{\theta}^{2}_{nT/2})-\gamma[H^{2}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})]^{-1}G^{2}_{nT/2}(\theta^{\dagger},\Pi^{\dagger}).$

Compute
$\tilde{\theta}_{b+1,nT}=2\theta_{b+1,nT}-[\theta^{1}_{b+1,nT/2}+\theta^{2}_{b+1,nT/2}]/2$,
we have the associated coupling:

$\displaystyle\tilde{\theta}^{\star}_{b+1,nT}$
$\displaystyle=\tilde{\theta}_{nT}+(1-\gamma)(\tilde{\theta}^{\star}_{b,nT}-\tilde{\theta}_{nT})-\gamma\Big{(}2[H_{nT}(\theta^{\dagger},\Pi^{\dagger})]^{-1}G_{nT}(\theta^{\dagger},\Pi^{\dagger})$
$\displaystyle-\Big{[}[H^{1}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})]^{-1}G^{1}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})+[H^{2}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})]^{-1}G^{2}_{nT/2}(\theta^{\dagger},\Pi^{\dagger})\Big{]}/2\Big{)},$

which has innovations that match the first-order expansion for
$\tilde{\theta}_{nT}$.

###  H.2 Example 2: A Probit with Many Regressors

Both rqn and MCMC are initialized at the same starting value $\theta_{0}$
which is computed as follows. A linear probability model is estimated by OLS
yielding coefficients $\tilde{\theta}_{n,\textsc{ols}}$. The starting value
$\theta_{0}$ is then computed using the linear approximation:
$x_{i}^{\prime}\tilde{\theta}_{n,\textsc{ols}}\simeq\Phi(x_{i}^{\prime}\theta_{0})\simeq\Phi(0)+\phi(0)x_{i}^{\prime}\theta_{0}$,
where $\phi,\Phi$ are the Gaussian pdf and cdf, $\Phi(0)=1/2$,
$\phi(0)=[2\pi]^{-1/2}$. For the intercept, this yields
$\theta_{0,\text{const}}=[\tilde{\theta}_{n,\text{const},\textsc{ols}}-\Phi(0)]/\phi(0)$.
For all other coefficients this yields
$\theta_{0,j}=\tilde{\theta}_{n,j,\textsc{ols}}/\phi(0)$.

#### Implementation of the MALA algorithm.

For Bayesian inference, let $Q_{n}$ be the sample negative log-likelihood and
$\pi$ the prior distribution. Let $G_{n}$ be the gradient of the negative log-
posterior distribution and $H_{n}$ its Hessian.222The negative log-likelihood
is minus one times the log-likelihood. The log-posterior distribution is the
log-likelihood plus the log-prior distribution. Given a draw $\theta_{b}$, the
baseline MALA algorithm produces the next draw $\theta_{b+1}$ using the
following steps:

$\displaystyle\text{Draw }\tilde{\theta}=\theta_{b}-\gamma
G_{n}(\theta_{b})+\sqrt{2\gamma}Z_{b},\,Z_{b}\sim\mathcal{N}(0,I),\text{ and
}u\sim\mathcal{U}_{[0,1]},$ $\displaystyle\text{If
}u\leq\frac{\exp[-Q_{n}(\tilde{\theta})]\pi(\tilde{\theta})}{\exp[-Q_{n}(\theta_{b})]\pi(\theta_{b})}\frac{q(\tilde{\theta}|\theta_{b})}{q(\theta_{b}|\tilde{\theta})},\text{
set }\theta_{b+1}=\tilde{\theta}.\text{ Otherwise, set
}\theta_{b+1}=\theta_{b},$

where $q(\tilde{\theta}|\theta_{b})=\phi([\tilde{\theta}-\theta_{b}+\gamma
G_{n}(\theta_{b})]/\sqrt{2\gamma})$ and
$q(\theta_{b}|\tilde{\theta})=\phi([\theta_{b}-\tilde{\theta}+\gamma
G_{n}(\tilde{\theta})]/\sqrt{2\gamma})$ are the transition probabilities,
$\phi$ is the pdf of the multivariate Gaussian. Compared to the random-walk
Metropolis-Hastings algorithm which uses the proposal
$\tilde{\theta}=\theta_{b}+\sqrt{2\gamma}Z_{b}$, MALA adds the gradient
descent term, $-\gamma G_{n}(\theta_{b})$, which directs the proposal towards
the maximum of the posterior distribution to get faster convergence and a
higher acceptance rate. Here, the proposal is modified for more direct
comparison with rnr, rqn. It requires computing $H_{n}(\hat{\theta}_{n})$:

$\tilde{\theta}=\theta_{b}-\gamma[H_{n}(\hat{\theta}_{n})]^{-1}G_{n}(\theta_{b})+\sqrt{2\gamma}Z_{b},\quad
Z_{b}\sim\mathcal{N}(0,[H_{n}(\hat{\theta}_{n})]^{-1}),$

the transition probability $q(\cdot|\cdot)$ is adjusted accordingly. Following
Roberts and Rosenthal (1998), $\gamma$ is tuned to get an acceptance rate of
0.57, after the initial convergence phase. This requires
$\gamma=O(1/d_{\theta})$ as illustrated below. The desired acceptance rate is
achieved for $\gamma\simeq 0.0034$ to be compared with $\gamma=0.1$ used for
rnr and rqn.

#### Choice of learning rate $\gamma$ for the MALA algorithm.

Recall the modified proposal:

$\tilde{\theta}=\theta_{b}-\gamma[H_{n}(\hat{\theta}_{n})]^{-1}G_{n}(\theta_{b})+\sqrt{2\gamma}Z_{b},\quad
Z_{b}\sim\mathcal{N}(0,[H_{n}(\hat{\theta}_{n})]^{-1}).$

To illustrate how the dimension of the parameters $d_{\theta}$ affects the
choice of $\gamma$, consider the case wher $\theta_{b}=\hat{\theta}_{n}$ so
that $G_{n}(\theta_{b})=0$. This implies that conditional on
$\theta_{b}=\hat{\theta}_{n}$,
$\tilde{\theta}=\hat{\theta}_{n}+\sqrt{2\gamma}Z_{b}\sim\mathcal{N}(\hat{\theta}_{n},2\gamma[H_{n}(\hat{\theta}_{n})]^{-1})$.
The log-posterior is approximately:

$Q_{n}(\tilde{\theta})-\log(\pi(\tilde{\theta}))\simeq
Q_{n}(\hat{\theta}_{n})-\log(\pi(\hat{\theta}_{n}))+\gamma
Z_{b}^{\prime}H_{n}(\hat{\theta}_{n})Z_{b},$

where $Z_{b}^{\prime}H_{n}(\hat{\theta}_{n})Z_{b}\sim\chi^{2}_{p}$. For
$d_{\theta}$ large, we have
$Z_{b}^{\prime}H_{n}(\hat{\theta}_{n})Z_{b}/d_{\theta}=1+o_{p}(1).$ Plug this
into the accept/reject step of the algorithm:

$\mathbb{P}\left(u\leq\frac{\exp[-Q_{n}(\tilde{\theta})]\pi(\tilde{\theta})}{\exp[-Q_{n}(\hat{\theta}_{n})]\pi(\hat{\theta}_{n})}\right)\simeq\mathbb{P}\left(u\leq\exp[-\gamma
d_{\theta}]\right)=\exp[-\gamma d_{\theta}],$

which equals $0.57$ for $\gamma=-\log(0.57)/d_{\theta}$. At
$\theta_{b}=\hat{\theta}_{n}$, the optimal choice of $\gamma$ is inversely
proportional to $d_{\theta}$. For $\theta_{b}\neq\hat{\theta}_{n}$, the
acceptance rate is higher because the gradient term sets $\tilde{\theta}$ in
directions of increasing posterior, on average. This implies that a larger
value $\gamma>-\log(0.57)/d_{\theta}$ should be used to get the desired
acceptance rate overall. In the empirical application, the learning rate was
set to $\gamma=-2\log(0.57)/d_{\theta}$ which resulted in the desired
acceptance rate.

#### Additional Empirical Results and Comparisons.

Table H7: Probit - Estimates and Standard Errors

| const | distance | border | island | landlock | legal | language | colonial | currency | FTA | religion
---|---|---|---|---|---|---|---|---|---|---|---
| Estimates |
mle | 2.537 | -0.618 | -0.380 | 0.355 | 0.220 | 0.072 | 0.275 | 0.288 | 0.530 | 1.854 | 0.249
rnr | 2.541 | -0.619 | -0.381 | 0.355 | 0.219 | 0.071 | 0.275 | 0.287 | 0.529 | 1.863 | 0.249
rqn | 2.539 | -0.618 | -0.380 | 0.355 | 0.219 | 0.071 | 0.275 | 0.286 | 0.529 | 1.860 | 0.249
mcmc | 2.533 | -0.617 | -0.380 | 0.355 | 0.215 | 0.072 | 0.274 | 0.289 | 0.534 | 1.854 | 0.250
| Standard Errors (iid) |
mle | 0.068 | 0.008 | 0.033 | 0.023 | 0.031 | 0.009 | 0.012 | 0.075 | 0.039 | 0.101 | 0.018
rnr | 0.066 | 0.008 | 0.035 | 0.022 | 0.030 | 0.009 | 0.011 | 0.077 | 0.038 | 0.118 | 0.017
rqn | 0.065 | 0.008 | 0.034 | 0.021 | 0.030 | 0.009 | 0.011 | 0.073 | 0.038 | 0.115 | 0.017
mcmc | 0.040 | 0.004 | 0.020 | 0.015 | 0.020 | 0.006 | 0.008 | 0.060 | 0.028 | 0.069 | 0.012
| Standard Errors (clustered) |
mle | 0.094 | 0.015 | 0.039 | 0.031 | 0.030 | 0.012 | 0.016 | 0.073 | 0.068 | 0.103 | 0.029
rnr | 0.087 | 0.015 | 0.038 | 0.032 | 0.030 | 0.011 | 0.014 | 0.073 | 0.069 | 0.102 | 0.028
rqn | 0.105 | 0.015 | 0.040 | 0.034 | 0.036 | 0.011 | 0.016 | 0.105 | 0.071 | 0.134 | 0.030

Legend: rnr/rqn computed using $\gamma=0.1$, Gaussian (iid) and exponential
(clustered) weights.
MCMC computed using $\gamma=0.0034$, Gaussian prior with mean zero, variance
10 for fixed effects, flat uniform prior on other parameters.

Figure H4: Stochastic Gradient Descent

Legend: x-axis $=k$ in millions,
$P_{b}=I_{d},m=1,\gamma_{k}=\gamma_{0}k^{-5/8}$, 5 million iterations.

Figure H5: Infeasible Stochastic Newton-Raphson

Legend: x-axis $=k$ in thousands,
$P_{b}=[H_{n}(\hat{\theta}_{n})+\gamma_{k}I_{d}]^{-1}$,
$m=1,\gamma_{k}=\gamma_{0}k^{-5/8}$, $100$ thousand iterations.

###  H.3 Example 3: NLS Estimation of Transportation Costs

#### Implementation of rnr.

The objective function can be written as:

$Q_{n}(\theta)=\frac{1}{n}\sum_{odt}[y_{odt}-f(x_{odt},\theta)]^{2},$

with $\theta=(\alpha,\delta,\beta)$, $\beta$ are the intercept and fixed
effect coefficients and $x_{odt}$ consists of $R_{t}$ and the intercept, fixed
effect dummies. $n$ is the total number of observations. The gradient and
hessian are:

$\displaystyle G_{n}(\theta)$
$\displaystyle=-\frac{2}{n}\sum_{odt}\partial_{\theta}f(x_{odt},\theta)[y_{odt}-f(x_{odt},\theta)],$
$\displaystyle H_{n}(\theta)$
$\displaystyle=\frac{2}{n}\sum_{odt}\partial_{\theta}f(x_{odt},\theta)\partial_{\theta}f(x_{odt},\theta)^{\prime}-\frac{2}{n}\sum_{odt}\partial^{2}_{\theta,\theta^{\prime}}f(x_{odt},\theta)[y_{odt}-f(x_{odt},\theta)].$

If $\mathbb{E}\left(y_{odt}-f(x_{odt},\theta^{\dagger})|x_{odt}\right)=0$,
then
$\hat{H}_{n}(\hat{\theta}_{n})=\frac{2}{n}\sum_{odt}\partial_{\theta}f(x_{odt},\hat{\theta}_{n})\partial_{\theta}f(x_{odt},\hat{\theta}_{n})^{\prime}$
is a consistent estimator of the population hessian
$H(\theta^{\dagger})=2\mathbb{E}[\partial_{\theta}f(x_{odt},\theta^{\dagger})\partial_{\theta}f(x_{odt},\theta^{\dagger})^{\prime}]$.

Using this estimator, a cost-effective implementation of rnr with $m=n$ and
multiplier weights clustered at the district level $d$ relies on:

$\displaystyle G_{m}^{(b)}(\theta)$
$\displaystyle=-\frac{2}{n}\sum_{odt}w^{(b)}_{d}\partial_{\theta}f(x_{odt},\theta)[y_{odt}-f(x_{odt},\theta)],$
$\displaystyle H_{m}^{(b)}(\theta)$
$\displaystyle=\phantom{-}\frac{2}{n}\sum_{odt}w^{(b)}_{d}\partial_{\theta}f(x_{odt},\theta)\partial_{\theta}f(x_{odt},\theta)^{\prime},$

where $w^{(b)}_{d}\sim\mathcal{N}(1,1)$ are iid over districts $d$ and
iterations $b$. The appeal of this approach is that $H_{m}^{(b)}(\theta)$ is
symmetric positive semidefinite for all $\theta$, even if there are values for
which $H_{n}(\theta)$ is non-definite. The added penalty ensures that
$H_{m}^{(b)}(\theta)+\lambda\partial^{2}_{\theta,\theta^{\prime}}\text{pen}(\theta)/n$
is symmetric positive definite. The derivative
$\partial_{\theta}f(x_{odt},\theta)$ is computed once for each iteration $b$,
and then used to evaluate both the gradient and hessian. More specifically
$\partial_{\alpha}f(x_{odt},\theta)$ is computed using finite differences and
$\partial_{(\delta,\beta)}f(x_{odt},\theta)$ is simply the vector of linear
regressors, i.e. $\log(\text{LCRED}(R_{t},\alpha)_{odt})$ and the intercept
and fixed effect dummies. The gradient and hessian of the penalty term are
then added to the re-weighted gradient and hessian above.

##  Appendix I Convergence in Some Non-Convex Settings

###  I.1 Recovering from a bad start

Lemmas 1 and 2 require the user-chosen learning rate $\gamma$ to satisfy
$A(\gamma)<1$ to get the desired contraction property, which leads to
convergence. In practice, feasible values of $\gamma$, depend on both the
choice of $P_{b}$ and the hessian $H_{n}$. The following discussion will focus
on nr in the context of Lemma 1 for simplicity. For values of $\gamma$ that
are too large, $A(\gamma)<1$ may not hold and convergence may fail. The
following illustrates how introducing a certain quadratic penalty in the first
few iterations can restore the contraction property.

From the starting value $\theta_{0}$, the firt nr iteration has the form:
$\theta_{1}=\theta_{0}-\gamma[H_{n}(\theta_{0})]^{-1}G_{n}(\theta_{0}),$ which
is only well-defined if $H_{n}(\theta_{0})$ is non-singular. Using the mean-
value theorem:
$\theta_{1}-\hat{\theta}_{n}=[I_{d}-\gamma[H_{n}(\theta_{0})]^{-1}H_{n}(\tilde{\theta}_{0})](\theta_{0}-\hat{\theta}_{n}),$
for some intermediate value $\tilde{\theta}_{0}$ between $\theta_{0}$ and
$\hat{\theta}_{n}$. A contraction only occurs if
$\sigma_{\text{max}}[I_{d}-\gamma[H_{n}(\theta_{0})]^{-1}H_{n}(\tilde{\theta}_{0})]<1$.
For $\|\theta_{0}-\hat{\theta}_{n}\|$ small, continuity of the hessian implies
$[H_{n}(\theta_{0})]^{-1}H_{n}(\tilde{\theta}_{0})]$ is close to the identity
matrix $I_{d}$ so that a small $\gamma<1$ always leads to a contraction closer
to $\hat{\theta}_{n}$.

Issues can arise for distant starting values $\theta_{0}$. Suppose
$H_{n}(\theta_{0})$ is close to singular, then
$[H_{n}(\theta_{0})]^{-1}H_{n}(\tilde{\theta}_{0})]$ can be very large and
$\gamma$ needs to be very small to get a contraction. In that situation, nr
iterations are poorly behaved if $\gamma$ is not sufficiently small. As a
solution, add a quadratic penalty:
$Q_{n}(\theta)+\frac{\lambda}{2n}\|\theta-\theta_{0}\|_{2}^{2}.$ Penalizing
towards $\theta_{0}$ implies the gradient is unchaged but the hessian becomes
$H_{n}(\theta_{0})+\frac{\lambda}{n}I_{d}$. The first nr iteration is now:
$\theta_{1}=\theta_{0}-\gamma[H_{n}(\theta_{0})+\lambda/nI_{d}]^{-1}G_{n}(\theta_{0}),$
which is always well-defined as long as $H_{n}(\theta_{0})$ is positive semi-
definite.333If $H_{n}(\theta_{0})$ has negative eigenvalues it is preferable
to use the positive semi-definite
$[H_{n}(\theta_{0})^{\prime}H_{n}(\theta_{0})]^{1/2}$ instead, as illustrated
in the next subsection. For the same intermediate value, we now have:
$\theta_{1}-\hat{\theta}_{n}=[I_{d}-\gamma[H_{n}(\theta_{0})+\lambda/nI_{d}]^{-1}H_{n}(\tilde{\theta}_{0})](\theta_{0}-\hat{\theta}_{n}),$
which even for $\gamma=1$ is a contraction, provided $\lambda$ is sufficiently
large. Note that the update above coincides with an iteration using the trust-
region algorithm.444The trust-region method is an alternative to line-search
which can improve on standard nr iterations, see Nocedal and Wright (2006,
Ch4.1) and Bertsekas (2016, Ch1.4.2).  The modification is also similar to the
$\tau_{b}$ used in the quasi-Newton update proposed in the paper. A simple
strategy to is then to pick $\gamma\in(0,1)$, and add a large penalty
$\lambda$ in the first couple of iterations. Valid inference requires the
penalty to be asymptotically negligible, so after the burn-in one can either
set $\lambda=0$, use a schedule where $\lambda_{b}\to 0$ sufficiently fast as
$b$ increases, or use a smaller $\lambda_{n}=o(1)$. Figure I6 below
illustrates the effect of different $\lambda$s during the initial iterations
for Example 1 in Section 5 where $H_{n}(\theta_{0})$ is singular. The
baseline, $\lambda=20$, is used in the Monte-Carlo simulations, see Appendix
H.1 for implementation details.

Figure I6: Illustration: initial iterations for different choice of $\lambda$

Legend: $\gamma=0.1$. Solid lines: black $\lambda=20$, orange $\lambda=0.1$,
blue $\lambda=100$. Dashed black line: $Q_{n}(\hat{\theta}_{n})$.

###  I.2 Suboptimal solutions

An important concern is when a classical optimizer converges to a point that
is suboptimal, i.e. is neither a local nor a global optimum. In these cases,
the user has to interve and restart an optimization at another value which can
be a time-consuming process. The following illustrates, using a stylized
example, the effect of the resampling noise in this scenario.
Consider a point $\theta_{\times}$ for which the Hessian is indefinite, i.e.
$H_{n}(\theta_{\times})$ has both strictly positive and negative
eigenvalues.555See Section 3.4 Nocedal and Wright (2006) for more details on
this issue. For simplicity, the following assumes that the objective is
quadratic:
$Q_{n}(\theta)=\frac{1}{2}(\theta-\theta_{\times})^{\prime}H_{n}(\theta-\theta_{\times})$,
with eigendecomposition $H_{n}=Q\Lambda Q^{\prime}$,
$\Lambda=\text{diag}(\lambda_{1},\dots,\lambda_{d_{\theta}})$ where
$0<\lambda_{S}\leq\lambda_{S-1}\leq\dots\leq\lambda_{1}$ and
$0>\lambda_{S+1}\geq\dots\geq\lambda_{d_{\theta}}$. There are $S$ directions
with positive curvature, and $d_{\theta}-S$ directions with negative
curvature. Let $q_{1},\dots,q_{d_{\theta}}$ be the corresponding eigenvectors.
Without constraints, the global minimum is $Q_{n}(\theta)=-\infty$ attained at
$\theta=+\infty\times q_{s}$, $s\geq S+1$, the optimizer should diverge:
$\|\theta_{b}\|\to\infty$.

#### Issues with classical Newton-Raphson:

Consider the Newton-Raphson iteration $\theta_{b+1}=\theta_{b}-\gamma
H_{n}^{-1}G_{n}(\theta_{b})$, here
$G_{n}(\theta_{b})=H_{n}(\theta_{b}-\theta_{\times})$ so that:
$\theta_{b+1}-\theta_{\times}=(1-\gamma)(\theta_{b}-\theta_{\times}),$ nr
converges exponentially fast to the suboptimal solution $\theta_{\times}$. It
is possible to improve the behaviour of nr by using a Hessian modification:
$\theta_{b+1}=\theta_{b}-\gamma(H_{n}^{2})^{-1/2}G_{n}(\theta_{b})$ which
enforces $(H_{n}^{2})^{-1/2}=Q(\Lambda^{2})^{-1/2}Q^{\prime}$ positive
definite. nr iterations are now such that:

$\theta_{b+1}-\theta_{\times}=(I+Q\Gamma
Q^{\prime})(\theta_{b}-\theta_{\times}),\quad\Gamma=\text{diag}(\underbrace{-\gamma,\dots,-\gamma}_{S\text{
times}},\underbrace{+\gamma,\dots,+\gamma}_{d_{\theta}-S\text{ times}}).$

If $\theta_{b}-\theta_{\times}$ equals zero or is orthogonal to
$q_{S+1},\dots,q_{d_{\theta}}$, then nr will converge to $\theta_{\times}$
exponentially fast. However, nr will diverge exponentially fast if
$\theta_{b}-\theta_{\times}=q_{S+1}$, for instance, because of the explosive
root $1+\gamma>1$ which is outside the unit circle.

#### The case of resampled Newton-Raphson:

The main difference between nr and rnr in this example is that rnr diverges
exponentially fast with high-probability even when
$\theta_{0}=\theta_{\times}$, which is always problematic for classical
optimizers such as gd, nr, or bfgs. To simplify notation, assume that $Q=I$
and $G_{m}^{(b)}(\theta)\sim\mathcal{N}\left(G_{n}(\theta),\Sigma/m\right)$,
$\Sigma$ non-singular, and $H_{m}^{(b)}=H_{n}$. For rnr:
$\theta_{b+1}=\theta_{b}-\gamma([H_{m}^{(b+1)}]^{2})^{-1/2}G_{m}^{(b+1)}(\theta_{b})$,
which implies:

$\theta_{b+1}-\theta_{\times}=(I+\Gamma)(\theta_{b}-\theta_{\times})+\Gamma
Z_{m}^{(b+1)},\quad
Z_{m}^{(b+1)}=G_{m}^{(b+1)}(\theta_{b})-G_{n}(\theta_{b}),$

where $\Gamma$ is the same as above. Notice that
$(\theta_{b+1}-\theta_{\times})$ follows a VAR(1) process with transition
matrix $I+\Gamma$ which has $d_{\theta}-S$ explosive roots outside the unit
circle, all equal to $1+\gamma$. Suppose $\theta_{0}=\theta_{\times}$ which is
problematic for nr. Let $\theta_{j,b+1}$ be the s-th row of $\theta_{b+1}$ and
$Z_{m,j}^{(b-j)}$ the $j$-th row of $Z_{m}^{(b-j)}$, denote its variance by
$\sigma^{2}_{j}>0$. For $j\geq S+1$ we have the AR(1) representation

$\theta_{j,b+1}-\theta_{j,\times}=\gamma(1+\gamma)^{b}\sum_{j=0}^{b}(1+\gamma)^{-j}Z_{m,j}^{(j)}\sim\mathcal{N}\left(0,\gamma^{2}(1+\gamma)^{2b}\frac{1-(1+\gamma)^{-2b-2}}{1-(1+\gamma)^{-2}}\sigma^{2}_{j}\right).$

Using Gaussianity and the Paley-Zygmund inequality we have the probability
bound:

$\mathbb{P}^{\star}\left(\frac{|\theta_{j,b+1}-\theta_{j,\times}|}{\gamma(1+\gamma)^{b}}\frac{[1-(1+\gamma)^{-2}]^{1/2}}{\sqrt{2/\pi}\sigma_{j}[1-(1+\gamma)^{-2b-2}]^{1/2}}\geq\varepsilon\right)\geq\frac{2(1-\varepsilon)^{2}}{\pi},$

for any $\varepsilon\in[0,1]$. This implies that, with high probability,
$\theta_{b+1}$ diverges exponentially fast.

#### Numerical Illustration:

To illustrate the above numerically, consider a two-dimensional optimization
problem where $\theta_{\times}=(0,0)$. We will set $\theta_{0}=q_{1}+cq_{2}$
and vary $c$ to compare the performance of R’s bfgs optimizer, which relies on
a line-search, with nr and rnr which use a fixed learning rate $\gamma=0.1$
and the hessian modification. For $c=0$, nr will converge to the saddle point
and $c\neq 0$ implies that nr will eventually diverge. Table I8 reports the
results for a range of $c$. Note that bfgs converges to $\theta_{\times}$ even
when $c\neq 0$. This is because, the bfgs update does not guarantee positive
definiteness of the conditioning matrix so that it behaves similarly to nr
without the hessian modification. nr diverges for $c\neq 0$ but slowly for
$c<1$. As predicted, rnr diverges quickly for all $c$.

Table I8: Behaviour of bfgs, nr, and rnr around a saddle point

| bfgs | nr | rnr
---|---|---|---
$c$ | $\theta_{1}$ | $\theta_{2}$ | $Q_{n}(\theta)$ | $\theta_{1}$ | $\theta_{2}$ | $Q_{n}(\theta)$ | $\theta_{1}$ | $\theta_{2}$ | $Q_{n}(\theta)$
0.0 | 0.00 | 0.00 | 0.00 | 0.01 | -0.01 | 0.00 | -0.01 | -90.42 | -91.35
0.1 | 0.05 | 0.05 | -0.00 | 4.78 | 5.30 | -0.25 | 5.30 | -85.65 | -81.95
0.5 | 0.23 | 0.25 | -0.00 | 23.88 | 26.53 | -6.37 | 26.53 | -66.54 | -49.47
1.0 | 0.45 | 0.50 | -0.00 | 47.76 | 53.06 | -25.48 | 53.06 | -42.66 | -20.33
5.0 | 2.26 | 2.51 | -0.06 | 238.80 | 265.32 | -637.09 | 265.32 | 148.37 | -245.96

Legend: bfgs: output of R’s optimizer (optim). nr and rnr: output after $50$
iterations.
True minimum is $Q_{n}=-\infty$. Saddle point is $\theta_{\times}=(0,0)$.