scripts/fooCopula.R

library(rugarch)
library(xts)
library(ADGofTest)
library(qqtest)
library(copula)
library(qrmtools)
library(scrAndFun)
library(Data)
library(dplyr)
library(ggplot2)

outdir = "C:/Users/Soren Schwartz/Dropbox/Egne dokumenter/Skole/master/opgave/Figures/Copula/"

coinList.xts <- list(
  BitCoinYahoo.xts = xts::xts(BitCoinYahoo[,2:7], order.by = BitCoinYahoo[,1]),
  BitCashYahoo.xts = xts::xts(BitCashYahoo[,2:7], order.by = BitCashYahoo[,1]),
  EthereumYahoo.xts = xts::xts(EthereumYahoo[,2:7], order.by = EthereumYahoo[,1]),
  RippleYahoo.xts = xts::xts(RippleYahoo[,2:7], order.by = RippleYahoo[,1])
)

MSCI.xts <- list(
  days = xts::xts(MSCI_DP[,2], order.by = MSCI_DP[,1]),
  weeks = xts::xts(MSCI_WP[,2], order.by = MSCI_WP[,1]),
  months = xts::xts(MSCI_MP[,2], order.by = MSCI_MP[,1])
)

dailyExcess <- Log_Period_Vol(listCoins = coinList.xts, listIndex = MSCI.xts,
                              period = "days",
                              delta = 60/61, annu = 365.25,
                              dateRange = list(
                                BitCoin = "2013/",
                                BitCash = "2017-08-01/",
                                Ethereum = "2016-02-01/",
                                Ripple = "2014-01-01/",
                                Cum = "2013/"
                              ))

crypto <- do.call(cbind, lapply(dailyExcess, function(x) {x$excessReturn}))[,1:4]
colnames(crypto) <- names(dailyExcess)[1:4]
crypto <- crypto[, c("BitCoin", "Ethereum", "Ripple")]



# Index portfolios without NA's
ind_BTC_ETH <- unique(which(is.na(crypto[,c("BitCoin","Ethereum")]), arr.ind = TRUE)[,1])
ind_BTC_XRP <- unique(which(is.na(crypto[,c("BitCoin","Ripple")]), arr.ind = TRUE)[,1])
ind_ETH_XRP <- unique(which(is.na(crypto[,c("Ethereum","Ripple")]), arr.ind = TRUE)[,1])
# Portfolio selection
crypto_BTC_ETH <- crypto[-ind_BTC_ETH,c("BitCoin","Ethereum")]
crypto_BTC_XRP <- crypto[-ind_BTC_XRP,c("BitCoin","Ripple")]
crypto_ETH_XRP <- crypto[-ind_ETH_XRP,c("Ethereum","Ripple")]
# BTC_ETH_train <- crypto_BTC_ETH["/2017-06"]
# BTC_XRP_train <- crypto_BTC_XRP["/2017-06"]
# BTC_ETH_test <- crypto_BTC_ETH["2017-07/"]
# BTC_XRP_test <- crypto_BTC_XRP["2017-07/"]
backtest_period <- 250
BTC_ETH_train <- crypto_BTC_ETH[1:(nrow(crypto_BTC_ETH)-backtest_period),]
BTC_XRP_train <- crypto_BTC_XRP[1:(nrow(crypto_BTC_XRP)-backtest_period),]
ETH_XRP_train <- crypto_ETH_XRP[1:(nrow(crypto_ETH_XRP)-backtest_period),]
# BTC_ETH_test <- crypto_BTC_ETH[(nrow(crypto_BTC_ETH)-backtest_period+1):nrow(crypto_BTC_ETH),]
## Fit ARMA
spec <- ugarchspec(variance.model = list(garchOrder = c(1, 1)),
                     mean.model = list(armaOrder = c(1, 1)),
                     distribution.model = "sstd")

# BTC_ETH
fitted_BTC_ETH <- fit_ARMA_GARCH(BTC_ETH_train, ugarchspec.list = rep(list(spec), 2))
Z_BTC_ETH <- as.matrix(do.call(merge, lapply(fitted_BTC_ETH$fit,
                                             rugarch::residuals, standardize = TRUE)))
colnames(Z_BTC_ETH) <- colnames(BTC_ETH_train)
U_BTC_ETH <- pobs(Z_BTC_ETH)

# BTC_XRP
fitted_BTC_XRP <- fit_ARMA_GARCH(BTC_XRP_train, ugarchspec.list = rep(list(spec), 2))
Z_BTC_XRP <- as.matrix(do.call(merge, lapply(fitted_BTC_XRP$fit,
                                             rugarch::residuals, standardize = TRUE)))
colnames(Z_BTC_XRP) <- colnames(BTC_XRP_train)
U_BTC_XRP <- pobs(Z_BTC_XRP)

# BTC_ETH
fitted_ETH_XRP <- fit_ARMA_GARCH(ETH_XRP_train, ugarchspec.list = rep(list(spec), 2))
Z_ETH_XRP <- as.matrix(do.call(merge, lapply(fitted_ETH_XRP$fit,
                                             rugarch::residuals, standardize = TRUE)))
colnames(Z_ETH_XRP) <- colnames(ETH_XRP_train)
U_ETH_XRP <- pobs(Z_ETH_XRP)

