I want to translate the "score bootstrap" that Kline and Santos (2011) propose into R code. The algorithm is described on page 6 as follows.
STEP 1: Obtain the full sample OLS estimate $\hat\beta$ and employ it to generate the fitted score contributions $\{X_i(Y_i-X^{'}_i\hat\beta)\}^n_{i=1}$.
STEP 2: Using an i.i.d. sample of random weights $\{W_i\}_{i=1}^n$ independent of $\{Y_i,X_i\}_{i=1}^n$ and satisfying $E[W_i]=0$ and $E[W^2_i]=1$, construct a new set of perturbed score contributions $\{X_i(Y_i-X^{'}_i\hat\beta)W_i\}_{i=1}^n$.
STEP 3: Multiply the constructed perturbed score by the inverse Hessian to obtain $H_n^{-1}n^{-\frac12}\sum_i(Y_i-X_i\hat\beta)X_iW_i$ and use its distribution conditional on $\{Y_i,X_i\}_ {i=1}^n$ as an estimate of the distribution of $\sqrt{n}(\hat\beta-\beta_0)$.
STEPS 1 and 2 seem to be straightforward.
set.seed(42)
n <- 100
x <- rnorm(n)
y <- x + rnorm(n)
fit <- lm(y ~ x)
X <- model.matrix(fit)
yhat <- fit$fitted.values # Xb
S <- X*(y - yhat) # score matrix
H <- - solve(vcov(fit)) # hessian
w <- ifelse(runif(n) < .5, 1, -1) # draws from Rademacher distr.
However, when I attempt to calculate STEP 3, I'm getting an error, because the number of columns of $H^{-1}$ won't match with the number of rows of $S$.
# sqrt(n) * solve(H) %*% S * w # STEP 3
## Error in solve(H) %*% S : non-conformable arguments
Does anyone see if I'm maybe misinterpreting the math somewhere?
##Edit
I obviously made an error with the multiplication in STEP 3, where I need to take account of the $\sum$, i.e. colSums
of the perturbed scores. Thus, after STEP 3 I get this result:
solve(H) * sqrt(n) * colSums(S * w) # STEP 3
# (Intercept) x
# (Intercept) 1.15362358 -0.03490467
# x 0.04228051 -1.30034588
To verify the result I set $w=1$, i.e. I don't apply any weights, and I get essentially zero.
w <- 1
solve(H) * sqrt(n) * colSums(S * w)
# (Intercept) x
# (Intercept) -4.149004e-16 1.255345e-17
# x -3.034329e-19 9.332144e-18
This feels right, because the sum of the scores is (by definition) zero, too:
colSums(S)
# (Intercept) x
# 5.023759e-15 -1.214306e-16
###Where I'm stuck
However I'm not sure what I'm actually getting here. According to the algorithm I have an "estimate of the distribution of $\sqrt{n}(\hat\beta-\beta_0)$". I know this has to do with maximum likelihood, but how to derive my variance estimator from this? If the result is already my variance matrix, then taking the square root of the diagonal should give me standard errors similar to those of the model summary, but they're pretty off:
sqrt(diag(solve(H) * colSums(t(S) %*% S) * solve(H)))
# (Intercept) x
# 0.06225221 0.07737588
summary(fit)$coe[, 2]
# (Intercept) x
# 0.09087774 0.08766507
Is this the wrong interpretation of the estimate of $\sqrt{n}(\hat\beta-\beta_0)$? Did I make an error in interpreting the equations? What I finally need is the test statistic of one bootstrap iteration, and I am not sure how to get it out of the equation of STEP 3.
Note: I'm asking for statistical help (distinct from programming help) to understand how to derive a variance estimator from the estimate of $\sqrt{n}(\hat\beta-\beta_0)$ of STEP 3.
For sake of completeness below the formula for the test statistic of the score bootstrap given in the paper on page 7.
$$T_n^{*score}:=\left(H_n^{-1}\sum_{n}^*(\hat\beta)H^{-1}_n\right)^{-\frac12}H^{-1}_n\frac1{\sqrt{n}}\sum_{i=1}^nX_i\epsilon_i^*,\quad where \quad\sum_n^*(\beta) := \frac1n\sum_iX_iX_i^T(Y_i^*-X_i^T\beta)^2$$
X*(y-yhat)
byapply(X, 1, `*`, y-yhat)
$\endgroup$dim
to debug your code. It's a good R skill. Your implementation isn't wrong from a stat side, you just need to sooth out your code. So probably off topic here. $\endgroup$