library(xts)
library(qrmdata)
library(rmgarch)
## Load some real data
tickers <- c("AAPL", "AMZN", "MSFT")
stocks <- lapply(tickers, quantmod::getSymbols , auto.assign = FALSE)
stocksAd <- do.call(cbind, lapply(stocks, function(x) {
quantmod::Ad(x)
}))
stocksdiff = diff(log(stocksAd))
colnames(stocksdiff) <- tickers
ReturnsCop = stocksdiff[-1,]
Cop_test_Par <- function(d, returnsCop = ReturnsCop, wL = windowLength, ARIMAfit = FALSE) {
ReturnsOffset = returnsCop[(1+d):(wL+d),]
if(ARIMAfit == TRUE) {
margins = list()
for(i in 1:ncol(ReturnsOffset)) {
Order = forecast::auto.arima(ReturnsOffset[,i])
Order <- as.character(Order)
Order <- stringr::str_extract_all(Order, "(?<=\\().*(?=\\))")
Order <- stringr::str_split_fixed(Order, ",", 3)
Order <- as.numeric(Order)
if(Order[2] !=0) {
for(i in 1:Order[2]) {
ReturnsOffset[,i] <- diff(ReturnsOffset[,i])
ReturnsOffset[1,i] <- 0
}
}
unifit = tryCatch(rugarch::ugarchspec(
variance.model=list(model = "sGARCH", garchOrder=c(1,1)),
mean.model=list(armaOrder=c(Order[1], Order[3]), include.mean=T),
distribution.model="std"),
error=function(e) e, warning=function(w) w)
if(is(unifit, "warning") | is(unifit, "error")) {
margins[[i]] = rugarch::ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
mean.model = list(armaOrder = c(2,1), include.mean = TRUE),
distribution.model = "std")
} else {
margins[[i]] = unifit
}
}
mulspec = tryCatch(rugarch::multispec(margins),
error=function(e) e, warning=function(w) w)
if(is(unifit, "warning") | is(unifit, "error")) {
## Specify univariate AR(1)-GARCH(1,1) for both marginal processes
uspec <- rugarch::ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
mean.model = list(armaOrder = c(2,1), include.mean = TRUE),
distribution.model = "std")
## Combine univariate specs to obtain spec for marginal models
marginspec <- rugarch::multispec(replicate(ncol(ReturnsOffset), uspec))
} else {
marginspec <- mulspec
}
} else {
## Specify univariate AR(1)-GARCH(1,1) for both marginal processes
uspec <- rugarch::ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
mean.model = list(armaOrder = c(2,1), include.mean = TRUE),
distribution.model = "std")
## Combine univariate specs to obtain spec for marginal models
marginspec <- rugarch::multispec(replicate(ncol(ReturnsOffset), uspec))
}
if("error" %in% class(marginspec)) {
## Specify univariate AR(1)-GARCH(1,1) for both marginal processes
uspec <- rugarch::ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
mean.model = list(armaOrder = c(2,1), include.mean = TRUE),
distribution.model = "std")
## Combine univariate specs to obtain spec for marginal models
marginspec <- rugarch::multispec(replicate(ncol(ReturnsOffset), uspec))
}
# A model with a changing copula
copspec <- rmgarch::cgarchspec(uspec = marginspec,
distribution.model = list(copula = "mvt", method = "ML",
time.varying = TRUE, transformation = "parametric"))
fit = tryCatch(rmgarch::cgarchfit(copspec, data = ReturnsOffset,
spd.control = list(lower = 0.1, upper = 0.9, type = "pwm",
kernel = "epanech"),
fit.control = list(eval.se = FALSE, trace = TRUE),
solver = "solnp"),
error=function(e) e, warning=function(w) w)
return(fit)
}
# CopSim_test_Par <- function(d, returnsCop = ReturnsCop,
# wL = windowLength, Copfit = forecastsCop) {
# ReturnsOffset = returnsCop[(1+d):(wL+d),]
# x = Copfit[[d+1]]
#
# if(any(class(x) != "cGARCHfit")) {
# return(rep(1, ncol(ReturnsOffset)))
# } else {
# preR = last(rmgarch::rcor(x),1)[,,1]
# diag(preR) = 1
# # maxGO = x@model$maxgarchOrder
# sim1 = rmgarch::cgarchsim(x, n.sim = 1, n.start = 0, m.sim = 1000,
# startMethod = "sample", preR = preR, rseed = 1:1000)
# sim = list()
# for(z in 1:ncol(sim1@msim$simX[[1]])) {
# sim[[z]] <- mean(sapply(sim1@msim$simX, FUN = function(x) x[,z]))
# }
# simmean = do.call(c, sim)
# return(as.numeric(ifelse(simmean < 0, -1, 1)))
# }
# }
windowLength = 100
foreLength = nrow(ReturnsCop["2018/"]) - windowLength
Period = ReturnsCop["2018/"]
Index = zoo::index(Period[(nrow(Period)-foreLength):nrow(Period),])
Name = NULL
# In parallel
require(parallel)
no_cores <- detectCores()
cl <- makeCluster(no_cores)
clusterEvalQ(cl, library(dplyr))
system.time(
forecastsCop <- parLapply(cl, 0:foreLength, Cop_test_Par,
returnsCop = Period, wL = windowLength,
ARIMAfit = FALSE)
)
# stopCluster(cl)
# cl <- makeCluster(no_cores)
# clusterEvalQ(cl, library(dplyr))
# system.time(
# foreSimCop <- parLapply(cl, 0:foreLength, CopSim_test_Par,
# returnsCop = Period, wL = windowLength,
# Copfit = forecastsCop)
# )
# stopCluster(cl)
# lapply(forecastsCop, function(x){any(class(x) == "cGARCHfit")})
# lapply(seq_along(x), function(y, n, i) { paste(n[[i]], y[[i]]) }, y=x, n=names(x))
system.time(
foreCopTest <- do.call(rbind, lapply(0:foreLength, function(returnsCop, wL, Copfit, d) {
ReturnsOffset = returnsCop[(1+d):(wL+d),]
x = Copfit[[d+1]]
if(any(class(x) != "cGARCHfit")) {
rep(1, ncol(Period))
} else {
preR = rmgarch::last(rmgarch::rcor(x),1)[,,1]
diag(preR) = 1
maxGO = x@model$maxgarchOrder
sim1 = rmgarch::cgarchsim(x, n.sim = 1, n.start = 0, m.sim = 1000, presigma = tail(rugarch::sigma(x), maxGO),
startMethod = "sample", preR = preR, prereturns = tail( as.matrix(Period), maxGO),
preresiduals = tail(rugarch::residuals(x), maxGO),rseed = 1:1000,
cluster = cl)
sim = list()
for(z in 1:ncol(sim1@msim$simX[[1]])) {
sim[[z]] <- mean(sapply(sim1@msim$simX, FUN = function(x) x[,z]))
}
simmean = do.call(c, sim)
as.numeric(ifelse(simmean < 0, -1, 1))
}}, returnsCop = Period, wL = windowLength,
Copfit = forecastsCop))
)
stopCluster(cl)
forCopTest <- foreCopTest %>%
dplyr::as_tibble()
readr::write_csv(forCopTest, path="./output/forecasts_Cop_mult.csv")
#### Univariate returns ####
Name = "AAPL"
fore_Cop <- readr::read_csv(file="./output/forecasts_Cop_mult.csv")
colnames(fore_Cop) <- paste0(tickers,".pos")
for_Cop <- fore_Cop %>%
dplyr::transmute(Date = lubridate::as_date(Index),
pos = as.numeric(pull(fore_Cop, paste0(Name, ".pos"))))
for_Copxts <- xts::xts(for_Cop$pos[1:(nrow(for_Cop)-1)],
order.by = for_Cop$Date[1:(nrow(for_Cop)-1)])
# Create the ARIMA+GARCH returns
Intersect_Cop = merge(for_Copxts, ReturnsCop[,Name], all = FALSE)
ArimaGarchReturns_Cop = Intersect_Cop[,1] * Intersect_Cop[,2]
# Create the backtests for ARIMA+GARCH and Buy & Hold
ArimaGarchCurve_Cop = exp( cumsum( ArimaGarchReturns_Cop )) - 1
BuyHoldCurve_Cop = exp( cumsum( Intersect_Cop[,2] )) - 1
CombinedCurve_Cop = merge( ArimaGarchCurve_Cop, BuyHoldCurve_Cop, all=F )
# Plot the equity curves
lattice::xyplot(
CombinedCurve_Cop,
superpose=T,
col=c("darkred", "darkblue"),
lwd=2,
key=list(
text=list(
c("Cop", "Buy & Hold")
),
lines=list(
lwd=2, col=c("darkred", "darkblue")
)
)
)
#### Multivariate returns ####
fore_Copmul <- readr::read_csv(file="./output/forecasts_Cop_mult.csv")
colnames(fore_Copmul) <- paste0(tickers,".pos")
for_Copmul <- fore_Copmul %>%
dplyr::mutate(Date = lubridate::as_date(Index))
for_Copmulxts <- xts::xts(for_Copmul[1:(nrow(for_Copmul)-1), 1:(ncol(for_Copmul)-1)],
order.by = for_Copmul$Date[1:(nrow(for_Copmul)-1)])
# Create the ARIMA+GARCH returns
ArimaGarchReturns_Copmul = xts::xts(rowSums(for_Copmulxts * ReturnsCop[zoo::index(for_Copmulxts)] * (1/ncol(for_Copmulxts))),
order.by = zoo::index(for_Copmulxts))
BuyHoldReturns_Copmul <- xts::xts(rowSums(ReturnsCop[zoo::index(for_Copmulxts)] * (1/ncol(for_Copmulxts))),
order.by = zoo::index(for_Copmulxts))
# Create the backtests for ARIMA+GARCH and Buy & Hold
ArimaGarchCurve_Copmul = exp( cumsum( ArimaGarchReturns_Copmul )) - 1
BuyHoldCurve_Copmul = exp( cumsum( BuyHoldReturns_Copmul )) - 1
CombinedCurve_Copmul = merge( ArimaGarchCurve_Copmul, BuyHoldCurve_Copmul, all=F )
# Plot the equity curves
lattice::xyplot(
CombinedCurve_Copmul,
superpose=T,
col=c("darkred", "darkblue"),
lwd=2,
key=list(
text=list(
c("Copmul", "Buy & Hold")
),
lines=list(
lwd=2, col=c("darkred", "darkblue")
)
)
)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.