I. Introduction

In this paper we apply the intervention analysis approach to the Retail and Food Services Sales data in the United States. Specifically, we want to incorporate the effect of the 2008 recession onto this time series. The remainder of the paper is as follows. Section II describes what intervention analysis is as well as the basic steps of the analysis. Section III describes the data. Section IV through Section VIII performs the five steps of the intervention analysis described in Section II. Section IX concludes.

II. Intervention Analysis

In this section we briefly describe the process of intervention analysis that will be used in the remainder of this paper. Intervention analysis is composed of primarily five steps:

  1. Specify the time-series model for data before the intervention
  2. Use the model from step (1) to predict the data after the intervention
  3. Examine the difference between the forecasted values and the actual values to specify the pulse and step functions
  4. Perform a joint estimation using all of the data
  5. Check the estimated model for adequacy

III. Downloading the Data

For this assignment we download monthly data for Retail and Food Services Sales from the St. Louis Federal Reserve Database (FRED). We download the data from FRED directly using the following code.

# Import the Data
Quandl.api_key("75HFthvbcxsrteNWupus")
data.RFSS <- Quandl("FRED/RSAFSNA", type = 'zoo')

To see why we want to perform intervention analysis on this data set we plot the original time series over the full sample period from January 1992 to February 2017.

From the above plot we can see that there is a clear drop in retail and food services sales after the recession of 2007 - 2009, which we estimate to have occured for this data set in July 2008. Thus, we split the data into a pre-intervention sample and a post-intervention sample. The pre-intervention sample starts in January 1992 and ends in June 2008. The post-intervention sample starts in July 2008 and ends Februray 2017. We give the code used to download and split our data below.

# Split the Data
data.RFSS.pre.intervention <- window(x = data.RFSS, end = 2008 + 5/12)
data.RFSS.post.intervention <- window(x = data.RFSS, start = 2008 + 6/12)

IV. Step 1 - Pre-Intervention Model

The first step of the intervention analysis approach is to specify a model over the pre-intervention sub-sample. To do so we follow a Box-Jenkins approach to select the best model.

Box-Jenkins Identification

In this section we examine the time series plot to determine if any transformations or differences need to made to produce an approximately weakly stationary time series. We also look at the auto-correlation function (ACF) and partial auto-correlation function (PACF) for each transformation to help us determine if any further transformations need to be made. Once all transformations are complete then we will use the final ACF and PACF to determine a set number of models and possible parameters that we will use in the estimation step of Box-Jenkins Methodology.

A. Logarithmic Transformation

We start by determining if we need to log transform our data. To do this we examine the simple time series of our original data as well as the time series of the log transformed data. Note, in this case, original data refers to the pre-intervention data, not the entire sample data.

original <- data.RFSS.pre.intervention
log.original <- log(data.RFSS.pre.intervention)

There is no evident reason why we should take the log from just examining the time series; however, the models we will be using can better use numerical methods approaches when we use the log-transformed data. We found that when using the original data (no log transformation) that many times the ARIMA models could not find numerical solutions. Furthermore, if we look at the original times series data it appears that there is more variation in the pre-intervention time period than in the post-intervention time period. Therefore, we will use the log-transformed data for the remainder of the paper.

B. Seasonal and Non-seasonal Differencing

For further examination of the data we look at the ACF and PACF of the log transformed data. Specifically we are looking for any evidence that might suggest the possibility of a unit root.

max.lag = 60    # Set to 60 because in months this is 5 years

Acf(log.original, type = 'correlation', lag = max.lag)
Acf(log.original, type = 'partial', lag = max.lag)

The slow decay of the ACF of the log transformed data suggests the presence of a unit root, and more importantly suggusts that the data is non-stationary. To fix this problem we need to take the difference of the data (first, second, third, etc.). To formally test which difference we should take we could run a series of ADF, KPSS, and ERS tests like we did in Homework Assignment #4. For the sake of time and space we will simply take the first difference of the data.

We should also point out that the time series of the log transformed data shows patterns of seasonality. We have monthly data and therefore we suspect that the seasonality can be fixed by differencing across the entire year (every twelve months).

We perform both of these differences on the log transformed data set separately and jointly as shown below.

first.difference <- diff(log.original, 1)
seasonal.difference <- diff(log.original, 12)
first.diff.seasonal.diff <- diff(diff(log.original, 12), 1)

We now plot the time series of the each of the above three data sets as well as the original log transformed data set. Ideally we want to find the time series that looks most like a weakly stationary time series (constant mean and constant co-variance with itself across time).