coef_arma_garch <- cbind(
sapply(1:2, function(j) fitted_BTC_ETH$fit[[j]]@fit$coef),
sapply(1:2, function(j) fitted_BTC_XRP$fit[[j]]@fit$coef),
sapply(1:2, function(j) fitted_ETH_XRP$fit[[j]]@fit$coef)) %>%
  round(digits = 4)
colnames(coef_arma_garch) <-  c("BTC_eth", "btc_ETH", "BTC_xrp",
                                "btc_XRP", "ETH_xrp", "eth_XRP")
coef_arma_garch
xtable::xtable(coef_arma_garch, digits = 4)

# Test if uniform margins
pdf(file = paste0(outdir, "pit_cop_models.pdf"))
par(mfrow= c(3,2))
PIT_1_ETH_XRP <- GAS::PIT_test(U_ETH_XRP[,1], G = 20, alpha = 0.05, plot = TRUE)
PIT_2_ETH_XRP <- GAS::PIT_test(U_ETH_XRP[,2], G = 20, alpha = 0.05, plot = TRUE)
#
PIT_1_BTC_XRP <- GAS::PIT_test(U_BTC_XRP[,1], G = 20, alpha = 0.05, plot = TRUE)
PIT_2_BTC_XRP <- GAS::PIT_test(U_BTC_XRP[,2], G = 20, alpha = 0.05, plot = TRUE)
#
PIT_1_BTC_ETH <- GAS::PIT_test(U_BTC_ETH[,1], G = 20, alpha = 0.05, plot = TRUE)
PIT_2_BTC_ETH <- GAS::PIT_test(U_BTC_ETH[,2], G = 20, alpha = 0.05, plot = TRUE)
par(mfrow = c(1,1))
dev.off()

hist_p_values <- cbind(PIT_1_ETH_XRP$Hist$pvalue,
PIT_2_ETH_XRP$Hist$pvalue,
PIT_1_BTC_XRP$Hist$pvalue,
PIT_2_BTC_XRP$Hist$pvalue,
PIT_1_BTC_ETH$Hist$pvalue,
PIT_2_BTC_ETH$Hist$pvalue)
colnames(hist_p_values) <-  c("BTC_eth", "btc_ETH", "BTC_xrp",
                                "btc_XRP", "ETH_xrp", "eth_XRP")
rownames(hist_p_values) <- c("p-value")
xtable::xtable(hist_p_values)

idd_p_values <-  cbind(PIT_1_ETH_XRP$IID$pvalue,
                       PIT_2_ETH_XRP$IID$pvalue,
                       PIT_1_BTC_XRP$IID$pvalue,
                       PIT_2_BTC_XRP$IID$pvalue,
                       PIT_1_BTC_ETH$IID$pvalue,
                       PIT_2_BTC_ETH$IID$pvalue)
colnames(idd_p_values) <-  c("$BTC_{ETH}$", "$ETH_{BTC}$", "$BTC_{XRP}$",
                             "$XRP_{BTC}$", "$ETH_{XRP}$", "$XRP_{ETH}$")
rownames(idd_p_values) <- c("$E[X]$", "$E[X^2]$", "$E[X^3]$", "$E[X^4]$")
xtable::xtable(idd_p_values)

# Fit copula
(LL_c_BTC_ETH <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta =0.5),
                                           function(par) {
                                             U_foo = U_BTC_ETH
                                             LL_clayton_cpp(par, U_foo)},
                                           lower = c(rep(-Inf,3)),
                               upper = c(rep(Inf, 3)),control = list(rel.tol = 1e-14,
                                                                     xf.tol = 2.2e-40)))
(LL_c_BTC_XRP <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta =0.5),
                                   function(par) {
                                     U_foo = U_BTC_XRP
                                     LL_clayton_cpp(par, U_foo)},
                                   lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
(LL_c_ETH_XRP <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta =0.5),
                                   function(par) {
                                     U_foo = U_ETH_XRP
                                     LL_clayton_cpp(par, U_foo)},
                                   lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
# Parameter estimates
(parameter <- cbind(BTC_ETH = LL_c_BTC_ETH$par,
                    BTC_XRP = LL_c_BTC_XRP$par,
                    ETH_XRP = LL_c_ETH_XRP$par))
