## Preamble

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

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