Out of the above time-series the data where we took the first difference model (third row) looks to be the most stationary out of all of the time-series. To more formally examine which time series is best we look at the ACF and PACF of each model.

From the second row marginally improves the slow decay of the ACF; however, there is still a slight decay and thus seaonsal differening alone does not alieviate the precense of a unit root. Seaonsal differencing alone cannot produce an approximately weakly stationary time series. The third row of the above figure shows us that first differencing appears to deal with the unit root problem but also shows us that seasonality is definitely present. It is important to note at this point that the ACF and PACF of the third row show good indications that our model is multiplicative seasonal AR model. We disregard the fourth row since the corresponding time series does not appear to have constant variance across time, and is therefore not weakly stationary.

C. Possible ARIMA Models

Using the ACF and PACFs above in Part B we now provide a list of possible ARIMA models that we think will be plausible. We write the models in the form:

\[ ARIMA(p, d, q)(P, D, Q)_{s} \]

where p denotes the non-seasonal AR order, d denotes the non-seasonal differencing order, q denotes the non-seasonal MA order, P denotes the seasonal AR order, D denotes the seasonal differencing order, Q denotes the seasonal MA order, and s denotes the time span of repeating seasonal pattern.

The Models we will be Estimating are:

  • Model 1: \(ARIMA(p, 1, 0)(1, 0, 0)_{12}\)
  • Model 2: \(ARIMA(p, 1, 0)(2, 0, 0)_{12}\)
  • Model 3: \(ARIMA(p, 1, 0)(1, 1, 0)_{12}\)
  • Model 4: \(ARIMA(p, 1, 0)(2, 1, 0)_{12}\)
  • Model 5: \(ARIMA(p, 1, 1)(1, 0, 0)_{12}\)
  • Model 6: \(ARIMA(p, 1, 1)(2, 0, 0)_{12}\)
  • Model 7: \(ARIMA(p, 1, 1)(1, 1, 0)_{12}\)
  • Model 8: \(ARIMA(p, 1, 1)(2, 1, 0)_{12}\)

We leave the parameter p as a character because we will actually be testing 24 models where we set p equal to 1, 4, 6, and 11.

Justification of the Model Choices:

We justify the above model choices using the ACF and PACF of the first difference of the log transformed data. In every model we take the first difference (d = 1) to correct for the unit root issue previously discussed.

The overall pattern of the ACF and PACF both suggest that our model is a multiplicative seasonal AR model. Therefore, we set both the non-seasonal and the seasonal MA component equal to zero (q = 0 and Q = 0). However, we do notice that the first lag of the ACF is highly statistically significant and therefore we also try setting the non-seasonal MA component equal to one (q = 0). Further, we see that the exponential decay of the ACF occurs every 12 lags starting at lag 12 hence we try setting the seasonal AR component to 1 (P = 1). Similarly, we see a smaller but similar pattern starting at lag 13, therefore we also attempt setting the seasonal AR component equal to 2 (P = 2). We also try all of our models with a first seasonal difference (D = 1).

To determine the non-seasonal AR component of the model we look at the PACF. From the PACF we see that the first four lags, and the sixth lags are statistically significant. Hence, we try setting the non-seasonal AR component of the model equal to 1, 4, and 6 (p = {1, 4, 6}). We also note that many of the lags of the PACF are statistically significant as expected in a multiplicative seasonal AR model, especially the eleventh lag, which appears to be abnormally significant. To ensure that these lags are not induced by the non-seasonl component we try up to lag 11 (p = 11). We do not go further than this due to the seasonal year effect for this being monthly data.

Box-Jenkins Estimation

In this section we estimate all 24 models suggested in the identificaiton section of Box-Jenkins methodology. To determine which model is best we look at the coefficients for statistical significance and look at the information criteria of each model. This is done in two rounds.

For the first round we only compare models in groups of eight as determined by their p value (the AR component). The best model of each p will be chosen based off of the AICc and BIC statistic; the lowest statistic determines the best model. Ideally we choose the model that has the lowest statistic for both. If no model has both the lowest AICc and BIC statistic then we will examine the residuals of the model to choose between the two.

For the second round we will the compare the residuals of each of the four models chosen from round one. The model that is most adequate when comparing the AICc and BIC statistics, the residuals, and the lowest in-sample esimtation errors will be the model chosen for forecasting.

A. Estimating the Model Parameters

Below is the code that we use to estimate the model parameters.

m1 <- Arima(log.original, order = c(p,1,0), seasonal = list(order = c(1,0,0), period = 12))
m2 <- Arima(log.original, order = c(p,1,0), seasonal = list(order = c(2,0,0), period = 12))
m3 <- Arima(log.original, order = c(p,1,0), seasonal = list(order = c(1,1,0), period = 12))
m4 <- Arima(log.original, order = c(p,1,0), seasonal = list(order = c(2,1,0), period = 12))

