#Setup
rm(list = ls(all = TRUE))
setwd('path.to/cv.ts')

#Load Packages
require(forecast)
require(doParallel)
source('R/cv.ts.R')
source('R/forecast functions.R')

#Download S&P 500 data and adjust from splits/dividends
library(quantmod)
getSymbols('^GSPC', from='1990-01-01')
GSPC <- adjustOHLC(GSPC, symbol.name='^GSPC')

#Calculate monthly returns
GSPC <- to.monthly(GSPC, indexAt='lastof')
GSPC <- Cl(GSPC)

#Convert from xts to ts
GSPC <- ts(GSPC, start=c(1990,1), frequency=12)

#Start a cluster to speed up cross validaiton
cl <- makeCluster(4, type='SOCK')
registerDoParallel(cl)

#Define cross validation parameters
myControl <- tseriescontrol(
                  minObs=60,
                  stepSize=1,
                  maxHorizon=12,
                  fixedWindow=TRUE,
                  preProcess=FALSE,
                  ppMethod='guerrero',
                  summaryFunc=tsSummary
                  )

#Forecast using several models
result_naive <- cv.ts(GSPC, naiveForecast, myControl)
myControl$preProcess <- TRUE
result_autoarima <- cv.ts(GSPC, auto.arimaForecast, myControl, ic='bic')
result_ets <- cv.ts(GSPC, etsForecast, myControl, ic='bic')

#Stop cluster
stopCluster(cl)

#Plot error
require(reshape2)
require(ggplot2)
plotData <- data.frame(
  horizon=1:12
  ,naive =result_naive$results$MAPE[1:12]
  ,arima=result_autoarima$results$MAPE[1:12]
  ,ets=result_ets$results$MAPE[1:12]
  )
plotData <- melt(plotData, id.vars='horizon', value.name='MAPE', variable.name='model')
ggplot(plotData, aes(horizon, MAPE, color=model)) + geom_line()