load("dataset_2.RData")
library(RiskPortfolios)
library(nlshrink)

#target return
Re <- 0.05
#estimation window
WIND <- 120
#sharinkage parameter
delta <- 0.2
#group size
group_size <- 3

#select the (excess) returns
NASDAQ100_EX_ret <- NASDAQ100_EX[,-c(1,2)]

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

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

#compute naive 1/N weights and portfolio returns
EXR_1N <- vector()
for (i in 1:length(MU)){
  w_1N <- rep(1/ncol(NASDAQ100_EX_ret),ncol(NASDAQ100_EX_ret))
  #compute out-of-sample excess returns
  EXR_1N[i] <- t(w_1N)%*%unlist(NASDAQ100_EX_ret_oos[i,])
}
#add the risk-free rate to the excess returns to get the returns
R_1N <- EXR_1N + RFree_oos
#compute mean and standard deviation
PERF_1N <- c(mean(R_1N), sd(R_1N), mean(EXR_1N)/sd(EXR_1N))
names(PERF_1N) <- c("mean","sd","SR")

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

#compute minimum variance weights and portfolio returns
V1 <- rep(1,ncol(NASDAQ100_EX_ret))
W_MINV <- list()
EXR_MINV <- vector()
for (i in 1:length(MU)){
  cov_i <- matrix(unlist(COV[i]),nrow = ncol(NASDAQ100_EX_ret), ncol = ncol(NASDAQ100_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_MINV[i] <- t(w_minv)%*%unlist(NASDAQ100_EX_ret_oos[i,])
}
#add the risk-free rate to the excess returns to get the returns
R_MINV <- EXR_MINV + RFree_oos
#compute mean and standard deviation
PERF_MINV <- c(mean(R_MINV), sd(R_MINV), mean(EXR_MINV)/sd(EXR_MINV))
names(PERF_MINV) <- c("mean","sd","SR")

#compute long-only minimum variance weights and portfolio returns
W_MINV_LONG <- list()
EXR_MINV_LONG <- vector()
for (i in 1:length(MU)){
  cov_i <- matrix(unlist(COV[i]),nrow = ncol(NASDAQ100_EX_ret), ncol = ncol(NASDAQ100_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_MINV_LONG[i] <- t(w_minv_long)%*%unlist(NASDAQ100_EX_ret_oos[i,])
}
#add the risk-free rate to the excess returns to get the returns
R_MINV_LONG <- EXR_MINV_LONG + RFree_oos
#compute mean and standard deviation
PERF_MINV_LONG <- c(mean(R_MINV_LONG), sd(R_MINV_LONG), mean(EXR_MINV_LONG)/sd(EXR_MINV_LONG))
names(PERF_MINV_LONG) <- c("mean","sd","SR")

#compute mean-variance weights and portfolio returns with target return using Ledoit & Wolf shrinkage covariance
W_MV_LW<- list()
EXR_MV_LW<- vector()
for (i in 1:length(MU)){
  cov_lw_i <- matrix(unlist(COV_LW[i]),nrow = ncol(NASDAQ100_EX_ret), ncol = ncol(NASDAQ100_EX_ret))
  mu_i <- unlist(MU[i])
  w_mv_lw<- as.vector((Re/(t(mu_i)%*%solve(cov_lw_i)%*%mu_i))%*%t(solve(cov_lw_i)%*%mu_i))
  W_MV_LW[[i]] <- w_mv_lw
  #compute out-of-sample excess returns
  EXR_MV_LW[i] <- t(w_mv_lw)%*%unlist(NASDAQ100_EX_ret_oos[i,])
}
#add the risk-free rate to the excess returns to get the returns
R_MV_LW<- EXR_MV_LW+ RFree_oos
#compute mean and standard deviation
PERF_MV_LW<- c(mean(R_MV_LW), sd(R_MV_LW), mean(EXR_MV_LW)/sd(EXR_MV_LW))
names(PERF_MV_LW) <- c("mean","sd","SR")

#compute minimum variance weights and portfolio returns using Ledoit & Wolf shrinkage covariance
V1 <- rep(1,ncol(NASDAQ100_EX_ret))
W_MINV_LW <- list()
EXR_MINV_LW <- vector()
for (i in 1:length(MU)){
  cov_lw_i <- matrix(unlist(COV_LW[i]),nrow = ncol(NASDAQ100_EX_ret), ncol = ncol(NASDAQ100_EX_ret))
  w_minv_lw <- as.vector((1/(t(V1)%*%solve(cov_lw_i)%*%V1))%*%t(solve(cov_lw_i)%*%V1))
  W_MINV_LW[[i]] <- w_minv_lw
  #compute out-of-sample excess returns
  EXR_MINV_LW[i] <- t(w_minv_lw)%*%unlist(NASDAQ100_EX_ret_oos[i,])
}
#add the risk-free rate to the excess returns to get the returns
R_MINV_LW <- EXR_MINV_LW + RFree_oos
#compute mean and standard deviation
PERF_MINV_LW <- c(mean(R_MINV_LW), sd(R_MINV_LW), mean(EXR_MINV_LW)/sd(EXR_MINV_LW))
names(PERF_MINV_LW) <- c("mean","sd","SR")

#compute long-only minimum variance weights and portfolio returns using Ledoit & Wolf shrinkage covariance
W_MINV_LW_LONG <- list()
EXR_MINV_LW_LONG <- vector()
for (i in 1:length(MU)){
  cov_lw_i <- matrix(unlist(COV_LW[i]),nrow = ncol(NASDAQ100_EX_ret), ncol = ncol(NASDAQ100_EX_ret))
  w_minv_lw_long <- optimalPortfolio(Sigma=cov_lw_i,mu=NULL,semiDev=NULL,control=list(type='minvol',constraint='lo'))
  W_MINV_LW_LONG[[i]] <- w_minv_lw_long
  #compute out-of-sample excess returns
  EXR_MINV_LW_LONG[i] <- t(w_minv_lw_long)%*%unlist(NASDAQ100_EX_ret_oos[i,])
}
#add the risk-free rate to the excess returns to get the returns
R_MINV_LW_LONG <- EXR_MINV_LW_LONG + RFree_oos
#compute mean and standard deviation
PERF_MINV_LW_LONG <- c(mean(R_MINV_LW_LONG), sd(R_MINV_LW_LONG), mean(EXR_MINV_LW_LONG)/sd(EXR_MINV_LW_LONG))
names(PERF_MINV_LW_LONG) <- c("mean","sd","SR")

#combine mean-variance and 1/N portfolios
W_MV_1N <- list()
EXR_MV_1N <- vector()
for (i in 1:length(MU)){
  w_mv_lw<- unlist(W_MV_LW[i])
  w_mv_1N <- delta*w_mv_lw+ (1-delta)*w_1N
  W_MV_1N[[i]] <- w_mv_1N
  #compute out-of-sample excess returns
  EXR_MV_1N[i] <- t(w_mv_1N)%*%unlist(NASDAQ100_EX_ret_oos[i,])
}
#add the risk-free rate to the excess returns to get the returns
R_MV_1N <- EXR_MV_1N + RFree_oos
#compute mean and standard deviation
PERF_MV_1N <- c(mean(R_MV_1N), sd(R_MV_1N), mean(EXR_MV_1N)/sd(EXR_MV_1N))
names(PERF_MV_1N) <- c("mean","sd","SR")

#compute mean-variance weights and portfolio returns with target return using the grouping strategy
index_matrix <- matrix(1:ncol(NASDAQ100_EX_ret), ncol = group_size, byrow = TRUE)

#create a dataset of stocks grouped by 3
NASDAQ100_EX_ret_grouped_list <- list()
for (i in 1:nrow(index_matrix)) {
  NASDAQ100_EX_ret_grouped_list[[i]] <- rowMeans(NASDAQ100_EX_ret[,index_matrix[i,]])
}
NASDAQ100_EX_ret_grouped <- do.call(cbind,NASDAQ100_EX_ret_grouped_list)

#estimate mean and covariance of the grouped dataset using a rolling window
MU_GROUP <- list()
COV_LW_GROUP <- list()
for (i in 1:(nrow(NASDAQ100_EX_ret_grouped)-WIND)) {
  MU_GROUP[[i]] <- colMeans(NASDAQ100_EX_ret_grouped[i:(WIND+i-1),])
  COV_LW_GROUP[[i]] <- linshrink_cov(as.matrix(NASDAQ100_EX_ret_grouped[i:(WIND+i-1),]))
}

#compute mean-variance weights and portfolio returns with target return using Ledoit & Wolf shrinkage covariance
W_MV_GROUP <- list()
EXR_MV_GROUP<- vector()
for (i in 1:length(MU_GROUP)){
  cov_lw_group_i <- matrix(unlist(COV_LW_GROUP[i]),nrow = ncol(NASDAQ100_EX_ret_grouped), ncol = ncol(NASDAQ100_EX_ret_grouped))
  mu_group_i <- unlist(MU_GROUP[i])
  w_mv_group_1 <- as.vector((Re/(t(mu_group_i)%*%solve(cov_lw_group_i)%*%mu_group_i))%*%t(solve(cov_lw_group_i)%*%mu_group_i))
  w_mv_group <- rep(w_mv_group_1, each = 3)
  W_MV_GROUP[[i]] <- w_mv_group
  #compute out-of-sample excess returns
  EXR_MV_GROUP[i] <- t(w_mv_group)%*%unlist(NASDAQ100_EX_ret_oos[i,])
}
#add the risk-free rate to the excess returns to get the returns
R_MV_GROUP<- EXR_MV_GROUP+ RFree_oos
#compute mean and standard deviation
PERF_MV_GROUP<- c(mean(R_MV_GROUP), sd(R_MV_GROUP), mean(EXR_MV_GROUP)/sd(EXR_MV_GROUP))
names(PERF_MV_GROUP) <- c("mean","sd","SR")

#print results
print(round(PERF_1N,4))
print(round(PERF_MV,4))
print(round(PERF_MV_LW,4))
print(round(PERF_MINV,4))
print(round(PERF_MINV_LW,4))
print(round(PERF_MINV_LONG,4))
print(round(PERF_MINV_LW_LONG,4))
print(round(PERF_MV_1N,4))
print(round(PERF_MV_GROUP,4))