rm(list=ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
# packages used
listofpackages <- c("tidyverse","dygraphs", "rugarch", "forecast","dplyr","ellipse","reshape2","ggplot2","xts","xlsx","readxl")

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

# setting the seed for replication
set.seed(77)

raw_data           = read_xlsx("../data/2023_monthly_stocks.xlsx") 
names(raw_data)[1] = 'Date'
typeof(raw_data)
typeof(raw_data$Date)
typeof(raw_data$AXP)
typeof(raw_data$CSCO)

dates <-seq(as.Date("1985-02-01"),length=462, by="months")
params <- c("Date","BA", "DJI")
data <- raw_data[, c(params)]
data<- na.omit(data)
data <- data %>% 
  mutate(Date = as.Date(Date, format = "%Y-%m-%d"))

params1 <- c("BA", "DJI")
tsdata <- xts(raw_data[, c(params1)], order.by=dates) # creates a time series object
tsdata <- na.omit(tsdata) # omitting the rows with NA presence
data<- na.omit(data)
## having created the database with all observation we generate a subset 
#tsdata1 <- tsdata["1992-02-01/1993-02-01"]
#data=subset(data,select=c(1:12))

## --------------------
# DATA TRANSFORMATIONS 
## --------------------
#1. from prices to returns
# exact monthly returns 
t1<-nrow(data)
data$BA_ret <- data$DJI_ret <- array(data = NA, dim = t1)
for (i in 2:t1) {
  data[i, "BA_ret"][[1]]=(data[i, "BA"][[1]]-data[i-1, "BA"][[1]])/data[i-1, "BA"][[1]]
  data[i, "DJI_ret"][[1]]=(data[i, "DJI"][[1]]-data[i-1, "DJI"][[1]])/data[i-1, "DJI"][[1]]
}


# same in .xts 
t1<-nrow(tsdata)

tsdata$BA_ret <- tsdata$DJI_ret<- array(data = NA, dim = t1)

for (i in 2:t1) {
  tsdata[i, "BA_ret"][[1]]=(tsdata[i, "BA"][[1]]-tsdata[i-1, "BA"][[1]])/data[i-1, "BA"][[1]]
  tsdata[i, "DJI_ret"][[1]]=(tsdata[i, "DJI"][[1]]-tsdata[i-1, "DJI"][[1]])/data[i-1, "DJI"][[1]]
 
}

## ----------------------------------
## VAR with CER-CAPM
## ----------------------------------


## -----------------------------------
## MODEL SPECIFICATION AND ESTIMATION
## -----------------------------------
start_date <- as.Date("1992-03-01")  # Replace with your start date
end_date <- as.Date("2005-12-01")    # Replace with your end date

# Extract observations between 'start_date' and 'end_date'

data_est <- subset(x = data,  Date >= start_date & Date <= end_date)

# estimation
cer_mkt <- lm(data_est$DJI_ret ~ 1)
capm_BA <- lm(data_est$BA_ret ~ data_est$DJI_ret)
summary(cer_mkt)
summary(capm_BA)

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

tt <- as.Date("2006-01-01")
tT <- as.Date("2015-12-01")
data_sim <- subset(x = data,  Date >= tt & Date <= tT)

# creating the containers
nrep <- 1000
BA_bt_2 <- mkt_bt_2 <- array(0, c(length(data_sim$DJI_ret), nrep))

# resampling the residuals
res_mkt_bt_2 <- matrix(sample(resid(cer_mkt), size = length(data_sim$DJI_ret) * nrep, replace = T),
                     nrow = length(data_sim$DJI_ret), ncol = nrep)
res_BA_bt_2 <- matrix(sample(resid(capm_BA), size = length(data_sim$DJI_ret) * nrep, replace = T), 
                    nrow = length(data_sim$DJI_ret), ncol = nrep)


# the loop
for (i in 1:nrep){
  for (j in 1:length(data_sim$DJI_ret)){
    mkt_bt_2[j, i] <- coef(cer_mkt)[1] + res_mkt_bt_2[j, i]
    BA_bt_2[j, i] <- coef(capm_BA)[1] + coef(capm_BA)[2] * mkt_bt_2[j, i] + res_BA_bt_2[j, i]
  }
}

# the quantiles
var_BA_capm <- array(0, length(data_sim$DJI_ret))
for (j in 1:length(data_sim$DJI_ret)){
  var_BA_capm[j] <- quantile(BA_bt_2[j, ], probs = 0.05)
}
data_sim$var_BA_capm<-var_BA_capm
# plotting
ggplot(data_sim, aes(x = Date)) +
  geom_line(aes(y = BA_ret, color = "BA"), size = 1, linetype = "solid") +
  geom_line(aes(y = var_BA_capm, color = "VaR"), size = 1, linetype = "solid") +
  labs(x = "Date", y = "Returns and VaR") +
  ylim(-0.15, 0.15) +
  theme_minimal() +
  theme(
    legend.position = c(0.15, 0.95),  # Set the legend position (top-left)
    legend.title = element_blank(),
    legend.text = element_text(size = 8),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 12, hjust = 0.5)
  ) +
  scale_color_manual(
    values = c("blue", "green"),
    labels = c("BA", "VaR")
  )



