# Wilmott CQF - Lecture 2.4 - R version
# Author: Martin Blais <blais@furius.ca>
#
# Save your SPX.xls Excel file as CSV format into SPX.csv.
# Run in batch mode like this:
#    R --quiet --no-save --args SPX.csv "S&P500 Vol" < vol.r >/dev/null
# You should obtain charts for the volatilities in output.
# Note: you can comment out library(moments) functions if you cannot install packages.

library(moments) # for skewness and kurtosis

args <- commandArgs(T)
fn_csv <- args[1]
title <- args[4]

fn_csv <- "SPX.csv"

# Read the data file and extract the price series.
data <- read.csv(fn_csv)
prices = rev(data$Adj.Close)
dates = rev(data$Date)

# Compute the returns, and get the mean and stddev.
lag <- prices[1:length(prices)-1]
returns <- diff(prices, 1)/lag

# Nb. of days in a year.
yeardays = 252

# Set the starting volatility and decay factor.
vol0 <- 0.2
lambda <- 0.99
## FIXME: todo, compute lambda from the nb. of days.

# Calculate EWMA/RiskMetrics volatility.
vol2a <- rep(0, length(returns)+1)
vol2a[1] = vol0*vol0
rpart = (1-lambda) * (returns**2) * 252 # precompute for speed
for (i in 2:length(vol2a)) {
  vol2a[i] = lambda * vol2a[i-1] + rpart[i-1]
}
vola = sqrt(vol2a)

# Create a plot for the EWMA vol.
ylim = c(0.05, 0.55)
png("vol_ewma.png", width=1200, height=800)
plot(as.integer(dates), vola, type="l", col="red", xlab="Time", ylab="Vol", ylim=ylim, main="S&P500 EWMA Vol")
dev.off()

# Calculate GARCH volatility.
alpha = 0.99
avgvol = sum(returns*returns)*252/length(returns)

vol2garch <- rep(0, length(returns)+1)
vol2garch[1] = vol0*vol0
for (i in 2:length(vol2a)) {
  ewma = lambda * vol2garch[i-1] + rpart[i-1]
  vol2garch[i] = alpha * ewma + (1-alpha) * avgvol
}
volgarch = sqrt(vol2garch)

# Create a plot for the GARCH vol.
png("vol_garch.png", width=1200, height=800)
plot(as.integer(dates), volgarch, type="l", col="red", xlab="Time", ylab="Vol", ylim=ylim, main="S&P500 GARCH Vol")
dev.off()