# Wilmott CQF - Lecture 1.1 - 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 SPX.png SPX.stats "S&P500" < SPX.r >/dev/null
# You should obtain a chart in SPX.png and numbers in file SPX.stats.
# 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]
fn_plot <- args[2]
fn_data <- args[3]
title <- args[4]

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

# Compute the returns, and get the mean and stddev.
lag <- prices[1:length(prices)-1]
returns <- diff(prices, 1)/lag
u <- mean(returns)
s <- sd(returns)
dt <- 1/252
#dt <- 5/(60*60*24*365)

# Normalize the returns and compute the range.
nreturns <- (returns - u)/s
rmin <- min(nreturns)
rmax <- max(nreturns)

# Create a histogram.
nbuckets <- 256
xlim <- c(-4,+4)
bucksize <- (xlim[2]-xlim[1])/nbuckets
h <- hist(nreturns, breaks=nbuckets, plot=F)

# Create a plot.
png(fn_plot, width=1200, height=800)
plot(h$mids, h$density, type="l", xlim=xlim, col="red",
     xlab="p-value", ylab="PDF", main=title)
lines(h$mids, dnorm(h$mids, mean=0, sd=1), col="green", lty="dashed")

# Compute some stats (you need the 'moments' package to make this work).
skew <- skewness(nreturns)
kurtosis <- kurtosis(nreturns)

# Compute the annualized drift and vol of the process.
drift <- u/dt
vol <- s/sqrt(dt)

# Output some values.
results <- list(name=title, mean=u, stddev=s,
                drift=drift,
                vol=vol,
                skew=skew,
                kurtosis=kurtosis)
write.csv(results, fn_data)