library("e1071") # kurtosis
library("forecast")
library("tseries") # garch
Daily Adjusted Close
data <- read.csv("http://ptrckprry.com/course/forecasting/data/nasdaq.csv")
date <- as.Date(data$date)
time <- seq_along(date)
NASDAQ <- data$nasdaq
log.NASDAQ <- log(NASDAQ)
plot(date, log.NASDAQ, type="l", col=2)
ARIMA Part
# Use p = 0, q = 0 for simplicity, like Black-Scholes
# (using values from AICc will work better in practice)
fit.mean <- Arima(log.NASDAQ, c(0, 1, 0), include.constant=TRUE, method="ML")
ARCH Part
resid <- as.numeric(residuals(fit.mean))
fit.var <- garch(resid, c(1,1), trace=FALSE)
ht <- fit.var$fit[,1]^2
drift <- coef(fit.mean)[["drift"]]
# Plot 21 most recent observations
n <- length(log.NASDAQ)
x <- log.NASDAQ[(n-20):n]
plot((-20):0, x, type="l", col=2, xlim=c(-20, 20),
xlab="time", ylab="log.NASDAQ")
# x0 is the most recent value of log.NASDAQ
x0 <- log.NASDAQ[n]
time.future <- 0:20
trend <- x0 + drift * time.future
lines(time.future, trend, col=1, lty=2)
GARCH Parameters
omega <- coef(fit.var)[["a0"]]
alpha <- coef(fit.var)[["a1"]]
beta <- coef(fit.var)[["b1"]]
Today’s Information
h0 <- ht[n] # Conditional variance
eps0 <- resid[n] # Shock
Simulate Code
garch.sim <- function(nsim, eps0, h0, omega, alpha, beta)
{
eps <- c(eps0, rep(NA, nsim))
h <- c(h0, rep(NA, nsim))
e <- c(NA, rnorm(nsim))
for (t in 2:(nsim+1)) {
h[t] <- omega + alpha * eps[t-1]^2 + beta * h[t-1]
eps[t] <- sqrt(h[t]) * e[t]
}
eps[-1]
}
Simulate 5000 Future Paths
nsim <- 20
nrep <- 5000
xsim <- matrix(NA, nrep, nsim + 1)
for (i in 1:nrep) {
eps <- garch.sim(nsim, eps0, h0, omega, alpha, beta)
xsim[i,] <- trend + diffinv(eps)
}
ylim <- range(c(x, as.vector(xsim[1:5,])))
plot((-20):0, x, type="l", col=2, xlim=c(-20, 20), ylim=ylim,
xlab="time", ylab="log.NASDAQ")
lines(time.future, trend, col=1, lty=2)
for (i in 1:5) {
lines(time.future, xsim[i,], col=1)
}