library(quantmod)
library(lattice)
library(timeSeries)
library(rugarch)
# Obtain the S&P500 returns and truncate the NA value
getSymbols("^GSPC", from="1950-01-01")
spReturns = diff(log(Cl(GSPC)))
spReturns[as.character(head(index(Cl(GSPC)),1))] = 0
# Create the forecasts vector to store the predictions
windowLength = 500
foreLength = length(spReturns) - windowLength
forecasts <- vector(mode="character", length=foreLength)
forenum <- numeric(foreLength)
foreDate <- character(foreLength)
for (d in 0:foreLength) {
# Obtain the S&P500 rolling window for this day
spReturnsOffset = spReturns[(1+d):(windowLength+d)]
# Fit the ARIMA model
final.aic <- Inf
final.order <- c(0,0,0)
for (p in 0:5) for (q in 0:5) {
if ( p == 0 && q == 0) {
next
}
arimaFit = tryCatch( arima(spReturnsOffset, order=c(p, 0, q)),
error=function( err ) FALSE,
warning=function( err ) FALSE )
if( !is.logical( arimaFit ) ) {
current.aic <- AIC(arimaFit)
if (current.aic < final.aic) {
final.aic <- current.aic
final.order <- c(p, 0, q)
final.arima <- arima(spReturnsOffset, order=final.order)
}
} else {
next
}
}
# Specify and fit the GARCH model
spec = ugarchspec(
variance.model=list(garchOrder=c(1,1)),
mean.model=list(armaOrder=c(final.order[1], final.order[3]), include.mean=T),
distribution.model="sged"
)
fit = tryCatch(
ugarchfit(
spec, spReturnsOffset, solver = 'hybrid'
), error=function(e) e, warning=function(w) w
)
# If the GARCH model does not converge, set the direction to "long" else
# choose the correct forecast direction based on the returns prediction
# Output the results to the screen and the forecasts vector
if(is(fit, "warning")) {
forecasts[d+1] = paste(index(spReturnsOffset[windowLength]), 1, sep=",")
forenum[d+1] <- 1
foreDate[d+1] <- as.character(index(spReturnsOffset[windowLength]))
print(paste(index(spReturnsOffset[windowLength]), 1, sep=","))
} else {
fore = ugarchforecast(fit, n.ahead=1)
ind = fore@forecast$seriesFor
forecasts[d+1] = paste(colnames(ind), ifelse(ind[1] < 0, -1, 1), sep=",")
forenum[d+1] <- ifelse(ind[1] < 0, -1, 1)
foreDate[d+1] <- colnames(ind)
print(paste(colnames(ind), ifelse(ind[1] < 0, -1, 1), sep=","))
}
}
# Output the CSV file to "forecasts.csv"
write.csv(forecasts, file="forecasts.csv", row.names=FALSE)
forexts <- xts::xts(forenum[1:(length(forenum)-1)],
order.by = lubridate::as_date(foreDate[1:(length(foreDate)-1)]))
# Create the ARIMA+GARCH returns
spIntersect = merge(forexts, spReturns, all = FALSE)
spArimaGarchReturns = spIntersect[,1] * spIntersect[,2]
# Create the backtests for ARIMA+GARCH and Buy & Hold
spArimaGarchCurve = log( cumprod( 1 + spArimaGarchReturns ) )
spBuyHoldCurve = log( cumprod( 1 + spIntersect[,2] ) )
spCombinedCurve = merge( spArimaGarchCurve, spBuyHoldCurve, all=F )
# Plot the equity curves
xyplot(
spCombinedCurve,
superpose=T,
col=c("darkred", "darkblue"),
lwd=2,
key=list(
text=list(
c("ARIMA+GARCH", "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.