m5 <- Arima(log.original, order = c(p,1,1), seasonal = list(order = c(1,0,0), period = 12))
m6 <- Arima(log.original, order = c(p,1,1), seasonal = list(order = c(2,0,0), period = 12))
m7 <- Arima(log.original, order = c(p,1,1), seasonal = list(order = c(1,1,0), period = 12))
m8 <- Arima(log.original, order = c(p,1,1), seasonal = list(order = c(2,1,0), period = 12))

As mentioned before we replace the above p values with 1, 4, 6, and 11.

B. First Round - Top Four Models

Below we have constructed a table that shows us each model and their respective AICc and BIC statistic. For convience we have also sorted the fourth through sixth column on the AICc statistic, and the seventh through ninth column on the BIC statistic. We perform this table for each value of p, the AR component.

The first set of models we look at are the models with p = 1.

Model AICc BIC Sort on AICc AICc Sort on BIC BIC
m1 ARIMA(1,1,0)(1,0,0) -886.14 -876.41 m6 ARIMA(1,1,1)(2,0,0) -933.27 m6 ARIMA(1,1,1)(2,0,0) -917.17
m2 ARIMA(1,1,0)(2,0,0) -888.96 -876.03 m5 ARIMA(1,1,1)(1,0,0) -928.31 m5 ARIMA(1,1,1)(1,0,0) -915.39
m3 ARIMA(1,1,0)(1,1,0) -866.31 -856.78 m8 ARIMA(1,1,1)(2,1,0) -919.32 m8 ARIMA(1,1,1)(2,1,0) -903.55
m4 ARIMA(1,1,0)(2,1,0) -879.76 -867.10 m7 ARIMA(1,1,1)(1,1,0) -906.42 m7 ARIMA(1,1,1)(1,1,0) -893.76
m5 ARIMA(1,1,1)(1,0,0) -928.31 -915.39 m2 ARIMA(1,1,0)(2,0,0) -888.96 m1 ARIMA(1,1,0)(1,0,0) -876.41
m6 ARIMA(1,1,1)(2,0,0) -933.27 -917.17 m1 ARIMA(1,1,0)(1,0,0) -886.14 m2 ARIMA(1,1,0)(2,0,0) -876.03
m7 ARIMA(1,1,1)(1,1,0) -906.42 -893.76 m4 ARIMA(1,1,0)(2,1,0) -879.76 m4 ARIMA(1,1,0)(2,1,0) -867.10
m8 ARIMA(1,1,1)(2,1,0) -919.32 -903.55 m3 ARIMA(1,1,0)(1,1,0) -866.31 m3 ARIMA(1,1,0)(1,1,0) -856.78

From the above analysis we should pick model 6 - \(ARIMA(1,1,1)(2,0,0)\) - since, model 6 has both the lowest AICc and BIC statistics.

The second set of models we look at are the models with p = 4.

Model AICc BIC Sort on AICc AICc Sort on BIC BIC
m1 ARIMA(4,1,0)(1,0,0) -930.08 -910.82 m6 ARIMA(4,1,1)(2,0,0) -945.67 m2 ARIMA(4,1,0)(2,0,0) -921.81
m2 ARIMA(4,1,0)(2,0,0) -944.20 -921.81 m2 ARIMA(4,1,0)(2,0,0) -944.20 m6 ARIMA(4,1,1)(2,0,0) -920.17
m3 ARIMA(4,1,0)(1,1,0) -918.94 -900.09 m8 ARIMA(4,1,1)(2,1,0) -935.39 m4 ARIMA(4,1,0)(2,1,0) -913.13
m4 ARIMA(4,1,0)(2,1,0) -935.04 -913.13 m4 ARIMA(4,1,0)(2,1,0) -935.04 m1 ARIMA(4,1,0)(1,0,0) -910.82
m5 ARIMA(4,1,1)(1,0,0) -932.47 -910.08 m5 ARIMA(4,1,1)(1,0,0) -932.47 m8 ARIMA(4,1,1)(2,1,0) -910.45
m6 ARIMA(4,1,1)(2,0,0) -945.67 -920.17 m1 ARIMA(4,1,0)(1,0,0) -930.08 m5 ARIMA(4,1,1)(1,0,0) -910.08
m7 ARIMA(4,1,1)(1,1,0) -920.11 -898.20 m7 ARIMA(4,1,1)(1,1,0) -920.11 m3 ARIMA(4,1,0)(1,1,0) -900.09
m8 ARIMA(4,1,1)(2,1,0) -935.39 -910.45 m3 ARIMA(4,1,0)(1,1,0) -918.94 m7 ARIMA(4,1,1)(1,1,0) -898.20

