# Wilmott CQF - Lecture 2.4 - R version # Author: Martin Blais # # 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()