The model includes a non-seasonal MA(1) term, a seasonal MA(1) term, no differencing, no AR terms and the seasonal period is S = 12.

The non-seasonal MA(1) polynomial is \(\theta(B) = 1 + \theta_1 B \). 

The seasonal MA(1) polynomial is \(\Theta(B^{12}) = 1 + \Theta_1B^{12}\). 

The model is \((x_t - \mu) = \Theta(B^{12})\theta(B)w_t = (1+\Theta_1B^{12})(1+\theta_1 B)w_t \).  

When we multiply the two polynomials on the right side, we get

\( (x_t - \mu) = (1 + \theta_1 B + \Theta_1 B^{12} + \theta_1 \Theta_1 B^{13})w_t \)

\( = w_t + \theta_1 w_{t-1} + \Theta_1 w_{t-12} + \theta_1 \Theta_1 w_{t-13} \).

Thus the model has MA terms at lags 1, 12, and 13. This leads many to think that the identifying ACF for the model will have non-zero autocorrelations only at lags 1, 12, and 13. There’s a slight surprise here. There will also be a non-zero autocorrelation at lag 11. We supply a proof in the Appendix below.

We simulated n = 1000 values from an ARIMA \(( 0,0,1 ) \times ( 0,0,1 ) _ { 12 }\). The non-seasonal MA (1) coefficient was \(\theta_1=0.7\). The seasonal MA(1) coefficient was \(\Theta_1=0.6 \). The sample ACF for the simulated series was as follows:

The model includes a non-seasonal AR(1) term, a seasonal AR(1) term, no differencing, no MA terms and the seasonal period is S = 12.

The non-seasonal AR(1) polynomial is \(\phi(B) = 1-\phi_1B \). 

The seasonal AR(1) polynomial is \( \Phi(B^{12}) = 1-\Phi_1B^{12} \).

The model is \( (1-\phi_1B)(1-\Phi_1B^{12})(x_t-\mu) = w_t \).

If we let \(z_t = x_t -\mu\) (for simplicity), multiply the two AR components and push all but z_t to the right side we get \(z_{t} = \phi_1 z_{t-1} + \Phi_1 z_{t-12} + \left( - \phi_1\Phi_1 \right) z_{t-13} + w_t\).

This is an AR model with predictors at lags 1, 12, and 13.

R can be used to determine and plot the PACF for this model, with \(\phi_1\)=.6 and \(\Phi_1\)=.5. That PACF (partial autocorrelation function) is:

It’s not quite what you might expect for an AR model, but it almost is. There are distinct spikes at lags 1, 12, and 13 with a bit of action coming before lag 12. Then, it cuts off after lag 13.

R commands were

thepacf=ARMAacf (ar = c(.6,0,0,0,0,0,0,0,0,0,0,.5,-.30),lag.max=30,pacf=T)
plot (thepacf,type="h")

The data series are a monthly series of a measure of the flow of the Colorado River, at a particular site, for n = 600 consecutive months.

Step 1

A time series plot is

With so many data points, it’s difficult to judge whether there is seasonality. If it was your job to work on data like this, you probably would know that river flow is seasonal – perhaps likely to be higher in the late spring and early summer, due to snow runoff.

Without this knowledge, we might determine means by month of the year. Below is a plot of means for the 12 months of the year. It’s clear that there are monthly differences (seasonality).

Looking back at the time series plot, it’s hard to judge whether there’s any long run trend. If there is, it’s slight.

Steps 2 and 3

We might try the idea that there is seasonality, but no trend. To do this, we can create a variable that gives the 12^th differences (seasonal differences), calculated as \(x _ { t } - x _ { t - 12}\). Then, we look at the ACF and the PACF for the 12^th difference series (not the original data). Here they are:

• Non-seasonal behavior: The PACF shows a clear spike at lag 1 and not much else until about lag 11. This is accompanied by a tapering pattern in the early lags of the ACF. A non-seasonal AR(1) may be a useful part of the model.
• Seasonal behavior: We look at what’s going on around lags 12, 24, and so on. In the ACF, there’s a cluster of (negative) spikes around lag 12 and then not much else. The PACF tapers in multiples of S; that is the PACF has significant lags at 12, 24, 36 and so on. This is similar to what we saw for a seasonal MA(1) component in Example 1 of this lesson.

Remembering that we’re looking at 12^th differences, the model we might try for the original series is ARIMA \(( 1,0,0 ) \times ( 0,1,1 ) _ { 12 }\).

Step 4

R results for the ARIMA \(( 1,0,0 ) \times ( 0,1,1 ) _ { 12 }\):

Step 5 (diagnostics)

The normality and Box-Pierce test results are shown in Lesson 4.2.  Besides normality, things look good.  The Box-Pierce statistics are all non-significant and the estimated ARIMA coefficients are statistically significant.

The ACF of the residuals looks good too:

What doesn’t look perfect is a plot of residuals versus fits. There’s non-constant variance.

We’ve got three choices for what to do about the non-constant variance: (1) ignore it, (2) go back to step 1 and try a variance stabilizing transformation like log or square root, or (3) use an ARCH model that includes a component for changing variances. We’ll get to ARCH models later in the course.

