HW5
Mai AlSwayan
2nd Question
-Obtain the monthly data for Retail and Food Services Sales FRED/RSAFSNA *Use intervention analysis approach to incorporate explicitly the effect of 2008 recession on this time series: (1) Identify and estimate a seasonal ARIMA model for pre-intervention period up until June 2008. (2)Identify and estimate a transfer function model for the period up until December 2014, using July 2008 as the date of intervention. (3) Use the estimated transfer function model to produce and plot a multistep forecast, and a sequence of one month ahead forecasts until the end of 2016. Comment on the accuracy of the two forecasts.
Ansewr:
Part 1
Data and data overview
library(xts)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(Quandl)
RFSS<- Quandl("FRED/RSAFSNA", type="ts")str(RFSS)## Time-Series [1:303] from 1992 to 2017: 146376 147079 159336 163669 170068 ...summary(RFSS)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 146400 238400 318500 314100 383100 541200head(RFSS)## [1] 146376 147079 159336 163669 170068 168663tail(RFSS)## [1] 453530 468304 541179 422799 419730 482257plot(RFSS, xlab="Years", ylab="", main="monthly data for Retail and Food Services Sales ")
Part 2
rece <- window(RFSS, end=c(2008,6))
lrece <- log(rece)
dlrece1 <-diff(lrece)
dlrece12 <- diff(lrece, lag=12)
d2lrece1_12 <- diff(diff(lrece), lag=12) par(mfrow=c(2,3))
plot(RFSS, main=expression(RFSS))
plot(lrece, main=expression(log(rece)))
plot.new()
plot(dlrece1, main=expression(paste(Delta, "log(rece)")))
plot(dlrece12, main=expression(paste(Delta[12], "log(rece)")))
plot(d2lrece1_12, main=expression(paste(Delta, Delta[12], "log(rece)")))
library(forecast)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 parameterAcf(lrece, type='correlation', lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("ACF for lag(rece)")))
Acf(dlrece1, type='correlation', lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("ACF for",Delta," lag(rece)")))
Acf(dlrece12, type='correlation',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("ACF for",Delta[12]," lag(rece)")))
Acf(d2lrece1_12, type='correlation',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("ACF for",Delta, Delta[12],"lag(rece)")))
Acf(lrece, type='partial',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("PACF for lag(rece)")))
Acf(dlrece1, type='partial',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("PACF for",Delta," lag(rece)")))
Acf(dlrece12, type='partial',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("PACF for",Delta[12]," lag(rece)")))
Acf(d2lrece1_12, type='partial',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("PACF for",Delta,Delta[12]," lag(rece)")))
m1 <- Arima(dlrece1, order = c(1,1,0), seasonal = list(order = c(1,0,0), period = 12))
m1## Series: dlrece1
## ARIMA(1,1,0)(1,0,0)[12]
##
## Coefficients:
## ar1 sar1
## -0.6250 0.9666
## s.e. 0.0557 0.0113
##
## sigma^2 estimated as 0.001257: log likelihood=360.86
## AIC=-715.73 AICc=-715.6 BIC=-705.89tsdiag(m1, gof.lag=36)
m2 <- Arima(dlrece1, order = c(1,1,0), seasonal = list(order = c(2,0,0), period = 12))
m2## Series: dlrece1
## 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.81tsdiag(m2, gof.lag=36)
m3 <- Arima(dlrece1, order = c(2,0,0), seasonal = list(order = c(1,0,0), period = 12))
m3## Series: dlrece1
## 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.72tsdiag(m3, gof.lag=36)
m4 <- Arima(dlrece1, order = c(2,0,0), seasonal = list(order = c(1,0,0), period = 12))
m4## Series: dlrece1
## 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.72tsdiag(m4, gof.lag=36)
m5 <- Arima(dlrece1, order = c(4,1,0), seasonal = list(order = c(1,0,0), period = 12))
m5## Series: dlrece1
## ARIMA(4,1,0)(1,0,0)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 sar1
## -1.3001 -1.2152 -0.6653 -0.2957 0.9717
## s.e. 0.0683 0.1050 0.1044 0.0685 0.0101
##
## sigma^2 estimated as 0.0006354: log likelihood=427.42
## AIC=-842.84 AICc=-842.4 BIC=-823.18tsdiag(m5, gof.lag=36)
m6 <- Arima(dlrece1, order = c(2,0,1), seasonal = list(order = c(2,0,0), period = 12))
m6## Series: dlrece1
## 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.66tsdiag(m6, gof.lag=36)
m7 <- Arima(dlrece1, order = c(3,0,0), seasonal = list(order = c(1,0,0), period = 12))
m7## Series: dlrece1
## 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.57tsdiag(m7, gof.lag=36)
m8 <- Arima(dlrece1, order = c(9,1,0), seasonal = list(order = c(1,0,0), period = 12))
m8## Series: dlrece1
## 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.22tsdiag(m8, gof.lag=36)
m9 <- Arima(dlrece1, order = c(9,1,0), seasonal = list(order = c(0,1,0), period = 12))
m9## Series: dlrece1
## ARIMA(9,1,0)(0,1,0)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7
## -1.6045 -1.9247 -1.9245 -1.9111 -1.6666 -1.3007 -0.9877
## s.e. 0.0728 0.1308 0.1808 0.2108 0.2220 0.2099 0.1810
## ar8 ar9
## -0.6656 -0.1638
## s.e. 0.1327 0.0745
##
## sigma^2 estimated as 0.0004338: log likelihood=453.37
## AIC=-886.73 AICc=-885.46 BIC=-854.58tsdiag(m9, gof.lag=36)
m10<- Arima(dlrece1, order = c(9,1,0), seasonal = list(order = c(2,0,0), period = 12))
m10## Series: dlrece1
## 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.62tsdiag(m10, gof.lag=36)
rece <- RFSS
m <- m10
m.f <- forecast(m, h=48)
m.f.err <- lrece - m.f$mean## Warning in .cbind.ts(list(e1, e2), c(deparse(substitute(e1))[1L],
## deparse(substitute(e2))[1L]), : non-intersecting seriesForcast
end.month <- 2008 + 6/12
rece1 <- window(rece, end = end.month)
rece2 <- window(rece, start = end.month + 1/12)Multistep Forecast
multiple<- forecast(m10, 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(lrece, 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(lrece, type="l", pch=16, lty="dashed")
Conclusion
we can see the both m10 actual data and the foretasted outcomes has fairly similar results.