From the above analysis it seems that we should pick either the sixth model - \(ARIMA(4,1,1)(2,0,0)\) - as chosen by the AICc statistic or the second model - \(ARIMA(4,1,0)(2,0,0)\) - as chosen by the BIC statistic. To pick between these two models we look at the residuals of each model.

tsdiag(m2, gof.lag = max.lag)

tsdiag(m6, gof.lag = max.lag)

From the residuals it appears that model 6 appears to perform better. Hence, we pick the \(ARIMA(4,1,1)(2,0,0)\) model.

The third set of models we look at are the models with p = 6.

Model AICc BIC Sort on AICc AICc Sort on BIC BIC
m1 ARIMA(6,1,0)(1,0,0) -931.20 -905.70 m2 ARIMA(6,1,0)(2,0,0) -945.77 m2 ARIMA(6,1,0)(2,0,0) -917.18
m2 ARIMA(6,1,0)(2,0,0) -945.77 -917.18 m6 ARIMA(6,1,1)(2,0,0) -943.80 m6 ARIMA(6,1,1)(2,0,0) -912.15
m3 ARIMA(6,1,0)(1,1,0) -921.00 -896.05 m4 ARIMA(6,1,0)(2,1,0) -936.50 m4 ARIMA(6,1,0)(2,1,0) -908.54
m4 ARIMA(6,1,0)(2,1,0) -936.50 -908.54 m8 ARIMA(6,1,1)(2,1,0) -934.84 m1 ARIMA(6,1,0)(1,0,0) -905.70
m5 ARIMA(6,1,1)(1,0,0) -929.06 -900.48 m1 ARIMA(6,1,0)(1,0,0) -931.20 m8 ARIMA(6,1,1)(2,1,0) -903.90
m6 ARIMA(6,1,1)(2,0,0) -943.80 -912.15 m5 ARIMA(6,1,1)(1,0,0) -929.06 m5 ARIMA(6,1,1)(1,0,0) -900.48
m7 ARIMA(6,1,1)(1,1,0) -919.22 -891.27 m3 ARIMA(6,1,0)(1,1,0) -921.00 m3 ARIMA(6,1,0)(1,1,0) -896.05
m8 ARIMA(6,1,1)(2,1,0) -934.84 -903.90 m7 ARIMA(6,1,1)(1,1,0) -919.22 m7 ARIMA(6,1,1)(1,1,0) -891.27

From the above analysis we should pick model 2 - \(ARIMA(6,1,0)(2,0,0)\) - since, model 2 has both the lowest AICc and BIC statistics.

The fourth set of models we look at are the models with p = 11.

Model AICc BIC Sort on AICc AICc Sort on BIC BIC
m1 ARIMA(11,1,0)(1,0,0) -950.81 -910.12 m2 ARIMA(11,1,0)(2,0,0) -962.10 m2 ARIMA(11,1,0)(2,0,0) -918.45
m2 ARIMA(11,1,0)(2,0,0) -962.10 -918.45 m6 ARIMA(11,1,1)(2,0,0) -960.15 m6 ARIMA(11,1,1)(2,0,0) -913.56
m3 ARIMA(11,1,0)(1,1,0) -938.78 -899.04 m4 ARIMA(11,1,0)(2,1,0) -952.73 m1 ARIMA(11,1,0)(1,0,0) -910.12
m4 ARIMA(11,1,0)(2,1,0) -952.73 -910.12 m8 ARIMA(11,1,1)(2,1,0) -951.30 m4 ARIMA(11,1,0)(2,1,0) -910.12
m5 ARIMA(11,1,1)(1,0,0) -948.57 -904.91 m1 ARIMA(11,1,0)(1,0,0) -950.81 m8 ARIMA(11,1,1)(2,1,0) -905.84
m6 ARIMA(11,1,1)(2,0,0) -960.15 -913.56 m5 ARIMA(11,1,1)(1,0,0) -948.57 m5 ARIMA(11,1,1)(1,0,0) -904.91
m7 ARIMA(11,1,1)(1,1,0) -936.72 -894.11 m3 ARIMA(11,1,0)(1,1,0) -938.78 m3 ARIMA(11,1,0)(1,1,0) -899.04
m8 ARIMA(11,1,1)(2,1,0) -951.30 -905.84 m7 ARIMA(11,1,1)(1,1,0) -936.72 m7 ARIMA(11,1,1)(1,1,0) -894.11