Lesson 4.2 for this week will give R guidance and an additional example or two.

R code for the Colorado River Analysis

The data are in the file coloradoflow.dat. We used scripts in the library(Stoffer’s routines as described in Lessons 1 & 3), so the first two lines of code are:

library(astsa)
flow <- ts(scan("coloradoflow.dat"))

The first step in the analysis is to plot the time series. The command is:

plot(flow, type="b")

The plot showed clear monthly effects and no obvious trend, so we examined the ACF and PACF of the 12^th differences (seasonal differencing). Commands are:

diff12 = diff(flow,12)
acf2(diff12, 48)

The acf2 command asks for information about 48 lags. On the basis of the ACF and PACF of the 12^th differences, we identified an ARIMA(1,0,0)×(0,1,1)₁₂ model as a possibility (See Lesson 4.1). The command for fitting this model is

sarima(flow, 1,0,0,0,1,1,12)

The parameters of the command just given are the data series, the non-seasonal specification of AR, differencing, and MA, and then the seasonal specification of seasonal AR, seasonal differencing, seasonal MA, and period or span for the seasonality.

Output from the sarima command is:

The output included these residual plots. The only difficulty we see is in the normal probability plot. The extreme standardized sample residuals (on both ends) are larger than they would be for normally distributed data. This could be related to the non-constant variance noted in Lesson 4.1.

We didn’t generate any forecasts in Lesson 4.1, but if we were to use R to generate forecasts for the next 24 months in R, one command that could be used is

sarima.for(flow, 24, 1,0,0,0,1,1,12)

Partial output for the forecast command follows (We skipped giving the standard errors.)

Note that the lower limits of some prediction intervals are negative - impossible for the flow of a river. In practice, we might truncate these lower limits to 0 when presenting them.

If you were to use R’s native commands to do the fit and forecasts, the commands might be:

themodel = arima(flow, order = c(1,0,0), seasonal = list(order = c(0,1,1), period = 12))
themodel
predict(themodel, n.ahead=24) 

The first command does the arima and stores results in an “object” called “themodel.” The second command, which is simply themodel, lists the results and the final command generates forecasts for 24 times ahead.

In Lesson 4.1, we presented a graph of the monthly means to make the case that the data are seasonal. The R commands are:

flowm = matrix(flow, ncol=12,byrow=TRUE)
col.means=apply(flowm,2,mean)
plot(col.means,type="b", main="Monthly Means Plot for Flow", xlab="Month", ylab="Mean")

In Lesson 1.1, this plot for quarterly beer production in Australia was given.

There is an upward trend and seasonality. The regular pattern of ups and downs is an indicator of seasonality. If you count points, you’ll see that the 4^th quarter of a year is the high point, the 2^nd and 3^rd quarters are low points, and the 1^st quarter value falls in between.

With trend and quarterly seasonality we will likely need both a first and a fourth difference. The plot above was produced in Minitab, but we’ll use R for the rest of this example. The commands for creating the differences and graphing the ACF and PACF of these differences are:

diff4 = diff(beer, 4)
diff1and4 = diff(diff4,1)
acf2(diff1and4,24)

Let's examine the TS plot and ACF/PACF of seasonal differences, diff4, to determine whether or not first differences are also necessary.

The TS plot still shows that the upward trend remains, however the ACF/PACF do not suggest the need to difference further. In this case, we could detrend by subtracting an estimated linear trend from each observation as follows:

dtb=residuals(lm (diff4~time(diff4)))
acf2(dtb)

The ACF and PACF of the detrended seasonally differenced data follow. The interpretation:

• Non-seasonal: Looking at just the first 2 or 3 lags, either a MA(1) or AR(1) might work based on the similar single spike in the ACF and PACF, if at all. Both terms are also possible with an ARMA(1,1), but with both cutting off immediately, this is less likely than a single order model. With S=4, the nonseasonal aspect can sometimes be difficult to interpret in such a narrow window.
• Seasonal: Look at lags that are multiples of 4 (we have quarterly data). Not much is going on there, although there is a (barely) significant spike in the ACF at lag 4 and a somewhat confusing spike at lag 9 (in ACF). Nothing significant is happening at the higher lags. Maybe a seasonal MA(1) or MA(2) might work.

We tried a few models. An initial guess was ARIMA(0,0,1)×(0,1,1)₄, which wasn’t a bad guess. We also tried ARIMA(1,0,0)×(0,1,1)₄, ARIMA(1,0,1)×(0,1,1)₄, and ARIMA(1,0,0)×(0,1,2)_4.

The commands for the 4 models are:

sarima (dtb, 0,0,1,0,0,1,4)
sarima (dtb, 1,0,0,0,0,1,4)
sarima (dtb, 0,0,0,0,0,1,4)
sarima (dtb, 1,0,0,0,0,2,4)

A summary of the results:

All models are quite close, though the second model is the best in terms of AIC and residual autocorrelation and might be worth keeping the marginally significant AR(1) term. The fourth model surprisingly comes through with a slightly lower MSE than the second, though the AIC and simplicity of the second model seems the best choice.