# Asset Allocation with CER 
# elaboration on the original code produced by E.Zivot by C. Favero 
# author: Carlo Favero 
# created: August, 2023
#
# comments: Original Examples are taken from  chapter 11 in Zivot and Wang (2006)

rm(list=ls()) #Removes all items in Environment!
#setwd(path)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#install.packages("fEcofin", repos="http://R-Forge.R-project.org")
library(fEcofin) 
# load required packages
listofpackages <- c("ellipse","dygraphs","ggplot2","reshape2")

for (j in listofpackages){
  if(sum(installed.packages()[, 1] == j) == 0) {
    install.packages(j)
  }
  library(j, character.only = T)
}

install.packages(c("cluster","mvoutlier","pastecs","fPortfolio"),repos="http://cran.r-project.org")
# load required packages
library(cluster)
library(mvoutlier)
library(pastecs)
library(fPortfolio)


# create data frame with dates as rownames
berndt.df = berndtInvest[, -1]
rownames(berndt.df) = as.character(as.Date(berndtInvest[, 1]))


################################################################################
# Derive the optimal portfolio weights (i.e. the weights in the tangency portfolio)
# using the CER for (i) the Minimun Variance Portfolio , (ii) the tangency portfolio. 
###############################################################################
returns.df=berndt.df[, c(1:9,11:16)]
#returns.df = berndt.df[, c(-10, -17)
exreturns.df=returns.df-berndt.df$RKFREE
returns.mat = as.matrix(exreturns.df)
# using ggplot to plot series in returns 
berndt.df$date <- as.Date(row.names(berndt.df))

# Create the time series plot using ggplot
ggplot(data = berndt.df, aes(x = date, y = WEYER)) +
  geom_line() +                    # Add a line plot
  labs(x = "Date", y = "WEYER")  # Label the axes


#
# compute global min variance portfolio
#
# use CER model: estimate the relevant unknown parameters with the sample covariances
returns.mat = as.matrix(exreturns.df)
n.obs = nrow(returns.mat)
cov.sample=var(returns.mat)
mu = matrix(colMeans(returns.mat), nrow = ncol(returns.mat), ncol = 1)
e = matrix(1, nrow = nrow(cov.sample), ncol = 1) # unitary column vector e
#
# compute GMIN portfolio
#
w.gmin.sample = solve(var(returns.mat))%*%rep(1,nrow(cov.sample))
w.gmin.sample = w.gmin.sample/sum(w.gmin.sample)
berndt.df$GMIN<-returns.mat%*%w.gmin.sample

barplot(t(w.gmin.sample), horiz=F, main="Weights", col="blue", cex.names = 0.75, las=2)

ggplot(data = berndt.df, aes(x = date, y = GMIN)) +
  geom_line() +                    # Add a line plot
  labs(x = "Date", y = "GMIN")  # Label the axes

#
# compute tangency portfolio
#
w.tan.sample = (solve(cov.sample)%*%as.numeric(mu))
w.tan.sample =w.tan.sample/as.numeric(t(e)%*%(solve(cov.sample)%*%(mu)))

berndt.df$TAN<-returns.mat%*%w.tan.sample

# visualize the differences
par(mfrow=c(1,2))
barplot(t(w.tan.sample), horiz=T, main="Tangency Port CER", col="blue", cex.names = 0.75, las=1)
barplot(t(w.gmin.sample), horiz=T, main="Min Var Port CER", col="red", cex.names = 0.75, las=1)
par(mfrow=c(1,1))

plot <- ggplot(data= berndt.df, aes(x = date)) +
  geom_line(aes(y = TAN, color = "TAN"), size = 1) +
  geom_line(aes(y = GMIN, color = "GMIN"), size = 1) +
  labs(title = "Returns",
       x = "Time", y = "Monthly Returns") +
  scale_color_manual(values = c("TAN" = "red", "GMIN" = "blue")) +
  theme_minimal() +
  theme(axis.line = element_line(color = "black"))

print(plot)


################################################################################
# Graphs the value over-time of 1 dollar invested in 1978:1 until the end of the 
# available sample in the two alternative tangency portfolios and in the market 
###############################################################################
berndt.df$Port_mkt <- berndt.df$Port_TAN<- berndt.df$Port_GMIN <- array(data = NA, dim = nrow(berndt.df))

