Preamble

library("e1071")     # kurtosis
library("forecast")
library("tseries")   # garch

NADAQ Series

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 (Adj. Close)

log.NASDAQ <- log(NASDAQ)
plot(date, log.NASDAQ, type="l", col=2)

ARIMA-ARCH Model

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

Future Trend

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)

Simulating Future Paths

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)
}

First 5 Simulated Paths

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)
}