# Manually bootstrapping linear regression in R

I am asking you for help as I am stucked with bootstrapping...

The task is: Use the nonparametric bootstrap to compute bootstrap standard error of CAPM beta estimate based on 1000 bootstrap replications and bootstrap sample size equal to the size of the original sample.

If I understand it correctly, I am supposed to run my regression model 1000 times to estimate different estimates of the beta and its standard error. However, I am not able to put my thoughts into an actual R code.

My code:

#1)fetch data from Yahoo
#AAPL prices
apple08 <- getSymbols('AAPL', auto.assign = FALSE, from = '2008-1-1', to =
"2008-12-31")[,6]
#market proxy
rm08<-getSymbols('^ixic', auto.assign = FALSE, from = '2008-1-1', to =
"2008-12-31")[,6]

#log returns of AAPL and market
logapple08<- na.omit(ROC(apple08)*100)
logrm08<-na.omit(ROC(rm08)*100)

#OLS for beta estimation
beta_AAPL_08<-summary(lm(logapple08~logrm08))$coefficients[2,1]  OK, I've got the coefficient estimate of AAPL beta for '08. Now, I would like to run bootstrap on the beta and its standard error 1000 times with the sample size same as the original. I have this code, it works but somehow it the resample is always the same, thus the bootstraping gives actualy the same estimate of st. dev as the original regression... Any ideas? #create df from AAPL returns and market returns df08<-cbind(logapple08,logrm08) set.seed(666) Boot_times=1000 sd.boot=rep(0,Boot) for(i in 1:Boot){ # nonparametric bootstrap bootdata=df08[sample(nrow(df08), size = 251, replace = TRUE),] sd.boot[i]= coef(summary(lm(logapple08~logrm08, data = bootdata)))[2,2] }  PS: In bootstrapping regression, I am supposed to sample/bootstrap both the predictor and dependent variable, right? Or am I only supposed to resample dependent and keep predictor sample unchanged? Thank you,Adam • For np bootstrap you want to sample whole observations. – generic_user Nov 30 '17 at 18:28 • You could also just sample the residuals: add the resampled residuals to the fitted values, and then use these new values to refit the model (in this way you can keep the predictor unchanged) – matteo Nov 30 '17 at 18:44 ## 1 Answer The following are notes from my Udemy course on MCMC methods. Disregard what is not relevant to you. However, you can follow along using the mtcars data set in R to get the general idea of using Bootstrap for linear regression analysis. ## Bootstrap Bootstrap methods are a class of Monte Carlo methods known as nonparametric Monte Carlo. Bootstrap methods in simple terms are methods of resampling observed data to estimate the CDF from which the observed data is supposed to have originate from. Suppose we observe independent samples$x_1, ..., x_n$from pdf/pmf$f$, and whose CDF$F$is unobservable. Well, given that$X = (x_1, ..., x_n)^T$originates from$F$, we can use$X$to generate the empirical CDF$F_n$which is itself an estimate of$F$. $$F_n \to F \text{ as } n \to \infty$$ If we sample (with replacement) another set of$n$observations from$F_n$, we will have$X^* = (x_1^*, ..., x_n^*)^T$. This new sample$X^*$can then generate another empirical CDF,$F^*_n$which is another estimate of$F$. That is,$F^*_n$is a bootstrap estimator of$F$. We can continue this process of resampling with replacement to obtain samples$X^*_1,X^*_2, ..., X^*_B$and$F^*_{n,1}, F^*_{n,2}, ..., F^*_{n,B}$. ## Bootstrap In addition to estimating the theoretical CDF$F$, there may be a statistic of interest$\theta$(e.g. mean). We can use bootstrap methods to calculate an empirical distribution of$\theta$. From our original sample$X$we can calculate estimate$\hat{\theta}$. Similarly, using the bootstrap samples we can also calcualte estimates for$\theta$:$\hat{\theta}^*_1, ..., \hat{\theta}^*_B$. We can also calculate Bias and make confidence intervals for our estimates. ## Bootstrap Algorithm A simple bootstrap algorithm for independent samples$X = (x_1, ..., x_n)^T$is: To generate B bootstrap samples, for b in 1, ..., B do 1. Sample$x_1, ..., x_n$with replacement to create sample set$X^*_b$. Each observation$x_i$has a probability of 1/n of being in the new sample. 2. For$X^*_b$calculate$\hat{\theta}^*_b$## Bootstrap Example We will use the mtcars data set to illustrate a simple implementation. data("mtcars") mpg = mtcars$mpg
n = length(mpg)
print(mean(mpg))
hist(x = mpg, probability = TRUE, xlab = "MPG", main = "Histogram of MPG")