xtable::xtable(parameter, digits = 3)
# Simulate
backtest_cop_c <- function(return = X_all, backtest_period = 200,
                           alpha = 0.99, m = 200, par) {
  future <- matrix(ncol = 1, nrow = backtest_period)
  CI_out <- matrix(ncol = 2, nrow = backtest_period)
  VaR_out <- matrix(ncol = 1, nrow = backtest_period)

  depen <- numeric(backtest_period)

  sim_list = list()

  for (i in 0:(backtest_period - 1)) {
    nY <- nrow(return) - (backtest_period - i)
    X <- return[1:nY, ]
    ## Fit marginal time series
    spec <- ugarchspec(variance.model = list(garchOrder = c(1, 1)),
                       mean.model = list(armaOrder = c(1, 1)),
                       distribution.model = "sstd")
    fit.ARMA.GARCH <- fit_ARMA_GARCH(X, ugarchspec.list = rep(list(spec), 2))
    stopifnot(sapply(fit.ARMA.GARCH$error, is.null)) # NULL = no error
    if(FALSE)
      fit.ARMA.GARCH$warning
    ## => Warning comes from finding initial values and can be ignored here
    fits <- fit.ARMA.GARCH$fit # fitted models
    Z <- as.matrix(do.call(merge, lapply(fits, residuals, standardize = TRUE))) # grab out standardized residuals
    colnames(Z) <- colnames(X)

    nu_mar <- vapply(fits, function(x) { x@fit$coef[["shape"]]}, NA_real_)
    gamma_mar <- vapply(fits, function(x) { x@fit$coef[["skew"]]}, NA_real_)

    d <- ncol(X) # dimension
    ## Compute pseudo-observations from the skewed standardized t residuals
    U <- pobs(Z)
    depen[i+1] <- tail(LL_clayton_sim_cpp(par, U), 1)
    ### Simulating paths from the full model ##
    m <- m
    U. <- rCopula(m, copula = claytonCopula(param = depen[i+1], dim = d))
    Z. <- sapply(1:d, function(j) {skewt::qskt(U.[,j], df = nu_mar[j],
                                               gamma = gamma_mar[j])})
    sim <- lapply(1:d, function(j)
      ugarchsim(fits[[j]], n.sim = 1, m.sim = m,
                startMethod="sample",
                custom.dist = list(name = "sample", distfit = t(Z.[,j, drop = FALSE]))))
    fitted_sim <- sapply(sim, function(x) fitted(x))
    colnames(fitted_sim) <- colnames(return)
    Xs. <- rowSums(fitted_sim)

    Xs.mean <- mean(Xs.)
    Xs.CI <- quantile(Xs., probs = c(0.025, 0.975))
    alpha <- alpha
    VaR <- quantile(Xs., probs = alpha)

    future[i+1, 1] <- Xs.mean
    CI_out[i+1, ] <- Xs.CI
    VaR_out[i+1, 1] <- VaR

    sim_list[[i+1]] = fitted_sim
  }
  return(list(
    forecast = future,
    CI = CI_out,
    VaR = VaR_out,
    depen = depen,
    sim_list = sim_list
  ))
}
# Back test using Claytong copula
set.seed(123)
back_BTC_ETH <- backtest_cop_c(crypto_BTC_ETH, backtest_period = backtest_period,
                           alpha = 0.99, m = 10000, par = LL_c_BTC_ETH$par)
back_BTC_XRP <- backtest_cop_c(crypto_BTC_XRP, backtest_period = backtest_period,
                               alpha = 0.99, m = 10000, par = LL_c_BTC_XRP$par)
back_ETH_XRP <- backtest_cop_c(crypto_ETH_XRP, backtest_period = backtest_period,
                               alpha = 0.99, m = 10000, par = LL_c_ETH_XRP$par)
# Future index
test_idx <- zoo::index(crypto_BTC_ETH[((nrow(crypto_BTC_ETH) -
                                          backtest_period)+1):nrow(crypto_BTC_ETH), ])
cbind(back_BTC_ETH$depen, back_BTC_XRP$depen, back_ETH_XRP$depen) %>%
  tibble::as.tibble() %>%
  transmute(BTC_ETH = V1, BTC_XRP = V2, ETH_XRP = V3, Date = test_idx) %>%
  tidyr::gather(key = "Portfolio", value = "theta", -Date) %>%
  ggplot(aes(x = Date, y = theta, col = Portfolio)) +
  geom_line()

back_test_results <- lapply(1:backtest_period, function(i){
  sims <- cbind(
  (back_BTC_ETH$sim_list[[i]][,c("BitCoin")] +
     back_BTC_XRP$sim_list[[i]][,c("BitCoin")]) / 2,
  (back_BTC_ETH$sim_list[[i]][,c("Ethereum")] +
    back_ETH_XRP$sim_list[[i]][,c("Ethereum")]) / 2,
  (back_BTC_XRP$sim_list[[i]][,c("Ripple")] +
     back_ETH_XRP$sim_list[[i]][,c("Ripple")]) / 2)
  colnames(sims) <- c("BTC", "ETH", "XRP")

  Xs. <- rowSums(sims)
  Xs.mean <- mean(Xs.)
  Xs.CI <- quantile(Xs., probs = c(0.025, 0.975))
  alpha <- 0.99
  VaR <- quantile(Xs., probs = c(1-alpha,alpha))
  X_col = colSums(sims)

  return(list(sims = sims,
              forecast = Xs.mean,
              future = X_col,
              CI = Xs.CI,
              VaR = VaR))

})
future_mean <- do.call(rbind, lapply(back_test_results, function(x){x$future}))
forecast <- do.call(rbind, lapply(back_test_results, function(x){x$forecast}))
CI = t(do.call(rbind, lapply(back_test_results, function(x){x$CI})))
VaR = t(do.call(rbind, lapply(back_test_results, function(x){x$VaR})))

X <- crypto[(1:(nrow(crypto) - backtest_period))[!(1:(nrow(crypto) - backtest_period) %in% unique(which(is.na(crypto), arr.ind = TRUE)[,1]))], ]
X <- X[(nrow(X)-20):nrow(X), ]
X_future <- crypto[((nrow(crypto) - backtest_period)+1):nrow(crypto), ]
past <- zoo::index(X)
future <- zoo::index(X_future)

