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

Normal Plot For 1-Day Ahead log NASDAQ

qqnorm(xsim[,2], col=2,
       main=paste("Ex. Kurtosis =", round(kurtosis(xsim[,2]), 1)))
qqline(xsim[,2], col=1, lty=2)

Normal Plot For 20-Day Ahead log NASDAQ

qqnorm(xsim[,21], col=2,
       main=paste("Ex. Kurtosis =", round(kurtosis(xsim[,21]), 1)))
qqline(xsim[,21], col=1, lty=2)

Black-Scholes

# price European call option
black.scholes <- function(S, K, r, h, sigma2)
{
    d1 <- ((log(S/K) + (r + sigma2/2) * h)
           / sqrt(sigma2 * h))
    d2 <- d1 - sqrt(sigma2 * h)
    pnorm(d1) * S - pnorm(d2) * K * exp(-r * h)
}

One Day At the Money Price

Black Scholes

bs.price <- black.scholes(S=exp(x0), K=exp(x0), r=drift, h=1, sigma2=fit.mean$sigma2)

cat("Black-Scholes Price: ", round(bs.price, 2), "\n", sep="")
Black-Scholes Price: 30.07

Simulation

payoff <- exp(xsim[,2]) - exp(x0)
in.money <- payoff > 0
payoff[!in.money] <- 0

hist(payoff[payoff > 0], col=2, xlab="Dollars",
     main = paste("Nonzero Payoffs (", round(mean(in.money) * 100, 2), "%)", sep=""))

cat("Monte-Carlo Estimate from GARCH Model: ",
    round(mean(payoff), 2),
    " (", round(sd(payoff) / sqrt(length(payoff)), 2), ")\n",
    sep="")
Monte-Carlo Estimate from GARCH Model: 14.64 (0.29)

20 Day At the Money Price

Black Scholes

bs.price <-  black.scholes(S=exp(x0), K=exp(x0), r=drift, h=20, sigma2=fit.mean$sigma2)
cat("Black-Scholes Price: ", round(bs.price, 2), "\n", sep="")
Black-Scholes Price: 150.96

Simulation

payoff <- exp(xsim[,21]) - exp(x0)
in.money <- payoff > 0
payoff[!in.money] <- 0

hist(payoff[payoff > 0], col=2, xlab="Dollars",
     main = paste("Nonzero Payoffs (", round(mean(in.money) * 100, 2), "%)", sep=""))

cat("Monte-Carlo Estimate from GARCH Model: ",
    round(mean(payoff), 2),
    " (", round(sd(payoff) / sqrt(length(payoff)), 2), ")\n",
    sep="")
Monte-Carlo Estimate from GARCH Model: 92.61 (1.72)

One Day 1% Out of Money Price

Black Scholes

bs.price <- black.scholes(S=exp(x0), K=1.01 * exp(x0), r=drift, h=1, sigma2=fit.mean$sigma2)

cat("Black-Scholes Price: ", round(bs.price, 2), "\n", sep="")
Black-Scholes Price: 9.14

Simulation

payoff <- exp(xsim[,2]) - 1.01 * exp(x0)
in.money <- payoff > 0
payoff[!in.money] <- 0

hist(payoff[payoff > 0], col=2, xlab="Dollars",
     main = paste("Nonzero Payoffs (", round(mean(in.money) * 100, 2), "%)", sep=""))

cat("Monte-Carlo Estimate from GARCH Model: ",
    round(mean(payoff), 2),
    " (", round(sd(payoff) / sqrt(length(payoff)), 2), ")\n",
    sep="")
Monte-Carlo Estimate from GARCH Model: 0.7 (0.06)

20 Day 1% Out of Money Price

Black Scholes

bs.price <-  black.scholes(S=exp(x0), K=1.01 * exp(x0), r=drift, h=20, sigma2=fit.mean$sigma2)
cat("Black-Scholes Price: ", round(bs.price, 2), "\n", sep="")
Black-Scholes Price: 121.48

Simulation

payoff <- exp(xsim[,21]) - 1.01 * exp(x0)
in.money <- payoff > 0
payoff[!in.money] <- 0

hist(payoff[payoff > 0], col=2, xlab="Dollars",
     main = paste("Nonzero Payoffs (", round(mean(in.money) * 100, 2), "%)", sep=""))

cat("Monte-Carlo Estimate from GARCH Model: ",
    round(mean(payoff), 2),
    " (", round(sd(payoff) / sqrt(length(payoff)), 2), ")\n",
    sep="")
Monte-Carlo Estimate from GARCH Model: 61.24 (1.46)

All Exercise Horizons

make.plot <- function(strike, main)
{
    h <- 1:20
    bs <- rep(NA, 20)
    garch <- rep(NA, 20)
    garch.se <- rep(NA, 20)

    for (i in 1:length(h)) {
        bs[i] <- black.scholes(S=exp(x0), K=strike, r=drift, h=h[i], sigma2=fit.mean$sigma2)

        payoff <- exp(xsim[,1 + h[i]]) - strike
        in.money <- payoff > 0
        payoff[!in.money] <- 0
        garch[i] <- mean(payoff)
        garch.se[i] <- sd(payoff) / sqrt(length(payoff))
    }

    plot(h, bs, col=1, main=main, xlab="Horizon", ylab="Price", type="l",
         ylim=range(c(bs, garch - 1.96 * garch.se, garch + 1.96 * garch.se)))
    points(h, garch, col=2, cex=0.5)
    lines(h, garch, col=2)
    segments(h, garch - 1.96 * garch.se, h, garch + 1.96 * garch.se, col=2)
    segments(h - 0.1, garch - 1.96 * garch.se, h + 0.1, garch - 1.96 * garch.se, col=2)
    segments(h - 0.1, garch + 1.96 * garch.se, h + 0.1, garch + 1.96 * garch.se, col=2)
    legend("topleft", inset=0.05,
           legend=c("Black-Scholes", "GARCH(1,1)"),
           lty=1, col=c(1,2))
}

make.plot(exp(x0), main="At the Money")

make.plot(1.01 * exp(x0), main="1% Out of the Money")