When N is small, the theory of portfolio optimization can be applied using sample means and the sample covariance matrix However, the expectation of estimation error, especially with larger values of N, can result in portfolios that only appear e±cient.
Assume that we want to find the tangency portfolio that maximizes Sharpe's ratio. However, we should suspect that 0.3681 must be an overestimate since this portfolio only maximizes Sharpe's ratio using estimated parameters, not the true means and covariance matrix. To evaluate the possible amount of overestimation, one can use the bootstrap. In bootstrap simulation , the sample is the \true population" so that the sample mean and covariance matrix are the \true parameters,“ and the resamples mimic the sampling process. Actual Sharpe's ratios are calculated with the sample means and covariance matrix, while estimated Sharpe's ratio use the means and covariance matrix of the resamples.First, 250 resamples were taken and for each the tangency portfolio was estimated. Resampling was done by sampling rows of the data matrix . For each of the 250 tangency portfolios estimated from the resamples, the actual and estimated Sharpe's ratios were calculated. Box-plots of the 250 actual and 250 estimated Sharpe's ratios One can see that all 250 estimated tangency portfolios have actual Sharpe's ratios below this value, as they must since the actual Sharpe's ratio is maximized by the true tangency portfolio, not the estimated tangency portfolios.
opts_chunk$set(comment = "test", results = "hide", warning = FALSE, error = FALSE,
message = FALSE)
library(xts)
library(MASS)
library(timeDate)
library(timeSeries)
library(fGarch)
library(zoo)
library("skewt")
library(Ecdat)
library(quadprog)
setwd("~/Desktop/Courses_backup/Stat_Method_Finance/Final/CSV_data")
data <- read.csv("stocks.csv", header = T)
date <- as.Date(data[, 1], format = "%m/%d/%Y")
zoodata <- zoo(data[, -1], date) # zoo price data for time-series plotting
xtsdata <- as.xts(zoodata) #xts price data for time-series plotting
data_aj <- apply(data[, -1], 2, function(x) diff(log(x))) # log return
date_aj <- date[-1] #cut off first date
zoo.r_aj <- zoo(data_aj, date_aj) #zoo log return data for time-series plotting
xts.r_aj <- as.xts(data_aj, date_aj) #xts log return data for time-series plotting
data.r_aj <- as.data.frame(xts.r_aj) # data.frame log return data for analysis
d <- data.r_aj
### PORTOFOLIO THEORY & ASSET ALLOCATION
#### MINIMUM VARIANCE PORTOFOLIO WITHOUT SHORT SALES
#shrinkage estimation
data<-data[,2:7]
n<-nrow(data)
return0<-apply(data,2,function(x) diff(x)/x[-n])
mean_vect_orig = apply(return0,2,mean)
cov_mat_orig = cov(return0)
sd_vect_orig = sqrt(diag(cov_mat_orig))
cm<-diag(1,6,6)
Amat = cbind(rep(1,6),mean_vect_orig,cm) # set the constraints matrix
muP = seq(min(mean_vect_orig),max(mean_vect_orig)-1e-4,length=300) # set of 300 possible target values
# for the expect portfolio return
sdP = muP # set up storage for standard deviations of portfolio returns
weights = matrix(0,nrow=300,ncol=6) # storage for portfolio weights
for (i in 1:length(muP)) # find the optimal portfolios for each target expected return
{
bvec = c(1,muP[i],rep(0,6)) # constraint vector
result =
solve.QP(Dmat=2*cov_mat_orig,dvec=rep(0,6),Amat=Amat,bvec=bvec,meq=2)
sdP[i] = sqrt(result$value)
weights[i,] = result$solution
}
mufree = 0.01/365 # input value of risk-free interest rate
ratio =( muP-mufree)/sdP # compute Sharpe ratios
ind = (ratio == max(ratio))
sharpeorig<-max(ratio)
w<-weights[ind,]
a<-0
sr_estimate<-list()
sr_actual<-list()
for(k in 1:8){
a=a+.1
sharpestimate<-c()
sharpeactual<-c()
for(j in 1:250){
returns<-apply(return0,2,function(x) sample(x,replace=T))
n<-nrow(returns)
mean_vect = apply(returns,2,mean)
mean_vect_a=a*mean_vect+(1-a)*mean(mean_vect)
cov_mat_a = cov(returns)
sd_vect_a = sqrt(diag(cov_mat_a))
cm<-diag(1,6,6)
Amat = cbind(rep(1,6),mean_vect_a,cm) # set the constraints matrix
if(min(mean_vect_a)<0) {
m=.9}else{
m=1.1
}
if(max(mean_vect_a)<0) {
mm=1.1}else{
mm=.9
}
muP = seq(min(mean_vect_a)*m,max(mean_vect_a)*mm,length=300) # set of 300 possible target values
# for the expect portfolio return
sdP = muP # set up storage for standard deviations of portfolio returns
weights = matrix(0,nrow=300,ncol=6) # storage for portfolio weights
for (i in 1:length(muP)) # find the optimal portfolios for each target expected return
{
bvec = c(1,muP[i],rep(0,6)) # constraint vector
result =
solve.QP(Dmat=2*cov_mat_a,dvec=rep(0,6),Amat=Amat,bvec=bvec,meq=2)
sdP[i] = sqrt(result$value)
weights[i,] = result$solution
}
mufree = 0.01/365 # input value of risk-free interest rate
ratio =( muP-mufree)/sdP # compute Sharpe ratios
ind = (ratio == max(ratio))
sharpestimate[j]<-ratio[ind]
w<-weights[ind,]
mean_vect_oa=a*mean_vect_orig+(1-a)*mean(mean_vect_orig)
sharpeactual[j]<-as.numeric((w%*%mean_vect_oa-mufree)/sqrt(t(w)%*%cov_mat_orig%*%w))
}
sr_estimate[[k]]<- sharpestimate
sr_actual[[k]]<- sharpeactual
}
boxplot(sr_actual[[1]], sr_estimate[[1]], main = "Bootstrapping Comparation",
sub = "under shrinkage adjustment", ylim = c(-0.04, 0.12))
axis(side = 1, at = c(1, 2), labels = c("actual", "estimate"))
abline(h = sharpeorig, lwd = 3, col = 2)
boxplot(sr_actual[[2]], sr_estimate[[2]], main = "Bootstrapping Comparation",
sub = "under shrinkage adjustment", ylim = c(-0.04, 0.12))
axis(side = 1, at = c(1, 2), labels = c("actual", "estimate"))
abline(h = sharpeorig, lwd = 3, col = 2)
boxplot(sr_actual[[3]], sr_estimate[[3]], main = "Bootstrapping Comparation",
sub = "under shrinkage adjustment", ylim = c(-0.04, 0.12))
axis(side = 1, at = c(1, 2), labels = c("actual", "estimate"))
abline(h = sharpeorig, lwd = 3, col = 2)
boxplot(sr_actual[[4]], sr_estimate[[4]], main = "Bootstrapping Comparation",
sub = "under shrinkage adjustment", ylim = c(-0.04, 0.14))
axis(side = 1, at = c(1, 2), labels = c("actual", "estimate"))
abline(h = sharpeorig, lwd = 3, col = 2)
boxplot(sr_actual[[5]], sr_estimate[[5]], main = "Bootstrapping Comparation",
sub = "under shrinkage adjustment", ylim = c(-0.04, 0.14))
axis(side = 1, at = c(1, 2), labels = c("actual", "estimate"))
abline(h = sharpeorig, lwd = 3, col = 2)
boxplot(sr_actual[[6]], sr_estimate[[6]], main = "Bootstrapping Comparation",
sub = "under shrinkage adjustment", ylim = c(-0.04, 0.14))
axis(side = 1, at = c(1, 2), labels = c("actual", "estimate"))
abline(h = sharpeorig, lwd = 3, col = 2)
boxplot(sr_actual[[7]], sr_estimate[[7]], ann = F, main = "Bootstrapping Comparation",
sub = "under shrinkage adjustment", ylim = c(-0.04, 0.14))
axis(side = 1, at = c(1, 2), labels = c("actual", "estimate"))
abline(h = sharpeorig, lwd = 3, col = 2)
boxplot(sr_actual[[8]], sr_estimate[[8]], main = "Bootstrapping Comparation",
sub = "under shrinkage adjustment", ylim = c(-0.04, 0.16))
axis(side = 1, at = c(1, 2), labels = c("actual", "estimate"))
abline(h = sharpeorig, lwd = 3, col = 2)
conclusion: when shrinkage parameter alpha decreases to zero, estimation turns to be less volatile but more bias. When choosing alpha, we should consider trade off between variance and bias to get a best estimation.