load("dataset.RData")
library(RiskPortfolios)

#risk-aversion parameter
gam_1 <- 2
gam_2 <- 5
#target return
Re <- 0.05
#estimation window
WIND <- 120

#select the (excess) returns
DJIA_ex_ret <- DJIA_ex[,-c(1,2)]

#asset returns to be used to compute out-of-sample portfolio returns
DJIA_ex_ret_oos <- DJIA_ex_ret[-(1:WIND),]
#risk free rate to be used to compute out-of-sample portfolio returns
RFree_oos <- DJIA_ex[-(1:WIND),2]

#estimate mean and covariance using a rolling window
MU <- list()
COV <- list()
for (i in 1:(nrow(DJIA_ex_ret)-WIND)) {
  MU[[i]] <- colMeans(DJIA_ex_ret[i:(WIND+i-1),])
  COV[[i]] <- cov(DJIA_ex_ret[i:(WIND+i-1),])
}

#compute optimal weights and portfolio returns for both values of gamma
W_MV_GAM_1 <- list()
W_MV_GAM_2 <- list()
EXR_MV_GAM_1 <- vector()
EXR_MV_GAM_2 <- vector()
for (i in 1:length(MU)){
  cov_i <- matrix(unlist(COV[i]),nrow = ncol(DJIA_ex_ret), ncol = ncol(DJIA_ex_ret))
  mu_i <- unlist(MU[i])
  #compute optimal weights with the first gamma value
  w_mv_gam_1 <- as.vector((1/gam_1)*(solve(cov_i))%*%mu_i)
  W_MV_GAM_1[[i]] <- w_mv_gam_1
  #compute optimal weights with the second gamma value
  w_mv_gam_2 <- as.vector((1/gam_2)*(solve(cov_i))%*%mu_i)
  W_MV_GAM_2[[i]] <- w_mv_gam_2
  #compute out-of-sample excess returns
  EXR_MV_GAM_1[i] <- t(w_mv_gam_1)%*%unlist(DJIA_ex_ret_oos[i,])
  EXR_MV_GAM_2[i] <- t(w_mv_gam_2)%*%unlist(DJIA_ex_ret_oos[i,])
}
#add the risk-free rate to the excess returns to get the returns
R_MV_GAM_1 <- EXR_MV_GAM_1 + RFree_oos
R_MV_GAM_2 <- EXR_MV_GAM_2 + RFree_oos
#compute mean, standard deviation, and utility
PERF_MV_GAM_1 <- c(mean(R_MV_GAM_1), sd(R_MV_GAM_1), mean(R_MV_GAM_1) - (gam_1/2)*var(R_MV_GAM_1))
names(PERF_MV_GAM_1) <- c("mean","sd","utility")
PERF_MV_GAM_2 <- c(mean(R_MV_GAM_2), sd(R_MV_GAM_2), mean(R_MV_GAM_2) - (gam_2/2)*var(R_MV_GAM_2))
names(PERF_MV_GAM_2) <- c("mean","sd","utility")

#compute 1/N weights and portfolio returns for both values of gamma
V1 <- rep(1,ncol(DJIA_ex_ret))
W_1N_GAM_1 <- list()
W_1N_GAM_2 <- list()
EXR_1N_GAM_1 <- vector()
EXR_1N_GAM_2 <- vector()
for (i in 1:length(MU)){
  cov_i <- matrix(unlist(COV[i]),nrow = ncol(DJIA_ex_ret), ncol = ncol(DJIA_ex_ret))
  mu_i <- unlist(MU[i])
  #compute optimal weights with the first gamma value
  w_1n_gam_1 <- as.vector((1/gam_1)*((t(V1)%*%mu_i)/(t(V1)%*%cov_i%*%V1))%*%V1)
  W_1N_GAM_1[[i]] <- w_1n_gam_1
  #compute optimal weights with the second gamma value
  w_1n_gam_2 <- as.vector((1/gam_2)*((t(V1)%*%mu_i)/(t(V1)%*%cov_i%*%V1))%*%V1)
  W_1N_GAM_2[[i]] <- w_1n_gam_2
  #compute out-of-sample excess returns
  EXR_1N_GAM_1[i] <- t(w_1n_gam_1)%*%unlist(DJIA_ex_ret_oos[i,])
  EXR_1N_GAM_2[i] <- t(w_1n_gam_2)%*%unlist(DJIA_ex_ret_oos[i,])
}
#add the risk-free rate to the excess returns to get the returns
R_1N_GAM_1 <- EXR_1N_GAM_1 + RFree_oos
R_1N_GAM_2 <- EXR_1N_GAM_2 + RFree_oos
#compute mean, standard deviation, and utility
PERF_1N_GAM_1 <- c(mean(R_1N_GAM_1), sd(R_1N_GAM_1), mean(R_1N_GAM_1) - (gam_1/2)*var(R_1N_GAM_1))
names(PERF_1N_GAM_1) <- c("mean","sd","utility")
PERF_1N_GAM_2 <- c(mean(R_1N_GAM_2), sd(R_1N_GAM_2), mean(R_1N_GAM_2) - (gam_2/2)*var(R_1N_GAM_2))
names(PERF_1N_GAM_2) <- c("mean","sd","utility")

