Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
# R code re: CapitalSpecator.com post for analyzing combination forecasts in R
# "Projecting The Payrolls Trend With Combination Forecasts"
# http://www.capitalspectator.com/projecting-the-payrolls-trend-with-combination-forecasts/#more-8309
# 22 Dec 2016
# By James Picerno
# http://www.capitalspectator.com/
# (c) 2016 by Beta Publishing LLC
# load packages
library(TTR)
library(quantmod)
library(vars)
library(forecast)
# download time series
fred.tickers <-c("USPRIV", "INDPRO", "DSPIC96", "PCEC96")
getSymbols(fred.tickers,src="FRED")
# process data
dat <-na.omit(cbind(USPRIV, INDPRO, DSPIC96, PCEC96))
dat.train <-ROC(dat['2002-04-01::2012-08-01'],12,"discrete",na.pad=FALSE)
colnames(dat.train) <-c("USPRIV", "INDPRO", "PI", "PCE")
# step-ahead instruction
sa=1 # number of steps ahead for forecasts
# VAR return forecasts
var.f.fun <-function(x) {
a <-VARselect(x)
b <-VAR(y=x, p=a$selection[3])
c <-predict(b, n.ahead=sa)
c$fcst$USPRIV[,1]
}
var.f <-sapply(seq(60,nrow(dat.train)), function(x) var.f.fun(dat.train[1:x,]))
# arima forecasts
aa.fun <-function(x) {
a <-auto.arima(x)
b <-forecast(a, h=sa)
c <-as.numeric(b$mean)
}
aa.f <-sapply(seq(60, length(dat.train$USPRIV)), function(x) aa.fun(dat.train$USPRIV[1:x]))
# linear regression forecasts
lr.fun <-function(x) {
a <-lm(x[,1] ~ ., data=x)
b <-fitted(a)
c <-forecast(b, h=sa)
d <-as.numeric(c$mean)
}
lr.f <-sapply(seq(60, length(dat.train$USPRIV)), function(x) lr.fun(dat.train$USPRIV[1:x]))
# neural net forecasts
nn.fun <-function(x) {
a <-nnetar(x)
b <-na.omit(fitted(a))
c <-forecast(b, h=sa)
d <-as.numeric(c$mean)
}
nn.f <-sapply(seq(60, length(dat.train$USPRIV)), function(x) nn.fun(dat.train$USPRIV[1:x]))
# naive forecasts (use last actual data point as 1-step-ahead forecast)
n.fun <-function(x) {
a <-naive(x, h=sa)
b <-a$mean
}
naive.f <-sapply(seq(60, length(dat.train$USPRIV)), function(x) n.fun(dat.train$USPRIV[1:x]))
# aggregate forecasts
actual.forecasts <-cbind(dat.train[61:nrow(dat.train),1],
head(var.f, -1),
head(aa.f, -1),
head(lr.f, -1),
head(nn.f, -1),
head(naive.f, -1)) # delete last entry: 1-step ahead forecast doesn't have actual data
# generate forecast based on mean equal weight of forecasts
ew.f <-xts(apply(actual.forecasts[,2:5],1,mean),as.Date(index(actual.forecasts)))
actual.forecasts$ew <-ew.f
colnames(actual.forecasts) <-c("uspriv", "var", "aa", "lr", "nn", "naive", "ew")
# root mean square error function
rmse <- function(error)
{sqrt(mean(error^2))}
# mean absolute error function
mae <- function(error)
{mean(abs(error))}
# generate error analytics
rmse.error <-round(sapply(2:7, function(x) rmse(actual.forecasts[,1] - actual.forecasts[,x])),5)
mae.error <-round(sapply(2:7, function(x) mae(actual.forecasts[,1] - actual.forecasts[,x])),5)
error.all <-rbind(rmse.error, mae.error)
colnames(error.all) <-c("var", "aa", "lr", "nn", "naive", "ew")
# simple charts
# RMSE barplot
barplot(error.all[1,], beside=T, main="Root Mean Square Error (RMSE)")
# Forecast error plot: EW vs. Naive, rolling 1yr % change
errors.1 <-cbind(actual.forecasts[,1]-actual.forecasts[,6],
actual.forecasts[,1]-actual.forecasts[,7])*100
par(mar = c(5.2,5,3,5),oma=c(0,0,0.2,0))
matplot(errors.1,type="l")
colnames(errors.1) <-c("naive.errors", "ew.errors")
legend("topleft",c("EW Combination", "Naive"),fill=c("red","black"))
# EW Combination Forecast vs. actual private payrolls, rolling 1yr % change
dat.a <-cbind(actual.forecasts[,1],actual.forecasts[,7])*100
par(mar = c(5.2,5,3,5),oma=c(0,0,0.2,0))
matplot(dat.a, type="l")
legend("topleft",c("Actual", "EW Combination"),fill=c("black","red"))
# END
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.