install.packages("ISLR")

library("ISLR")

data(Portfolio)

ts.plot(Portfolio,col=1:2)

plot(Portfolio)

alpha = function (data ,index){
  X=data$X [index]
  Y=data$Y [index]
  return ((var(Y)-cov (X,Y))/(var(X)+var(Y) -2* cov(X,Y)))
}

# Estimating alpha and correlation of X and Y
a = alpha(Portfolio,1:100)
r = cor(Portfolio)[1,2]
c(a,r)

set.seed(1)
B = 1000
a.bootstrap = rep(0,B)
r.bootstrap = rep(0,B)
for (i in 1:B){
  ind = sample(1:100,100,replace=TRUE)
  Portfolio.b = Portfolio[ind,]
  a.bootstrap[i] = alpha(Portfolio.b,1:100)
  r.bootstrap[i] = cor(Portfolio.b)[1,2]
}

par(mfrow=c(1,2))
hist(a.bootstrap)
abline(v=a,col=2)
hist(r.bootstrap)
abline(v=r,col=2)


# Regressao linear simples
y = Portfolio[,1]
x = Portfolio[,2]

reg = lm(y~x)
bhat = reg$coef
summary(reg)
X = cbind(1,x)
iXtX = solve(t(X)%*%X)
Xty = t(X)%*%y
bhat = iXtX%*%Xty
s    = sqrt(sum((y-X%*%bhat)^2)/(100-2))
s.bhat = s*sqrt(diag(iXtX))

set.seed(1)
B = 5000
b.bootstrap = matrix(0,B,2)
for (i in 1:B){
  ind = sample(1:100,100,replace=TRUE)
  reg = lm(y[ind]~x[ind])
  b.bootstrap[i,] = reg$coef
}

par(mfrow=c(1,2))
hist(b.bootstrap[,1],prob=TRUE)
abline(v=bhat[1],col=2)
xxx = seq(-0.4,0.4,length=100)
lines(xxx,dnorm(xxx,bhat[1],s.bhat[1]),col=4,lwd=2)

hist(b.bootstrap[,2],prob=TRUE)
abline(v=bhat[2],col=2)
xxx = seq(0.2,0.8,length=100)
lines(xxx,dnorm(xxx,bhat[2],s.bhat[2]),col=4,lwd=2)



# Parametric bootstrap
set.seed(1)
B = 5000
b1.bootstrap = matrix(0,B,2)
for (i in 1:B){
  ind = sample(1:100,100,replace=TRUE)
  X1 = X[ind,]
  y1 = rnorm(100,X1%*%bhat,s)
  b1.bootstrap[i,] = lm(y1~X1-1)$coef
}

par(mfrow=c(1,2))
hist(b1.bootstrap[,1],prob=TRUE)
abline(v=bhat[1],col=2)
xxx = seq(-0.4,0.4,length=100)
lines(xxx,dnorm(xxx,bhat[1],s.bhat[1]),col=4,lwd=2)

hist(b1.bootstrap[,2],prob=TRUE)
abline(v=bhat[2],col=2)
xxx = seq(0.2,0.8,length=100)
lines(xxx,dnorm(xxx,bhat[2],s.bhat[2]),col=4,lwd=2)