berndt.df[1, c("Port_mkt", "Port_TAN","Port_GMIN")] <- 1 
t1<-nrow(berndt.df)
for (i in 2:t1) {
  berndt.df[i, "Port_mkt"][[1]]=berndt.df[i-1, "Port_mkt"][[1]]*(1+berndt.df[i, "MARKET"][[1]])
  berndt.df[i, "Port_TAN"][[1]]=berndt.df[i-1, "Port_TAN"][[1]]*(1+berndt.df[i, "TAN"][[1]])
  berndt.df[i, "Port_GMIN"][[1]]=berndt.df[i-1, "Port_GMIN"][[1]]*(1+berndt.df[i, "GMIN"][[1]])
}


# time series Plot of the three Portfolios

plot <- ggplot(data= berndt.df, aes(x = date)) +
  geom_line(aes(y = Port_mkt, color = "Port_mkt"), size = 1) +
  geom_line(aes(y = Port_GMIN, color = "Port_GMIN"), size = 1) +
  geom_line(aes(y = Port_TAN, color = "Port_TAN"), size = 1) +
  labs(title = "Returns",
       x = "Time", y = "Monthly Returns") +
  scale_color_manual(values = c("Port_mkt" = "red", "Port_GMIN" = "blue","Port_TAN" = "green")) +
  theme_minimal() +
  theme(axis.line = element_line(color = "black"))


# compare means and sd values on global min variance portfolios

mu.gmin.sample = as.numeric(colMeans(berndt.df$GMIN))
mu.tan.sample = as.numeric(colMeans(berndt.df$TAN))
sd.gmin.sample = as.numeric(apply(berndt.df$GMIN,2,sd))
sd.tan.sample = as.numeric(apply(berndt.df$TAN,2,sd))
cbind(mu.tan.sample,mu.gmin.sample, sd.tan.sample, sd.gmin.sample)

## ------------------------------
#  BACKTESTING   with fPortfolio
## ------------------------------
companies <- colnames(berndt.df)
#getting the data in ts format 
tsdata <- ts(berndt.df, start = c(1978, 1), frequency = 12, names = companies)
data01ts <- as.timeSeries(tsdata)
ddown<-drawdowns(data01ts)
ddowndata <- ts(ddown, start = c(1978, 1), frequency = 12, names = companies)
s1 <- window(ddowndata[, "TAN"], start = c(1978, 1), end = c(1987, 12))
dygraph(s1, ylab = "TAN", main = "drawdowns")
drawdownsStats(data01ts[, "TAN"])
## ------------------------------
#  out-of-sample BACKTESTING   
## ------------------------------

Data <- data01ts
Spec <- portfolioSpec()
Constraints <- "LongOnly"
Backtest <- portfolioBacktest()
setWindowsHorizon(Backtest) <- "60m"
equidistWindows(data = Data, backtest = Backtest)


#Specify assets for backtesting
#Formula <- MARKET ~ CITCRP + CONED + CONTIL + DATGEN + DEC + DELTA +
# +     GENMIL + GERBER +IBM+MOBIL+PANAM+PSNH+TANDY+TEXACO+WEYER
Formula <- MARKET ~ CITCRP + CONED + CONTIL + DATGEN + DEC + DELTA +
  GENMIL + GERBER + IBM + MOBIL + PANAM + PSNH + TANDY + TEXACO + WEYER

#Optimize rolling portfolios and run backtests
#btportfolios <- portfolioBacktesting(formula = Formula,
#                                       +data = data01ts, spec = Spec, constraints = Constraints,
#                                      + backtest = Backtest, trace = FALSE)

btportfolios <- portfolioBacktesting(formula = Formula,
                                     data = data01ts, spec = Spec, constraints = Constraints,
                                     backtest = Backtest, trace = FALSE)


#Weights are rebalanced on a monthly basis
Weights <- round(100 * btportfolios$weights, 2)[1:60, ]
Weights

