How to implement bootstrap efficiently in a matrix-based language?

Question

I am working on a Monte Carlo study of bootstrapping in an AR(1) model for a homework assignment (I'm using Matlab). The goal is to say something about the empirical rejection probabilities of the bootstrap in this specific context. However, I've been running into some computational problems.

Some context first. Suppose one has a regression $y_t = \gamma + \rho y_{t-1} + \varepsilon_t$, and wants to test the two-tailed hypothesis $\rho = \rho_0$ by the bootstrap. To do this, one first estimates the regression by, say, OLS, and then uses the obtained residuals $\hat \varepsilon_t$ and coefficients $\hat \gamma, \; \hat \rho$ for constructing the bootstrap samples. However, in this case, the bootstrap sample has to be constructed recursively. What is even worse for the computations, I have to repeat these calculations many times using different simulated data sets and different values of $\rho$. (The reason is that this allows one to get some sense of the empirical level of the test at different values of $\rho$.)

So far I have been doing it in the naive way of just estimating many regressions, generating the bootstrap samples, and so on. I was able to eliminate quite a few loops by using the filter function when generating bootstrap samples but still computational times are quite huge. I compared it with the built-in bootstrap in gretl, and it is many times slower. (I'm aware that gretl is written in C, but I sense that there is more to this than that.)

Hence, my question would be: what general advice could you give to implement bootstrapping efficiently in a matrix based language? Are there some smart ways to avoid estimation of many OLS regressions and generation of the bootstrap samples? Although I'm using Matlab for the computations, but if you have some general advice on say, doing this with R, all the help would be very much appreciated. Also, since this is a homework assignment I have to do, using existing libraries is not really an option. :)

Answer 1

As a daily Matlab user, I suggest eliminating all for-loops. Design each function, from the start to work on a bunch of problems in parallel. For example, if you are going to bootstrap a length $n$ vector $k$ times, you would be better off writing code to return a $n \times k$ matrix, rather than writing a vector bootstrap which you then call $k$ times in a loop. Here is the simplest bootstrap function of this form:

function [xboot] = simple_boot(x,k);
%simple bootstrap function
if (nargin < 2)
  k = 1;  % default to single bootstrap
end    
xboot = x(ceil(numel(x) * rand(numel(x),k)));
end  % function

Similarly, if you are going to perform the lagged variable regression multiple times on different vectors, it would be better to assemble them in a matrix, column by column, and solve them all at once. This would be fairly difficult in a general setting, but for your regression, it should be pretty easy:

% let Y be the matrix whose columns are samples of y.
f_demean = @(x)(bsxfun(@minus,x,mean(x,1)));   % function which subtracts the mean.
rho_hat = f_demean(Y(2:end,:)) ./ f_demean(Y(1:end-1,:));  % gives a row of rho
% if you needed estimates of gamma, it would be better not to discard
% the sample means, as I do above. Also, the two mean computations
% can be consolidated if efficiency is a premium

This just gives a taste of what kind of funny tricks must be done to eliminate for-loops. Without know more specifics of your problem, I cannot do much more.