From the above analysis we should pick model 2 - \(ARIMA(11,1,0)(2,0,0)\) - since, model 2 has both the lowest AICc and BIC statistics.

We should note that every time model 2 and model 6 consistently performed well for every value of p. We should also point out that as we increased the number of AR lags the MA lag (from model 6) became less useful as shown by the fact that model 2 began to outpeform model 6 when the number of lags were higher than 4. The consistency should also give us confidence that we are headed in the right direction for choosing the best model.

C. Second Round - Final Model

From the previous section we have narrowed down our forecasting model to four choices. The four models are shown below.

m.p1 <- Arima(log.original, order = c(1,1,1), seasonal = list(order = c(2,0,0), period = 12))
m.p4 <- Arima(log.original, order = c(4,1,1), seasonal = list(order = c(2,0,0), period = 12))
m.p6 <- Arima(log.original, order = c(6,1,0), seasonal = list(order = c(2,0,0), period = 12))
m.p11 <- Arima(log.original, order = c(11,1,0), seasonal = list(order = c(2,0,0), period = 12))

We now look at each model’s coefficients and their standard errors. Specifically we are looking to see which coefficients are statistically significant.

m.p1
## Series: log.original 
## ARIMA(1,1,1)(2,0,0)[12]                    
## 
## Coefficients:
##           ar1      ma1    sar1    sar2
##       -0.1251  -0.6886  0.7686  0.2040
## s.e.   0.0893   0.0586  0.0723  0.0751
## 
## sigma^2 estimated as 0.0004186:  log likelihood=471.79
## AIC=-933.59   AICc=-933.27   BIC=-917.17

Every coefficient is statistically significant in the \(ARIMA(1,1,1)(2,0,0)\) model.

m.p4
## Series: log.original 
## ARIMA(4,1,1)(2,0,0)[12]                    
## 
## Coefficients:
##           ar1      ar2     ar3      ar4      ma1    sar1    sar2
##       -0.3582  -0.3296  0.0179  -0.2304  -0.4260  0.6698  0.3086
## s.e.   0.1801   0.1541  0.1245   0.0744   0.1763  0.0736  0.0751
## 
## sigma^2 estimated as 0.000381:  log likelihood=481.22
## AIC=-946.43   AICc=-945.67   BIC=-920.17

All but one coefficient (the third lag) are statistically significant in the \(ARIMA(4,1,1)(2,0,0)\) model. Hence, we will try a restricted model as well.

m.p6
## Series: log.original 
## ARIMA(6,1,0)(2,0,0)[12]                    
## 
## Coefficients:
##           ar1      ar2      ar3      ar4     ar5     ar6    sar1    sar2
##       -0.7861  -0.6348  -0.2623  -0.3165  -0.133  0.0449  0.6599  0.3197
## s.e.   0.0716   0.0927   0.1024   0.1011   0.093  0.0738  0.0728  0.0742
## 
## sigma^2 estimated as 0.0003768:  log likelihood=482.37
## AIC=-946.73   AICc=-945.77   BIC=-917.18

Surprisingly all but two coefficients are statistically significant in the \(ARIMA(6,1,0)(2,0,0)\) model, namely the fifth and sixth lags. This implies that this model simplifies to an \(ARIMA(4,1,0)(2,0,0)\) model if we restrict the last two lags to equal zero. Therefore, since we have already shown that the \(ARIMA(4,1,1)(2,0,0)\) model appears to perform better than the \(ARIMA(4,1,0)(2,0,0)\) model then we will no longer regard this model.

m.p11
## Series: log.original 
## ARIMA(11,1,0)(2,0,0)[12]                    
## 
## Coefficients:
##           ar1      ar2      ar3      ar4      ar5     ar6     ar7     ar8
##       -0.7431  -0.6344  -0.2808  -0.3407  -0.0896  0.0224  0.0717  0.0887
## s.e.   0.0705   0.0902   0.0985   0.0993   0.1033  0.1015  0.1013  0.0980
##          ar9    ar10    ar11    sar1    sar2
##       0.3939  0.1345  0.1621  0.6590  0.3271
## s.e.  0.0970  0.0947  0.0818  0.0832  0.0835
## 
## sigma^2 estimated as 0.0003286:  log likelihood=496.21
## AIC=-964.41   AICc=-962.1   BIC=-918.45

In the above model the fifth, sixth, seventh, and eight lags are not significant. The tenth and eleventh lags also may not be significant. Hence, we will try two different restrictions for this model (especially since we have already eliminated one model).

Below we give the code used to estimate our restricted models.