pdf(file = paste0(outdir, "copula_GARCH_sim.pdf"),
width = 8, height = 6)
plot(past, rowSums(X), type = "l", xlim = range(c(past, future)),
     xlab = "Time", ylab = expression(r[t]),
     main = c("copula-GARCH simulations"),
     ylim = c(-(max(range(VaR))+0.1), (max(range(VaR))+0.1)))
polygon(c(future, rev(future)), c(CI[1,], rev(CI[2,])),
        border = NA, col = "grey80") # CI region
lines(future, rowSums(X_future),  col = "orange")
lines(future, forecast, col = "royalblue3") # predicted aggregated return
lines(future, CI[1,], col = "grey50") # lower CI
lines(future, CI[2,], col = "grey50") # upper CI
lines(future, VaR[1,], col = "maroon3") # lower VaR_alpha
lines(future, VaR[2,], col = "maroon3") # upper VaR_alpha
legend("bottomright", bty = "n", lty = rep(1, 4),
       col = c("black", "orange", "royalblue3", "grey50", "maroon3"),
       legend = c("(Aggregated) returns", "True returns",
                  "(Simulated) predicted return",
                  "95% CIs",
                  as.expression(substitute("Simulated"~VaR[a], list(a = alpha)))))
dev.off()
mse = mean((rowSums(X_future) - forecast)^2)

## TSMOM extended
sigma_TS <- do.call(cbind, lapply(dailyExcess, function(x) {x$sigma}))[,1:4]
colnames(sigma_TS) <- names(dailyExcess)[1:4]
sigma_TS <- sigma_TS[, c("BitCoin", "Ethereum", "Ripple")]

sig_BTC_ETH <- sigma_TS[-ind_BTC_ETH, c("BitCoin", "Ethereum")]
sig_BTC_XRP <- sigma_TS[-ind_BTC_XRP, c("BitCoin", "Ripple")]
sig_ETH_XRP <- sigma_TS[-ind_ETH_XRP, c("Ethereum", "Ripple")]

BTC_XRP_train_TS <- as.matrix(crypto_BTC_XRP[5:(nrow(crypto_BTC_XRP)-backtest_period),])
# BTC_XRP_train_TS_lag <- as.matrix(crypto_BTC_XRP[2:(nrow(crypto_BTC_XRP)-backtest_period-12),])
sigma_BTC_XRP_1 <- as.matrix(sig_BTC_XRP[4:(nrow(sig_BTC_XRP)-backtest_period-1),])
# sigma_BTC_XRP_13 <- as.matrix(sig_BTC_XRP[1:(nrow(sig_BTC_XRP)-backtest_period-13),])
BTC_XRP_train_TS_lag <- as.matrix(crypto_BTC_XRP[2:(nrow(crypto_BTC_XRP)-backtest_period-3),])
sigma_BTC_XRP_lag <- as.matrix(sig_BTC_XRP[1:(nrow(sig_BTC_XRP)-backtest_period-4),])

ETH_XRP_train_TS <- as.matrix(crypto_ETH_XRP[5:(nrow(crypto_ETH_XRP)-backtest_period),])
ETH_XRP_train_TS_lag <- as.matrix(crypto_ETH_XRP[2:(nrow(crypto_ETH_XRP)-backtest_period-3),])
sigma_ETH_XRP_1 <- as.matrix(sig_ETH_XRP[4:(nrow(sig_ETH_XRP)-backtest_period-1),])
sigma_ETH_XRP_lag <- as.matrix(sig_ETH_XRP[1:(nrow(sig_ETH_XRP)-backtest_period-4),])

BTC_ETH_train_TS <- as.matrix(crypto_BTC_ETH[5:(nrow(crypto_BTC_ETH)-backtest_period),])
BTC_ETH_train_TS_lag <- as.matrix(crypto_BTC_ETH[2:(nrow(crypto_BTC_ETH)-backtest_period-3),])
sigma_BTC_ETH_1 <- as.matrix(sig_BTC_ETH[4:(nrow(sig_BTC_ETH)-backtest_period-1),])
sigma_BTC_ETH_lag <- as.matrix(sig_BTC_ETH[1:(nrow(sig_BTC_ETH)-backtest_period-4),])

BTC_XRP_1 <- lm((BTC_XRP_train_TS/sigma_BTC_XRP_1)[,1] ~ (BTC_XRP_train_TS_lag/sigma_BTC_XRP_lag)[,1])
BTC_XRP_2 <- lm((BTC_XRP_train_TS/sigma_BTC_XRP_1)[,2] ~ (BTC_XRP_train_TS_lag/sigma_BTC_XRP_lag)[,2])

ETH_XRP_1 <- lm((ETH_XRP_train_TS/sigma_ETH_XRP_1)[,1] ~ (ETH_XRP_train_TS_lag/sigma_ETH_XRP_lag)[,1])
ETH_XRP_2 <- lm((ETH_XRP_train_TS/sigma_ETH_XRP_1)[,2] ~ (ETH_XRP_train_TS_lag/sigma_ETH_XRP_lag)[,2])

