Constructing ARMA Models

Sección 6.4–6.7 material impreso

This example uses the monthly UK house price series which was already previously incorporated
(‘ukhp.RData’). So first we re-load the workfile. There is a total of 326 monthly observations running
from February 1991 (recall that the January observation is ‘NA’ because we are looking at returns) to
March 2018 for the percentage change in house price series. Although in general it is not recommended to
drop observations unless necessary, it will be more convenient to drop the first row of UKHP as the ‘NA’
value can cause trouble in several functions. To do so, execute the following line.
UKHP = UKHP [ - 1 ,]

The objective of this exercise is to build an ARMA model for the house price changes. Recall that there
are three stages involved: identification, estimation and diagnostic checking. The first stage is carried
out by looking at the autocorrelation and partial autocorrelation coefficients to identify any structure
in the data.

Estimating Autocorrelation Coefficients

To generate autocorrelations and partial correlations, the stats package provides the functions acf and
pacf. Setting the argument lag.max to twelve, you obtain a plot of the first 12 autocorrelations and
partial autocorrelations as shown in Figure 1.

Figure 1: Autocorrelation and Partial Autocorrelation Functions

It is clearly evident that the series is quite persistent given that it is already in percentage change
form: the autocorrelation function dies away rather slowly. Only the first two partial autocorrelation
coefficients appear strongly significant.
Remember that as a rule of thumb, a given autocorrelation coefficient is classed as significant if it is
outside a ±1.96 × 1/(T )1/2 band, where T is the number of observations. In this case, it would imply
that a correlation coefficient is classed as significant if it is bigger than approximately 0.11 or smaller than
–0.11. This band is also drawn in Figure 1 with blue colour. The band is of course wider when the sampling
frequency is monthly, as it is here, rather than daily, where there would be more observations. It can be
deduced that the first six autocorrelation coefficients (then nine through twelve) and the first

43
two partial autocorrelation coefficients (then nine, eleven and twelve) are significant under this rule.
Since the first acf coefficient is highly significant, the joint test statistic presented in column 4 rejects
the null hypothesis of no autocorrelation at the 1% level for all numbers of lags considered. It could be
concluded that a mixed ARMA process might be appropriate, although it is hard to precisely determine
the appropriate order given these results. In order to investigate this issue further, information criteria
are now employed.

Using Information Criteria to Decide on Model Orders

An important point to note is that books and statistical packages often differ in their construction of the
test statistic. For example, the formulae given in Brooks (2019) for Akaike’s and Schwarz’s Information
Criteria are
2k
AIC = ln(σ̂ 2 ) + (2)
T
k
SBIC = ln(σ̂ 2 ) + ln(T ) (3)
T
where σ̂ 2 is the estimator of the variance of regressions disturbances ut , k is the number of parameters
and T is the sample size. When using the criterion based on the estimated standard errors, the model
with the lowest value of AIC and SBIC should be chosen. The AIC() function of the stats package in
R uses a formulation of the test statistic based on maximum likelihood estimation. The corresponding
formulae are

AICL = −2 · ln(L) + 2 · n (4)

SBICL = −2 · ln(L) + ln(T ) · n (5)

where ln(L) is the maximised log-likelihood of the model. Unfortunately, this modification is not benign,
since it affects the relative strength of the penalty term compared with the error variance, sometimes
leading different packages to select different model orders for the same data and criterion!
Suppose that it is thought that ARMA models from order (0,0) to (5,5) are plausible for the house
price changes. This would entail considering 36 models (ARMA(0,0), ARMA(1,0), ARMA(2,0), . . . ,
ARMA(5,5)), i.e., up to 5 lags in both the autoregressive and moving average terms.
This can be done by separately estimating each of the models and noting down the value of the
information criteria in each case. We can do this in the following way. Using the function arima() for
an ARMA(1,1) model, we specify the argument order as c(1,0,1). An equivalent formulation would be

∆HPt = µ + ϕ∆HPt−1 + θt−1 + t , (6)

with mean µ, AR coefficient ϕ and MA coefficient θ. We type

arima ( UKHP$dhp , order = c ( 1 ,0 ,1 ) )

into the Console and the following estimation output appears.

44
 In theory, the output would be discussed in a similar way to the simple linear regression model previously
