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

Load the data

data(greeklandings, package = "atsalibrary")
landings <- greeklandings
data(chinook, package = "atsalibrary")
chinook <- chinook.month

Problem 1

  1. What is the null hypothesis for the Dickey-Fuller and Augmented Dickey-Fuller tests?

    The null hypothesis is that the time-series has a random walk (unit root) in it and is thus non-stationary.

  2. How do the Dickey-Fuller and Augmented Dickey-Fuller tests differ?

    They differ in terms of the number of lags considered in the stationarity tests, so the length (in lags) of the autocorrelation. The Dickey-Fuller just has \(x_{t-1}\) in the test (specified as lag=0 or k=0). The augmented Dickey-Fuller allows for longer lagged correlations.

  3. For adf.test(), does the test allow the data to have a non-zero level? Does the test allow the data to be stationarity around a trend (a linear slope)?

    Yes it includes both. If you type ?adf.test you will see The general regression equation which incorporates a constant and a linear trend is used. The constant is the non-zero level. So this test uses the most general (but also less powerful) test. If you have removed the level or trend already (via differencing say), then using the constant + trend test would be overly strict.

  4. For ur.df(), what does type = “none”, “drift”, and “trend” mean? Which one gives you the same result as adf.test()? What do you have to set the lags equal to get the default lags in adf.test()?

    “none” means mean zero data with no trend. “drift” means include a non-zero mean (or level). “trend” means include both mean level and a linear trend (trend against time). If you look at ?adf.test, you’ll see that the default lags are trunc((length(x)-1)^(1/3)) so you would need to set lags in ur.df() to that and set type="trend". Alternatively set k in adf.test() and lags in ur.df() to the same values (like to 0 to do the Dickey-Fuller test).

  5. For ur.df(), how do you determine if the null hypothesis (that the data are non-stationary) is rejected?

    There are 2 ways. Say test <- ur.df(x), then summary(test) has the test statistics and critical values at the bottom and you are looking for tau (which will be called tau, tau2 or tau3 depending on what you set type to). However summary() will not tell you which test statistic corresponds to tau. So you can also use attributes(test) and look for teststats and cval which will have tau labeled.

    The null hypothesis of non-stationarity (unit root) is rejected if the tau test statistic value is lower than the tau critical value at the alpha level that you choose. The alpha level is called 1pct, 5pct and 10pct in the output.

  6. For ur.df(), how do you determine if there is a significant trend in the data? How do you determine if the intercept is different than zero?

    You can get an idea of this by looking at the linear regression in the top of the summary(test) output. The test is based on \(y_t = \alpha + \beta t + \phi y_{t-1} + e_t\) but tested as \((y_t - y_{t-1}) = \alpha + \beta t + (\phi-1)y_{t-1} + e_t\) (shown at the top of the summary output). So if the intercept in the output is non-zero and/or the trend is non-zero, that helps you determine if the mean or trend is non-zero. Note, the p-values are not right since that is for non-temporally correlated data but it will give you an idea.

    The tests are in the phi shown at the bottom of the summary output. If type="trend", 3 tests are done and the output is

    • tau3 test: \(\gamma = 0\) and \(\alpha\) and \(\beta\) may or may not be zero.
    • phi2 test: \(\gamma = \beta = 0\) and \(\alpha\) may or may not be.
    • phi3 test: \(\gamma = \beta = \alpha =0\)

Problem 2 KPSS tests in R.

  1. What is the null hypothesis for the KPSS test

    The null hypothesis is stationarity (around a mean or linear trend), so opposite of Dickey-Fuller tests. This test allows more general forms of non-stationarity.

  2. For kpss.test(), what does setting null equal to “Level” versus “Trend” change?

    Setting null="Trend" allows for stationarity around a linear trend.

Problem 3

Repeat the stationarity tests for sardine 1964-1987 in the landings data set. Here is how to set up the data.

```r
datdf <- subset(landings, Species=="Sardine")
dat <- ts(datdf$log.metric.tons, start=1964)
dat <- window(dat, start=1964, end=1987)
```

  1. Do a Dickey-Fuller (DF) test using ur.df() and adf.test(). You will have to set the lags. What does the result tell you? Note for ur.df() use summary(ur.df(...)) and look at the bottom of the summary information for the test statistics and critical values. The first test statistic is the one you want, labeled tau (or tau3).

    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)
    ## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
    ## 
##  Augmented Dickey-Fuller Test
## 
## data:  dat
## Dickey-Fuller = -2.854, Lag order = 0, p-value = 0.247
## alternative hypothesis: stationary

    The 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.19056 -0.07610  0.03307  0.08018  0.16923 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  5.233710   1.839900   2.845  0.01002 * 