BTC_ETH_1 <- lm((BTC_ETH_train_TS/sigma_BTC_ETH_1)[,1] ~ (BTC_ETH_train_TS_lag/sigma_BTC_ETH_lag)[,1])
BTC_ETH_2 <- lm((BTC_ETH_train_TS/sigma_BTC_ETH_1)[,2] ~ (BTC_ETH_train_TS_lag/sigma_BTC_ETH_lag)[,2])

lm_param <- cbind(BTC_ETH_1$coefficients,
                  BTC_ETH_2$coefficients,
                  BTC_XRP_1$coefficients,
                  BTC_XRP_2$coefficients,
                  ETH_XRP_1$coefficients,
                  ETH_XRP_2$coefficients)
colnames(lm_param) <-  c("BTC_eth", "btc_ETH", "BTC_xrp",
                              "btc_XRP", "ETH_xrp", "eth_XRP")
rownames(lm_param) <- c("alpha_OLS", "beta_OLS")
xtable::xtable(lm_param, digits = 4)

res_BTC_XRP <- cbind(resid(BTC_XRP_1), resid(BTC_XRP_2))
res_ETH_XRP <- cbind(resid(ETH_XRP_1), resid(ETH_XRP_2))
res_BTC_ETH <- cbind(resid(BTC_ETH_1), resid(BTC_ETH_2))
# res_BTC_XRP <- BTC_XRP_train_TS/sigma_BTC_XRP_1 - BTC_XRP_train_TS_lag/sigma_BTC_XRP_lag
# res_ETH_XRP <- ETH_XRP_train_TS/sigma_ETH_XRP_1 - ETH_XRP_train_TS_lag/sigma_ETH_XRP_lag
# res_BTC_ETH <- BTC_ETH_train_TS/sigma_BTC_ETH_1 - BTC_ETH_train_TS_lag/sigma_BTC_ETH_lag

# fit skew-t distributions
st_fit_BTC_XRP_1 = sn::st.mple(y=res_BTC_XRP[,1])
st_fit_BTC_XRP_2 = sn::st.mple(y=res_BTC_XRP[,2])

st_fit_ETH_XRP_1 = sn::st.mple(y=res_ETH_XRP[,1])
st_fit_ETH_XRP_2 = sn::st.mple(y=res_ETH_XRP[,2])

st_fit_BTC_ETH_1 = sn::st.mple(y=res_BTC_ETH[,1])
st_fit_BTC_ETH_2 = sn::st.mple(y=res_BTC_ETH[,2])


skewed_t_param <- cbind(st_fit_BTC_ETH_1$dp,
                        st_fit_BTC_ETH_2$dp,
                        st_fit_BTC_XRP_1$dp,
                        st_fit_BTC_XRP_2$dp,
                        st_fit_ETH_XRP_1$dp,
                        st_fit_ETH_XRP_2$dp)
colnames(skewed_t_param) <-  c("BTC_eth", "btc_ETH", "BTC_xrp",
                         "btc_XRP", "ETH_xrp", "eth_XRP")
xtable::xtable(skewed_t_param, digits = 3)


U_BTC_XRP_TS <- cbind(sn::pst(res_BTC_XRP[,1], dp = st_fit_BTC_XRP_1$dp),
                   sn::pst(res_BTC_XRP[,2], dp = st_fit_BTC_XRP_2$dp))
colnames(U_BTC_XRP_TS) <- colnames(BTC_XRP_train_TS)
U_ETH_XRP_TS <- cbind(sn::pst(res_ETH_XRP[,1], dp = st_fit_ETH_XRP_1$dp),
                      sn::pst(res_ETH_XRP[,2], dp = st_fit_ETH_XRP_2$dp))
colnames(U_ETH_XRP_TS) <- colnames(ETH_XRP_train_TS)
U_BTC_ETH_TS <- cbind(sn::pst(res_BTC_ETH[,1], dp = st_fit_BTC_ETH_1$dp),
                      sn::pst(res_BTC_ETH[,2], dp = st_fit_BTC_ETH_2$dp))
colnames(U_BTC_ETH_TS) <- colnames(BTC_ETH_train_TS)

# Check if uniform and independent
pdf(file = paste0(outdir, "pit_cop_models_TS.pdf"))
par(mfrow= c(3,2))
PIT_1_ETH_XRP_TS <- GAS::PIT_test(U_ETH_XRP_TS[,1], G = 20, alpha = 0.05, plot = TRUE)
PIT_2_ETH_XRP_TS <- GAS::PIT_test(U_ETH_XRP_TS[,2], G = 20, alpha = 0.05, plot = TRUE)
#
PIT_1_BTC_XRP_TS <- GAS::PIT_test(U_BTC_XRP_TS[,1], G = 20, alpha = 0.05, plot = TRUE)
PIT_2_BTC_XRP_TS <- GAS::PIT_test(U_BTC_XRP_TS[,2], G = 20, alpha = 0.05, plot = TRUE)
#
PIT_1_BTC_ETH_TS <- GAS::PIT_test(U_BTC_ETH_TS[,1], G = 20, alpha = 0.05, plot = TRUE)
PIT_2_BTC_ETH_TS <- GAS::PIT_test(U_BTC_ETH_TS[,2], G = 20, alpha = 0.05, plot = TRUE)
par(mfrow = c(1,1))
dev.off()

