In this assignment we incorporate the effect of the 2008 recession onto the monthly data for Total Private Residential Construction Spending FRED/PRRESCON. We are going to identify and estimate a seasonal ARIMA model, also, identify and estimate a transfer function model.
we imported the date from Quandl website (FRED/RSAFSNA), the monthly data for Retail and Food Services Sales FRED/RSAFSNA. The data until the end of 2013.
library(Quandl)
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(forecast)
library(xts)
library(forecast)
Quandl.api_key("3fnnc4EE2uzFb1a44mAQ")
rfss <- Quandl("FRED/RSAFSNA", type = "ts")
y <- window(rfss, end=c(2008,6))
ly <- log(y)
dly1 <-diff(ly)
dly12 <- diff(ly, lag=12)
d2ly1_12 <- diff(diff(ly), lag=12)
par(mfrow=c(2,3))
plot(rfss, main=expression(rfss))
plot(ly, main=expression(log(y)))
plot.new()
plot(dly1, main=expression(paste(Delta, "log(y)")))
plot(dly12, main=expression(paste(Delta[12], "log(y)")))
plot(d2ly1_12, main=expression(paste(Delta, Delta[12], "log(y)")))
maxlag <- 20
par(mfrow=c(2,4), max=c(3,3,4,2))
## Warning in par(mfrow = c(2, 4), max = c(3, 3, 4, 2)): "max" is not a
## graphical parameter
Acf(ly, type='correlation', lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("ACF for lag(y)")))
Acf(dly1, type='correlation', lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("ACF for",Delta," lag(y)")))
Acf(dly12, type='correlation',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("ACF for",Delta[12]," lag(y)")))
Acf(d2ly1_12, type='correlation',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("ACF for",Delta, Delta[12],"lag(y)")))
Acf(ly, type='partial',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("PACF for lag(y)")))
Acf(dly1, type='partial',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("PACF for",Delta," lag(y)")))
Acf(dly12, type='partial',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("PACF for",Delta[12]," lag(y)")))
Acf(d2ly1_12, type='partial',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("PACF for",Delta,Delta[12]," lag(y)")))
m1 <- Arima(dly1, order = c(1,0,0), seasonal = list(order = c(1,0,0), period = 12))
m2 <- Arima(dly1, order = c(1,1,0), seasonal = list(order = c(2,0,0), period = 12))
m3 <- Arima(dly1, order = c(2,0,0), seasonal = list(order = c(1,0,0), period = 12))
m4 <- Arima(dly1, order = c(2,1,0), seasonal = list(order = c(1,0,0), period = 12))
m5 <- Arima(dly1, order = c(2,0,1), seasonal = list(order = c(2,0,0), period = 12))
m6 <- Arima(dly1, order = c(3,0,0), seasonal = list(order = c(1,0,0), period = 12))
m7 <- Arima(dly1, order = c(3,1,0), seasonal = list(order = c(1,0,0), period = 12))
m8 <- Arima(dly1, order = c(6,0,1), seasonal = list(order = c(1,0,0), period = 12))
m9 <- Arima(dly1, order = c(6,1,0), seasonal = list(order = c(2,0,0), period = 12))
m10 <- Arima(dly1, order = c(9,1,0), seasonal = list(order = c(1,0,0), period = 12))
m11 <- Arima(dly1, order = c(9,1,0), seasonal = list(order = c(2,0,0), period = 12))
m1
## Series: dly1
## ARIMA(1,0,0)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 sar1 mean
## -0.4746 0.9687 0.0026
## s.e. 0.0627 0.0107 0.0163
##
## sigma^2 estimated as 0.0005405: log likelihood=446.15
## AIC=-884.3 AICc=-884.09 BIC=-871.17
m2
## Series: dly1
## ARIMA(1,1,0)(2,0,0)[12]
##
## Coefficients:
## ar1 sar1 sar2
## -0.6253 0.8858 0.0841
## s.e. 0.0557 0.0748 0.0767
##
## sigma^2 estimated as 0.001253: log likelihood=361.46
## AIC=-714.92 AICc=-714.71 BIC=-701.81
m3
## Series: dly1
## ARIMA(2,0,0)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 ar2 sar1 mean
## -0.6777 -0.4288 0.9706 0.0028
## s.e. 0.0645 0.0642 0.0102 0.0107
##
## sigma^2 estimated as 0.0004412: log likelihood=466.07
## AIC=-922.14 AICc=-921.82 BIC=-905.72
m4
## Series: dly1
## ARIMA(2,1,0)(1,0,0)[12]
##
## Coefficients:
## ar1 ar2 sar1
## -1.0214 -0.6273 0.9755
## s.e. 0.0564 0.0560 0.0088
##
## sigma^2 estimated as 0.0007566: log likelihood=408.84
## AIC=-809.68 AICc=-809.47 BIC=-796.57
m5
## Series: dly1
## ARIMA(2,0,1)(2,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 ar2 ma1 sar1 sar2 mean
## -0.3899 -0.2998 -0.4412 0.7088 0.2667 0.0033
## s.e. 0.2078 0.1430 0.2218 0.0810 0.0840 0.0071
##
## sigma^2 estimated as 0.0004048: log likelihood=475.32
## AIC=-936.64 AICc=-936.05 BIC=-913.66
m6
## Series: dly1
## ARIMA(3,0,0)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 sar1 mean
## -0.7098 -0.4814 -0.0768 0.9692 0.0029
## s.e. 0.0710 0.0809 0.0721 0.0107 0.0096
##
## sigma^2 estimated as 0.0004422: log likelihood=466.63
## AIC=-921.27 AICc=-920.83 BIC=-901.57
m7
## Series: dly1
## ARIMA(3,1,0)(1,0,0)[12]
##
## Coefficients:
## ar1 ar2 ar3 sar1
## -1.2107 -0.9395 -0.3130 0.9751
## s.e. 0.0680 0.0870 0.0692 0.0090
##
## sigma^2 estimated as 0.0006876: log likelihood=418.55
## AIC=-827.1 AICc=-826.79 BIC=-810.71
m8
## Series: dly1
## ARIMA(6,0,1)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ma1 sar1
## -0.2909 -0.2716 -0.0416 -0.1979 -0.0191 0.0581 -0.4788 0.9662
## s.e. 0.5972 0.4621 0.3833 0.2178 0.2241 0.1292 0.5918 0.0116
## mean
## 0.0033
## s.e. 0.0059
##
## sigma^2 estimated as 0.0004197: log likelihood=474.18
## AIC=-928.35 AICc=-927.17 BIC=-895.52
m9
## Series: dly1
## ARIMA(6,1,0)(2,0,0)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 sar1 sar2
## -1.5512 -1.7380 -1.4222 -1.1746 -0.7348 -0.2394 0.6122 0.3711
## s.e. 0.0700 0.1214 0.1534 0.1533 0.1208 0.0695 0.0715 0.0725
##
## sigma^2 estimated as 0.0004582: log likelihood=457.99
## AIC=-897.99 AICc=-897.02 BIC=-868.48
m10
## Series: dly1
## ARIMA(9,1,0)(1,0,0)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7
## -1.6051 -1.9386 -1.9608 -1.9519 -1.6934 -1.3168 -0.9821
## s.e. 0.0706 0.1271 0.1766 0.2062 0.2175 0.2059 0.1773
## ar8 ar9 sar1
## -0.6367 -0.1426 0.9721
## s.e. 0.1289 0.0711 0.0101
##
## sigma^2 estimated as 0.0004297: log likelihood=466.64
## AIC=-911.28 AICc=-909.84 BIC=-875.22
m11
## Series: dly1
## ARIMA(9,1,0)(2,0,0)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7
## -1.6369 -2.0074 -1.9864 -1.9970 -1.7098 -1.3333 -0.9817
## s.e. 0.0705 0.1284 0.1836 0.2143 0.2273 0.2144 0.1850
## ar8 ar9 sar1 sar2
## -0.6449 -0.1156 0.5703 0.4140
## s.e. 0.1331 0.0728 0.0692 0.0707
##
## sigma^2 estimated as 0.0003634: log likelihood=480.98
## AIC=-937.96 AICc=-936.25 BIC=-898.62
par(mar=c(1,1,1,1))
tsdiag(m1, gof.lag=36)
tsdiag(m2, gof.lag=36)
tsdiag(m3, gof.lag=36)
tsdiag(m4, gof.lag=36)
tsdiag(m5, gof.lag=36)
tsdiag(m6, gof.lag=36)
tsdiag(m7, gof.lag=36)
tsdiag(m8, gof.lag=36)
tsdiag(m9, gof.lag=36)
tsdiag(m10, gof.lag=36)
tsdiag(m11, gof.lag=36)
## according to the results, the best model over all 11 models that we is ARIMA(9,1,0),(2,0,0) tested, so we will use it for the forcasting.
RFSS <- rfss
m <- m10
m.f <- forecast(m, h=48)
end.month <- 2008 + 6/12
rece1 <- window(RFSS, end = end.month)
rece2 <- window(RFSS, start = end.month + 1/12)
fst <- 1992-1 # 1992-01-01
lst <- 2017-2 # 2017-02-01
multiple<- forecast(m11, length(rece2))
plot(multiple, type="o", pch=16, xlim=c(2008,2014), ylim=c(-0.03,0.03),
main="Multistep Forecasts: Retail and Food Services Sales")
lines(multiple$mean, type="p", pch=16, lty="dashed", col="blue")
lines(ly, type="o", pch=16, lty="dashed")
m1.fcast <- forecast(m9, h=25)
plot(m1.fcast, type="o", pch=16, xlim=c(2008,2014), ylim=c(-0.03,0.03),
main=" Forecasts: Retail and Food Services Sales")
lines(multiple$mean, type="l", pch=16, lty="dashed", col="blue")
lines(ly, type="l", pch=16, lty="dashed")
Since ARIMA(9,1,0)(2,1,0) has similar orginal and forcasting data, So ARIMA(9,1,0)(2,1,0) is the best model for the Retail and Food Services Sales and 1-month Ahead Recursive Forecasting works better performence than multistep ahead Forecasting.