## -------------------
## GARCH MODELLING 
## -------------------

# the market GARCH regression
## specification
mkt_garch <- ugarchspec(variance.model = list(garchOrder = c(1, 1)),
                      mean.model = list(armaOrder = c(0, 0)))
## estimation
mkt_garchfit <- ugarchfit(mkt_garch, data = data_est$DJI_ret)
mkt_garchfit


# forecasting and plotting the results
horizon <- 10*12 # ten years
mygarchforecast <- ugarchforecast(mkt_garchfit, n.ahead = 10*12)

plotdata <- cbind(mygarchforecast@forecast$seriesFor,
                  mygarchforecast@forecast$seriesFor + mygarchforecast@forecast$sigmaFor*1.96,
                  mygarchforecast@forecast$seriesFor - mygarchforecast@forecast$sigmaFor*1.96)
colnames(plotdata) <- c("mean", "upper", "lower")
dygraph(ts(plotdata, start = c(2006,1), frequency = 12), main = "Forecast of the mean") %>%
  dySeries(c("lower", "mean", "upper"))

plotdata2 <- as.matrix(mygarchforecast@forecast$sigmaFor^2)
colnames(plotdata2) <- "var"
dygraph(ts(plotdata2, start = c(2006, 1), frequency = 12), main = "Forecast of the variance")


## ---------------------
## GARCH SIMULATION 
## ---------------------

## coefficients
gamma0 <- coef(mkt_garchfit)[1]
omega0 <- coef(mkt_garchfit)[2]
omega1 <- coef(mkt_garchfit)[3]
omega2 <- coef(mkt_garchfit)[4]
sigma2 <- sigma(mkt_garchfit) # this constructs the series of standard deviations conditional on information at "t-1". Is thus a vector.

# the CAPM 
## estimation
capm_BA <- lm(data_est$BA_ret ~ data_est$DJI_ret)
summary(capm_BA)

beta0 <- coef(capm_BA)[1]
beta1 <- coef(capm_BA)[2]



## -----------------
# simulation
# output containers
nrep <- 1000
BA_bt <- mkt_bt<- sigma<- array(0, c(length(data_sim$DJI_ret), nrep))

# extracting the errors and resampling 

res_mkt <- as.numeric(residuals(mkt_garchfit, standardize = T)) # the standardized residuals from the market equation

res_mkt_bt <- matrix(sample(res_mkt, size = length(data_sim$DJI_ret) * nrep, replace = T),
                     nrow = length(data_sim$DJI_ret), ncol = nrep)
res_BA_bt <- matrix(sample(resid(capm_BA), size = length(data_sim$DJI_ret) * nrep, replace = T), 
                    nrow = length(data_sim$DJI_ret), ncol = nrep)

# initial values
mkt_bt[1, ] <- data_sim$DJI_ret[1] 
BA_bt[1, ] <- beta0 + beta1*mkt_bt[1, ]
sigma[1, ]<- ugarchforecast(mkt_garchfit)@forecast$sigmaFor[1] #takes the first value (the one step ahead)
# the loop
for (i in 1:nrep){
  for (j in 2:length(data_sim$DJI_ret)){
    sigma[j, i] <- sqrt(omega0+omega1*( data_sim$DJI_ret[j-1]- gamma0)^2 + omega2*(sigma[j-1, i])^2) 
    mkt_bt[j, i] <- gamma0 + res_mkt_bt[j, i] * sigma[j,i]
    BA_bt[j, i] <- beta0 + beta1 * mkt_bt[j, i] + res_BA_bt[j, i]
  }
}

# getting the quantiles
var_BA_garch <- array(0, length(data_sim$DJI_ret))
for (j in 1:length(data_sim$DJI_ret)){
  var_BA_garch[j] <- quantile(BA_bt[j, ], probs = 0.05)
}

data_sim$var_BA_garch<-var_BA_garch

# plotting
tt <- as.Date("2006-02-01")
tT <- as.Date("2015-12-01")
data_simplot <- subset(x = data_sim,  Date >= tt & Date <= tT)
ggplot(data_simplot, aes(x = Date)) +
  geom_line(aes(y = BA_ret, color = "BA"), size = 1, linetype = "solid") +
  geom_line(aes(y = var_BA_capm, color = "Var CAPM"), size = 1, linetype = "solid") +
  geom_line(aes(y = var_BA_garch, color = "VaR GARCH"), size = 1, linetype = "solid") +
  labs(x = "Date", y = "Returns and VaR") +
  ylim(-0.30, 0.20) +
  theme_minimal() +
  theme(
    legend.position = c(0.15, 0.95),  # Set the legend position (top-left)
    legend.title = element_blank(),
    legend.text = element_text(size = 8),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 12, hjust = 0.5)
  ) +
  scale_color_manual(
    values = c("blue", "green", "red"),
    labels = c("BA", "VaR CAPM", "Var GARCH")
  )

save(data_simplot,file="VaRdata.Rdata")