m.p4.restricted <- Arima(log.original, order = c(4,1,1), seasonal = list(order = c(2,0,0), period = 12), 
                       fixed = c(NA, NA, 0, 0, NA, NA, NA))

m.p11.restricted.a <- Arima(log.original, order = c(11,1,0), seasonal = list(order = c(2,0,0), period = 12), 
                        fixed = c(NA, NA, NA, NA, 0, 0, 0, 0, NA, NA, NA, NA, NA))

m.p11.restricted.b <- Arima(log.original, order = c(11,1,0), seasonal = list(order = c(2,0,0), period = 12), 
                        fixed = c(NA, NA, NA, NA, 0, 0, 0, 0, NA, 0, 0, NA, NA))

We again check to see which of the restricted models are statistically significant.

m.p4.restricted
## Series: log.original 
## ARIMA(4,1,1)(2,0,0)[12]                    
## 
## Coefficients:
##           ar1      ar2  ar3  ar4      ma1    sar1    sar2
##       -0.4030  -0.3076    0    0  -0.4253  0.7071  0.2691
## s.e.   0.2205   0.1498    0    0   0.2380  0.0819  0.0849
## 
## sigma^2 estimated as 0.0004075:  log likelihood=475.21
## AIC=-938.41   AICc=-937.97   BIC=-918.71

As expected all lags are now significant.

m.p11.restricted.a
## Series: log.original 
## ARIMA(11,1,0)(2,0,0)[12]                    
## 
## Coefficients:
##           ar1      ar2      ar3      ar4  ar5  ar6  ar7  ar8     ar9
##       -0.7127  -0.5991  -0.2215  -0.2728    0    0    0    0  0.3390
## s.e.   0.0682   0.0851   0.0835   0.0677    0    0    0    0  0.0687
##         ar10    ar11    sar1    sar2
##       0.0761  0.1266  0.6336  0.3525
## s.e.  0.0839  0.0778  0.0788  0.0792
## 
## sigma^2 estimated as 0.0003332:  log likelihood=494.78
## AIC=-969.57   AICc=-968.39   BIC=-936.74

All lags are significant except for the tenth and elventh lag as we mentioned might be the case. Hence we now look to the second restricted model.

m.p11.restricted.b
## Series: log.original 
## ARIMA(11,1,0)(2,0,0)[12]                    
## 
## Coefficients:
##           ar1      ar2      ar3      ar4  ar5  ar6  ar7  ar8     ar9  ar10
##       -0.7198  -0.5670  -0.2162  -0.2785    0    0    0    0  0.2975     0
## s.e.   0.0653   0.0834   0.0845   0.0673    0    0    0    0  0.0531     0
##       ar11    sar1    sar2
##          0  0.5876  0.3959
## s.e.     0  0.0707  0.0718
## 
## sigma^2 estimated as 0.0003396:  log likelihood=493.44
## AIC=-970.88   AICc=-970.12   BIC=-944.62

As expected all lags are statistically significant.

We now compare the infomration criterian of each of the statistically significant models.

Model AICc BIC
m.p1 ARIMA(1,1,1)(2,0,0) -933.27 -917.17
m.p4.restricted ARIMA(4,1,1)(2,0,0) - Restricted -937.65 -912.15
m.p11.restricted.b ARIMA(11,1,0)(2,0,0) - Restricted -968.57 -924.92

Even with as many parameters as the model has the restricted \(ARIMA(11,1,0)(2,0,0)\) model seems to perform the best because it has both a lower AICc statistic and a lower BIC statistic. To further check we also look at the residual plots of each of the three models above.

tsdiag(m.p1, gof.lag = max.lag)

tsdiag(m.p4.restricted, gof.lag = max.lag)

tsdiag(m.p11.restricted.b, gof.lag = max.lag)

The residual plots also conclude that the restricted \(ARIMA(11,1,0)(2,0,0)\) far outperforms the other two models. To further the strength of each model we also look at the in-sample accuracy measures of each model below.

ME RMSE MAE MPE MAPE MASE ACF1
ARIMA(1,1,1)(2,0,0) 0.00055 0.02020 0.01606 0.00469 0.12907 0.30607 -0.00641
ARIMA(4,1,1)(2,0,0) - Restricted 0.00033 0.01978 0.01558 0.00285 0.12517 0.29684 0.03172
ARIMA(11,1,0)(2,0,0) - Restricted 0.00011 0.01776 0.01398 0.00108 0.11231 0.26628 -0.00902

For every single type of accuracy measure the restricted \(ARIMA(11,1,0)(2,0,0)\) outperforms the other two models.