#compute mean-variance weights and portfolio returns with target return
W_MV_RE <- list()
EXR_MV_RE <- vector()
for (i in 1:length(MU)){
  cov_i <- matrix(unlist(COV[i]),nrow = ncol(DJIA_ex_ret), ncol = ncol(DJIA_ex_ret))
  mu_i <- unlist(MU[i])
  w_mv_re <- as.vector((Re/(t(mu_i)%*%solve(cov_i)%*%mu_i))%*%t(solve(cov_i)%*%mu_i))
  W_MV_RE[[i]] <- w_mv_re
  #compute out-of-sample excess returns
  EXR_MV_RE[i] <- t(w_mv_re)%*%unlist(DJIA_ex_ret_oos[i,])
}
#add the risk-free rate to the excess returns to get the returns
R_MV_RE <- EXR_MV_RE + RFree_oos
#compute mean and standard deviation
PERF_MV_RE <- c(mean(R_MV_RE), sd(R_MV_RE), mean(EXR_MV_RE)/sd(EXR_MV_RE))
names(PERF_MV_RE) <- c("mean","sd","SR")

#compute minimum variance weights and portfolio returns
W_MINV <- list()
EXR_MINW <- vector()
for (i in 1:length(MU)){
  cov_i <- matrix(unlist(COV[i]),nrow = ncol(DJIA_ex_ret), ncol = ncol(DJIA_ex_ret))
  w_minv <- as.vector((1/(t(V1)%*%solve(cov_i)%*%V1))%*%t(solve(cov_i)%*%V1))
  W_MINV[[i]] <- w_minv
  #compute out-of-sample excess returns
  EXR_MINW[i] <- t(w_minv)%*%unlist(DJIA_ex_ret_oos[i,])
}
#add the risk-free rate to the excess returns to get the returns
R_MINW <- EXR_MINW + RFree_oos
#compute mean and standard deviation
PERF_MINW <- c(mean(R_MINW), sd(R_MINW), mean(EXR_MINW)/sd(EXR_MINW))
names(PERF_MINW) <- c("mean","sd","SR")

#compute long-only minimum variance weights and portfolio returns
W_MINV_LONG <- list()
EXR_MINW_LONG <- vector()
for (i in 1:length(MU)){
  cov_i <- matrix(unlist(COV[i]),nrow = ncol(DJIA_ex_ret), ncol = ncol(DJIA_ex_ret))
  w_minv_long <- optimalPortfolio(Sigma=cov_i,mu=NULL,semiDev=NULL,control=list(type='minvol',constraint='lo'))
  W_MINV_LONG[[i]] <- w_minv_long
  #compute out-of-sample excess returns
  EXR_MINW_LONG[i] <- t(w_minv_long)%*%unlist(DJIA_ex_ret_oos[i,])
}
#add the risk-free rate to the excess returns to get the returns
R_MINW_LONG <- EXR_MINW_LONG + RFree_oos
#compute mean and standard deviation
PERF_MINW_LONG <- c(mean(R_MINW_LONG), sd(R_MINW_LONG), mean(EXR_MINW_LONG)/sd(EXR_MINW_LONG))
names(PERF_MINW_LONG) <- c("mean","sd","SR")

#print results
print(round(PERF_MV_GAM_1,4))
print(round(PERF_MV_GAM_2,4))
print(round(PERF_1N_GAM_1,4))
print(round(PERF_1N_GAM_2,4))
print(round(PERF_MV_RE,4))
print(round(PERF_MINW,4))
print(round(PERF_MINW_LONG,4))