scripts/Copula_Crypto.R
In 3schwartz/SpecialeScrAndFun: Script and functions used for master thesis.
# outdir = "C:/Users/sch/Dropbox/Egne dokumenter/Skole/master/opgave/Figures/Copula/"
outdir = "C:/Users/Soren Schwartz/Dropbox/Egne dokumenter/Skole/master/opgave/Figures/Copula/"
library(xts)
library(copula)
library(qrmdata)
library(qrmtools)
library(Data)
library(scrAndFun)
library(ggplot2)
#### Data manupulation ####
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/"
))
#### 1 Working with the data ####
SX <- returns(stocks.xts[,-4]) # or diff(log(x))[-1,]
IX <- returns(indx.xts)
FX <- returns(cur.xts)
COMX <- returns(com.xts)
ZCB <- ZCB.xts[,c("THREEFY1","THREEFY3", "THREEFY10")]
colnames(ZCB) <- c("Y1", "Y3", "Y10")
ZCBX <- returns(ZCB)
crypto <- do.call(cbind, lapply(dailyExcess, function(x) {x$excessReturn}))[,1:4]
colnames(crypto) <- names(dailyExcess)[1:4]
crypto <- crypto[, c("BitCoin", "Ethereum", "Ripple")]
cryptoX <- crypto[-which(is.na(crypto), arr.ind = TRUE)[,1],]
periods <- lapply(list(SX, IX, FX, COMX, ZCBX, cryptoX), xts::periodicity)
names(periods) <- c("SX", "IX", "FX", "COMX", "ZCBX", "cryptoX"); periods
SX_sub <- SX[paste0(range(zoo::index(cryptoX)), collapse = "/")]
IX_sub <- IX[paste0(range(zoo::index(cryptoX)), collapse = "/")]
FX_sub <- FX[paste0(range(zoo::index(cryptoX)), collapse = "/")]
COMX_sub <- COMX[paste0(range(zoo::index(cryptoX)), collapse = "/")]
ZCBX_sub <- ZCBX[paste0(range(zoo::index(cryptoX)), collapse = "/")]
periods_sub <- lapply(list(SX_sub, IX_sub, FX_sub, COMX_sub, ZCBX_sub, cryptoX), xts::periodicity)
names(periods_sub) <- c("SX_sub", "IX_sub", "FX_sub", "COMX_sub", "ZCBX_sub", "cryptoX"); periods_sub
## Compute pseudo-observations
U_SX <- as.matrix(pobs(SX_sub))
plot(U_SX, xlab = expression(U[1]^"hej"), ylab = expression(U[2]))
### Fitting statis copulas ####
foo <- list(crypto = cryptoX, Stocks = SX_sub, Index = IX_sub,
Forex = FX_sub, Com = COMX_sub, ZCB = ZCBX_sub)
get_log_like <- function(prices, type = c("ex", "un")) {
out <- matrix(ncol = 4, nrow = length(prices))
if (type == "ex") {
rho_matrix <- matrix(ncol = 4, nrow = length(prices))
}
for (i in 1:length(prices)) {
U <- as.matrix((pobs(prices[[i]])))
fit.N <- fitCopula(normalCopula(dim = dim(U)[2], dispstr = type), data = U)
fit.t <- fitCopula(tCopula(dim = dim(U)[2], dispstr = type), data = U) # df of freedom are estimated, too
fit.C <- fitCopula(claytonCopula(dim = dim(U)[2]), data = U)
fit.G <- fitCopula(gumbelCopula(dim = dim(U)[2]), data = U)
out[i, ] <- c(N = fit.N@loglik, t = fit.t@loglik, C = fit.C@loglik, G = fit.G@loglik)
if (type == "ex") {
rho_matrix[i, ] <- c((2/pi) * asin(fit.N@estimate), (2/pi) * asin(fit.t@estimate[1]),
fit.C@estimate / (fit.C@estimate + 2), (fit.G@estimate - 1) / fit.G@estimate)
}
}
colnames(out) <- c("N", "t", "C", "G")
rownames(out) <- names(prices)
if (type == "ex") {
colnames(rho_matrix) <- c("N", "t", "C", "G")
rownames(rho_matrix) <- names(prices)
}
if (type == "ex") {
return(list(
log_like = out,
rho = rho_matrix))
} else {
return(list(
log_like = out))
}
}
system.time(
log_like_est_un <- get_log_like(prices = foo, type = "un")
)
system.time(
log_like_est_ex <- get_log_like(prices = foo, type = "ex")
)
xtable::xtable(log_like_est_ex$log_like)
xtable::xtable(log_like_est_un$log_like)
xtable::xtable(log_like_est_ex$rho)
xtable::xtable(
cbind(log_like_est_ex$log_like, log_like_est_un$log_like[,c("N", "t")]))
get_pairs_plot <- function(prices, outdir) {
for (i in 1:length(prices)) {
pdf(file = paste0(outdir, names(foo)[i], ".pdf"),
height = 6, width = 8)
as.matrix((pobs(foo[[i]]))) %>%
pairs()
dev.off()
}
}
get_pairs_plot(foo, outdir)
#### Timevarying copula ####
crypto <- do.call(cbind, lapply(dailyExcess, function(x) {x$excessReturn}))[,1:4]
colnames(crypto) <- names(dailyExcess)[1:4]
crypto <- crypto[, c("BitCoin", "Ethereum", "Ripple")]
cryptoX <- crypto[-which(is.na(crypto), arr.ind = TRUE)[,1],]
# Index portfolios without NA's
ind_BTC_ETH <- which(is.na(crypto[,c("BitCoin","Ethereum")]), arr.ind = TRUE)[,1]
ind_BTC_XRP <- which(is.na(crypto[,c("BitCoin","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")]
# Pseudo observations
U_BTC_ETH <- as.matrix(pobs(crypto_BTC_ETH))
U_BTC_XRP <- as.matrix(pobs(crypto_BTC_XRP))
# Fit of BTC_ETH
fit_N_St_BTC_ETH <- fitCopula(normalCopula(dim = dim(U_BTC_ETH)[2], dispstr = "un"), data = U_BTC_ETH)
fit_C_St_BTC_ETH <- fitCopula(claytonCopula(dim = dim(U_BTC_ETH)[2]), data = U_BTC_ETH)
fit_Crot_St_BTC_ETH <- fitCopula(claytonCopula(dim = dim(U_BTC_ETH)[2]), data = 1-U_BTC_ETH)
(maxLL_gaus_TV_BTC_ETH <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta = 0.5),
function(par) {
U_foo = U_BTC_ETH
LL_gaus_cpp(par, U_foo, 10)},
lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
(maxLL_clayton_TV_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))))
(maxLL_clayton_surv_TV_BTC_ETH <- stats::nlminb(start = c(w = 0.1, alpha = 0.1, beta = 0.1),
function(par) {
U_foo = U_BTC_ETH
LL_clayton_surv_cpp(par, U_foo)},
lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
# Fit of BTC_XRP
fit_N_St_BTC_XRP <- fitCopula(normalCopula(dim = dim(U_BTC_XRP)[2], dispstr = "un"), data = U_BTC_XRP)
(maxLL_clayton_surv_BTC_XRP <- stats::nlminb(start = c(th = 0.5),
function(th) {
U_foo = 1-U_BTC_XRP
LL_clayton_static_cpp(th, U_foo)},
lower = c(-1), upper = c(Inf)))
(maxLL_clayton_BTC_XRP <- stats::nlminb(start = c(th = 0.5),
function(th) {
U_foo = U_BTC_XRP
LL_clayton_static_cpp(th, U_foo)},
lower = c(-1), upper = c(Inf)))
(maxLL_gaus_TV_BTC_XRP <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta = 0.5),
function(par) {
U_foo = U_BTC_XRP
LL_gaus_cpp(par, U_foo, 10)},
lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
(maxLL_clayton_TV_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))))
(maxLL_clayton_surv_TV_BTC_XRP <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta = 0.5),
function(par) {
U_foo = U_BTC_XRP
LL_clayton_surv_cpp(par, U_foo)},
lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
# Log-likelihood
loglike_TV <- (rbind(c(N = fit_N_St_BTC_ETH@loglik,
C = fit_C_St_BTC_ETH@loglik,
C_rot = fit_Crot_St_BTC_ETH@loglik,
N_TV = -maxLL_gaus_TV_BTC_ETH$objective,
C_TV = -maxLL_clayton_TV_BTC_ETH$objective,
C_TV_rot = -maxLL_clayton_surv_TV_BTC_ETH$objective),
c(N = fit_N_St_BTC_XRP@loglik,
C = -maxLL_clayton_BTC_XRP$objective,
C_rot = -maxLL_clayton_surv_BTC_XRP$objective,
N_TV = -maxLL_gaus_TV_BTC_XRP$objective,
C_TV = -maxLL_clayton_TV_BTC_XRP$objective,
C_TV_rot = -maxLL_clayton_surv_TV_BTC_XRP$objective)))
rownames(loglike_TV) <- c("BTC-ETH", "BTC-XRP")
xtable::xtable(loglike_TV)
# Plot of pseudo-observations
pdf(file = paste0(outdir, "pseudoplot_bivar_cops.pdf"),
width = 8, height = 6)
par(mfrow = c(1,2), mar = c(4,5,2,1) + 0.1)
plot(U_BTC_ETH, xlab = expression(U[1]^"BTC"), ylab = expression(U[2]^"ETH"),
main = "Pseudo-observation of BTC-ETH")
plot(U_BTC_XRP[,1],U_BTC_XRP[,2], xlab = expression(U[1]^"BTC"), ylab = expression(U[2]^"XRP"),
main = "Pseudo-observation of BTC-XRP")
par(mfrow = c(1,1))
dev.off()
# Parameter estimates
parameter_table <- cbind(N_BTC_ETH = maxLL_gaus_TV_BTC_ETH$par,
C_BTC_ETH = maxLL_clayton_TV_BTC_ETH$par,
C_BTC_ETH = maxLL_clayton_surv_TV_BTC_ETH$par,
N_BTC_XRP = maxLL_gaus_TV_BTC_XRP$par,
C_BTC_XRP = maxLL_clayton_TV_BTC_XRP$par,
C_BTC_XRP = maxLL_clayton_surv_TV_BTC_XRP$par)
xtable::xtable(parameter_table, digits = 3)
# Plot of time-varying parameter
par_BTC_ETH <- cbind(Gaussian = LL_gaus_sim_cpp(maxLL_gaus_TV_BTC_ETH$par, U_BTC_ETH, 10),
Clayton = LL_clayton_sim_cpp(maxLL_clayton_TV_BTC_ETH$par, U_BTC_ETH),
`Clayton rotated` = LL_clayton_sim_cpp(maxLL_clayton_surv_TV_BTC_ETH$par, U_BTC_ETH))
par_BTC_XRP <- cbind(Gaussian = LL_gaus_sim_cpp(maxLL_gaus_TV_BTC_XRP$par, U_BTC_XRP, 10),
Clayton = LL_clayton_sim_cpp(maxLL_clayton_TV_BTC_XRP$par, U_BTC_XRP),
`Clayton rotated` = LL_clayton_sim_cpp(maxLL_clayton_surv_TV_BTC_XRP$par, U_BTC_XRP))
param_BTC_ETH <- par_BTC_ETH %>%
tibble::as_tibble() %>%
dplyr::mutate(Date = zoo::index(crypto_BTC_ETH))
rho_BTC_ETH_tidy <- param_BTC_ETH[-c(1:10),] %>%
dplyr::transmute(Gaussian = (2/pi)*asin(Gaussian),
Clayton = Clayton/(Clayton + 2),
`Clayton rotated` = `Clayton rotated`/(`Clayton rotated` + 2),
Date = Date) %>%
tidyr::gather(key = "Model", value = "rho_t", - Date)
param_BTC_XRP <- par_BTC_XRP %>%
tibble::as_tibble() %>%
dplyr::mutate(Date = zoo::index(crypto_BTC_XRP))
rho_BTC_XRP_tidy <- param_BTC_XRP[-c(1:10),] %>%
dplyr::transmute(Gaussian = (2/pi)*asin(Gaussian),
Clayton = Clayton/(Clayton + 2),
`Clayton rotated` = `Clayton rotated`/(`Clayton rotated` + 2),
Date = Date) %>%
tidyr::gather(key = "Model", value = "rho_t", - Date)
#
pdf(file = paste0(outdir, "btc_ETH_dyn_tau.pdf"),
width = 8, height = 6)
ggplot(rho_BTC_ETH_tidy, aes(x = Date, y = rho_t, col = Model)) +
geom_line(size = 0.7) +
theme_classic() +
labs(title = "Dynamics of time-varying parameter",
subtitle = "BTC-ETH",
caption = expression(hat(rho)[t]~" transformed to Kendall's "~tau),
y = expression(hat(rho)[t]^tau)) +
ggplot2::scale_color_manual(values = c(RColorBrewer::brewer.pal(3,"Set1"))) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.caption = element_text(size=12))
dev.off()
pdf(file = paste0(outdir, "btc_XRP_dyn_tau.pdf"),
width = 8, height = 6)
ggplot(rho_BTC_XRP_tidy, aes(x = Date, y = rho_t, col = Model)) +
geom_line(size = 0.7) +
theme_classic() +
labs(title = "Dynamics of time-varying parameter",
subtitle = "BTC-XRP",
caption = expression(hat(rho)[t]~" transformed to Kendall's "~tau),
y = expression(hat(rho)[t]^tau)) +
ggplot2::scale_color_manual(values = c(RColorBrewer::brewer.pal(3,"Set1"))) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.caption = element_text(size=12))
dev.off()
#### Copula estimation ####
library(rugarch)
library(xts)
library(ADGofTest)
library(qqtest)
library(copula)
library(qrmtools)
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 <- which(is.na(crypto[,c("BitCoin","Ethereum")]), arr.ind = TRUE)[,1]
ind_BTC_XRP <- which(is.na(crypto[,c("BitCoin","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")]
# Pseudo observations
U_BTC_ETH <- as.matrix(pobs(crypto_BTC_ETH))
U_BTC_XRP <- as.matrix(pobs(crypto_BTC_XRP))
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/"]
## Fit ARMA
fit_ARMA_BTC_ETH_1 <- forecast::auto.arima(BTC_ETH_train[,1], d = 0)
fit_ARMA_BTC_ETH_2 <- forecast::auto.arima(BTC_ETH_train[,2], d = 0)
aic_arma_1 <- aic_arma_2 <- Inf
a1_order <- a2_order <- c(0,0)
for (p in 0:4) for (q in 0:4) {
fit_1 <- arima(BTC_ETH_train[,1], order = c(p, 0, q))
fit_2 <- arima(BTC_ETH_train[,2], order = c(p, 0, q))
if(fit_1$aic < aic_arma_1) {
aic_arma_1 <- fit_1$aic
a1_order <- c(p, q)
}
if(fit_2$aic < aic_arma_2) {
aic_arma_2 <- fit_2$aic
a2_order <- c(p, q)
}
}
aic_arma_1; aic_arma_2;
a1_order; a2_order;
## Fit GARCH
aic1 <- aic2 <- Inf
g1_order <- g2_order <- c(0,0)
for (P in 0:3) for (Q in 0:3) {
uspec1 <- ugarchspec(variance.model = list(garchOrder = c(P, Q)),
mean.model = list(armaOrder = a1_order),#forecast::arimaorder(fit_ARMA_BTC_ETH_1)[c(1,3)]),
distribution.model = "sstd")
uspec2 <- ugarchspec(variance.model = list(garchOrder = c(P, Q)),
mean.model = list(armaOrder = a2_order),#forecast::arimaorder(fit_ARMA_BTC_ETH_2)[c(1,3)]),
distribution.model = "sstd")
fit_both <- fit_ARMA_GARCH(BTC_ETH_train, ugarchspec.list = list(uspec1, uspec2))
result_1 <- tryCatch({
current1 <- rugarch::infocriteria(fit_both$fit[[1]])[1]
if(current1 < aic1) {
aic1 <- current1; g1_order <- c(P,Q)
}
}, error = function(cond){
message(cond)
})
result_2 <- tryCatch({
current2 <- rugarch::infocriteria(fit_both$fit[[2]])[1]
if(current2 < aic2) {
aic2 <- current2; g2_order <- c(P,Q)
}
}, error = function(cond){
message(cond)
})
}
aic1 ; aic2 ;
g1_order ; g2_order ;
## Fit final
fit_final_BTC_ETH = fit_ARMA_GARCH(BTC_ETH_train, ugarchspec.list = list(
ugarchspec(variance.model = list(garchOrder = g1_order),
mean.model = list(armaOrder = a1_order),#forecast::arimaorder(fit_ARMA_BTC_ETH_1)[c(1,3)]),
distribution.model = "sstd"),
ugarchspec(variance.model = list(garchOrder = g2_order),
mean.model = list(armaOrder = a2_order),#forecast::arimaorder(fit_ARMA_BTC_ETH_2)[c(1,3)]),
distribution.model = "sstd")))
Z_BTC_ETH <- as.matrix(do.call(merge, lapply(fit_final_BTC_ETH$fit, rugarch::residuals, standardize = TRUE)))
colnames(Z_BTC_ETH) <- colnames(BTC_ETH_train)
# Extract coefficients of skewed t and getprobability transformations
(nu_mar_BTC_ETH <- vapply(fit_final_BTC_ETH$fit, function(x) { x@fit$coef[["shape"]]}, NA_real_))
(gamma_mar_BTC_ETH <- vapply(fit_final_BTC_ETH$fit, function(x) { x@fit$coef[["skew"]]}, NA_real_))
n_BTC_ETH <- nrow(Z_BTC_ETH); d_BTC_ETH <- ncol(Z_BTC_ETH); U_BTC_ETH <- pobs(Z_BTC_ETH)
# Test if uniform margins
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)
# Fit copula
(maxLL_gaus_TV_BTC_ETH <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta = 0.5),
function(par) {
U_foo = U_BTC_ETH
LL_gaus_cpp(par, U_foo, 10)},
lower = c(rep(-1,3)), upper = c(rep(1, 3))))
(maxLL_clayton_TV_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(-1,3)), upper = c(rep(1, 3))))
(maxLL_clayton_surv_TV_BTC_ETH <- stats::nlminb(start = c(w = 0, alpha = 0, beta = 0),
function(par) {
U_foo = U_BTC_ETH
LL_clayton_surv_cpp(par, U_foo)},
lower = c(rep(-1,3)), upper = c(rep(1, 3))))
# Fit of BTC_XRP
(maxLL_gaus_TV_BTC_XRP <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta = 0.5),
function(par) {
U_foo = U_BTC_XRP
LL_gaus_cpp(par, U_foo, 10)},
lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
(maxLL_clayton_TV_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))))
(maxLL_clayton_surv_TV_BTC_XRP <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta = 0.5),
function(par) {
U_foo = U_BTC_XRP
LL_clayton_surv_cpp(par, U_foo)},
lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
# Parameter estimates
parameter_BTC_ETH <- cbind(N_BTC_ETH = maxLL_gaus_TV_BTC_ETH$par,
C_BTC_ETH = maxLL_clayton_TV_BTC_ETH$par,
C_BTC_ETH = maxLL_clayton_surv_TV_BTC_ETH$par)
parameter_BTC_XRP <- cbind(N_BTC_XRP = maxLL_gaus_TV_BTC_XRP$par,
C_BTC_XRP = maxLL_clayton_TV_BTC_XRP$par,
C_BTC_XRP = maxLL_clayton_surv_TV_BTC_XRP$par)
# Simulate
test = BTC_ETH_test; parameters = parameter_BTC_ETH[,2];
fits = fit_final_BTC_ETH$fit; m = 1000; train = U_BTC_ETH;
nu_mar = nu_mar_BTC_ETH; gamma_mar = gamma_mar;
alpha = 0.99
# 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/"]
d = ncol(test)
Z_BTC_ETH
U_sim <- as.matrix(rbind(U_BTC_ETH,
matrix(0, ncol = d, nrow = nrow(BTC_ETH_test))))
Z_foo <- as.matrix(rbind(Z_BTC_ETH,
matrix(0, ncol = d, nrow = nrow(BTC_ETH_test))))
CI_out <- matrix(ncol = d, nrow = nrow(test))
future <- VaR_out <- numeric(nrow(test))
data_all <- as.matrix(rbind(BTC_ETH_train, BTC_ETH_test))
u_1 <- c(as.numeric(fitted(fits[[1]])), numeric(length(test)))
u_2 <- c(as.numeric(fitted(fits[[1]])), numeric(length(test)))
sigma_1 <- c(as.numeric(sigma(fits[[1]])), numeric(length(test)))
sigma_2 <- c(as.numeric(sigma(fits[[2]])), numeric(length(test)))
dependece <- c(LL_clayton_sim_cpp(parameters, as.matrix(train)),
numeric(length(test)))
mu_1 <- coef(fits[[1]])[[1]]; omega_1 <- coef(fits[[1]])[[2]];
alpha_1 <- coef(fits[[1]])[[4]]; beta_1 <- coef(fits[[1]])[[4]];
mu_2 <- coef(fits[[2]])[[1]]; ar1_2 <- coef(fits[[2]])[[2]];
ar2_2 <- coef(fits[[2]])[[3]]; ar3_2 <- coef(fits[[2]])[[4]];
ar4_2 <- coef(fits[[2]])[[5]]; ma1_2 <- coef(fits[[2]])[[6]];
ma2_2 <- coef(fits[[2]])[[7]]; ma3_2 <- coef(fits[[2]])[[8]];
omega_2 <- coef(fits[[2]])[[9]]; alpha_2 <- coef(fits[[2]])[[10]];
beta_2 <- coef(fits[[2]])[[11]];
end = nrow(train) + nrow(test)
start = nrow(train) + 1
for (i in start:end) {
U. <- rCopula(m, copula = claytonCopula(param = 1000, dim = d))
U.[U.[,2] == 0,2] = U.[U.[,2] == 0,1]
Z. <- sapply(1:d, function(i) {skewt::qskt(U.[,i], df = nu_mar[i],
gamma = gamma_mar[i])})
u_1[i] <- mu_1
sigma_1[i] <- omega_1 + alpha_1*(data_all[i-1,1] - u_1[i-1])^2 + beta_1*sigma_1[i-1]
u_2[i] <- mu_2 + ar1_2 * (data_all[i-1,2] - mu_2) + ar2_2 * (data_all[i-2,2] - mu_2) +
ar3_2 * (data_all[i-3,2] - mu_2) + ar4_2 * (data_all[i-4,2] - mu_2) +
ma1_2 * (data_all[i-1,2] - u_2[i-1]) + ma2_2 * (data_all[i-2,2] - u_2[i-2]) +
ma3_2 * (data_all[i-3,2] - u_2[i-3])
sigma_2[i] <- omega_2 + alpha_2*(data_all[i-1,1] - u_2[i-1])^2 + beta_2*sigma_2[i-1]
X_sim = cbind(u_1[i] + sigma_1[i] * Z.[,1],
u_2[i] + sigma_2[i] * Z.[,2])
Xs. <- rowSums(X_sim)
Xs.mean <- mean(Xs.)
Xs.CI <- quantile(Xs., probs = c(0.025, 0.975))
VaR <- quantile(Xs., probs = alpha)
Z_foo[i,] <- (data_all[i, ] - c(u_1[i], u_2[i])) / c(sigma_1[i], sigma_2[i])
U_sim[i, ] <- tail(pobs(Z_foo[1:i,]),1)
dependece[i] <- clayton_theta_cpp(parameters, 1e-10, U_sim[i-1,])
# dependece[i] <- ifelse(dependece[i] < 1000,
# ifelse(dependece[i] != 0, dependece[i], dependece[i-1]),
# dependece[i-1])
future[i-nrow(train)] <- Xs.mean
CI_out[i-nrow(train),] <- Xs.CI
VaR_out[i-nrow(train)] <- VaR
}
backtest_cop_t <- function(return = X_all, backtest_period = 200,
alpha = 0.99, m = 200) {
future <- matrix(ncol = 1, nrow = backtest_period)
CI_out <- matrix(ncol = 2, nrow = backtest_period)
VaR_out <- matrix(ncol = 1, nrow = backtest_period)
for (i in 0:(backtest_period - 1)) {
nY <- nrow(return) - (backtest_period - i)
X <- return[1:nY, ]
## Fit marginal time series
uspec <- rep(list(ugarchspec(distribution.model = "std")), ncol(X))
fit.ARMA.GARCH <- fit_ARMA_GARCH(X, ugarchspec.list = uspec)
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_) # vector of estimated df
d <- ncol(X) # dimension
### Fitting copulas ##
## Compute pseudo-observations from the standardized t residuals
U <- pobs(Z)
## Fitting a t copula
fit.tc <- fitCopula(tCopula(dim = d, dispstr = "un"),
data = U)#, method = "itau.mpl")
tc <- fit.tc@copula # fitted copula
### Simulating paths from the full model ##
## Simulate paths
m <- m
U. <- rCopula(m, copula = tc)
Z. <- sapply(1:d, function(j) sqrt((nu.mar[j] - 2) / nu.mar[j]) * qt(U.[,j], df = nu.mar[j]))
sim <- lapply(1:d, function(j)
ugarchsim(fits[[j]], n.sim = 1, m.sim = m,
custom.dist = list(name = "sample", distfit = t(Z.[,j, drop = FALSE]))))
fitted_sim <- sapply(sim, function(x) fitted(x))
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
}
return(list(
forecast = future,
CI = CI_out,
VaR = VaR_out
))
}
backtest_period <- 50
back_out <- backtest_cop_t(X_all, backtest_period = backtest_period)
nX = nrow(X_all) - backtest_period
X <- X_all[1:nX, ]
X_future <- X_all[(nX+1):nrow(X_all), ]
Xs <- rowSums(X) # aggregated return; n-vector
Xs.mean <- back_out$forecast # predicted aggregated return; m-vector
Xs.CI <- t(back_out$CI)# CIs; (2, m)-matrix
VaR <- back_out$VaR# VaR_alpha; m-vector
## Plot
past <- zoo::index(X)
future <- zoo::index(X_future)
plot(past, Xs, type = "l", xlim = range(c(past, future)), xlab = "", ylab = "",
ylim = c(-(max(range(VaR))+0.1), (max(range(VaR))+0.1))) # actual (past) returnes
polygon(c(future, rev(future)), c(Xs.CI[1,], rev(Xs.CI[2,])),
border = NA, col = "grey80") # CI region
lines(future, rowSums(X_future), col = "orange")
lines(future, Xs.mean, col = "royalblue3") # predicted aggregated return
lines(future, Xs.CI[1,], col = "grey50") # lower CI
lines(future, Xs.CI[2,], col = "grey50") # upper CI
lines(future, VaR, col = "maroon3") # 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)))))
mse = mean((rowSums(X_future) - Xs.mean)^2)
3schwartz/SpecialeScrAndFun documentation built on May 4, 2019, 6:29 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# outdir = "C:/Users/sch/Dropbox/Egne dokumenter/Skole/master/opgave/Figures/Copula/"
outdir = "C:/Users/Soren Schwartz/Dropbox/Egne dokumenter/Skole/master/opgave/Figures/Copula/"
library(xts)
library(copula)
library(qrmdata)
library(qrmtools)
library(Data)
library(scrAndFun)
library(ggplot2)
#### Data manupulation ####
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/"
))
#### 1 Working with the data ####
SX <- returns(stocks.xts[,-4]) # or diff(log(x))[-1,]
IX <- returns(indx.xts)
FX <- returns(cur.xts)
COMX <- returns(com.xts)
ZCB <- ZCB.xts[,c("THREEFY1","THREEFY3", "THREEFY10")]
colnames(ZCB) <- c("Y1", "Y3", "Y10")
ZCBX <- returns(ZCB)
crypto <- do.call(cbind, lapply(dailyExcess, function(x) {x$excessReturn}))[,1:4]
colnames(crypto) <- names(dailyExcess)[1:4]
crypto <- crypto[, c("BitCoin", "Ethereum", "Ripple")]
cryptoX <- crypto[-which(is.na(crypto), arr.ind = TRUE)[,1],]
periods <- lapply(list(SX, IX, FX, COMX, ZCBX, cryptoX), xts::periodicity)
names(periods) <- c("SX", "IX", "FX", "COMX", "ZCBX", "cryptoX"); periods
SX_sub <- SX[paste0(range(zoo::index(cryptoX)), collapse = "/")]
IX_sub <- IX[paste0(range(zoo::index(cryptoX)), collapse = "/")]
FX_sub <- FX[paste0(range(zoo::index(cryptoX)), collapse = "/")]
COMX_sub <- COMX[paste0(range(zoo::index(cryptoX)), collapse = "/")]
ZCBX_sub <- ZCBX[paste0(range(zoo::index(cryptoX)), collapse = "/")]
periods_sub <- lapply(list(SX_sub, IX_sub, FX_sub, COMX_sub, ZCBX_sub, cryptoX), xts::periodicity)
names(periods_sub) <- c("SX_sub", "IX_sub", "FX_sub", "COMX_sub", "ZCBX_sub", "cryptoX"); periods_sub
## Compute pseudo-observations
U_SX <- as.matrix(pobs(SX_sub))
plot(U_SX, xlab = expression(U[1]^"hej"), ylab = expression(U[2]))
### Fitting statis copulas ####
foo <- list(crypto = cryptoX, Stocks = SX_sub, Index = IX_sub,
Forex = FX_sub, Com = COMX_sub, ZCB = ZCBX_sub)
get_log_like <- function(prices, type = c("ex", "un")) {
out <- matrix(ncol = 4, nrow = length(prices))
if (type == "ex") {
rho_matrix <- matrix(ncol = 4, nrow = length(prices))
}
for (i in 1:length(prices)) {
U <- as.matrix((pobs(prices[[i]])))
fit.N <- fitCopula(normalCopula(dim = dim(U)[2], dispstr = type), data = U)
fit.t <- fitCopula(tCopula(dim = dim(U)[2], dispstr = type), data = U) # df of freedom are estimated, too
fit.C <- fitCopula(claytonCopula(dim = dim(U)[2]), data = U)
fit.G <- fitCopula(gumbelCopula(dim = dim(U)[2]), data = U)
out[i, ] <- c(N = fit.N@loglik, t = fit.t@loglik, C = fit.C@loglik, G = fit.G@loglik)
if (type == "ex") {
rho_matrix[i, ] <- c((2/pi) * asin(fit.N@estimate), (2/pi) * asin(fit.t@estimate[1]),
fit.C@estimate / (fit.C@estimate + 2), (fit.G@estimate - 1) / fit.G@estimate)
}
}
colnames(out) <- c("N", "t", "C", "G")
rownames(out) <- names(prices)
if (type == "ex") {
colnames(rho_matrix) <- c("N", "t", "C", "G")
rownames(rho_matrix) <- names(prices)
}
if (type == "ex") {
return(list(
log_like = out,
rho = rho_matrix))
} else {
return(list(
log_like = out))
}
}
system.time(
log_like_est_un <- get_log_like(prices = foo, type = "un")
)
system.time(
log_like_est_ex <- get_log_like(prices = foo, type = "ex")
)
xtable::xtable(log_like_est_ex$log_like)
xtable::xtable(log_like_est_un$log_like)
xtable::xtable(log_like_est_ex$rho)
xtable::xtable(
cbind(log_like_est_ex$log_like, log_like_est_un$log_like[,c("N", "t")]))
get_pairs_plot <- function(prices, outdir) {
for (i in 1:length(prices)) {
pdf(file = paste0(outdir, names(foo)[i], ".pdf"),
height = 6, width = 8)
as.matrix((pobs(foo[[i]]))) %>%
pairs()
dev.off()
}
}
get_pairs_plot(foo, outdir)
#### Timevarying copula ####
crypto <- do.call(cbind, lapply(dailyExcess, function(x) {x$excessReturn}))[,1:4]
colnames(crypto) <- names(dailyExcess)[1:4]
crypto <- crypto[, c("BitCoin", "Ethereum", "Ripple")]
cryptoX <- crypto[-which(is.na(crypto), arr.ind = TRUE)[,1],]
# Index portfolios without NA's
ind_BTC_ETH <- which(is.na(crypto[,c("BitCoin","Ethereum")]), arr.ind = TRUE)[,1]
ind_BTC_XRP <- which(is.na(crypto[,c("BitCoin","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")]
# Pseudo observations
U_BTC_ETH <- as.matrix(pobs(crypto_BTC_ETH))
U_BTC_XRP <- as.matrix(pobs(crypto_BTC_XRP))
# Fit of BTC_ETH
fit_N_St_BTC_ETH <- fitCopula(normalCopula(dim = dim(U_BTC_ETH)[2], dispstr = "un"), data = U_BTC_ETH)
fit_C_St_BTC_ETH <- fitCopula(claytonCopula(dim = dim(U_BTC_ETH)[2]), data = U_BTC_ETH)
fit_Crot_St_BTC_ETH <- fitCopula(claytonCopula(dim = dim(U_BTC_ETH)[2]), data = 1-U_BTC_ETH)
(maxLL_gaus_TV_BTC_ETH <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta = 0.5),
function(par) {
U_foo = U_BTC_ETH
LL_gaus_cpp(par, U_foo, 10)},
lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
(maxLL_clayton_TV_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))))
(maxLL_clayton_surv_TV_BTC_ETH <- stats::nlminb(start = c(w = 0.1, alpha = 0.1, beta = 0.1),
function(par) {
U_foo = U_BTC_ETH
LL_clayton_surv_cpp(par, U_foo)},
lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
# Fit of BTC_XRP
fit_N_St_BTC_XRP <- fitCopula(normalCopula(dim = dim(U_BTC_XRP)[2], dispstr = "un"), data = U_BTC_XRP)
(maxLL_clayton_surv_BTC_XRP <- stats::nlminb(start = c(th = 0.5),
function(th) {
U_foo = 1-U_BTC_XRP
LL_clayton_static_cpp(th, U_foo)},
lower = c(-1), upper = c(Inf)))
(maxLL_clayton_BTC_XRP <- stats::nlminb(start = c(th = 0.5),
function(th) {
U_foo = U_BTC_XRP
LL_clayton_static_cpp(th, U_foo)},
lower = c(-1), upper = c(Inf)))
(maxLL_gaus_TV_BTC_XRP <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta = 0.5),
function(par) {
U_foo = U_BTC_XRP
LL_gaus_cpp(par, U_foo, 10)},
lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
(maxLL_clayton_TV_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))))
(maxLL_clayton_surv_TV_BTC_XRP <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta = 0.5),
function(par) {
U_foo = U_BTC_XRP
LL_clayton_surv_cpp(par, U_foo)},
lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
# Log-likelihood
loglike_TV <- (rbind(c(N = fit_N_St_BTC_ETH@loglik,
C = fit_C_St_BTC_ETH@loglik,
C_rot = fit_Crot_St_BTC_ETH@loglik,
N_TV = -maxLL_gaus_TV_BTC_ETH$objective,
C_TV = -maxLL_clayton_TV_BTC_ETH$objective,
C_TV_rot = -maxLL_clayton_surv_TV_BTC_ETH$objective),
c(N = fit_N_St_BTC_XRP@loglik,
C = -maxLL_clayton_BTC_XRP$objective,
C_rot = -maxLL_clayton_surv_BTC_XRP$objective,
N_TV = -maxLL_gaus_TV_BTC_XRP$objective,
C_TV = -maxLL_clayton_TV_BTC_XRP$objective,
C_TV_rot = -maxLL_clayton_surv_TV_BTC_XRP$objective)))
rownames(loglike_TV) <- c("BTC-ETH", "BTC-XRP")
xtable::xtable(loglike_TV)
# Plot of pseudo-observations
pdf(file = paste0(outdir, "pseudoplot_bivar_cops.pdf"),
width = 8, height = 6)
par(mfrow = c(1,2), mar = c(4,5,2,1) + 0.1)
plot(U_BTC_ETH, xlab = expression(U[1]^"BTC"), ylab = expression(U[2]^"ETH"),
main = "Pseudo-observation of BTC-ETH")
plot(U_BTC_XRP[,1],U_BTC_XRP[,2], xlab = expression(U[1]^"BTC"), ylab = expression(U[2]^"XRP"),
main = "Pseudo-observation of BTC-XRP")
par(mfrow = c(1,1))
dev.off()
# Parameter estimates
parameter_table <- cbind(N_BTC_ETH = maxLL_gaus_TV_BTC_ETH$par,
C_BTC_ETH = maxLL_clayton_TV_BTC_ETH$par,
C_BTC_ETH = maxLL_clayton_surv_TV_BTC_ETH$par,
N_BTC_XRP = maxLL_gaus_TV_BTC_XRP$par,
C_BTC_XRP = maxLL_clayton_TV_BTC_XRP$par,
C_BTC_XRP = maxLL_clayton_surv_TV_BTC_XRP$par)
xtable::xtable(parameter_table, digits = 3)
# Plot of time-varying parameter
par_BTC_ETH <- cbind(Gaussian = LL_gaus_sim_cpp(maxLL_gaus_TV_BTC_ETH$par, U_BTC_ETH, 10),
Clayton = LL_clayton_sim_cpp(maxLL_clayton_TV_BTC_ETH$par, U_BTC_ETH),
`Clayton rotated` = LL_clayton_sim_cpp(maxLL_clayton_surv_TV_BTC_ETH$par, U_BTC_ETH))
par_BTC_XRP <- cbind(Gaussian = LL_gaus_sim_cpp(maxLL_gaus_TV_BTC_XRP$par, U_BTC_XRP, 10),
Clayton = LL_clayton_sim_cpp(maxLL_clayton_TV_BTC_XRP$par, U_BTC_XRP),
`Clayton rotated` = LL_clayton_sim_cpp(maxLL_clayton_surv_TV_BTC_XRP$par, U_BTC_XRP))
param_BTC_ETH <- par_BTC_ETH %>%
tibble::as_tibble() %>%
dplyr::mutate(Date = zoo::index(crypto_BTC_ETH))
rho_BTC_ETH_tidy <- param_BTC_ETH[-c(1:10),] %>%
dplyr::transmute(Gaussian = (2/pi)*asin(Gaussian),
Clayton = Clayton/(Clayton + 2),
`Clayton rotated` = `Clayton rotated`/(`Clayton rotated` + 2),
Date = Date) %>%
tidyr::gather(key = "Model", value = "rho_t", - Date)
param_BTC_XRP <- par_BTC_XRP %>%
tibble::as_tibble() %>%
dplyr::mutate(Date = zoo::index(crypto_BTC_XRP))
rho_BTC_XRP_tidy <- param_BTC_XRP[-c(1:10),] %>%
dplyr::transmute(Gaussian = (2/pi)*asin(Gaussian),
Clayton = Clayton/(Clayton + 2),
`Clayton rotated` = `Clayton rotated`/(`Clayton rotated` + 2),
Date = Date) %>%
tidyr::gather(key = "Model", value = "rho_t", - Date)
#
pdf(file = paste0(outdir, "btc_ETH_dyn_tau.pdf"),
width = 8, height = 6)
ggplot(rho_BTC_ETH_tidy, aes(x = Date, y = rho_t, col = Model)) +
geom_line(size = 0.7) +
theme_classic() +
labs(title = "Dynamics of time-varying parameter",
subtitle = "BTC-ETH",
caption = expression(hat(rho)[t]~" transformed to Kendall's "~tau),
y = expression(hat(rho)[t]^tau)) +
ggplot2::scale_color_manual(values = c(RColorBrewer::brewer.pal(3,"Set1"))) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.caption = element_text(size=12))
dev.off()
pdf(file = paste0(outdir, "btc_XRP_dyn_tau.pdf"),
width = 8, height = 6)
ggplot(rho_BTC_XRP_tidy, aes(x = Date, y = rho_t, col = Model)) +
geom_line(size = 0.7) +
theme_classic() +
labs(title = "Dynamics of time-varying parameter",
subtitle = "BTC-XRP",
caption = expression(hat(rho)[t]~" transformed to Kendall's "~tau),
y = expression(hat(rho)[t]^tau)) +
ggplot2::scale_color_manual(values = c(RColorBrewer::brewer.pal(3,"Set1"))) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.caption = element_text(size=12))
dev.off()
#### Copula estimation ####
library(rugarch)
library(xts)
library(ADGofTest)
library(qqtest)
library(copula)
library(qrmtools)
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 <- which(is.na(crypto[,c("BitCoin","Ethereum")]), arr.ind = TRUE)[,1]
ind_BTC_XRP <- which(is.na(crypto[,c("BitCoin","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")]
# Pseudo observations
U_BTC_ETH <- as.matrix(pobs(crypto_BTC_ETH))
U_BTC_XRP <- as.matrix(pobs(crypto_BTC_XRP))
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/"]
## Fit ARMA
fit_ARMA_BTC_ETH_1 <- forecast::auto.arima(BTC_ETH_train[,1], d = 0)
fit_ARMA_BTC_ETH_2 <- forecast::auto.arima(BTC_ETH_train[,2], d = 0)
aic_arma_1 <- aic_arma_2 <- Inf
a1_order <- a2_order <- c(0,0)
for (p in 0:4) for (q in 0:4) {
fit_1 <- arima(BTC_ETH_train[,1], order = c(p, 0, q))
fit_2 <- arima(BTC_ETH_train[,2], order = c(p, 0, q))
if(fit_1$aic < aic_arma_1) {
aic_arma_1 <- fit_1$aic
a1_order <- c(p, q)
}
if(fit_2$aic < aic_arma_2) {
aic_arma_2 <- fit_2$aic
a2_order <- c(p, q)
}
}
aic_arma_1; aic_arma_2;
a1_order; a2_order;
## Fit GARCH
aic1 <- aic2 <- Inf
g1_order <- g2_order <- c(0,0)
for (P in 0:3) for (Q in 0:3) {
uspec1 <- ugarchspec(variance.model = list(garchOrder = c(P, Q)),
mean.model = list(armaOrder = a1_order),#forecast::arimaorder(fit_ARMA_BTC_ETH_1)[c(1,3)]),
distribution.model = "sstd")
uspec2 <- ugarchspec(variance.model = list(garchOrder = c(P, Q)),
mean.model = list(armaOrder = a2_order),#forecast::arimaorder(fit_ARMA_BTC_ETH_2)[c(1,3)]),
distribution.model = "sstd")
fit_both <- fit_ARMA_GARCH(BTC_ETH_train, ugarchspec.list = list(uspec1, uspec2))
result_1 <- tryCatch({
current1 <- rugarch::infocriteria(fit_both$fit[[1]])[1]
if(current1 < aic1) {
aic1 <- current1; g1_order <- c(P,Q)
}
}, error = function(cond){
message(cond)
})
result_2 <- tryCatch({
current2 <- rugarch::infocriteria(fit_both$fit[[2]])[1]
if(current2 < aic2) {
aic2 <- current2; g2_order <- c(P,Q)
}
}, error = function(cond){
message(cond)
})
}
aic1 ; aic2 ;
g1_order ; g2_order ;
## Fit final
fit_final_BTC_ETH = fit_ARMA_GARCH(BTC_ETH_train, ugarchspec.list = list(
ugarchspec(variance.model = list(garchOrder = g1_order),
mean.model = list(armaOrder = a1_order),#forecast::arimaorder(fit_ARMA_BTC_ETH_1)[c(1,3)]),
distribution.model = "sstd"),
ugarchspec(variance.model = list(garchOrder = g2_order),
mean.model = list(armaOrder = a2_order),#forecast::arimaorder(fit_ARMA_BTC_ETH_2)[c(1,3)]),
distribution.model = "sstd")))
Z_BTC_ETH <- as.matrix(do.call(merge, lapply(fit_final_BTC_ETH$fit, rugarch::residuals, standardize = TRUE)))
colnames(Z_BTC_ETH) <- colnames(BTC_ETH_train)
# Extract coefficients of skewed t and getprobability transformations
(nu_mar_BTC_ETH <- vapply(fit_final_BTC_ETH$fit, function(x) { x@fit$coef[["shape"]]}, NA_real_))
(gamma_mar_BTC_ETH <- vapply(fit_final_BTC_ETH$fit, function(x) { x@fit$coef[["skew"]]}, NA_real_))
n_BTC_ETH <- nrow(Z_BTC_ETH); d_BTC_ETH <- ncol(Z_BTC_ETH); U_BTC_ETH <- pobs(Z_BTC_ETH)
# Test if uniform margins
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)
# Fit copula
(maxLL_gaus_TV_BTC_ETH <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta = 0.5),
function(par) {
U_foo = U_BTC_ETH
LL_gaus_cpp(par, U_foo, 10)},
lower = c(rep(-1,3)), upper = c(rep(1, 3))))
(maxLL_clayton_TV_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(-1,3)), upper = c(rep(1, 3))))
(maxLL_clayton_surv_TV_BTC_ETH <- stats::nlminb(start = c(w = 0, alpha = 0, beta = 0),
function(par) {
U_foo = U_BTC_ETH
LL_clayton_surv_cpp(par, U_foo)},
lower = c(rep(-1,3)), upper = c(rep(1, 3))))
# Fit of BTC_XRP
(maxLL_gaus_TV_BTC_XRP <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta = 0.5),
function(par) {
U_foo = U_BTC_XRP
LL_gaus_cpp(par, U_foo, 10)},
lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
(maxLL_clayton_TV_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))))
(maxLL_clayton_surv_TV_BTC_XRP <- stats::nlminb(start = c(w = 0.5, alpha = 0.5, beta = 0.5),
function(par) {
U_foo = U_BTC_XRP
LL_clayton_surv_cpp(par, U_foo)},
lower = c(rep(-Inf,3)), upper = c(rep(Inf, 3))))
# Parameter estimates
parameter_BTC_ETH <- cbind(N_BTC_ETH = maxLL_gaus_TV_BTC_ETH$par,
C_BTC_ETH = maxLL_clayton_TV_BTC_ETH$par,
C_BTC_ETH = maxLL_clayton_surv_TV_BTC_ETH$par)
parameter_BTC_XRP <- cbind(N_BTC_XRP = maxLL_gaus_TV_BTC_XRP$par,
C_BTC_XRP = maxLL_clayton_TV_BTC_XRP$par,
C_BTC_XRP = maxLL_clayton_surv_TV_BTC_XRP$par)
# Simulate
test = BTC_ETH_test; parameters = parameter_BTC_ETH[,2];
fits = fit_final_BTC_ETH$fit; m = 1000; train = U_BTC_ETH;
nu_mar = nu_mar_BTC_ETH; gamma_mar = gamma_mar;
alpha = 0.99
# 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/"]
d = ncol(test)
Z_BTC_ETH
U_sim <- as.matrix(rbind(U_BTC_ETH,
matrix(0, ncol = d, nrow = nrow(BTC_ETH_test))))
Z_foo <- as.matrix(rbind(Z_BTC_ETH,
matrix(0, ncol = d, nrow = nrow(BTC_ETH_test))))
CI_out <- matrix(ncol = d, nrow = nrow(test))
future <- VaR_out <- numeric(nrow(test))
data_all <- as.matrix(rbind(BTC_ETH_train, BTC_ETH_test))
u_1 <- c(as.numeric(fitted(fits[[1]])), numeric(length(test)))
u_2 <- c(as.numeric(fitted(fits[[1]])), numeric(length(test)))
sigma_1 <- c(as.numeric(sigma(fits[[1]])), numeric(length(test)))
sigma_2 <- c(as.numeric(sigma(fits[[2]])), numeric(length(test)))
dependece <- c(LL_clayton_sim_cpp(parameters, as.matrix(train)),
numeric(length(test)))
mu_1 <- coef(fits[[1]])[[1]]; omega_1 <- coef(fits[[1]])[[2]];
alpha_1 <- coef(fits[[1]])[[4]]; beta_1 <- coef(fits[[1]])[[4]];
mu_2 <- coef(fits[[2]])[[1]]; ar1_2 <- coef(fits[[2]])[[2]];
ar2_2 <- coef(fits[[2]])[[3]]; ar3_2 <- coef(fits[[2]])[[4]];
ar4_2 <- coef(fits[[2]])[[5]]; ma1_2 <- coef(fits[[2]])[[6]];
ma2_2 <- coef(fits[[2]])[[7]]; ma3_2 <- coef(fits[[2]])[[8]];
omega_2 <- coef(fits[[2]])[[9]]; alpha_2 <- coef(fits[[2]])[[10]];
beta_2 <- coef(fits[[2]])[[11]];
end = nrow(train) + nrow(test)
start = nrow(train) + 1
for (i in start:end) {
U. <- rCopula(m, copula = claytonCopula(param = 1000, dim = d))
U.[U.[,2] == 0,2] = U.[U.[,2] == 0,1]
Z. <- sapply(1:d, function(i) {skewt::qskt(U.[,i], df = nu_mar[i],
gamma = gamma_mar[i])})
u_1[i] <- mu_1
sigma_1[i] <- omega_1 + alpha_1*(data_all[i-1,1] - u_1[i-1])^2 + beta_1*sigma_1[i-1]
u_2[i] <- mu_2 + ar1_2 * (data_all[i-1,2] - mu_2) + ar2_2 * (data_all[i-2,2] - mu_2) +
ar3_2 * (data_all[i-3,2] - mu_2) + ar4_2 * (data_all[i-4,2] - mu_2) +
ma1_2 * (data_all[i-1,2] - u_2[i-1]) + ma2_2 * (data_all[i-2,2] - u_2[i-2]) +
ma3_2 * (data_all[i-3,2] - u_2[i-3])
sigma_2[i] <- omega_2 + alpha_2*(data_all[i-1,1] - u_2[i-1])^2 + beta_2*sigma_2[i-1]
X_sim = cbind(u_1[i] + sigma_1[i] * Z.[,1],
u_2[i] + sigma_2[i] * Z.[,2])
Xs. <- rowSums(X_sim)
Xs.mean <- mean(Xs.)
Xs.CI <- quantile(Xs., probs = c(0.025, 0.975))
VaR <- quantile(Xs., probs = alpha)
Z_foo[i,] <- (data_all[i, ] - c(u_1[i], u_2[i])) / c(sigma_1[i], sigma_2[i])
U_sim[i, ] <- tail(pobs(Z_foo[1:i,]),1)
dependece[i] <- clayton_theta_cpp(parameters, 1e-10, U_sim[i-1,])
# dependece[i] <- ifelse(dependece[i] < 1000,
# ifelse(dependece[i] != 0, dependece[i], dependece[i-1]),
# dependece[i-1])
future[i-nrow(train)] <- Xs.mean
CI_out[i-nrow(train),] <- Xs.CI
VaR_out[i-nrow(train)] <- VaR
}
backtest_cop_t <- function(return = X_all, backtest_period = 200,
alpha = 0.99, m = 200) {
future <- matrix(ncol = 1, nrow = backtest_period)
CI_out <- matrix(ncol = 2, nrow = backtest_period)
VaR_out <- matrix(ncol = 1, nrow = backtest_period)
for (i in 0:(backtest_period - 1)) {
nY <- nrow(return) - (backtest_period - i)
X <- return[1:nY, ]
## Fit marginal time series
uspec <- rep(list(ugarchspec(distribution.model = "std")), ncol(X))
fit.ARMA.GARCH <- fit_ARMA_GARCH(X, ugarchspec.list = uspec)
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_) # vector of estimated df
d <- ncol(X) # dimension
### Fitting copulas ##
## Compute pseudo-observations from the standardized t residuals
U <- pobs(Z)
## Fitting a t copula
fit.tc <- fitCopula(tCopula(dim = d, dispstr = "un"),
data = U)#, method = "itau.mpl")
tc <- fit.tc@copula # fitted copula
### Simulating paths from the full model ##
## Simulate paths
m <- m
U. <- rCopula(m, copula = tc)
Z. <- sapply(1:d, function(j) sqrt((nu.mar[j] - 2) / nu.mar[j]) * qt(U.[,j], df = nu.mar[j]))
sim <- lapply(1:d, function(j)
ugarchsim(fits[[j]], n.sim = 1, m.sim = m,
custom.dist = list(name = "sample", distfit = t(Z.[,j, drop = FALSE]))))
fitted_sim <- sapply(sim, function(x) fitted(x))
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
}
return(list(
forecast = future,
CI = CI_out,
VaR = VaR_out
))
}
backtest_period <- 50
back_out <- backtest_cop_t(X_all, backtest_period = backtest_period)
nX = nrow(X_all) - backtest_period
X <- X_all[1:nX, ]
X_future <- X_all[(nX+1):nrow(X_all), ]
Xs <- rowSums(X) # aggregated return; n-vector
Xs.mean <- back_out$forecast # predicted aggregated return; m-vector
Xs.CI <- t(back_out$CI)# CIs; (2, m)-matrix
VaR <- back_out$VaR# VaR_alpha; m-vector
## Plot
past <- zoo::index(X)
future <- zoo::index(X_future)
plot(past, Xs, type = "l", xlim = range(c(past, future)), xlab = "", ylab = "",
ylim = c(-(max(range(VaR))+0.1), (max(range(VaR))+0.1))) # actual (past) returnes
polygon(c(future, rev(future)), c(Xs.CI[1,], rev(Xs.CI[2,])),
border = NA, col = "grey80") # CI region
lines(future, rowSums(X_future), col = "orange")
lines(future, Xs.mean, col = "royalblue3") # predicted aggregated return
lines(future, Xs.CI[1,], col = "grey50") # lower CI
lines(future, Xs.CI[2,], col = "grey50") # upper CI
lines(future, VaR, col = "maroon3") # 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)))))
mse = mean((rowSums(X_future) - Xs.mean)^2)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.