library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.3.3
library(ggfortify)
## Warning: package 'ggfortify' was built under R version 3.3.3
library(forecast)
## Warning: package 'forecast' was built under R version 3.3.3
library(Quandl)
## Warning: package 'Quandl' was built under R version 3.3.3
## Loading required package: xts
## Warning: package 'xts' was built under R version 3.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 3.3.3
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
# Read The Data
RPCE <- Quandl("FRED/PCECC96", type ="ts", collapse = "quarterly", start_date = "1955-01-01", end_date = "2017-12-04")
str(RPCE)
## Time-Series [1:252] from 1955 to 2018: 1599 1630 1650 1671 1673 ...
autoplot(RPCE, main = "The Quarterly Real Personal Consumption Expenditures")
# Construct the log changes in the Quarterly Real Personal Consumption Expenditures
dlRPCE <- diff(log(RPCE))
autoplot(dlRPCE, main = "The log changes in the Quarterly Real Personal Consumption Expenditures")
# split Sample - estimation sub-sample dates
dlRPCE1 <- window(dlRPCE, end = c(2008,4)) # Series in Sub-Sample One
autoplot(dlRPCE1, main = "The log changes in the Quarterly Real Personal Consumption Expenditures in Sub-Sample One")
dlRPCE2 <- window(dlRPCE, start = c(2009,1)) # Series in Sub-Sample Two
autoplot(dlRPCE2, main = "The log changes in the Quarterly Real Personal Consumption Expenditures in Sub-Sample TWO")
# find the best model using ARIMA model
# for Sub-Sample One
# Series: dlRPCE1
# ARIMA(3,0,2) with non-zero mean
m1 <- auto.arima(dlRPCE1, ic="aic", seasonal=FALSE, stationary=TRUE, stepwise=FALSE, approximation=FALSE)
m1
## Series: dlRPCE1
## ARIMA(3,0,2) with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 mean
## 0.1545 -0.4081 0.3333 0.0927 0.6006 0.0084
## s.e. 0.2643 0.2062 0.0837 0.2688 0.1796 0.0008
##
## sigma^2 estimated as 4.325e-05: log likelihood=777.99
## AIC=-1541.99 AICc=-1541.45 BIC=-1518.39
# for Sub-Sample Two
# Series: dlRPCE2
# ARIMA(2,0,0) with zero mean
m2 <- auto.arima(dlRPCE2, ic="aic", seasonal=FALSE, stationary=TRUE, stepwise=FALSE, approximation=FALSE)
m2
## Series: dlRPCE2
## ARIMA(2,0,0) with zero mean
##
## Coefficients:
## ar1 ar2
## 0.4588 0.4631
## s.e. 0.1463 0.1479
##
## sigma^2 estimated as 9.406e-06: log likelihood=157.38
## AIC=-308.77 AICc=-308.02 BIC=-304.02
# Accuracy in both Sub-Samples
# Accuracy in Sub-Sample One
accuracy(m1)
## ME RMSE MAE MPE MAPE
## Training set -3.687619e-05 0.006484194 0.004884715 30.00586 192.1146
## MASE ACF1
## Training set 0.6696203 -0.008823866
# Accuracy in Sub-Sample Two
accuracy(m2)
## ME RMSE MAE MPE MAPE
## Training set 0.000827199 0.002980545 0.002282254 24.04451 64.57332
## MASE ACF1
## Training set 0.7136509 -0.08991558
# Plot Inverted AR and MA roots in both Sub-Samples
# Inverted AR and MA roots in Sub-Sample One
plot(m1)
# Inverted AR and MA roots in Sub-Sample Two
plot(m2)
# Diagnose residulas using ggtsdiag in both Sub-Samples
nlags <- 40
# Diagnose residulas using ggtsdiag in Sub-Sample One
ggtsdiag(m1, gof.lag = nlags)
## Warning: package 'bindrcpp' was built under R version 3.3.3
# Diagnose residulas using ggtsdiag in Sub-Sample Two
ggtsdiag(m2, gof.lag = nlags)
## Warning: Removed 5 rows containing missing values (geom_point).
# forecaste to generate 1 to 36 step ahead forecast for the prediction Sub-Sample Two, 2009Q1-2017Q4
# create 1, 2, ..., h=36 step ahead forecasts, with 2008Q4 as forecast origin
# Multi-Step Forecast
m1.Mul.F <- forecast(m1, length(dlRPCE2))
plot(m1.Mul.F)
autoplot(m1.Mul.F)+geom_hline(yintercept = 0, linetype = "dashed")
## Rolling Scheme Forecast
m1.Rol.F <- zoo()
fstQ <- 1955.25
lstQ <- 2008.75
for(i in 1:length(dlRPCE2)) {
temp <- window(dlRPCE, start=fstQ + (i-1)/4, end = lstQ + (i-1)/4)
m1.update <- arima(temp, order=c(3,0,2))
m1.Rol.F <- c(m1.Rol.F, forecast(m1.update, 1)$mean)
}
m1.Rol.F <- as.ts(m1.Rol.F)
plot(m1.Rol.F)
# from 1955
plot(m1.Mul.F, type="o", pch=20, xlim=c(1955,2017), ylim=c(-0.03,0.03),
main = " Multi-step and Rolling Forecasts")
lines(m1.Mul.F$mean, type="p", pch=20, col="blue")
lines(m1.Rol.F, type="o", pch=20, col="red")
lines(dlRPCE, type="o", pch=20)
legend(1955, -0.014, c("The Log Changes in the Quarterly Real Personal Consumption Expenditures", "Multi-Step Forecast", "Rolling Scheme forecast"),
col = c(1, 2, 4), lty = c(1, 1, 1), merge = TRUE, bg = "gray90")
# from 2000
plot(m1.Mul.F, type="o", pch=20, xlim=c(2000,2017), ylim=c(-0.03,0.03),
main = " Multi-Step and Rolling Forecasts")
lines(m1.Mul.F$mean, type="p", pch=20, col="blue")
lines(m1.Rol.F, type="o", pch=20, col="red")
lines(dlRPCE, type="o", pch=20)
legend(2000, -0.014, c("The Log Changes in the Quarterly Real Personal Consumption Expenditures", "Multi-Step Forecast", "Rolling Scheme forecast"),
col = c(1, 2, 4), lty = c(1, 1, 1), merge = TRUE, bg = "gray90")
## Accuracy
# Multi-Step Forecast Accuracy
accuracy(m1.Mul.F$mean, dlRPCE2)
## ME RMSE MAE MPE MAPE ACF1
## Test set -0.002376687 0.003727399 0.003067461 150.5172 285.7756 0.1923808
## Theil's U
## Test set 0.5141596
# Rolling Scheme Forecast Accuracy
accuracy(m1.Rol.F, dlRPCE2)
## ME RMSE MAE MPE MAPE ACF1
## Test set -0.001035403 0.002979235 0.002444962 143.2454 226.7712 -0.1121612
## Theil's U
## Test set 0.1001448
# The Rolling scheme forecast method has less error than the Multi-Step Forecast method.
# Therefore, the Rolling scheme Forecast method have more accuracy than the Multi-Step
# Forecast method in predicting the Log Changes in the Quarterly Real Personal Consumption
# Expenditures.
# Summary
# As a result, the ARIMA(3,0,2) Rolling Scheme Forecast method performed the best than the
# ARIMA(3,0,2) Multi-Step Forecast Method.