setSmootherLambda(btportfolios$backtest) <- "1m"
SmoothPortfolios <- portfolioSmoothing(object = btportfolios,trace = FALSE)
smoothWeights <- round(100 * SmoothPortfolios$smoothWeights,2)[1:60, ]
smoothWeights

backtestPlot(SmoothPortfolios, cex = 0.6, font = 1, family = "mono")

netPerformance(SmoothPortfolios)

## --------------------
#  MODEL-SIMULATION   
## --------------------

## -----------------------------------
# model specification and estimation 
mod_TAN <- lm(berndt.df$TAN ~ 1)
summary(mod_TAN) 


## -----------------------------------
# parameter calibration and choice of the number of replications in the simulation and of the sample size for simulated data 
vol <- sd(mod_TAN$residuals)
alpha <- mod_TAN$coefficients[[1]]
nrep <- 1000
TT <- nrow(berndt.df)
# here I create the containers to be filled with the generated data.
y_bt <- y_mc <- array(1, c(TT, nrep))
x_bt <- x_mc <- array(alpha, c(TT, nrep))

# now, the loop

for (i in 1:nrep){
  u <- rnorm(TT)
  res <- sample(mod_TAN$residuals, replace = T) # this (re)samples from the data
  
  x_mc[, i] <- alpha + vol * u # the Monte Carlo way
  x_bt[, i] <- alpha+res             # the bootstrap way
  
  # now we simply construct and store the bootstrapped and MC cumulative returns
  for (j in 2:TT){
    y_mc[j, i] <- y_mc[j-1, i] * (1 + x_mc[j, i])
    y_bt[j, i] <- y_bt[j-1, i] * (1 + x_bt[j, i])
  }
}

# now we want to construct the series of means and quantiles of the resulting collection of drawn series
for (i in 1:TT){
  # obtaining the means
  berndt.df$y_bt_mean[i] <- mean(y_bt[i, ])
  berndt.df$x_bt_mean[i] <- mean(x_bt[i, ])
  berndt.df$y_mc_mean[i] <- mean(y_mc[i, ])
  berndt.df$x_mc_mean[i] <- mean(x_mc[i, ])
  
  # and the quantiles
  berndt.df$y_bt_q05[i] <- quantile(y_bt[i, ], 0.05)
  berndt.df$x_bt_q05[i] <- quantile(x_bt[i, ], 0.05)
  berndt.df$y_mc_q05[i] <- quantile(y_mc[i, ], 0.05)
  berndt.df$x_mc_q05[i] <- quantile(x_mc[i, ], 0.05)
  
  berndt.df$y_bt_q95[i] <- quantile(y_bt[i, ], 0.95)
  berndt.df$x_bt_q95[i] <- quantile(x_bt[i, ], 0.95)
  berndt.df$y_mc_q95[i] <- quantile(y_mc[i, ], 0.95)
  berndt.df$x_mc_q95[i] <- quantile(x_mc[i, ], 0.95)
  
}


## --------------------------------------------------------------------------------------------------
# plotting 
plot <- ggplot(data= berndt.df, aes(x = date)) +
  geom_line(aes(y = TAN, color = "TAN"), size = 1) +
  geom_line(aes(y = x_mc_mean, color = "x_mc_mean"), size = 1) +
  geom_line(aes(y = x_mc_q05, color = "x_mc_q05"), size = 1) +
  geom_line(aes(y = x_mc_q95, color = "x_mc_q95"), size = 1) +
  labs(title = "Simulation",
       x = "Time", y = "Monthly Returns") +
  scale_color_manual(values = c("TAN" = "blue", "x_mc_q05" = "red","x_mc_q95" = "red","x_mc_mean" = "green")) +
  theme_minimal() +
  theme(axis.line = element_line(color = "black"))

## --------------------------------------------
# Value at Risk  via Monte Carlo simulation
## --------------------------------------------
s1_mc=x_mc[2,]
hist(s1_mc, breaks = seq(min(s1_mc), max(s1_mc), l = 20+1),prob=TRUE, main = "histogram of monthly returns")
curve(dnorm(x,mean=mean(s1_mc),sd=sd(s1_mc)),col='darkblue',lwd=2,add=TRUE)
VaR_mc <- quantile(s1_mc, 0.05)
VaR_mc