discussed. However, in reality it is very difficult to interpret the parameter estimates in the sense of, for
example, saying ‘a 1 unit increase in x leads to a β unit increase in y’. In part because the construction of
ARMA models is not based on any economic or financial theory, it is often best not to even try to interpret
the individual parameter estimates, but rather to examine the plausibility of the model as a whole, and to
determine whether it describes the data well and produces accurate forecasts (if this is the objective of the
exercise, which it often is).
In order to generate the information criteria corresponding to the ARMA(1,1) model, we use the
function AIC(). Although the AIC criterion is already given in the output above, it can also be
computed using this function. For the SBIC, the only change necessary is to set k to ln(T). Having
saved the model from before into a new variable ar11, we apply the functions and obtain the output
below.

We see that the AIC has a value of 928.30 and the BIC a value of 943.45. However, by themselves these
two statistics are relatively meaningless for our decision as to which ARMA model to choose. Instead,
we need to generate these statistics for the competing ARMA models and then select the model with
the lowest information criterion. Repeating these steps for the other ARMA models would give all of
the required values for the information criteria. The quickest way to this is within two loops running
over the two variable the order of the AR part and the order of the MA part. A possible set of code
could look like the one below.
1 aic_table = array ( NA , c ( 6 ,6 ,2 ) )
2 for ( ar in 0 : 5 ) {
3 for ( ma in 0 : 5 ) {
4 arma = arima ( UKHP$dhp , order = c ( ar , 0 , ma ) )
5 aic_table [ ar + 1 , ma + 1 ,1 ] = AIC ( arma )
6 aic_table [ ar + 1 , ma + 1 ,2 ] = AIC ( arma , k = log ( nrow ( UKHP ) ) )
7 }
8 }

45
As this section is dealing with ARMA models and the way to choose the right order, the above code will
not be discussed in detail. However, note that it makes use of the possibility to define three dimensional
arrays in line 1. The table created by the code above is presented below with the AIC values in the
first table and the SBIC values in the second table. Note that due to the fact that row and column
indices start with 1, the respective AR and MA order represented is to be reduced by 1, i.e., entry (5,3)
represents an ARMA(4,2) model.

So which model actually minimises the two information criteria? A quick way of finding the index of
the minimal entry in a matrix is to use the function which.min(). Note that the index is not given in
two dimensions, but as if the matrix would be a list of concatenated columns. In this case, the criteria
choose different models: AIC selects an ARMA(4,2), while SBIC selects the smaller ARMA(2,0) model.
It will always be the case that SBIC selects a model that is at least as small (i.e., with fewer or the same
number of parameters) as AIC, because the former criterion has a stricter penalty term. This means
that SBIC penalises the incorporation of additional terms more heavily. Many different models provide
almost identical values of the information criteria, suggesting that the chosen models do not provide
particularly sharp characterisations of the data and that a number of other specifications would fit the
data almost as well.

46
Forecasting Using ARMA Models
Sección 6.8 material impreso

Suppose that an AR(2) model selected for the house price percentage changes series was estimated using
observations February 1991–December 2015, leaving 27 remaining observations to construct forecasts for
and to test forecast accuracy (for the period January 2016–March 2018). Let us first estimate the
ARMA(2,0) model for the time period from February 1991–December 2015. To account for this reduced
dataset, simply put the restrictions into the square brackets for the dataset. The regression output is
presented below.

Now that we have fitted the model, we can produce the forecasts for the period January 2016 to March
2018. There are two methods for constructing forecasts: dynamic and static. Dynamic forecasts are
multi-step forecasts starting from the first period in the forecast sample. Static forecasts imply a
sequence of one-step-ahead forecasts, rolling the sample forwards one observation after each forecast. In
R, the simplest way to forecast is using the function predict which is part of the system library stats.1
Although a static forecast only uses one-step ahead predictions, it is not as easy to implement, since the
data has to be updated and predict does not allow for this. A dynamic forecast instead does not need
any further input as it builds on the forecasted values to compute the next forecasts.
Thus, we can compute the static forecast directly after estimating the AR(2) model by typing
ar 2 = arima ( UKHP$dhp [ UKHP$Month <= " 2 0 1 5 -1 2 -0 1 " ] , order = c ( 2 ,0 ,0 ) )
dynamic_fc = predict ( ar 2 ,n . ahead = 2 7 )

