# 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))

# set output options
options(width = 70, digits=4)

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

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)

####################################################
# Data Loadings and Transform: Descriptive Analysis 
####################################################

# create data frame with dates as rownames
berndt.df = berndtInvest[, -1]
berndt.df$date <- as.Date(berndtInvest[, 1])
rownames(berndt.df) = as.character(berndtInvest[, 1])
colnames(berndt.df)
dimnames(berndt.df)[[2]] #command alternative to the previous one 

# transform the data and compute cumulative returns 

t0 <- which(berndt.df$date == "1978-01-01")
t1 <- which(berndt.df$date == "1987-12-01")

series_names <- c("CITCRP","CONED","CONTIL","DATGEN","DEC","DELTA","GENMIL","GERBER","IBM","MARKET","MOBIL","PANAM","PSNH","TANDY","TEXACO","WEYER","RKFREE")

for (name in series_names) {
  P_col_name <- paste0(name,"_P")
  LP_col_name<- paste0("L",P_col_name)
  berndt.df[t0, P_col_name] <- 1
  for (i in (t0+1):(t1)) {
    berndt.df[i, P_col_name][[1]] <- berndt.df[i-1, P_col_name][[1]] * (1+berndt.df[i, name][[1]] )
  } 
  berndt.df[, LP_col_name] <- log(berndt.df[, P_col_name]) 
}
# add a trend to the database 
berndt.df$TREND<- array(data = NA, dim = nrow(berndt.df))
berndt.df[t0, c("TREND")] <- 1 # don't need to repeat the value to make the array being assigned be of the same length. be careful though as it is one of the few cases of exception

############################
# Descriptive Analysis 
############################
  
  #We can now plot, please note the difference with plotting from a time-series object
  
plot(berndt.df$date[t0:t1],berndt.df$TEXACO[t0:t1],ylab = "TEXACO",xlab="year", main = "Monthly Returns", col = "blue", lwd = 2,type="l")
  
 
plot(berndt.df$date[t0:t1],berndt.df$TEXACO[t0:t1], col = 'blue', type = "l", 
     ylab = "returns TEXACO and MKT", xlab = "date",lwd = 2)
lines(y = rep(mean(berndt.df$TEXACO[t0:t1], na.rm = T), length(berndt.df$TEXACO[t0:t1])), x = berndt.df$date[t0:t1], col = "red")
lines(y = berndt.df$MARKET[t0:t1], x = berndt.df$date[t0:t1], col = "green",lwd = 2)
legend("topleft", legend = c("TEXACO", "MKT"),
       col = c("blue", "green"), lty = 1)

plot(berndt.df$date[t0:t1],berndt.df$LTEXACO_P[t0:t1], col = 'blue', type = "l", 
     ylab = "portfolios TEXACO and MKT", xlab = "date",ylim = c(-0.5, 2),lwd = 2)
lines(y = berndt.df$LMARKET_P[t0:t1], x = berndt.df$date[t0:t1], col = "green",lwd = 2)
legend("topleft", legend = c("TEXACO", "MKT"),
       col = c("blue", "green"), lty = 1)

# Create the plot using ggplot, as generated by Chat GPT 
ggplot(berndt.df, aes(x = date)) +
  geom_line(aes(y = LTEXACO_P), color = "blue", size = 2, linetype = "solid") +
  geom_line(aes(y = LMARKET_P), color = "green", size = 2, linetype = "solid") +
  labs(x = "Date", y = "Portfolios TEXACO and MKT") +
  ylim(-0.5, 2) +
  theme_minimal() +
  theme(
    legend.position = "topleft",
    legend.title = element_blank(),
    legend.text = element_text(size = 12),
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 14),
    plot.title = element_text(size = 16, hjust = 0.5)
  ) +
  scale_color_manual(
    values = c("blue", "green"),
    guide = guide_legend(override.aes = list(size = 2, linetype = "solid"))
  ) +
  guides(fill = guide_legend(override.aes = list(size = 2)))

############################
# Asset Allocation with CER 
############################
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)
n.obs = nrow(returns.mat)

#Estimation of CER model parameters
cov.sample=var(returns.mat)
mu = matrix(colMeans(returns.mat), nrow = ncol(returns.mat), ncol = 1)

#
# compute tangency portfolio
#

e = matrix(1, nrow = nrow(cov.sample), ncol = 1) # unitary column vector e
w.tan.sample = (solve(cov.sample)%*%(mu))/as.numeric(t(e)%*%(solve(cov.sample)%*%(mu)))

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


################################
# Using the fportfolio package  
################################