hist_p_values_TS <- cbind(PIT_1_ETH_XRP_TS$Hist$pvalue,
                       PIT_2_ETH_XRP_TS$Hist$pvalue,
                       PIT_1_BTC_XRP_TS$Hist$pvalue,
                       PIT_2_BTC_XRP_TS$Hist$pvalue,
                       PIT_1_BTC_ETH_TS$Hist$pvalue,
                       PIT_2_BTC_ETH_TS$Hist$pvalue)
colnames(hist_p_values_TS) <-  c("BTC_eth_TS", "btc_ETH_TS", "BTC_xrp_TS",
                              "btc_XRP_TS", "ETH_xrp_TS", "eth_XRP_TS")
rownames(hist_p_values_TS) <- c("p-value")
xtable::xtable(hist_p_values_TS)

idd_p_values_TS <-  cbind(PIT_1_ETH_XRP_TS$IID$pvalue,
                       PIT_2_ETH_XRP_TS$IID$pvalue,
                       PIT_1_BTC_XRP_TS$IID$pvalue,
                       PIT_2_BTC_XRP_TS$IID$pvalue,
                       PIT_1_BTC_ETH_TS$IID$pvalue,
                       PIT_2_BTC_ETH_TS$IID$pvalue)
colnames(idd_p_values_TS) <-  c("BTC_eth_TS", "btc_ETH_TS", "BTC_xrp_TS",
                                "btc_XRP_TS", "ETH_xrp_TS", "eth_XRP_TS")
rownames(idd_p_values_TS) <- c("$E[X]$", "$E[X^2]$", "$E[X^3]$", "$E[X^4]$")
xtable::xtable(idd_p_values_TS)


# Fit copula
(LL_BTC_XRP_TS <- stats::nlminb(start = c(w = 0.5, alpha = 1, beta =0.5),
                               function(par) {
                                 U_foo = U_BTC_XRP_TS
                                 LL_clayton_cpp(par, U_foo)},
                               lower = c(-Inf, 0, 0),
                               upper = c(rep(Inf, 3))))
(LL_ETH_XRP_TS <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta =0.5),
                               function(par) {
                                 U_foo = U_ETH_XRP_TS
                                 LL_clayton_cpp(par, U_foo)},
                               lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
(LL_BTC_ETH_TS <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta =0.5),
                               function(par) {
                                 U_foo = U_BTC_ETH_TS
                                 LL_clayton_cpp(par, U_foo)},
                               lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
# Parameter estimates
(parameter_TS <- cbind(BTC_XRP = LL_BTC_XRP_TS$par,
                       ETH_XRP = LL_ETH_XRP_TS$par,
                    BTC_ETH = LL_BTC_ETH_TS$par))

xtable::xtable(parameter_TS)
# Simulate
backtest_TSMOM_cop <- function(return = X_all, backtest_period = 200,
                           alpha = 0.99, m = 200, par_cop, par_dp = cbind(NA, NA),
                           sigma = sig, par_coef = cbind(NA, NA),
                           lag = 3) {
  depen <- numeric(backtest_period)

  future <- matrix(ncol = 1, nrow = backtest_period)
  CI_out <- matrix(ncol = 2, nrow = backtest_period)
  VaR_out <- matrix(ncol = 1, nrow = backtest_period)

  sim_list = list()

  for (i in 0:(backtest_period - 1)) {

    TS <- as.matrix(return[(lag+2):(nrow(return)-backtest_period - i),])
    TS_lag <- as.matrix(return[2:(nrow(return)-backtest_period-lag - i),])
    sigma_1 <- as.matrix(sigma[(lag+1):(nrow(sigma)-backtest_period-1 - i),])
    sigma_lag <- as.matrix(sigma[1:(nrow(sigma)-backtest_period-lag - 1 - i),])

    res <- TS/sigma_1 - par_coef[1,] - par_coef[2,] * TS_lag/sigma_lag

    U_TS <- cbind(sn::pst(res[,1], dp = par_dp[,1]),
                          sn::pst(res[,2], dp = par_dp[,2]))
    colnames(U_TS) <- colnames(TS)

    depen[i+1] <- tail(LL_clayton_sim_cpp(par_cop, U_TS), 1)
    U. <- rCopula(m, copula = claytonCopula(param = depen[i+1], dim = 2))
    # if(any(U. > 0.99997)) {
    #   U.[which(U. > 0.99997)] = 0.99997
    # }
    # if(any(U. < 1 - 0.99997)) {
    #   U.[which(U. < 1 - 0.99997)] = 1 - 0.99997
    # }
    Z. <- sapply(1:2, function(j) {sn::qst(U.[,j], dp=par_dp[,j],
                                           tol = 1e-04)})

    ini_fore <- (as.matrix(return[(nrow(return)-backtest_period-lag - 1 - i),]) /
                   as.matrix(sigma[(nrow(sigma)-backtest_period-lag - i),]))

    sim <- Z. + as.numeric(ini_fore) * par_coef[2,]# + par_coef[1,]
    colnames(sim) <- colnames(return)
    Xs. <- rowSums(sim)

    Xs.mean <- mean(Xs.)
    Xs.CI <- quantile(Xs., probs = c(0.025, 0.975))
    alpha <- alpha
    VaR <- quantile(Xs., probs = alpha)

    future[i+1, 1] <- Xs.mean
    CI_out[i+1, ] <- Xs.CI
    VaR_out[i+1, 1] <- VaR

    sim_list[[i+1]] = sim

    print(i)
  }
  return(list(
    forecast = future,
    CI = CI_out,
    VaR = VaR_out,
    depen = depen,
    sim_list = sim_list
  ))
}
set.seed(123)
system.time(
back_TS_BTC_XRP <- backtest_TSMOM_cop(return = crypto_BTC_XRP, backtest_period = backtest_period,
                                      alpha = 0.99, m = 5000, par_cop = LL_BTC_XRP_TS$par,
                                      par_dp = cbind(st_fit_BTC_XRP_1$dp, st_fit_BTC_XRP_2$dp),
                                      sigma = sig_BTC_XRP,
                                      par_coef = cbind(coef(BTC_XRP_1),
                                                        coef(BTC_XRP_2)), lag = 3)
)
system.time(
back_TS_BTC_ETH <- backtest_TSMOM_cop(return = crypto_BTC_ETH, backtest_period = backtest_period,
                                      alpha = 0.99, m = 5000, par_cop = LL_BTC_ETH_TS$par,
                                      par_dp = cbind(st_fit_BTC_ETH_1$dp, st_fit_BTC_ETH_2$dp),
                                      sigma = sig_BTC_ETH,
                                      par_coef = cbind(coef(BTC_ETH_1),
                                                        coef(BTC_ETH_2)), lag = 3)
)
system.time(
back_TS_ETH_XRP <- backtest_TSMOM_cop(return = crypto_ETH_XRP, backtest_period = backtest_period,
                                      alpha = 0.99, m = 5000, par_cop = LL_ETH_XRP_TS$par,
                                      par_dp = cbind(st_fit_ETH_XRP_1$dp, st_fit_ETH_XRP_2$dp),
                                      sigma = sig_ETH_XRP,
                                      par_coef = cbind(coef(ETH_XRP_1),
                                                        coef(ETH_XRP_2)), lag = 3)
)

