> 
> n <- length(Y)
> 
> toluca.reg1 <- lm(Y ~ X)
> summary(toluca.reg1)

Call:
lm(formula = Y ~ X)

Residuals:
    Min      1Q  Median      3Q     Max 
-83.876 -34.088  -5.982  38.826 103.528 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   62.366     26.177   2.382   0.0259 *  
X              3.570      0.347  10.290 4.45e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 48.82 on 23 degrees of freedom
Multiple R-squared: 0.8215,     Adjusted R-squared: 0.8138 
F-statistic: 105.9 on 1 and 23 DF,  p-value: 4.449e-10 

> Yhat <- as.matrix(predict(toluca.reg1))
> e <- as.matrix(residuals(toluca.reg1))
> b1 <- coef(toluca.reg1)[2]
> 
> 
> 
> X <- as.matrix(cbind(rep(1,n),X))
> 
> 
> num.boot <- 1000
> b1.boot <- rep(0,num.boot)
> 
> for (i in 1:num.boot) {
+ e.boot <- as.matrix(sample(e,size=n,replace=TRUE))
+ Y.boot <- Yhat + e.boot
+ 
+ b.boot <- solve(t(X) %*% X) %*% t(X) %*% Y.boot
+ b1.boot[i] <- b.boot[2,1]
+ }
> 
> hist(b1.boot,breaks=seq(2.40,4.80,0.12))
> 
> (b1.boot_025 <- quantile(b1.boot,.025))
    2.5% 
2.938487 
> (b1.boot_975 <- quantile(b1.boot,.975))
   97.5% 
4.189233 
> 
> (b1.boot.sd <- sd(b1.boot))
[1] 0.3272935
> 
> (d1 <- b1-b1.boot_025)
       X 
0.631715 
> (d2 <- b1.boot_975-b1)
    97.5% 
0.6190311 
> 
> (beta1.95CI <- c(b1-d2,b1+d1))
       X        X 
2.951171 4.201917 
> 
>