We will use the dataset elecequip from the package fpp for this Lab. These data show the number of new orders for electrical equipment (computer, electronic and optical products) in the Eurozone area (16 countries). The data have already been adjusted by the number of working days and normalized so a value of 100 corresponds to 2005.

require(fpp)

Let us see if there is a trend in the dataset elecequip over the years:

data(elecequip)
head(elecequip) # see what is in the dataset
##        Jan   Feb   Mar   Apr   May   Jun
## 1996 79.43 75.86 86.40 72.67 74.93 83.88
plot(elecequip,main="Time series plot for Electrical Equipment Manufacturing", xlab="", ylab="New orders index")

Note that elecequip is a ts object.

class(elecequip)
## [1] "ts"

Let us decompose the variable in elecequip into the time series components: Cycle, Trend, Seasonal and Irregularity. The additive model is:

\[ Y_t = S_t + T_t + \epsilon_t \]

fit <- decompose(elecequip, type="additive") #additive decomposition
plot(fit)

fitm <- decompose(elecequip, type="multiplicative") #multiplicative
plot(fitm)

Let us now see how the seasonal data looks for each month in the series:

fits <- stl(elecequip, s.window=7)
plot(elecequip, col="gray",
 main="Electrical equipment manufacturing",
 ylab="New orders index", xlab="")
lines(fits$time.series[, "trend"], col="red", ylab="Trend")

monthplot(fits, choice = "seasonal")

monthplot(fits, choice = "trend")

monthplot(fits, choice = "remainder")

Let us see what the seasonally adjusted plot would look like:

plot(elecequip, col="grey",
 main="Electrical equipment manufacturing",
  xlab="", ylab="New orders index")
lines(seasadj(fit),col="red",ylab="Seasonally adjusted")

Seasonal Trend Lowess approach

How about if we detrend the series using STL (seasonal trend using lowess)? Note the trend window is set to 15 months and seasonal is set to periodic.

eeadj <- seasadj(stl(elecequip, s.window="periodic"))
plot(eeadj) # plot with only seasonal adjustment

fitloess <- stl(elecequip, t.window=15, s.window="periodic", robust=TRUE)
plot(fitloess)#plot with both trend and seasonal adjusment

Autoregressive models

An autoregressive model uses previous instances of each measurement as a predictor in a linear model for future instances.

\[ Y_t = \beta_0 + \sum_{i=1}^p \beta_i Y_{t-i} + \epsilon_t \,, \] where the parameters to be fit are the \(\beta_i\)’s. As usual, the error terms are assumed to follow a normal distribution with mean 0 (e.g. \(\epsilon \sim N(0, \sigma)\)). We can fit the model using OLS regression, but here we will use a special function in R.

The autoregressive model \(AR(p)\) uses \(p\) lagged terms.

ar.mod <- ar(eeadj, order.max = 5)
ar.mod
## 
## Call:
## ar(x = eeadj, order.max = 5)
## 
## Coefficients:
##       1        2        3        4  
##  0.6705   0.2928   0.2867  -0.2999  
## 
## Order selected 4  sigma^2 estimated as  11.04
ar.hats <- predict(ar.mod, n.ahead=24)
plot(eeadj, xlim=c(1996, 2014))
lines(ar.hats$pred, col="blue")
lines(ar.hats$pred + ar.hats$se, col="red")
lines(ar.hats$pred - ar.hats$se, col="red")

Moving Average

Moving average models are similar, but use the error terms rather than the previous values. That is, here the measurements are based on the lagged errors, rather than the lagged observed values.

The moving average model of order \(q\) is denoted \(MA(q)\):

\[ Y_t = \mu + \sum_{i=1}^q \theta_i \epsilon_{t-i} + \epsilon_t \,, \] where the parameters to be fit are the \(\theta_i\)’s.

Note that the above plot is not adjusted for trend or cycles so is not very useful. Let us look at a plot after computing a 2x12 moving average which is suitable for monthly data.

ma.mod <- ma(eeadj, order = 12)
plot(eeadj, ylab="New orders index", col="gray",
 main="Electrical equipment manufacturing (Euro area)")
lines(ma.mod, col="red") #moving average for m=12

ARIMA model

An ARMA model combines the ideas of autoregressive model with moving average models.

\[ Y_t = \mu + \sum_{i=1}^p \beta_i Y_{t-i} + \sum_{i=1}^q \theta_i \epsilon_{t-i} + \epsilon_t \,. \]

Let us detrend using the difference method

tsdisplay(diff(eeadj), main="")

The ARIMA model is a generalization of the ARMA model that includes a third parameter for taking the differences of consecutive terms. This the parameter \(d\) in the \(ARIMA(p,d,q)\) model.

Now let us fit the ARIMA model to find the lowest AIC. Choices include \(ARIMA (2, 1, 0)\), \(ARIMA(3, 1, 0)\), \(ARIMA(3, 1, 1)\). Try a few choices of your own and see if you can get a lower AIC.

fitarima <- arima(eeadj, order=c(3,1,1)) # gave the lowest AIC from above choices
summary(fitarima)
## 
## Call:
## arima(x = eeadj, order = c(3, 1, 1))
## 
## Coefficients:
##          ar1     ar2     ar3      ma1
##       0.0519  0.1191  0.3730  -0.4542
## s.e.  0.1840  0.0888  0.0679   0.1993
## 
## sigma^2 estimated as 9.532:  log likelihood = -484.08,  aic = 978.17
## 
## Training set error measures:
##                        ME     RMSE      MAE         MPE     MAPE      MASE
## Training set -0.001227744 3.079373 2.389267 -0.04290849 2.517748 0.8676442
##                     ACF1
## Training set 0.008928479

After you decide on a model you need to check the residual plots for autocorrelation.

Acf(residuals(fitarima))

Box.test(residuals(fitarima), lag=24, fitdf=4, type="Ljung")
## 
##  Box-Ljung test
## 
## data:  residuals(fitarima)
## X-squared = 20.496, df = 20, p-value = 0.4273

Much of this material was taken from Forecasting: Prinicples and Practice by Hyndman and Athanaosopoulos.