In conclusion, we will use the restricted \(ARIMA(11,1,0)(2,0,0)\) model to forecast our data.

V. Step 2 - Forecasting Pre-Intervention Model

Using the restricted \(ARIMA(11,1,0)(2,0,0)\) model we now forecast this model onto the post-intervention data using the code below.

multistep.forecast <- forecast(m.p11.restricted.b, length(data.RFSS.post.intervention))
multistep.forecast.error <- as.ts(log(data.RFSS.post.intervention)) - multistep.forecast$mean

We now graph the forecast next to the error from the forecast, as well as the ACF and PACF of the residuals from the forecast.

It appears that there is a gradual decaying effect with a zero effect in the long-run. This now helps us to form the pulse and step functions.

VI. Step 3 - Pulse and Step Functions

Before we create our pulse and step functions we first split our sample into an estimation sub-sample and a prediction sub-sample.

start.month <- 1992
end.month <- 2013 + (11/12)
data.RFSS.estimation <- window(data.RFSS, end = end.month)
data.RFSS.prediction <- window(data.RFSS, start = end.month + (1/12))

From Step 2 we believe that the intervention is temporary and decaying across time. Therefore, we think that the change in the retail and food servicse data is from an intervention with a pulse function. We model the pulse function in R below.

pulse.RFSS <- 1*(index(data.RFSS.estimation)==2008+(7-1)/12)

VII. Step 4 - Estimating the Model

We then estimate our model with a transfer function that assumes that the effect of the 2008 recession has a gradually disappearing effect on the retail and food services sales. The code to produce this model is given below.

m <- arimax(log(data.RFSS.estimation), order = c(11,1,0), seasonal = list(order = c(2,0,0), period = 12), 
            fixed = c(NA, NA, NA, NA, 0, 0, 0, 0, NA, 0, 0, NA, NA, NA, NA),  
            xtransf = data.frame(pulse.RFSS), transfer=list(c(1,0)))

VIII. Step 5 - Adequacy of the Model

A. Residuals and Significance

To examine the adequacy of the model we first examine the residuals below.

tsdiag(m, gof.lag = max.lag)

As we can see the residuals perform fairly well. We next check the statistical significance of each coefficient. We will not adjust the previous coefficients since these were determine signifcant in the pre and post intervention analysis. We only examine the coefficeints of the two additional terms from the pulse function.

m
## 
## Call:
## arimax(x = log(data.RFSS.estimation), order = c(11, 1, 0), seasonal = list(order = c(2, 
##     0, 0), period = 12), fixed = c(NA, NA, NA, NA, 0, 0, 0, 0, NA, 0, 0, NA, 
##     NA, NA, NA), xtransf = data.frame(pulse.RFSS), transfer = list(c(1, 0)))
## 
## Coefficients:
##           ar1      ar2     ar3      ar4  ar5  ar6  ar7  ar8     ar9  ar10
##       -0.5429  -0.2927  0.0817  -0.0956    0    0    0    0  0.2405     0
## s.e.   0.0607   0.0715  0.0702   0.0609    0    0    0    0  0.0514     0
##       ar11    sar1    sar2  pulse.RFSS-AR1  pulse.RFSS-MA0
##          0  0.5914  0.3901          0.5939          0.0292
## s.e.     0  0.0600  0.0609          0.1542          0.0149
## 
## sigma^2 estimated as 0.0004021:  log likelihood = 635.11,  aic = -1252.22

Both coefficients for the pulse functions are marginally significant.

B. Forecasting

The real evaluation of the adequacy of the model is to perform a forecast of the model. We will run two forecasts to check how accurate our model is:

  1. Multistep Forecast Model
  2. One Month ahead Forecasts

Both forecasts are forecasted over the period between January 2014 and December 2016. To perform our forecasts we extend our transfer function across time to include our prediction sample period using the following code.

pulse.RFSS <- c(pulse.RFSS, rep(0,length(data.RFSS.prediction)))

transfer.function <- m$coef["pulse.RFSS-MA0"]*filter(pulse.RFSS, filter = m$coef["pulse.RFSS-AR1"], method = 'recursive')

model.x <- Arima(log(data.RFSS.estimation), order = c(11,1,0), 
                 seasonal = list(order = c(2,0,0), period = 12), 
                 fixed = c(NA, NA, NA, NA, 0, 0, 0, 0, NA, 0, 0, NA, NA, NA),
                 xreg = transfer.function[1:length(data.RFSS.estimation)])

Forecast 1 - Multistep Forecast