#returns.df=berndt.df[, c(1:9,11:16)]
#exreturns.df=returns.df-berndt.df$RKFREE
companies <- colnames(exreturns.df)
#ts
tsdata <- ts(exreturns.df, start = c(1978, 1), frequency = 12, names = companies)
s1 <- window(tsdata[, "TEXACO"], start = c(1978, 1), end = c(1987, 12))
dygraph(s1, ylab = "TEXACO", main = "monthly excess returns")
data01ts <- as.timeSeries(tsdata)
# financial data description 
ddown<-drawdowns(data01ts)
ddowndata <- ts(ddown, start = c(1978, 1), frequency = 12, names = companies)
s1 <- window(ddowndata[, "TEXACO"], start = c(1978, 1), end = c(1987, 12))
dygraph(s1, ylab = "TEXACO", main = "drawdowns")
drawdownsStats(data01ts[, "TEXACO"])
#-----------------------
# Portfolio Allocation 
#-----------------------

# Step 1 define the data in our case 15 excess returns data in data01ts
showClass("fPFOLIODATA")

lppData <- portfolioData(data = data01ts, spec = portfolioSpec())
# once the data have been defined we can get info on them 
str(lppData, width = 65, strict.width = "cut")
print(lppData)
getData(portfolioData(lppData))[-1]
getStatistics(portfolioData(lppData))

# Step 2 Set Portfolio Constraints

showClass("fPFOLIOCON")
#default constraints: long-only 
Data<-data01ts
Spec <- portfolioSpec()
setTargetReturn(Spec) <- mean(Data)
Constraints <- "LongOnly"
defaultConstraints <- portfolioConstraints(Data, Spec, Constraints)
str(defaultConstraints, width = 65, strict.width = "cut")
print(defaultConstraints)

# short constraints 
shortConstraints <- "Short"
portfolioConstraints(Data, Spec, shortConstraints)

# box constraints 
box.1 <- "minW[1:15] = 0.1"
box.2 <- "maxW[1:15] = 1" # you can have more boxes before combining them 
boxConstraints <- c(box.1, box.2)
boxConstraints
portfolioConstraints(Data, Spec, boxConstraints)


# Step 3 Computing Optimal Portfolios 

#3.0 A benchmark: equal weight portfolio 
ewSpec <- portfolioSpec()
nAssets <- ncol(data01ts)
setWeights(ewSpec) <- rep(1/nAssets, times = nAssets)
ewPortfolio <- feasiblePortfolio(
  data = data01ts,
  spec = ewSpec,
  constraints = "LongOnly")
print(ewPortfolio)

# Efficient Frontier plot 
setNFrontierPoints(ewSpec) <- 25
eff_ew_frontier <- portfolioFrontier(data = data01ts, spec = ewSpec, constraints = "LongOnly")
tailoredFrontierPlot(object = eff_ew_frontier)

#3.1 Long-Only 
tgSpec <- portfolioSpec()
setRiskFreeRate(tgSpec) <- 0
constraints <- "longOnly"
tgPortfolio <- tangencyPortfolio(
  data = data01ts,
  spec = tgSpec, constraints = constraints) 
print(tgPortfolio)

#printing the results 
col <- seqPalette(ncol(data01ts), "BuPu")
weightsPie(tgPortfolio, box = FALSE, col = col)
mtext(text = "Tangency MV Portfolio", side = 3, line = 1.5,
      font = 2, cex = 0.7, adj = 0)
weightedReturnsPie(tgPortfolio, box = FALSE, col = col)
mtext(text = "Tangency MV Portfolio", side = 3, line = 1.5,
      font = 2, cex = 0.7, adj = 0)
covRiskBudgetsPie(tgPortfolio, box = FALSE, col = col)
mtext(text = "Tangency MV Portfolio", side = 3, line = 1.5,
      font = 2, cex = 0.7, adj = 0)    

efficient_frontier <- portfolioFrontier(data = data01ts, spec = tgSpec, constraints = constraints)
print(efficient_frontier)
# Efficient Frontier plot 
setNFrontierPoints(tgSpec) <- 25
efficient_frontier <- portfolioFrontier(data = data01ts, spec = tgSpec, constraints = constraints)
tailoredFrontierPlot(object = efficient_frontier)



#---------------------
#3.2 Box-Constraints  
#---------------------
boxSpec <- portfolioSpec()
setRiskFreeRate(boxSpec) <- 0
boxConstraints <- c("minW[1:15]=0.05", "maxW[1:15]=0.5")
tgPortfolio1 <- tangencyPortfolio(
  data = data01ts,
  spec = boxSpec, constraints = boxConstraints) 
print(tgPortfolio1)

#printing the results 
col <- seqPalette(ncol(data01ts), "BuPu")
weightsPie(tgPortfolio1, box = FALSE, col = col)
mtext(text = "Tangency MV Portfolio", side = 3, line = 1.5,
      font = 2, cex = 0.7, adj = 0)
weightedReturnsPie(tgPortfolio, box = FALSE, col = col)
mtext(text = "Tangency MV Portfolio", side = 3, line = 1.5,
      font = 2, cex = 0.7, adj = 0)
covRiskBudgetsPie(tgPortfolio, box = FALSE, col = col)
mtext(text = "Tangency MV Portfolio", side = 3, line = 1.5,
      font = 2, cex = 0.7, adj = 0)