Skip to content

Instantly share code, notes, and snippets.

What would you like to do?
# R code re: post for analyzing combination forecasts in R
# "Projecting The Payrolls Trend With Combination Forecasts"
# 22 Dec 2016
# By James Picerno
# (c) 2016 by Beta Publishing LLC
# load packages
# download time series
fred.tickers <-c("USPRIV", "INDPRO", "DSPIC96", "PCEC96")
# 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 <-function(x) {
a <-VARselect(x)
b <-VAR(y=x, p=a$selection[3])
c <-predict(b, n.ahead=sa)
var.f <-sapply(seq(60,nrow(dat.train)), function(x)[1:x,]))
# arima forecasts <-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)$USPRIV[1:x]))
# linear regression forecasts <-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)$USPRIV[1:x]))
# neural net forecasts <-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)$USPRIV[1:x]))
# naive forecasts (use last actual data point as 1-step-ahead forecast) <-function(x) {
a <-naive(x, h=sa)
b <-a$mean
naive.f <-sapply(seq(60, length(dat.train$USPRIV)), function(x)$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)
# mean absolute error function
mae <- function(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],
par(mar = c(5.2,5,3,5),oma=c(0,0,0.2,0))
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"))
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.