In a multi-step forecast the model parameters for the restricted \(ARIMA(11,1,0)(2,0,0)\) are estimated once over the estimation sub-sample. The parameters estimated for the estimation sub-sample are then used across the entire prediction sample. The code used to produce the multi-step forecast is given below. We also show a plot of the multi-step forecast, along with the actual values over the prediction period.

multistep.forecast <- forecast(model.x, length(data.RFSS.prediction), 
                               xreg = transfer.function[length(data.RFSS.estimation+1):
                                                          (length(data.RFSS.estimation) + 
                                                             length(data.RFSS.prediction))])

From the plot alone we can see that the multi-step forecast performs extremely well since the estimated data are close to the actual data. In other words, the restricted and transfered \(ARIMA(11,1,0)(2,0,0)\) model produces data that looks like the actual data.

Forecast 2 - Rolling Forecast

In a rolling forecast the estimation sample moves across time as time progresses. The estimation sample updates to contain the same amount of data that were in the original model, and each update contains new parameters for the \(ARIMA(3,1,0)(3,0,0)\) model to forecast the one month ahead Total Private Residential Construction Spending. The model is re-estimated again and again until we have forecasted the entirety of the prediction sample. The code used to produce the rolling forecast is given below. We also show a plot of the rolling forecast, along with the actual values over the prediction period.

My attempt for the rolling scheme forecast is shown below. However, I was unable to get the code working properly. As a result this forecast will not be performed.

Rolling.Forecast <- function(DATA, first.month, last.month){
  library(TSA)
  pulse.RFSS <- 1*(index(data.RFSS.estimation)==2008+(7-1)/12)
  
  rolling.forecast <- list()
  log.full.data <- log(DATA)
  
  for(i in 1:length(data.RFSS.prediction)) {
    temp <- window(log.full.data, start = first.month + (i-1)/12, end = last.month + (i-1) / 12)
    model.update <- arimax(temp, order = c(11,1,0), seasonal = list(order = c(2,0,0), period = 12), 
                           fixed = c(NA, NA, NA, NA, 0, 0, 0, 0, NA, 0, 0, NA, NA, NA, NA),  
                           xtransf = data.frame(pulse.RFSS), transfer=list(c(1,0)))
    pulse.RFSS <- c(pulse.RFSS, 0)
    transfer.function <- m$coef["pulse.RFSS-MA0"]*filter(pulse.RFSS, filter = m$coef["pulse.RFSS-AR1"], method = 'recursive')
    model.x.update <- Arima(temp, order = c(11,1,0), 
                            seasonal = list(order = c(2,0,0), period = 12), 
                            fixed = c(NA, NA, NA, NA, 0, 0, 0, 0, NA, 0, 0, NA, NA, NA),
                            xreg = transfer.function[1:length(data.RFSS.estimation)])
    
    rolling.forecast.update <- forecast(model.x.update, 1, xreg = transfer.function[(length(data.RFSS.estimation)+1)])
    rolling.forecast$mean <- rbind(rolling.forecast$mean , as.zoo(rolling.forecast.update$mean))
    rolling.forecast$lower <- rbind(rolling.forecast$lower, rolling.forecast.update$lower)
    rolling.forecast$upper <- rbind(rolling.forecast$upper, rolling.forecast.update$upper)
    
  }
  
  rolling.forecast <- as.ts(rolling.forecast)
  
  return(rolling.forecast)
}

rolling.forecast <- Rolling.Forecast(data.RFSS, 1992, 2013 + (11/12))

Comparing the Accuracy of the Forecasts

Since I was unable to get one of the forecasting methods to work, I will only show the accuracy of the multi-step forecast using the following code.

accuracy(f = multistep.forecast, x = log(data.RFSS.prediction))

In the table below ME stands for Mean Error, RMSE stands for Relative Mean Squared Error, MAE stands for Mean Absolute Error, MPE stands for Mean Percent Error, and lastly MAPE stands for Mean Absolute Percentage Error. All of the above measures quantify the out-of-sample accuracy of each forecast.

ME RMSE MAE MPE MAPE MASE
Multi-Step Forecast - In-Sample 0.0002133 0.0200281 0.0153804 0.0017527 0.1225455 0.2844621
Multi-Step Forecast - Out-of-Sample -0.0061198 0.0198774 0.0172319 -0.0475118 0.1327047 0.3187044

Overall, the multi-step forecast of the restricted \(ARIMA(11,1,0)(2,0,0)\) model appears to work extermely well in forecasting RFSS data even with the recession occuring.

IX. Conclusion

We recommend that a restricted \(ARIMA(11,0,0)(2,0,0)\) model with a decaying transfer function with zero permanent effect be used when modeling or forecasting retail and food service sales data. This model specifically accounts for the affect of the 2008 recession.