Since there are 27 out-of-sample observations, set the argument n.ahead to 27 to produce a 27-steps
ahead forecast. For a one-step ahead forecast, n.ahead would only need to be altered to 1. But as the
static forecast uses the out-of-sample observations, it is less straight-forward to compute. However, we
can just use the formula for the one-step ahead forecast of AR models and update the data manually.
Remember that for an AR(2) model like

∆HPt = µ + ϕ1 ∆HPt−1 + ϕ2 ∆HPt−2 + t , (7)

the one-step ahead forecast is given by

∆HP
\ t+1 = µ + ϕ1 ∆HPt + ϕ2 ∆HPt−1 + t , (8)
1
Note that there is also the package forecast which extends some forecasting techniques that are not necessary for
this section.

47
and for any step h the static forecast is given as
∆HP
\ t+h = µ + ϕ1 ∆HPt+h−1 + ϕ2 ∆HPt+h−2 + t , (9)
Putting this into code is very efficiently done by typing
static_fc = ar 2 $coef [ 3 ]+ ar 2 $coef [ 1 ]* UKHP$dhp [ 2 9 9 : 3 2 5 ]+ ar 2 $coef [ 2 ]* UKHP$dhp [ 2
98:324]

which again makes use of the way R interprets combinations of vectors and scalars. The variable coef
contains the estimated coefficients ϕ1 , ϕ2 and µ. Multiplied with the vectors of the past observations of
dhp, the code resembles Equation (9) for h=1 to 27.
To spot differences between the two forecasts and to compare them to the actual values of the changes
in house prices that were realised over this period, it is useful to create a graph of the three series. The
code below is sufficient to do so. We have added some graphical options to show how to make graphs
look more professional, while leaving the editing at a minimum.
par ( lwd = 2 , cex . axis = 2 )
plot ( UKHP$Month [ 3 0 0 : 3 2 6 ] , UKHP$dhp [ 3 0 0 : 3 2 6 ] , type = " l " , xlab = " " , ylab = " " )
lines ( UKHP$Month [ 3 0 0 : 3 2 6 ] , dynamic_fc$mean , col = " blue " )
lines ( UKHP$Month [ 3 0 0 : 3 2 6 ] , static_fc , col = " red " )
legend ( " topright " , legend = c ( " Actual " , " Dynamic " , " Static " ) , col = c ( " black " ,"
blue " ," red " ) , lty = 1 )

The resulting graph should resemble Figure 2.

Figure 2: Graph Comparing Static and Dynamic Forecasts with the Actual Series

Let us have a closer look at the graph. For the dynamic forecasts, it is clearly evident that the forecasts
quickly converge upon the long-term unconditional mean value as the horizon increases. Of course, this
does not occur with the series of 1-step-ahead forecasts which seem to more closely resemble the actual
‘dhp’ series.
A robust forecasting exercise would of course employ a longer out-of-sample period than the two years
or so used here, would perhaps employ several competing models in parallel, and would also compare
the accuracy of the predictions by examining the forecast error measures, such as the square root of the
mean squared error (RMSE), the MAE, the MAPE, and Theil’s U-statistic.

48
Estimating Exponential Smoothing Models
Sección 6.9 material impreso

Exponential smoothing models can be estimated in R using the package smooth. There is a variety of
smoothing methods available, including single and double, or various methods to allow for seasonality
and trends in the data. We will focus on single-exponential smoothing method.
Using the same in-sample data as before and forecasting 27 steps ahead can be done using the
function es() with the following input.

es ( data = UKHP$dhp [ 1 : 2 9 9 ] , h = 2 7 )

Here we use the simpler notation for the observations before January 2016, as we now know these are
the first 299 observations. The argument h, as before, defines the out-of-sample observations. After
pressing Enter, the following output is obtained.

The output includes the value of the estimated smoothing coefficient α (0.235 in this case), together
with the root mean squared error (RMSE) or residual standard deviation for the whole forecast. The
final in-sample smoothed value will be the forecast for those 27 observations (which in this case would
be 0.2489906). You can find these forecasts by saving the forecast model to a variable, e.g., smooth fc
and query the value forecast.

49