In this exercise, I compare a Montecarlo estimator with a bootstrap one.

Suppose that we wish to invest a fixed sum of money in two financial assets that yield returns of $$X$$ and $$Y$$, respectively, where $$X$$ and $$Y$$ are random quantities. We will invest a fraction $$\alpha$$ of our money in $$X$$, and will invest the remaining $$1-\alpha$$ in $$Y$$ . Since there is variability associated with the returns on these two assets, we wish to choose $$\alpha$$ to minimize the total risk, or variance, of our investment. In other words, we want to minimize $$Var(\alpha X + (1 − \alpha)Y )$$. One can show that the value that minimizes the risk is given by

$$$\alpha=\frac{\sigma^2_Y-\sigma_{XY}}{\sigma^2_X+\sigma^2_Y-2\sigma_{XY}}\label{eq::1}$$$

where $$\sigma^2_X = Var(X)$$, $$\sigma^2_Y = Var(Y)$$, and $$\sigma^2_{XY} = Cov(X,Y)$$. The interest is to estimate $$\alpha$$, based on the bootstrap technique.

1. Simulate $$100$$ observations with set.seed(123), such that

$X,Y\sim N_2\left[ \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 1&0.5\\ 0.5 & 1.25 \end{pmatrix} \right]$

Call this data set original. Compute $$\hat{\alpha}$$ with $$s^2_X$$, $$s^2_Y$$ and $$s_{XY}$$
library(mnormt)
mean=c(0,0)
varcov=matrix(c(1,0.5,0.5,1.25),2,2)

set.seed(123)
original=rmnorm(n=100,mean=mean,varcov=varcov)

vX=var(original[,1])
vY=var(original[,2])
covXY=cov(original[,1],original[,2])

alpha.or=(vY-covXY)/(vX+vY-2*covXY)
vX
## [1] 0.8090975
vY
## [1] 1.030543
covXY
## [1] 0.2542461
alpha.or
## [1] 0.5831784
1. Using the set.seed(Sys.time()) function, remove set.seed(123) . By using the lapply function, make a list called data.simulation and keep $$1000$$ data sets with the same previous characteristics.
set.seed(Sys.time())
n.dataset=1000
data.simulation=lapply(1:n.dataset,
function(a) rmnorm(n=100,mean=mean,varcov=varcov))
class(data.simulation[[1]])
## [1] "matrix"
1. Using both the lapply and the apply functions compute $$s^2_X$$, $$s^2_Y$$ and $$s_{XY}$$ from each element of data.simulation. Now, compute $$\hat{\alpha}$$ from each element of data.simulation. Keep them in only one data set of $$\hat{\alpha}$$.
varXY.sim=lapply(1:n.dataset,
function(a) apply(data.simulation[[a]],2,var))

## [[1]]
## [1] 1.078819 1.363634
##
## [[2]]
## [1] 1.039702 1.331717
##
## [[3]]
## [1] 0.920765 1.033606
##
## [[4]]
## [1] 1.097335 1.259797
##
## [[5]]
## [1] 1.091336 1.250979
##
## [[6]]
## [1] 1.236733 1.240761
covXY.sim=lapply(1:n.dataset,
function(a) cov(data.simulation[[a]][,1],data.simulation[[a]][,2]))

## [[1]]
## [1] 0.5511685
##
## [[2]]
## [1] 0.4647191
##
## [[3]]
## [1] 0.2827962
##
## [[4]]
## [1] 0.2983015
##
## [[5]]
## [1] 0.3985677
##
## [[6]]
## [1] 0.5235005
alpha.sim=lapply(1:n.dataset,
function(a)
(varXY.sim[[a]][2]-covXY.sim[[a]])/(varXY.sim[[a]][1]+varXY.sim[[a]][2]-2*covXY.sim[[a]]))

## [[1]]
## [1] 0.6062652
##
## [[2]]
## [1] 0.6012547
##
## [[3]]
## [1] 0.5406261
##
## [[4]]
## [1] 0.5461401
##
## [[5]]
## [1] 0.5516581
##
## [[6]]
## [1] 0.5014081
alpha=do.call('rbind',alpha.sim)
##           [,1]
## [1,] 0.6062652
## [2,] 0.6012547
## [3,] 0.5406261
## [4,] 0.5461401
## [5,] 0.5516581
## [6,] 0.5014081
1. Make a histogram of $$\hat{\alpha}$$ with a vertical red line in $$\hat{\alpha}=0.6$$ (the true value of $$\alpha$$). Compute the mean and s.d of $$\hat{\alpha}$$. This is the $$\alpha$$ Montecarlo Estimator and it standard error.
hist(alpha)
abline(v=0.6,col="red")

mean(alpha);sd(alpha)
## [1] 0.5995154
## [1] 0.07850607
1. Bootstrap function: Make a function called alpha.function that returns $$\hat{\alpha}$$ in any data set based on applying the $$\alpha$$ function to the observations indexed by the argument index. Apply your function in order to compute $$\hat{\alpha}$$ from data set using all the observations. Now, repeat this $$1000$$ times using the sample() function to randomly select $$100$$ observations from the range $$1$$ to $$100$$, with replacement. Make a histogram of $$\hat{\alpha}$$ with a vertical red line in $$\hat{\alpha}=0.6$$ (the true value of $$\alpha$$). Compute the mean and s.d of $$\hat{\alpha}$$. Compare the results with them obtained in 4.
alpha.function=function(data,index)
{
X=data[index,1]
Y=data[index,2]
return((var(Y)-cov(X,Y))/(var(X)+var(Y)-2*cov(X,Y)))
}
alpha.function(original,1:dim(original)[1])
## [1] 0.5831784
## This value must equal to alpha.or

alpha.sample=do.call('rbind',
lapply(1:1000,
function(a) alpha.function(original,sample(1:100,100,replace=TRUE))))
##           [,1]
## [1,] 0.5967613
## [2,] 0.6441368
## [3,] 0.4665356
## [4,] 0.4673642
## [5,] 0.5233725
## [6,] 0.5335732
par(mfrow=c(1,2))
hist(alpha.sample)
abline(v=0.6,col="red")
hist(alpha)
abline(v=0.6,col="red")

mean(alpha);sd(alpha) ## results from 4.)
## [1] 0.5995154
## [1] 0.07850607
mean(alpha.sample);sd(alpha.sample) ## bootstrap results
## [1] 0.5844732
## [1] 0.06534835
1. From the boot library use the boot function with $$R=1000$$ (The number of bootstrap replicates) bootstrap estimates for $$\alpha$$. Compare on the results obtained in 5.
library(boot)
boot(original,alpha.function,R=1000)
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = original, statistic = alpha.function, R = 1000)
##
##
## Bootstrap Statistics :
##      original       bias    std. error
## t1* 0.5831784 -0.001275389  0.06517734
mean(alpha.sample);sd(alpha.sample) ## results from 5.)
## [1] 0.5844732
## [1] 0.06534835