## z.lag.1     -0.566363   0.198445  -2.854  0.00981 **
## tt           0.002727   0.003686   0.740  0.46805   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1139 on 20 degrees of freedom
## Multiple R-squared:  0.2895, Adjusted R-squared:  0.2185 
## F-statistic: 4.075 on 2 and 20 DF,  p-value: 0.03278
## 
## 
## Value of test-statistic is: -2.854 2.8125 4.0747 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -4.38 -3.60 -3.24
## phi2  8.21  5.68  4.67
## phi3 10.61  7.24  5.91

    Ignore the part in the summary with the p-value. You need to look at the test-statistic at the bottom of the summary output and see if the tau statistic is less than the critical value at 5pct; if greater it means you do not reject the null of non-stationarity. The first test statistic is tau3, the second phi2 and the third is phi3. The test statistic for tau3 is is greater than the critical value for tau3. Thus the null of \(\gamma=0\) (random walk or unit root) is not rejected at \(\alpha=0.05\).

    For our purposes, we just need to focus on the tau. The phi are for other null hypotheses involving the presence or absence of the drift and trend terms.

  2. Do an Augmented Dickey-Fuller (ADF) test using ur.df(). How did you choose to set the lags? How is the ADF test different than the DF test?

    The ADF test allows longer lags (so longer temporal correlation) in the data compared to the Dickey-Fuller which only includes lag-1 correlation.

    There are a few different approaches to use to set the lags. You could use the rule of thumb that adf.test() uses which is trunc((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.3781 0.7047 0.9769

    You could also let ur.df() select the best lag based on AIC or BIC. It chooses 2 also. Look for the number of z.diff.lags in the summary table.

  3. Do a KPSS test using kpss.test(). What does the result tell you?

    tseries::kpss.test(dat, null="Trend")
    ## 
##  KPSS Test for Trend Stationarity
## 
## data:  dat
## KPSS Trend = 0.13266, Truncation lag parameter = 2, p-value = 0.07471

    The null of trend stationary is not rejected at \(\alpha=0.05\). The null hypothesis for KPSS is stationarity so the data do pass this stationarity test. It is important to set null="Trend" to allow for stationarity around a linear trend.

    So in this case, the KPSS and ADF tests give somewhat different results. This suggests that that non-stationarity cannot be ruled out but the data are not obviously non-stationary.

Problem 4

Use the anchovy 1964-2007 data [Corrected 1/20. If you did the HW with 1964-1987, that’s fine but part b won’t have any models within 2 of the best for the shorter series.]. Fit this time series using auto.arima() with trace=TRUE.

    datdf <- subset(landings, Species=="Anchovy")
    dat <- ts(datdf$log.metric.tons, start=1964)
    anchovy <- window(dat, start=1964, end=2007)
    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.67

  1. Fit each of the models listed using Arima() and show that you can produce the same AICc value that is shown in the trace table.

    You need to pass in order and include.drift to Arima(). AICc is in $aicc in the output.

    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

  2. What models are within \(\Delta\)AICc of 2 of the best model (model with lowest AICc)? What is different about these models? The AICc is output in fit$aicc.

    The best model is ARIMA(0,1,1) no drift and ARIMA(0,1,0) no drift and ARIMA(0,1,1) with drift are within 2. These models are all random walks (no AR, one difference) but different in whether they have autocorrelated error (MA component) and drift (trend) or not.

Problem 5

Repeat the stationarity tests and differencing tests for anchovy using the following two time ranges: 1964-1987 and 1988-2007.

    datdf <- subset(landings, Species=="Anchovy")
    dat <- ts(datdf$log.metric.tons, start=1964)
    dat64.87 <- window(dat, start=1964, end=1987)
    dat88.07 <- window(dat, start=1988, end=2007)

  1. Plot the time series for the two time periods. For the kpss.test(), which null is appropriate, “Level” or “Trend”?

    par(mfrow = c(1,2))
plot(dat64.87, ylab = "Log Metric Tons of Sardines", xlab = "Year")
plot(dat88.07, ylab = "Log Metric Tons of Sardines", xlab = "Year")

    We will definitely need “Trend” since there is a strong upward trend in the early data.

  2. Do the conclusions regarding stationarity and the amount of differencing needed change depending on which time period you analyze? For both time periods, use adf.test() with default values and kpss.test() with null=“Trend”.

    tseries::adf.test(dat88.07)
    ## 
##  Augmented Dickey-Fuller Test
## 
## data:  dat88.07
## Dickey-Fuller = -2.1734, Lag order = 2, p-value = 0.5063
## alternative hypothesis: stationary
    tseries::kpss.test(dat88.07, null="Trend")
    ## Warning in tseries::kpss.test(dat88.07, null = "Trend"): p-value greater than
## printed p-value
    ## 
##  KPSS Test for Trend Stationarity
## 
## data:  dat88.07
## KPSS Trend = 0.071148, Truncation lag parameter = 2, p-value = 0.1
    forecast::ndiffs(dat88.07, test="kpss")
    ## [1] 0
    forecast::ndiffs(dat88.07, test="adf")
    ## [1] 1
    tseries::adf.test(dat64.87)
    ## 
##  Augmented Dickey-Fuller Test
## 
## data:  dat64.87
## Dickey-Fuller = -0.57814, Lag order = 2, p-value = 0.9685
## alternative hypothesis: stationary
    tseries::kpss.test(dat64.87, null="Trend")
    ## 
##  KPSS Test for Trend Stationarity
## 
## data:  dat64.87
## KPSS Trend = 0.19182, Truncation lag parameter = 2, p-value = 0.01907
    forecast::ndiffs(dat64.87, test="kpss")
    ## [1] 1
    forecast::ndiffs(dat64.87, test="adf")
    ## [1] 1

    With 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 not reject that the data are non-stationary. Correspondingly, no differencing is required to pass the KPSS test 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.

    With the early data, both tests indicate non-stationarity and ndiffs() reports needing one difference.

  3. Fit each time period using auto.arima(). Do the selected models change? What do the coefficients mean? Coefficients means the mean and drifts terms and the AR and MA terms.

    forecast::auto.arima(dat64.87)
    ## Series: dat64.87 
## ARIMA(0,1,1) with drift 
## 
## Coefficients:
##           ma1   drift
##       -0.5731  0.0641
## s.e.   0.1610  0.0173
## 
## sigma^2 estimated as 0.03583:  log likelihood=6.5
## AIC=-6.99   AICc=-5.73   BIC=-3.58
    forecast::auto.arima(dat88.07)
    ## Series: dat88.07 
## ARIMA(0,0,1) with non-zero mean 
## 
## Coefficients:
##          ma1    mean
##       0.4776  9.5492
## s.e.  0.2587  0.0538
## 
## sigma^2 estimated as 0.03022:  log likelihood=7.54
## AIC=-9.07   AICc=-7.57   BIC=-6.09

    With the early data ARIMA(0,1,1) with drift is chosen while with the later data ARIMA(0,0,1) with non-zero mean is chosen.

    The first one is \(x_t = x_{t-1} + u + e_t\) where \(e_t\) is lag-1 correlated (MA lag 1). The second one is \(x_t = u + e_t\) where again \(e_t\) is MA-1, so lag 1 correlated.

  4. Discuss the best models for each time period. How are they different?

    This question and the next confused most of you. When thinking about ARIMA() models, I find it helpful to write out the equation, especially if I am comparing models with and without \(d\) (differencing). The equations are shown above. The first one is a random walk with drift (trend) and autocorrelated error (MA-1). A random walk DOES have an AR component; the coefficient on that is 1 but by differencing we get rid of that and make our data stationary. The second ARIMA model is stationary and is just MA-1 (autocorrelated) error around a mean. Given that the latter data is borderline non-stationary, I would force d=1 when comparing these two time periods.

  5. You cannot compare the AIC values for an Arima(0,1,0) and Arima(0,0,1). Why do you think that is? Hint when comparing AICs, the data being fit must be the same for each model.

    The data are different. In an ARIMA(0,1,0) the data are \(x_t-x_{t-1}\) and in ARIMA(0,0,1), the data are \(x_t\). Note, you can compare any ARIMA(p,d,q) models with AIC as long as the \(d\) part is identical between the ARIMA models.

Problem 6

For the anchovy 1964-2007 data, use auto.arima() with stepwise=FALSE to fit models.

    datdf <- subset(landings, Species=="Anchovy")
    anchovy <- ts(datdf$log.metric.tons, start=1964)

  1. find the set of models within \(\Delta AICc=2\) of the top model.

    forecast::auto.arima(anchovy, stepwise = FALSE, trace = TRUE)
    ## 
##  ARIMA(0,1,0)                    : -8.876343
##  ARIMA(0,1,0) with drift         : -7.274254
##  ARIMA(0,1,1)                    : -9.892016
##  ARIMA(0,1,1) with drift         : -9.062605
##  ARIMA(0,1,2)                    : -8.0402
##  ARIMA(0,1,2) with drift         : -7.292236
##  ARIMA(0,1,3)                    : -6.228872
##  ARIMA(0,1,3) with drift         : -5.08747
##  ARIMA(0,1,4)                    : -3.663422
##  ARIMA(0,1,4) with drift         : -2.402711
##  ARIMA(0,1,5)                    : -1.155922
##  ARIMA(0,1,5) with drift         : 0.1419822
##  ARIMA(1,1,0)                    : -8.861307
##  ARIMA(1,1,0) with drift         : -7.565745
##  ARIMA(1,1,1)                    : -7.81313
##  ARIMA(1,1,1) with drift         : -7.033147
##  ARIMA(1,1,2)                    : -5.943474
##  ARIMA(1,1,2) with drift         : -5.016628
##  ARIMA(1,1,3)                    : -3.661973
##  ARIMA(1,1,3) with drift         : -2.389519
##  ARIMA(1,1,4)                    : -0.9700597
##  ARIMA(1,1,4) with drift         : Inf
##  ARIMA(2,1,0)                    : -8.571668
##  ARIMA(2,1,0) with drift         : -7.660279
##  ARIMA(2,1,1)                    : -6.140817
##  ARIMA(2,1,1) with drift         : -5.101195
##  ARIMA(2,1,2)                    : Inf
##  ARIMA(2,1,2) with drift         : -3.451482
##  ARIMA(2,1,3)                    : Inf
##  ARIMA(2,1,3) with drift         : Inf
##  ARIMA(3,1,0)                    : -6.141033
##  ARIMA(3,1,0) with drift         : -5.099723
##  ARIMA(3,1,1)                    : Inf
##  ARIMA(3,1,1) with drift         : Inf
##  ARIMA(3,1,2)                    : Inf
##  ARIMA(3,1,2) with drift         : 0.4594539
##  ARIMA(4,1,0)                    : -3.573995
##  ARIMA(4,1,0) with drift         : -2.407022
##  ARIMA(4,1,1)                    : -0.8887037
##  ARIMA(4,1,1) with drift         : 0.4437702
##  ARIMA(5,1,0)                    : -1.055414
##  ARIMA(5,1,0) with drift         : 0.3626202
## 
## 
## 
##  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.67

    The chosen model is an ARIMA(0,1,1) model with AICc = -9.89. The models within 2 of this are:

    • ARIMA(0,1,0), AICc = -8.88

    • ARIMA(0,1,1) with drift, AICc = -9.06

    • ARIMA(0,1,2), AICc = -8.04

    • ARIMA(1,1,0), AICc = -8.86

  2. Use Arima() to fit the models with Inf or -Inf in the list. Does the set of models within \(\Delta AICc=2\) change?

    There are a number of models with an AICc value of inf. For example

    • ARIMA(1,1,4) with drift

    If we fit these with Arima(), we can get the AICc and you’ll see that they are not close to the best model.

    fit <- forecast::Arima(anchovy, order=c(1,1,4), include.drift=TRUE)

    If you plot this, you’ll see that it has a root close to the unit circle.

    forecast::autoplot(fit)

  3. Create a 5-year forecast for each of the top 3 models according to AICc.

    Here are the top 3 models. They are all random walks but vary in terms of having a trend or not or having autocorrelated errors or not.

    require(forecast)
    ## Loading required package: forecast
    ## Warning: package 'forecast' was built under R version 4.0.3
    mod1 <- forecast::Arima(anchovy, order=c(0,1,1), include.drift=FALSE)
mod2 <- forecast::Arima(anchovy, order=c(0,1,1), include.drift=TRUE)
mod3 <- forecast::Arima(anchovy, order=c(0,1,0), include.drift=FALSE)

    We can forecast with forecast() and plot. I put on the same y scale so you can compare the prediction intervals easier.

    par(mfrow=c(1,3))
plot(forecast(mod1, h=5), main="mod1", ylim=c(8,11))
plot(forecast(mod2, h=5), main="mod2", ylim=c(8,11))
plot(forecast(mod3, h=5), main="mod3", ylim=c(8,11))

  4. How do the forecasts differ in trend and size of prediction intervals?

    mod2 is the only one with a trend so that’s pretty different. mod3 is a pure random walk and has the largest prediction intervals. The models with MA errors have slightly smaller intervals.

Problem 7

Using the chinook data set,

  1. Set up a monthly time series object for the Chinook log metric tons catch for Jan 1990 to Dec 2015.

    load("chinook.RData")
chinookts <- ts(chinook$log.metric.tons, start=c(1990,1), 
            frequency=12)
chinookts <- window(chinookts, c(1990,1), c(2015,12))

  2. Fit a seasonal model to the Chinook Jan 1990 to Dec 1999 data using auto.arima().

    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.7523:  log likelihood=-126.74
## AIC=263.49   AICc=264.08   BIC=276.9

    The fit is ARIMA(1,0,0)(0,1,2)[12] with drift.

  3. Create a forecast through 2015 using the model in part b.

    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)

  4. 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)

  5. The model from part b has drift. Fit this model using Arima() without drift and compare the 2015 forecast with this model.

    Look at ?Arima to see how to specify a seasonal model. We want to fit ARIMA(1,0,0)(0,1,2)[12] without drift.

    fit2 <- forecast::Arima(traindat, order=c(1,0,0), seasonal=c(0,1,2), include.drift=FALSE)
fr2 <- forecast::forecast(fit2, h=16*12)

    Now plot:

    plot(fr2)
points(testdat)