B = 1000 ## number of bootstraps
results = numeric(B) ## vector to hold results
for(b in 1:B){
i = sample(x = 1:n, size = n, replace = TRUE) ## sample indices
bootSample = mpg[i] ## get data
thetaHat = mean(bootSample) ## calculate the mean for bootstrap sample
results[b] = thetaHat ## store results
}

hist(x = results, probability = TRUE,
main = "Bootstrapped Samples of Mean_mpg",
xlab = "theta estimates")


## Bootstrap Example | Precaution

Before enbarking on resampling methods we must ask what variables are iid in order to determine a correct bootstrapping approach.

Bootstrap methods are not a method of generating new data for, say, a regression setting when observed samples are low.

In the above example, it is assumed that each observation in the mpg data set is indpendent and identically distributed from an unknown distribution $f$.

However, if there were to have existed some autocorrelation structure (as exist in time-series data) then we would need to adjust our resampling methodology to account for this correlation.

When dealing with time-series data, we will use a method called block bootsrap.

## Paired Bootstrapping

Let's continue to work with the mtcars data set. Say we wanted to make inferences about the linear regression parameters.

library(ggplot2, quietly = TRUE) ## for graphics
mtcars$am <- as.factor(mtcars$am) ## Transmission (0 = automatic, 1 = manual
fit = lm(formula = mpg ~ wt + am, data = mtcars)
data.frame(coefficients = coefficients(fit), CI = confint(fit), check.names = FALSE)

qplot(x = as.factor(am), y = mpg, data = mtcars, geom = "boxplot",
main = "Boxplot: MPG ~ AM", ylab = "MPG", xlab = "AM",
colour = am)

qplot(x = wt, y = mpg, data = mtcars, geom = c("point", "smooth"),
main = "Boxplot: MPG ~ Weight", ylab = "MPG", xlab = "Weight",
method = "lm", formula = y~x)

## save coefficients
beta_int = coefficients(fit)[1]
beta_wt = coefficients(fit)[2]
beta_am = coefficients(fit)[3]

n = dim(mtcars)[1] ## number of obs in data
B = 1000 ## number of bootstrap samples

results = matrix(data = NA, nrow = B, ncol = 3,
dimnames = list(NULL, c("Intercept", "wt", "am")))

## begin bootstrap for-loop
for(b in 1:B){
i = sample(x = 1:n, size = n, replace = TRUE) ## sample indices
temp = mtcars[i,] ## temp data set
temp_model =  lm(formula = mpg ~ wt + am, data = temp) ## train model
coeff = matrix(data = coefficients(temp_model), ncol = 3) ## get coefficients
results[b,] = coeff ## save coefficients in matrix
}

results <- data.frame(results, check.names = FALSE)

summary(results) ## take a look at the samples

boot_int = results[,"Intercept"]
boot_wt = results[,"wt"]
boot_am = results[,"am"]

par(mfrow = c(2,2))
hist(boot_int, main = "Bootstrapped Coefficients for Intercept",
xlab = "Coefficients for Intercept", probability = TRUE)
abline(v = coefficients(fit)[1], col = "black", lty=2)

hist(boot_wt, main = "Bootstrapped Coefficients for Weight",
xlab = "Coefficients for Weight", probability = TRUE)
abline(v = coefficients(fit)[2], col = "blue", lty=2)

hist(boot_am, main = "Bootstrapped Coefficients for AM = 1",
xlab = "Coefficients for Automatic Transmission", probability = TRUE)
abline(v = coefficients(fit)[3], col = "green", lty=2)


Now we can estimate bias for each parameter estimate. Define Bias as $Bias(\theta) = E[\theta^*] - \theta$, where in our scenario we have $Bias(\hat\theta) = E[\hat\theta^*] - \hat\theta$. Our bootstrap bias corrected estimates are then $\hat\theta_{BC} = \hat\theta - Bias(\hat\theta)$.

bias_int = mean(boot_int - beta_int)
print(bias_int)

bias_wt = mean(boot_wt - beta_wt)
print(bias_wt)

bias_am = mean(boot_am - beta_am)
print(bias_am)


Now you can incorporate our bias into the coefficients. We now have bias corrected coefficients

intercept = beta_int - bias_int
print(intercept)

wt = beta_wt - bias_wt
print(wt)

am = beta_am - bias_am
print(am)

• you didn't just copy this from Udemy? – Aksakal Nov 30 '17 at 19:34
• I am the course author – Jon Nov 30 '17 at 19:37
• It would be nice if you could add bootstrapping of residuals. In particular with non-linear models you often want to preserve the distribution of x-values when bootstrapping. – Roland Dec 1 '17 at 8:26