# Future index
test_idx <- zoo::index(crypto_BTC_ETH[((nrow(crypto_BTC_ETH) -
                                          backtest_period)+1):nrow(crypto_BTC_ETH), ])
# cbind(back_TS_BTC_XRP$depen, back_TS_BTC_ETH$depen, back_TS_ETH_XRP$depen) %>%
#   tibble::as.tibble() %>%
#   transmute(BTC_ETH = V1, BTC_XRP = V2, ETH_XRP = V3, Date = test_idx) %>%
#   tidyr::gather(key = "Portfolio", value = "theta", -Date) %>%
#   ggplot(aes(x = Date, y = theta, col = Portfolio)) +
#   geom_line()

back_test_TS_results <- lapply(1:backtest_period, function(i){
  sims <- cbind(
    (back_TS_BTC_XRP$sim_list[[i]][,c("BitCoin")] +
       back_TS_BTC_ETH$sim_list[[i]][,c("BitCoin")]) / 2,
    (back_TS_BTC_ETH$sim_list[[i]][,c("Ethereum")] +
       back_TS_ETH_XRP$sim_list[[i]][,c("Ethereum")]) / 2,
    (back_TS_BTC_XRP$sim_list[[i]][,c("Ripple")] +
       back_TS_ETH_XRP$sim_list[[i]][,c("Ripple")]) / 2)
  colnames(sims) <- c("BTC", "ETH", "XRP")

  Xs. <- rowSums(sims)
  Xs.mean <- mean(Xs.)
  Xs.CI <- quantile(Xs., probs = c(0.025, 0.975))
  alpha <- 0.99
  VaR <- quantile(Xs., probs = c(1 - alpha, alpha))
  X_col = colSums(sims)

  return(list(sims = sims,
              forecast = Xs.mean,
              future = X_col,
              CI = Xs.CI,
              VaR = VaR))

})
future_mean_TS <- do.call(rbind, lapply(back_test_TS_results, function(x){x$future}))
forecast_TS <- do.call(rbind, lapply(back_test_TS_results, function(x){x$forecast}))
CI_TS = t(do.call(rbind, lapply(back_test_TS_results, function(x){x$CI})))
VaR_TS = t(do.call(rbind, lapply(back_test_TS_results, function(x){x$VaR})))

X_TS <- crypto[(1:(nrow(crypto) - backtest_period))[!(1:(nrow(crypto) - backtest_period) %in% unique(which(is.na(crypto), arr.ind = TRUE)[,1]))], ]
X_TS <- X_TS[(nrow(X_TS)-20):nrow(X_TS), ]
X_future_TS <- crypto[((nrow(crypto) - backtest_period)+1):nrow(crypto), ]
past_TS <- zoo::index(X_TS)
future_TS <- zoo::index(X_future_TS)

pdf(file = paste0(outdir, "TSMOM_extended_sim.pdf"),
    width = 8, height = 6)
