load("dataset_2.RData")
load("FF3.RData")
library(PerformanceAnalytics)
library(nlshrink)
library(lmtest)

#estimation window
WIND <- 120
#benchmark for the semivariance and Sortino ratio. When working with excess returns, 0=risk-free rate
bench <- 0
#confidence level used in VaR calculations
percentile <- 0.95
#transaction costs
tcost <- 0.001

#function that evaluates the results (mean,sd,dd,sk,exKurt,sharpe,sortino,AdjSharpe,ModAdjSharpe,turnover)
Eval <- function(exReturns,Returns,Returns_oos,benchmark,confidence,weights,MktReturns){
  exr_mean <- mean(exReturns)
  names(exr_mean) <- "ExR_mean"
  r_mean <- mean(Returns)
  names(r_mean) <- "R_mean"
  exr_sd <- sd(exReturns)
  names(exr_sd) <- "R_sd"
  r_sd <- sd(Returns)
  names(r_sd) <- "R_sd"
  r_dd <- DownsideDeviation(exReturns, MAR=benchmark, method="full")
  names(r_dd) <- "R_DownDev"
  r_sk <- skewness(Returns)
  names(r_sk) <- "R_skew"
  r_exKurt <- PerformanceAnalytics::kurtosis(Returns, method="excess")
  names(r_exKurt) <- "R_exKurt"
  sharpe <- exr_mean/exr_sd
  names(sharpe) <- "SharpeR"
  sortino <- SortinoRatio(exReturns, MAR=benchmark)
  names(sortino) <- "SortinoR"
  CVaR <- ETL(exReturns,p=confidence,method="historical")
  names(CVaR) <- "CVaR"
  alpha_reg <- lm(Returns ~ MktReturns)
  alpha <- alpha_reg$coefficients[1]
  names(alpha) <- "alpha"
  alpha_pr <- lmtest::coeftest(alpha_reg,vcov=sandwich::NeweyWest(alpha_reg))[1,4]
  names(alpha_pr) <- "alpha p-value"
  to <- vector()
  for(i in 2:length(weights)){
    weights_prev <- (unlist(weights[i-1])+unlist(weights[i-1])*Returns_oos[i-1,])/sum(unlist(weights[i-1])+unlist(weights[i-1])*Returns_oos[i-1,])
    w_diff <- abs(unlist(weights[i]) - weights_prev)
    to[i] <- sum(w_diff)
  }
  to <- mean(to[-1])
  names(to) <- "Turnover"
  results <- c(r_mean,r_sd,r_dd,r_sk,r_exKurt,sharpe,sortino,CVaR,alpha,alpha_pr,to)
  return(results)
}

#function that computes transaction costs at each out-of-sample period
TransactionCosts <- function(weights,tCost,Returns_oos){
  to_full <- vector()
  for(i in 2:nrow(weights)){
    weights_prev <- (weights[i-1,]+weights[i-1,]*Returns_oos[i-1,])/sum(weights[i-1,]+weights[i-1,]*Returns_oos[i-1,])
    w_Diff <- abs(weights[i,] - weights_prev)
    to_full[i] <- sum(w_Diff)
  }
  to_full <- to_full[-1]
  tCosts <- c(0,to_full*tCost)
  return(tCosts)
}

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

#asset excess returns to be used to compute out-of-sample portfolio excess 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]
#asset returns to be used to compute the portfolio turnover
NASDAQ100_ret_oos <- NASDAQ100_ret[-(1:WIND),]
#out-of-sample market excess returns used to compute the alpha
MKT_EX_oos <- FF3$MKT[-(1:WIND)] - RFree_oos

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

#compute naive 1/N weights and portfolio returns
W_1N <- list()
EXR_1N <- vector()
for (i in 1:length(COV)){
  w_1N <- rep(1/ncol(NASDAQ100_EX_ret),ncol(NASDAQ100_EX_ret))
  W_1N[[i]] <- w_1N
  #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
#evaluate results
PERF_1N <- Eval(exReturns=EXR_1N,Returns=R_1N,Returns_oos=NASDAQ100_ret_oos,benchmark=bench,
                confidence=percentile,weights=W_1N,MktReturns=MKT_EX_oos)
#compute wealth dynamics
wealth_1N <- c(1, cumprod(1 + R_1N))

#results net of transaction costs
TCOSTS_1N <- TransactionCosts(weights=do.call(rbind,W_1N),tCost=tcost,Returns_oos=NASDAQ100_ret_oos)
R_1N_NET <- R_1N - TCOSTS_1N
EXR_1N_NET <- R_1N_NET - RFree_oos
PERF_1N_NET <- Eval(exReturns=EXR_1N_NET,Returns=R_1N_NET,Returns_oos=NASDAQ100_ret_oos,benchmark=bench,
                    confidence=percentile,weights=W_1N,MktReturns=MKT_EX_oos)
#compute wealth dynamics net of transaction costs
wealth_1N_NET <- c(1, cumprod(1 + R_1N_NET))

#compute minimum variance weights and portfolio returns using Ledoit & Wolf shrinkage covariance
V1 <- rep(1,ncol(NASDAQ100_EX_ret))
W_MINV <- list()
EXR_MINV <- vector()
for (i in 1:length(COV)){
  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
#evaluate results
PERF_MINV <- Eval(exReturns=EXR_MINV,Returns=R_MINV,Returns_oos=NASDAQ100_ret_oos,benchmark=bench,
                  confidence=percentile,weights=W_MINV,MktReturns=MKT_EX_oos)
#compute wealth dynamics
wealth_MINV <- c(1, cumprod(1 + R_MINV))

#results net of transaction costs
TCOSTS_MINV <- TransactionCosts(weights=do.call(rbind,W_MINV),tCost=tcost,Returns_oos=NASDAQ100_ret_oos)
R_MINV_NET <- R_MINV - TCOSTS_MINV
EXR_MINV_NET <- R_MINV_NET - RFree_oos
PERF_MINV_NET <- Eval(exReturns=EXR_MINV_NET,Returns=R_MINV_NET,Returns_oos=NASDAQ100_ret_oos,benchmark=bench,
                      confidence=percentile,weights=W_MINV,MktReturns=MKT_EX_oos)
#compute wealth dynamics net of transaction costs
wealth_MINV_NET <- c(1, cumprod(1 + R_MINV_NET))

#print results
print(round(PERF_1N,4))
print(round(PERF_1N_NET,4))
print(round(PERF_MINV,4))
print(round(PERF_MINV_NET,4))

#plot wealth
DATE <- seq(as.Date("2009-01-01"),length=nrow(NASDAQ100_EX_ret_oos)+1,by="months")-1
plot(wealth_1N~DATE, type="l", lwd=2, col="red", ylab="Wealth", main="1/N portfolio")
lines(wealth_1N_NET~DATE, lwd=2, col="orange")
plot(wealth_MINV~DATE, type="l", lwd=2, col="blue", ylab="Wealth", main="Minimum variance portfolio")
lines(wealth_MINV_NET~DATE, lwd=2, col="green3")