HW 2 Box-Jenkins Chapter Key
EE Holmes
1/28/2019
Here are the answers for the homework problems based on the material in Box-Jenkins method in R. You can download the Rmd file for this key from here.
Solutions Chapter 5
Problem 1
Because the information on how to do a Dickey-Fuller test was buried in the material that I gave you and not in the help file for these functions, I did not grade problem 1.
Repeat the stationarity tests for another species in the landings data set. Sardine, Chub.mackerel, and Horse.mackerel have enough data. Here is how to set up the data for another species.
load("landings.RData")
datdf <- subset(landings, Species=="Chub.mackerel" & Year<=1989)
dat <- ts(datdf$log.metric.tons, start=1964)Do a Dickey-Fuller test using
ur.df()andadf.test(). What do you need to set the lags to to achieve that?You need to set the lags to 0 to do a Dickey-Fuller test. The Dickey-Fuller test assumes an AR(1) process for the stationary process while an Augmented Dickey-Fuller test allows AR(\(k+1\)) where \(k\) is the lag value that you pass into the test function. The \(k\) is referring to the \(\delta\) lags in the test equation: \[\Delta y_t = \alpha + \beta t + \gamma y_{t-1} + \delta_1 \Delta y_{t-1} + \delta_2 \Delta y_{t-2} + \dots\] Setting lags = 0 will be this model and we are testing if \(\gamma = 0\).
\[\Delta y_t = \alpha + \beta t + \gamma y_{t-1}\]
tseries::adf.test(dat, k=0)## ## Augmented Dickey-Fuller Test ## ## data: dat ## Dickey-Fuller = -2.8932, Lag order = 0, p-value = 0.2321 ## alternative hypothesis: stationaryThe null hypothesis of non-stationarity is not rejected, i.e. we cannot reject that \(\gamma=0\).
Getting the information on whether the null is rejected is harder for
ur.df().library(urca) test <- ur.df(dat, lags=0, type="trend") summary(test)## ## ############################################### ## # Augmented Dickey-Fuller Test Unit Root Test # ## ############################################### ## ## Test regression trend ## ## ## Call: ## lm(formula = z.diff ~ z.lag.1 + 1 + tt) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.35165 -0.13534 -0.09449 0.10013 0.67696 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3.34754 1.18997 2.813 0.01013 * ## z.lag.1 -0.49845 0.17228 -2.893 0.00844 ** ## tt 0.02991 0.01027 2.913 0.00806 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.2637 on 22 degrees of freedom ## Multiple R-squared: 0.3104, Adjusted R-squared: 0.2478 ## F-statistic: 4.952 on 2 and 22 DF, p-value: 0.01676 ## ## ## Value of test-statistic is: -2.8932 3.4645 4.9522 ## ## Critical values for test statistics: ## 1pct 5pct 10pct ## tau3 -4.15 -3.50 -3.18 ## phi2 7.02 5.13 4.31 ## phi3 9.31 6.73 5.61Ignore the part in the summary with the p-value. You need to look at the test-statistic and see if the
taustatistic is greater than the critical value at 5pct; if greater it means you do not reject the null of non-stationarity. The first test statistic istau3, the secondphi2and the third isphi3. The test statistic fortau3is -2.89 and is greater than -3.50. Thus the null of \(\gamma=0\) is not rejected at \(\alpha=0.05\).For our purposes, we just need to focus on the
tau. Thephiare for other null hypotheses involving the presence or absence of the drift and trend terms.Do an Augmented Dickey-Fuller test using
ur.df(). How did you choose to set the lags?There are a few different approaches to use to set the lags. You could use the rule of thumb that
adf.test()uses which istrunc((length(x)-1)^(1/3))which is 2 for 26 data points.urca::ur.df(dat, lags=2, type="trend")## ## ############################################################### ## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # ## ############################################################### ## ## The value of the test statistic is: -1.7704 1.8467 1.6961You could also let
ur.df()select the best lag based on AIC or BIC. It chooses 2 also. Look for the number ofz.diff.lags in the summary table.library(urca) test<-ur.df(dat, lags=8, type="trend", selectlags="AIC") summary(test)## ## ############################################### ## # Augmented Dickey-Fuller Test Unit Root Test # ## ############################################### ## ## Test regression trend ## ## ## Call: ## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.28502 -0.08179 0.01515 0.09477 0.26411 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 7.93107 1.97731 4.011 0.00173 ** ## z.lag.1 -1.27315 0.31986 -3.980 0.00183 ** ## tt 0.09858 0.02605 3.785 0.00260 ** ## z.diff.lag1 0.55866 0.26393 2.117 0.05586 . ## z.diff.lag2 0.59190 0.23347 2.535 0.02617 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.1744 on 12 degrees of freedom ## Multiple R-squared: 0.5818, Adjusted R-squared: 0.4424 ## F-statistic: 4.174 on 4 and 12 DF, p-value: 0.02402 ## ## ## Value of test-statistic is: -3.9804 6.7354 7.9258 ## ## Critical values for test statistics: ## 1pct 5pct 10pct ## tau3 -4.15 -3.50 -3.18 ## phi2 7.02 5.13 4.31 ## phi3 9.31 6.73 5.61Do a KPSS test using
kpss.test().tseries::kpss.test(dat, null="Trend")## ## KPSS Test for Trend Stationarity ## ## data: dat ## KPSS Trend = 0.15569, Truncation lag parameter = 2, p-value = ## 0.04193The null of trend stationary is rejected at \(\alpha=0.05\). Again the data do not pass the stationarity tests.
Problem 2
Repeat the stationarity tests and differencing tests for anchovy using the full data. The full data go to 2007.
load("landings.RData")
datdf <- subset(landings, Species=="Anchovy")
anchovy <- ts(datdf$log.metric.tons, start=1964)Do the conclusions regarding stationarity and the amount of differencing needed change?
tseries::adf.test(anchovy)## ## Augmented Dickey-Fuller Test ## ## data: anchovy ## Dickey-Fuller = -1.5069, Lag order = 3, p-value = 0.769 ## alternative hypothesis: stationarytseries::kpss.test(anchovy)## Warning in tseries::kpss.test(anchovy): p-value smaller than printed p- ## value## ## KPSS Test for Level Stationarity ## ## data: anchovy ## KPSS Level = 0.82215, Truncation lag parameter = 3, p-value = 0.01forecast::ndiffs(anchovy, test="kpss")## [1] 1forecast::ndiffs(anchovy, test="adf")## [1] 1In the chapter, we tested the anchovy catch data up to 1989 and needed 1 difference to pass the stationarity tests. With the longer time series, the null is not rejected in the ADF test, is rejected in the KPSS test, and one difference is required to pass the KPSS and ADF tests. So the conclusions are the same.
However notice that the time series is fairly different after 1987.
plot(anchovy)
Repeat with only the data after 1987. Are the conclusions different for only the recent data?
anchovy88.07 <- window(anchovy, start=1988)tseries::adf.test(anchovy88.07)## ## Augmented Dickey-Fuller Test ## ## data: anchovy88.07 ## Dickey-Fuller = -2.1734, Lag order = 2, p-value = 0.5063 ## alternative hypothesis: stationarytseries::kpss.test(anchovy88.07)## Warning in tseries::kpss.test(anchovy88.07): p-value greater than printed ## p-value## ## KPSS Test for Level Stationarity ## ## data: anchovy88.07 ## KPSS Level = 0.24418, Truncation lag parameter = 2, p-value = 0.1forecast::ndiffs(anchovy88.07, test="kpss")## [1] 0forecast::ndiffs(anchovy88.07, test="adf")## [1] 1With only the data after 1987, the null is not rejected in the ADF test nor in the KPSS test. So for the KPSS, we cannot reject (at \(\alpha=0.05\)) that the data are stationary while with ADF, we do reject that the data are non-stationary. Correspondingly, no difference is required to pass the KPSS while one difference is required to pass the ADF test. Basically the data are on the borderline of appearing stationary and a conservative approach would use one difference.
Fit the using
auto.arima()to the full data (1964-2007) and recent data (1987 to present). Do the selected models change?anchovy87.07 <- window(anchovy, start=1987) fit1 <- forecast::auto.arima(anchovy, stepwise=FALSE) fit2 <- forecast::auto.arima(anchovy87.07, stepwise=FALSE) fit1## Series: anchovy ## ARIMA(0,1,1) ## ## Coefficients: ## ma1 ## -0.3112 ## s.e. 0.1548 ## ## sigma^2 estimated as 0.04299: log likelihood=7.1 ## AIC=-10.19 AICc=-9.89 BIC=-6.67fit2## Series: anchovy87.07 ## ARIMA(0,0,1) with non-zero mean ## ## Coefficients: ## ma1 mean ## 0.6684 9.5810 ## s.e. 0.2284 0.0647 ## ## sigma^2 estimated as 0.0361: log likelihood=5.83 ## AIC=-5.66 AICc=-4.25 BIC=-2.53Both are MA(1), with different parameter values, however the differencing changes so we are fitting to different data. If we difference the 1987 to 2007 data and fit, the model choose is white noise (0,0,0) which would be a ARIMA(0,1,0) for the undifferenced data. ARIMA(0,1,0) is a random walk.
fit3 <- forecast::auto.arima(diff(anchovy87.07), stepwise=FALSE) fit3## Series: diff(anchovy87.07) ## ARIMA(0,0,0) with zero mean ## ## sigma^2 estimated as 0.04109: log likelihood=3.54 ## AIC=-5.08 AICc=-4.86 BIC=-4.09We cannot directly compare the AICc from fitting an ARIMA(0,1,0) to a ARIMA(0,0,1) with
Arima()because the data have changed. In one case we are fitting to the first differences and in the other we are not. Instead let’s fit an AR(1) and compare to the MA(1) which was the best fit to the undifferenced data.fitar1 <- forecast::Arima(anchovy87.07, order=c(1,0,0)) fitma1 <- forecast::Arima(anchovy87.07, order=c(0,0,1)) tab <- data.frame(model=c("AR(1)","MA(1) with mean"), AICc=c(fitar1$aicc, fitma1$aicc) ) tab## model AICc ## 1 AR(1) -3.519512 ## 2 MA(1) with mean -4.252483The AICc are very similar (<1 difference) though the estimated \(\theta\) for the AR(1) is nowhere close to 1. 1 would mean a random walk.
fitar1## Series: anchovy87.07 ## ARIMA(1,0,0) with non-zero mean ## ## Coefficients: ## ar1 mean ## 0.5764 9.6055 ## s.e. 0.2141 0.0943 ## ## sigma^2 estimated as 0.03772: log likelihood=5.47 ## AIC=-4.93 AICc=-3.52 BIC=-1.8
Problem 3
Fit each of the models listed under the auto.arima() trace using Arima() and show that you can produce the same AICc value that is shown in the trace table.
forecast::auto.arima(anchovy, trace=TRUE)##
## ARIMA(2,1,2) with drift : -3.451482
## ARIMA(0,1,0) with drift : -7.274254
## ARIMA(1,1,0) with drift : -7.565745
## ARIMA(0,1,1) with drift : -9.062605
## ARIMA(0,1,0) : -8.876343
## ARIMA(1,1,1) with drift : -7.033147
## ARIMA(0,1,2) with drift : -7.292236
## ARIMA(1,1,2) with drift : -5.016628
## ARIMA(0,1,1) : -9.892016
## ARIMA(1,1,1) : -7.81313
## ARIMA(0,1,2) : -8.0402
## ARIMA(1,1,0) : -8.861307
## ARIMA(1,1,2) : -5.943474
##
## Best model: ARIMA(0,1,1)## Series: anchovy
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## -0.3112
## s.e. 0.1548
##
## sigma^2 estimated as 0.04299: log likelihood=7.1
## AIC=-10.19 AICc=-9.89 BIC=-6.67The AICc is output in `fit$aicc`.
```r
aiccs <- c(
forecast::Arima(anchovy, order=c(2,1,2), include.drift=TRUE)$aicc,
forecast::Arima(anchovy, order=c(0,1,0), include.drift=TRUE)$aicc,
forecast::Arima(anchovy, order=c(1,1,0), include.drift=TRUE)$aicc,
forecast::Arima(anchovy, order=c(0,1,1), include.drift=TRUE)$aicc,
forecast::Arima(anchovy, order=c(0,1,0), include.drift=FALSE)$aicc,
forecast::Arima(anchovy, order=c(1,1,1), include.drift=TRUE)$aicc,
forecast::Arima(anchovy, order=c(0,1,2), include.drift=TRUE)$aicc,
forecast::Arima(anchovy, order=c(1,1,2), include.drift=TRUE)$aicc,
forecast::Arima(anchovy, order=c(0,1,1), include.drift=FALSE)$aicc,
forecast::Arima(anchovy, order=c(1,1,1), include.drift=FALSE)$aicc,
forecast::Arima(anchovy, order=c(0,1,2), include.drift=FALSE)$aicc,
forecast::Arima(anchovy, order=c(1,1,2), include.drift=FALSE)$aicc )
models <- c("(2,1,2) with drift","(0,1,0) with drift","(1,1,0) with drift",
"(0,1,1) with drift","(0,1,0) no drift","(1,1,1) with drift",
"(0,1,2) with drift","(1,1,2) with drift","(0,1,1) no drift",
"(1,1,1) no drift","(0,1,2) no drift","(1,1,2) no drift"
)
data.frame(model=models, AICc=aiccs)
```
```
## model AICc
## 1 (2,1,2) with drift -3.451482
## 2 (0,1,0) with drift -7.274254
## 3 (1,1,0) with drift -7.565745
## 4 (0,1,1) with drift -9.062605
## 5 (0,1,0) no drift -8.876343
## 6 (1,1,1) with drift -7.033147
## 7 (0,1,2) with drift -7.292236
## 8 (1,1,2) with drift -5.016628
## 9 (0,1,1) no drift -9.892016
## 10 (1,1,1) no drift -7.813130
## 11 (0,1,2) no drift -8.040200
## 12 (1,1,2) no drift -5.943474
```Problem 4
Use auto.arima() with stepwise=FALSE to find the top 3 models (according to AICc) for the anchovy.
We fit using trace=TRUE and scan to find the 3 models with the lowest AICc. They are ARIMA(0,1,1), ARIMA(0,1,1) with drift, and ARIMA(0,1,0). The all are within 1 of the lowest AICc so have similar data support.
forecast::auto.arima(anchovy, stepwise=FALSE, trace=TRUE)Use
Arima()to fit the top 3 models to the data.The top 3 models are all random walks. The first has no drift but autocorrelated errors, the 2nd and 3rd have independent errors but the 2nd has drift and the 3rd does not.
fit1 <- forecast::Arima(anchovy, order=c(0,1,1), include.drift=FALSE) fit2 <- forecast::Arima(anchovy, order=c(0,1,0), include.drift=TRUE) fit3 <- forecast::Arima(anchovy, order=c(0,1,0), include.drift=FALSE)Create a 5-year forecast for each of the top 3 models using the fits from (a).
fr1 <- forecast::forecast(fit1, h=5) fr2 <- forecast::forecast(fit2, h=5) fr3 <- forecast::forecast(fit3, h=5)Do the forecasts differ?
The forecast from the model with drift has a trend so that is definitely different.
plot(fr1)
plot(fr2)
plot(fr3)
The first and 3rd models differ in having (or not) the MA(1) term which is autocorrelation in the errors. It is hard to see but this leads to the forecasts from the first model having smaller prediction intervals. You can see this is you compare 20-year forecasts. Compare the high prediction intervals for 2027.
fr1 <- forecast::forecast(fit1, h=20) fr3 <- forecast::forecast(fit3, h=20) fr1## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 ## 2008 9.625964 9.360238 9.891690 9.219571 10.03236 ## 2009 9.625964 9.303303 9.948626 9.132496 10.11943 ## 2010 9.625964 9.255006 9.996923 9.058632 10.19330 ## 2011 9.625964 9.212310 10.039619 8.993334 10.25859 ## 2012 9.625964 9.173626 10.078303 8.934172 10.31776 ## 2013 9.625964 9.137999 10.113930 8.879686 10.37224 ## 2014 9.625964 9.104802 10.147127 8.828915 10.42301 ## 2015 9.625964 9.073597 10.178332 8.781191 10.47074 ## 2016 9.625964 9.044062 10.207866 8.736022 10.51591 ## 2017 9.625964 9.015956 10.235973 8.693037 10.55889 ## 2018 9.625964 8.989089 10.262839 8.651948 10.59998 ## 2019 9.625964 8.963311 10.288618 8.612523 10.63941 ## 2020 9.625964 8.938498 10.313431 8.574575 10.67735 ## 2021 9.625964 8.914550 10.337378 8.537950 10.71398 ## 2022 9.625964 8.891383 10.360546 8.502519 10.74941 ## 2023 9.625964 8.868924 10.383004 8.468171 10.78376 ## 2024 9.625964 8.847113 10.404816 8.434814 10.81712 ## 2025 9.625964 8.825896 10.426033 8.402365 10.84956 ## 2026 9.625964 8.805227 10.446702 8.370755 10.88117 ## 2027 9.625964 8.785066 10.466862 8.339922 10.91201fr3## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 ## 2008 9.681593 9.408639 9.954547 9.264146 10.09904 ## 2009 9.681593 9.295578 10.067609 9.091234 10.27195 ## 2010 9.681593 9.208823 10.154364 8.958553 10.40463 ## 2011 9.681593 9.135685 10.227501 8.846699 10.51649 ## 2012 9.681593 9.071249 10.291937 8.748153 10.61503 ## 2013 9.681593 9.012995 10.350192 8.659060 10.70413 ## 2014 9.681593 8.959425 10.403762 8.577132 10.78605 ## 2015 9.681593 8.909562 10.453624 8.500874 10.86231 ## 2016 9.681593 8.862731 10.500456 8.429251 10.93394 ## 2017 9.681593 8.818437 10.544750 8.361509 11.00168 ## 2018 9.681593 8.776307 10.586880 8.297077 11.06611 ## 2019 9.681593 8.736052 10.627134 8.235513 11.12767 ## 2020 9.681593 8.697443 10.665743 8.176466 11.18672 ## 2021 9.681593 8.660292 10.702894 8.119648 11.24354 ## 2022 9.681593 8.624447 10.738740 8.064827 11.29836 ## 2023 9.681593 8.589777 10.773410 8.011804 11.35138 ## 2024 9.681593 8.556175 10.807012 7.960414 11.40277 ## 2025 9.681593 8.523547 10.839639 7.910514 11.45267 ## 2026 9.681593 8.491814 10.871373 7.861983 11.50120 ## 2027 9.681593 8.460905 10.902281 7.814712 11.54847
Problem 5
Fit a seasonal model to the Chinook data up to 1999 using
auto.arima().load("chinook.RData") chinookts <- ts(chinook$log.metric.tons, start=c(1990,1), frequency=12) traindat <- window(chinookts, c(1990,1), c(1999,12)) fit <- forecast::auto.arima(traindat) fit## Series: traindat ## ARIMA(1,0,0)(0,1,2)[12] with drift ## ## Coefficients: ## ar1 sma1 sma2 drift ## 0.3222 -0.2499 -0.2151 -0.0246 ## s.e. 0.1077 0.1253 0.1348 0.0064 ## ## sigma^2 estimated as 0.6872: log likelihood=-126.74 ## AIC=263.49 AICc=264.08 BIC=276.9Then create a forecast through 2015.
We need to forecast 16 years and 12 seasons into the future so
h=16*12. The output is not shown since it is verbose.fr <- forecast::forecast(fit, h=16*12)Plot the forecast with the 2014 and 2015 actual landings added as data points.
plot(fr) testdat <- window(chinookts, c(2014,1), c(2015,12)) points(testdat)