plot(past_TS, rowSums(X_TS), type = "l", xlim = range(c(past_TS, future_TS)),
     xlab = "Time", ylab = expression(r[t]),
     main = c("TSMOM extended simulations"),
     ylim = c(-(max(range(VaR_TS))+0.1), (max(range(VaR_TS))+0.1)))
polygon(c(future_TS, rev(future_TS)), c(CI_TS[1,], rev(CI_TS[2,])),
        border = NA, col = "grey80") # CI region
lines(future_TS, rowSums(X_future_TS),  col = "orange")
lines(future_TS, forecast_TS, col = "royalblue3") # predicted aggregated return
lines(future_TS, CI_TS[1,], col = "grey50") # lower CI
lines(future_TS, CI_TS[2,], col = "grey50") # upper CI
lines(future_TS, VaR_TS[1,], col = "maroon3") # lower VaR_alpha
lines(future_TS, VaR_TS[2,], col = "maroon3") # upper VaR_alpha
legend("bottomright", bty = "n", lty = rep(1, 4),
       col = c("black", "orange", "royalblue3", "grey50", "maroon3"),
       legend = c("(Aggregated) returns", "True returns",
                  "(Simulated) predicted return",
                  "95% CIs",
                  as.expression(substitute("Simulated"~VaR[a], list(a = alpha)))))
dev.off()

# TSMOM
# X_sign_12 <- sign(as.matrix(crypto[(((nrow(crypto) - backtest_period)+1)-12):(nrow(crypto)-12), ]))
X_sign_lag <- sign(as.matrix(crypto[(((nrow(crypto) - backtest_period)+1)-1):(nrow(crypto)-1), ]))

sigma <- do.call(cbind, lapply(dailyExcess, function(x) {x$sigma}))[,1:4]
colnames(sigma) <- names(dailyExcess)[1:4]
sigma <- sigma[, c("BitCoin", "Ethereum", "Ripple")]
sigma_t_1 <- as.matrix(sigma[(((nrow(sigma) - backtest_period)+1)-1):(nrow(sigma)-1), ])

stra_TSMOM <- (1/3) * rowSums((0.4 / sigma_t_1) * X_sign_lag * as.matrix(X_future))
names(stra_TSMOM) <- NULL
cum_TSMOM <- c(0, cumsum(as.numeric(stra_TSMOM)))
# TSMOM with copula
stra_TSMOM_cop <- (1/3) * rowSums((0.4 / sigma_t_1) * sign(future_mean_TS) * as.matrix(X_future))
cum_TSMOM_cop <- c(0, cumsum(as.numeric(stra_TSMOM_cop)))
# copula-GARCH strategy
stra_cop_GARCH <- (1/3) * rowSums(sign(future_mean) * X_future * (0.4 / sigma_t_1))
cum_cop_GARCH <- c(0, cumsum(as.numeric(stra_cop_GARCH)))
# Buy and hold
stra_Buy_Hold <- (1/3) *  rowSums((0.4 / sigma_t_1) * X_future)
cum_Buy_Hold <- c( 0, cumsum(as.numeric(stra_Buy_Hold)))


Date <- lubridate::as_date(zoo::index(crypto[(nrow(crypto) - backtest_period):nrow(crypto), ]))

stra_wide <- matrix(ncol = 5, nrow = (backtest_period + 1)) %>%
  as_tibble() %>%
  transmute(`Buy and Hold` = cum_Buy_Hold,
            TSMOM = cum_TSMOM,
            `cop-GARCH` = cum_cop_GARCH,
            `TSMOM extended` = cum_TSMOM_cop,
            Date = Date)
stra_tidy <- stra_wide %>%
  tidyr::gather(key = "Strategy", value = "Cumulative return", -Date)
#
pdf(file = paste0(outdir, "cumulate_strategy_comparison.pdf"),
    height = 6, width = 8)
ggplot(stra_tidy, aes(x = Date, y = `Cumulative return`, col = Strategy)) +
  geom_line(size = 0.7, alpha = 0.7) +
  theme_classic() +
  labs(title = "Cumulative return over time of different strategies",
       caption = "2017-10-12 until 2018-09-26 (250 days)") +
  ggplot2::scale_color_manual(values = c(RColorBrewer::brewer.pal(4,"Set1"))) +
  theme(plot.title = element_text(hjust = 0.5))
dev.off()

sum_stats <- matrix(c(tseries::sharpe(cum_Buy_Hold, scale = sqrt(350)),
                      tseries::sharpe(cum_TSMOM, scale = sqrt(350)),
                      tseries::sharpe(cum_cop_GARCH, scale = sqrt(350)),
                      tseries::sharpe(cum_TSMOM_cop, scale = sqrt(350)),
                      tail(cum_Buy_Hold, 1),
                      tail(cum_TSMOM, 1),
                      tail(cum_cop_GARCH, 1),
                      tail(cum_TSMOM_cop, 1)), nrow = 2, byrow = TRUE)
colnames(sum_stats) <- c("Buy and Hold", "TSMOM", "cop-GARCH", "TSMOM extended")
rownames(sum_stats) <- c("Sharpe Ratio", "Cumulative return")

xtable::xtable(sum_stats)
3schwartz/SpecialeScrAndFun documentation built on May 4, 2019, 6:29 a.m.