library("quantmod")
library("ggfortify")

# Install BEKKs from source
# install.packages("~/BEKKs_1.4.4.tar.gz", repos = NULL, type="source")
# or from CRAN: 
# install.packages("BEKKs")

library("BEKKs")

# Bivariate BEKK(1,1) model with S&P 500 ETF and 20+ Treasury bond ETF

assets <- c('SPY', 'TLT')

getSymbols(assets, from = "2006-05-01", to = "2022-03-10",
           src = "yahoo")

s_data <- na.omit(ROC(merge(Ad(SPY), Ad(TLT)), type = 'discrete')) * 100

# alternative if retrieving data from yahoo finance is not working
# load("s_data.RData")

colnames(s_data) <- assets
autoplot(s_data) + theme_bw()

# Fitting a symmetric BEKK(1,1) model
objBEKK1 <- bekk_spec()
m1 <- bekk_fit(objBEKK1, s_data, QML_t_ratios = FALSE, max_iter = 50)
summary(m1)

library("tseries")
jarque.bera.test(residuals(m1)[,1])
jarque.bera.test(residuals(m1)[,2])

portmanteau.test(m1, lags = 5)
portmanteau.test(m1, lags = 15)
portmanteau.test(m1, lags = 30)

objBEKK1.1 <- bekk_spec()
m1.1 <- bekk_fit(objBEKK1.1, s_data, QML_t_ratios = TRUE, max_iter = 50)
summary(m1.1)


# Illustrating BHHH convergence
#plot(m1.1, diagnostic = TRUE)

# Illustrating estimated standard deviation and correlations (10x7)
plot(m1.1)

# Calculating Value at risk for a standard 60/40 stocks and bonds portfolio
m1.1_var <- VaR(m1.1, p = 0.99, portfolio_weights = c(0.6, 0.4), distribution = "empirical")
plot(m1.1_var)


# Large BEKK(1,1) model with S&P 500 ETF, 20+ Treasury bond ETF, Gold ETF, BTC and US Oil fund (tracking WTI prices)

assets <- c('SPY', 'TLT', 'GLD', 'USO')

getSymbols(assets, from = "2006-05-01", to = "2022-03-10",
           src = "yahoo")


s_data_4dim <- na.omit(ROC(merge(Ad(SPY), Ad(TLT), Ad(GLD), Ad(USO)), 
                           type = "discrete")) * 100

# alternative if retrieving data from yahoo finance is not working
# load("s_data_4dim.RData")
colnames(s_data_4dim) <- assets


# Fitting an asymmetric BEKK(1,1) model
objBEKK2 <- bekk_spec(
  model = list(type = "bekk", asymmetric = TRUE)
)
m2 <- bekk_fit(objBEKK2, s_data_4dim, QML_t_ratios = TRUE, max_iter = 150)

# Fitting an asymmetric BEKK(1,1) model with signs = c(-1,1,1,-1)
objBEKK2.1 <- bekk_spec(
  model = list(type = "bekk", asymmetric = TRUE), signs = c(-1, 1, 1, -1)
)
m2.1 <- bekk_fit(objBEKK2.1, s_data_4dim, QML_t_ratios = TRUE, max_iter = 150)
summary(m2.1)

objBEKK2.2 <- bekk_spec(
  model = list(type = "bekk", asymmetric = FALSE)
)
m2.2 <- bekk_fit(objBEKK2.2, s_data_4dim, QML_t_ratios = TRUE, max_iter = 150)

objBEKK2.3 <- bekk_spec(
  model = list(type = "dbekk", asymmetric = FALSE)
)
m2.3 <- bekk_fit(objBEKK2.3, s_data_4dim, QML_t_ratios = TRUE, max_iter = 150)

objBEKK2.4 <- bekk_spec(
  model = list(type = "dbekk", asymmetric = TRUE)
)
m2.4 <- bekk_fit(objBEKK2.4, s_data_4dim, QML_t_ratios = TRUE, max_iter = 150)

objBEKK2.5 <- bekk_spec(
  model = list(type = "sbekk", asymmetric = FALSE)
)
m2.5 <- bekk_fit(objBEKK2.5, s_data_4dim, QML_t_ratios = TRUE, max_iter = 150)

objBEKK2.6 <- bekk_spec(
  model = list(type = "sbekk", asymmetric = TRUE)
)
m2.6 <- bekk_fit(objBEKK2.6, s_data_4dim, QML_t_ratios = TRUE, max_iter = 150)

logLik(m2, m2.1, m2.2, m2.3, m2.4, m2.5, m2.6)


# Illustrating estimated standard deviation and correlations
plot(m2.1)

# Calculating Value at risk for a weighted portfolio

m2.1_var <- VaR(m2.1, p = 0.99, portfolio_weights = c(0.6, 0.2, 0.1, 0.1),
                distribution ="t")
plot(m2.1_var)

#backtesting VaR with best asymmetric model
m2.1_backtest <- backtest(m2.1, window_length = 2600, p = 0.99, 
                          portfolio_weights = c(0.6, 0.2, 0.1, 0.1),
                          n.ahead = 1, nc = 8,  distribution = "t")
plot(m2.1_backtest)
summary(m2.1_backtest)

m2.6_backtest <- backtest(m2.6, window_length = 2600, p = 0.99,
                          portfolio_weights = c(0.6, 0.2, 0.1, 0.1), 
                          n.ahead = 1, nc = 8,  distribution = "t")

summary(m2.6_backtest)
plot(m2.6_backtest)


# virfs for Lehman default

m2.2_virf <- virf(m2.2, time = "2008-09-15", q = 0.01, 
                  index_series = 1, n.ahead = 260, ci = 0.9, 
                  time_shock = FALSE)
plot(m2.2_virf)


## Code for benchmark ----

# computation time comparison between BEKKs, mgarchBEKK and bmgarch
library("mgarchBEKK")
library("bmgarch")


spec <- bekk_spec()
system.time(bekk_fit(spec, s_data_4dim, max_iter = 150))
system.time(BEKK(as.matrix(s_data_4dim)))
system.time(bmgarch(s_data_4dim, parameterization = "BEKK"))