ARIMA modelling in R: minimizing AIC criterion

In earlier posts I talked about the arima() function in R. In particular, I gave some examples showing that it does some funny things when calculating residuals and I also mentioned that three different estimation methods could be specified using the method=() argument. Again, the estimation methods are "CSS-ML", "ML" and "CSS".

One way of specifying an ARIMA model is to choose the specification that minimizes Aikake's information criterion. This process, which I'll refer to as specification search, is data-mining. As far as I am aware, there is an R function called step() that can automatically identify the ARIMA model with the smallest AIC value for a given series of data.

Personally, I am not a fan of this modelling procedure as I have preference for following the Box-Jenkins methodology when modelling univariate time-series. Part of the BJ methodology relies on investigating the estimated autocorrelation function and partial autocorrelation function in order to try to identify possible models worth entertaining as candidate models.

Let me try to justify my preference using a single argument.

One has to remember that there are both producers and consumers of econometric models. The producer builds the model with some purpose, such as for forecasting or to guide policy decisions, and the consumer of the model uses it if he or she believes that the model is valid.

In my opinion, a producer who follows the Box-Jenkins methodology, or some other sound systematic procedure that can be explained to and understood by the consumer, should be in a position with decent chances of convincing the consumer that the final model is valid.

On the other hand, a producer who builds models by simply clicking a button and waiting for the result to pop up, may have trouble convincing consumers that the model is valid. This second kind of producer may find that he or she has very few customers willing to do business (even if the model actually does a good job at representing the data generating mechanism)! I wouldn't argue that an R function like step() is useless or should be avoided, but it probably shouldn't be at the helm during the model building procedure.

Next, I want to discuss a little problem in R that is not related to step().

The arima() function reports an AIC value as part of its output alongside estimated coefficients and standard errors as is standard for any popular statistical package. Now, caution must be taken here. This AIC value is only valid if the estimation method is specified as "ML". So, this AIC value is only useful to the analyst if their model is of the form: arima(x, c(p,d,q), method="ML"). The arima help page tells us the following about the reported AIC value:

aic: the AIC value corresponding to the log-likelihood. Only valid for 'method = "ML"' fits.

There has been some confusion about this and it has been raised in a discussion on the R mailing list, so it's worth pointing it out. I have also come across these University Lecture Notes on Model Selection in R Using AIC which happen to contain errors; the author estimates ARIMA models using the default ("CSS-ML") method and then uses the reported AIC values to select a final model despite the fact that the AIC values are only valid when method="ML"!

There are two key points here. The first is that if a model builder is going to use R for data analysis, he or she should carefully read the documentation. Secondly, model builders should recognize that information about their models, and quite possibly their model building procedure, may have to be communicated to consumers of their models. It's worth noting that consumers - including managers and politicians - may not have a great deal of technical knowledge of model building or even request a technical explanation of the modelling procedure, but they may put more trust in the models if they hear more than "the computer told me... this one".

Bayesian VARs

There are various approaches to modelling economic time-series. One very popular univariate method is known as the Box-Jenkins methodology and it's based on the idea that the analyst should "let the data speak for itself". Moving from the single-equation framework to a vector of equations, we have vector autoregressions and these are associated with Christopher Sims which he introduced in his paper Macroeconomics and Reality. Other econometric methodologies are linked to David Hendry and Edward Leamer; both of whom take different approaches.

Although I've listed various methodologies and said that they differ, they do, however, contain some elements which are the same across the board. A paper by Adrian Pagan (Three Econometric Methodologies: A Critical Appraisal) gives an excellent overview of the Sims, Hendry and Leamer methodologies and also outlines their similarities and differences.

Atheoretic models, such as univariate Box-Jenkins models and the Unrestricted VAR, do not rely on economic theory as the basis for specifying the model. They let the data do the talking. On the other hand, big simultaneous equation models, often built by a team of economists, have their foundations based on economic theory. While, in the Bayesian approach, beliefs (also known as priors) are allowed to enter the model.

Richard Todd discusses combining two of the aforementioned approaches, namely, the Bayesian approach and Sim's VAR framework, in his paper Improving Economic Forecasting with Bayesian Vector Autoregression. Mixing these two methods results in an economic forecasting model known as the Bayesian vector autoregression (BVAR).

We can immediately see the difference between the data-driven Box-Jenkins methodology and BVARs (data + beliefs) when Todd says:

The Bayesian approach to statistics, a general method for combining beliefs with data, suggests an objective procedure for blending beliefs and data in economic forecasting models. This procedure provides a framework that forecasters can use to document and discuss their beliefs, which can help make economic forecasting more of a science and less of an art.

So, how does the BVAR improve the basic Unrestricted VAR? Well, for a first, a vector autoregression with $N$ variables and $p$ lags, involves $N^{2}p$ parameters. In other words, adding extra variables and lags to a VAR consumes many degrees of freedom. This is a well known feature of vector autoregressions. Consider an example given by Clive Granger:

If N=10, p=4 there are 400 parameters to be estimated from 2000 pieces of data, say, if 200 terms are available for each series. It is clear why there is a desire for simpler models.

Clive's comment about the desire for simpler models contains an appeal to the principle of parsimony which says that we should build models that use the smallest number of coefficients needed to explain the available data. In other words, we should not include any unnecessary coefficients in our models.

Structural econometric models use exclusion restrictions to reduce the number of coefficients to be estimated in a VAR model.

So, how does one go about BVAR modelling in practice? Well, according to these slides by Patrick Brandt at the University of Texas, RATS, Eviews, Stata, and Matlab, are not very conducive to modelling BVARs. However, we're in luck because there is a package that has been developed for modelling BVARs in R. The package is called BMR. Strangely enough, this package was developed by one of my old UCD classmates! Better yet, there appears to be a second package MSBVAR that should also do the trick.

Impulse Response Functions

I bought a copy of Ed Leamer's book Macroeconomic Patterns and Stories after hearing him speaking about it on EconTalk. Actually, Professor Leamer is probably my favourite guest on EconTalk, so purchasing his book was a no-brainer for me. If you've heard of him before it's probably due to his famous paper Let's Take the Con Out of Econometrics. He is also well known for advocating "extreme bounds analysis", which, in basic terms, involves finding a possible range of values for an estimate and investigating the fragility of the model by seeing if the value of the estimate is sensitive to small changes in model specification.

The example of EBA given in Let's Take the Con Out of Econometrics is to do with the effect of capital punishment on the murder rate. Using a number of regressions, with different, but somewhat similar specifications, Leamer shows that at one end of the spectrum, the estimate is negative, but at the other end, the estimate is positive!

In his book, Macroeconomic Patterns and Stories, Professor Leamer talks about four key macro variables in detail; namely, Growth, Unemployment, Inflation, and Interest Rates. In this post, I want to focus on these four variables and in particular to the topic of Chapter 6 of the book which is Extrapolative Forecasting.

In Chapter 6, a table is shown and it contains regression results. There are four separate regressions and each of them are autoregressive. Growth is modeled as an AR(2) process, so too is the unemployment rate and 3-month treasuries, while inflation is modeled as an AR(4) process. Admittedly, Professor Leamer says that a more careful data analysis could be done on the data, but the results from these regressions still permit us to get an insight into the dynamic nature of the four variables.

In R, I've estimated how each of the four variables respond to a one unit shock. I've plotted the results in the Figure 1 below. A similar graph appears on page x in Leamer's book and this gives me confidence that the code I wrote in R is correct.

Figure 1

To explain the dynamic behaviour of each of the four variables in response to a shock, Leamer introduces two concepts to the reader: persistence and momentum.

R Time Series Issue: Residuals (Revisited)

In the previous post, I investigated the arima() function and its method of calculating residuals. The issue that I raised was that the first residual associated with arima(x, c(1,0,0)) happens to be non-zero when it should probably either be treated as zero or classified as n.a. (I prefer the latter) when doing time series analysis.

Last time, I specified only two arguments in the arima() function; namely, x and order. For more information type ?arima in R.

In this post, I'd like to point out some more details about calculating residuals in R via arima() and to do so, I'll talk about a third argument of arima() - the method argument - which is used to specify which estimation method is to be used to estimate the parameters of the model. I briefly mentioned last time that maximum-likelihood and non-linear least squares were two possible estimation methods for ARIMA models.

The method argument has three options: "CSS-ML", "ML", and "CSS". The default is set to "CSS-ML", so that was the method used in the previous post. The help page for the arima() function states the following about the method argument:

Fitting method: maximum likelihood or minimize conditional sum-of-squares. The default (unless there are missing values) is to use conditional-sum-of-squares to find starting values, then maximum likelihood. 

I cannot confirm it yet, but I'm almost certain that "CSS-ML" (the default) refers to Marquardt's compromise, "ML" to maximum likelihood, and "CSS" to what I (and the documentation that I've read) call non-linear least squares.

To begin my investigation, I estimated an AR(1) model for the same data in the previous post using, in turn, each of the three estimation methods. That is, I estimated these three models in R:

1. arima(x, c(1,0,0), method = c("CSS-ML")),
2. arima(x, c(1,0,0), method = c("ML")),
3. arima(x, c(1,0,0), method = c("CSS")).

Now, here is the problem...

Upon inspection of the residuals associated with each model, I found that, in two of the three cases ("CSS-ML", "ML") the first residual is non-zero, while in the third case ("CSS"), the first residual is precisely equal to zero.

I've tabulated the residuals for each of the three models in Table 1 below...

If the analyst accepts these first residuals as appropriate then they can fire away and use x$\$$residuals as they please. On the other hand, if the analyst is of the opinion that the first residual should be omitted from their diagnostic checking procedures, they can do this in R by specifying x$residuals[-1]. This command will return all but the first residual as a vector object.

R Time Series Issue: Residuals

Professor David S. Stoffer at the Department of Statistics at the University of Pittsburgh has documented some issues with using R for time-series analysis. In this post, I would like to make an addition to the issues pointed out by Professor Stoffer.

The issue I'd like to bring to the fore is on the subject of residuals and how they are calculated after using the arima() function. Let's go through a numerical example - I suggest grabbing a pencil and paper - with an application in R:

Suppose we have a realization $ z_{t} = (z_{1}, z_{2}, \ldots, z_{t-1}, z_{t})$ where $t = 60$. The data is shown in Table 1 below, although I have omitted a number of data points because we won't be needing them here (don't worry, you can download all of the data, and an R script, at the bottom of the post in order to replicate all of the results presented here). Notice that I have included a row in the table showing that we do not have an observation for $z_{t}$ where $t = 0$. This will be important later on and will, in fact, show us the crux of the problem that I'm raising.

Table 1 \begin{array}{|c|c|} \hline t & z_{t} \\ \hline 0 & n.a. \\ 1 & 4.4 \\ 2 & 5.8 \\ 3 & 6.7 \\ 4 & 7.1 \\ 5 & 5.7 \\ \vdots & \vdots \\ 59 & 13.3 \\ 60 & 6.2 \\ \hline \end{array} Now suppose that given this data, we use the Box-Jenkins methodology and happen to identify an AR(1) process as a tentative model. That is, we investigated the acutocorrelation function and partial autocorrelation for the series and found that the autocorrelation decayed rapidly and there was a single significant spike in the partial autocorrelation function.

In backshift notation an AR(1) can be written as \[(1 - \phi_{1}B) \tilde{z}_{t} = a_{t} \] where $\tilde{z}_{t} = z_{t} - \mu$ is just the original time-series expressed in deviations from the mean, $a_{t}$ are random shocks, B is the backshift operator, and $(1 - \phi_{1}B)$ is an AR operator which can be written in a generalized AR(p) form as: $(1 - \phi_{1}B - \phi_{2}B^{2} - \cdots - \phi_{p}B^{p})$.

Rewriting the model in common algebraic form, because this will help with explanations later on, results in the following series of equations: \begin{aligned} (1 - \phi_{1}B) \tilde{z}_{t} &= a_{t}\\ (1 - \phi_{1}B) (z_{t} - \mu) &= a_{t}\\ z_{t} - \mu - \phi_{1}B z_{t} + \phi_{1}B \mu &= a_{t}\\ z_{t} - \mu - \phi_{1} z_{t-1} + \phi_{1} \mu &= a_{t}\\ z_{t} &= \mu + \phi_{1} z_{t-1} - \phi_{1} \mu + a_{t}\\ z_{t} &= \mu(1 - \phi_{1}) + \phi_{1} z_{t-1} + a_{t}\\ z_{t} &= C + \phi_{1} z_{t-1} + a_{t} \end{aligned} where $\mu(1 - \phi_{1}) = C$ is the constant term. (Tip: Chapter 2 of Hamilton's Time Series Analysis gives an excellent round-up of the algebra of the backshift (or lag) operator. Note, however, that I've used the letter B where Hamilton uses the letter L. Additionally, a concise list of properties of the lag operator can be found in Chapter 8 (p.168-9) of Harrison & Waldron's Mathematics for Economics and Finance.)

At the second stage of the Box-Jenkins methodology - the estimation stage - we fit the model and find estimates for the parameters in our model. In this case, we need to find estimates for $\phi_{1}$ and the constant, $C$. These estimates of the true underlying values are denoted $\hat{\phi}_{1}$ and $\hat{C}$, respectively. ARIMA estimation is usually carried out using nonlinear least squares, however, there is a trade off between using least squares and maximum likelihood which will not discussed further since it is outside the scope of this post.

In R, we can fit an AR(1) model using arima(x, order=c(1,0,0)) and in doing so using the data from Table 1, we get the following results: \[ \hat{\phi}_{1} = 0.6803 \text{ and } \hat{C} = 1.931723.\] Actually, R reports $\hat{\mu}$ as the intercept, so it was necessary to calculate $\hat{C}$ "by hand". That is, the arima() output tells us that the intercept is $6.0423$, but when it refers to the intercept it means $\hat{\mu}$ rather than $\hat{C}$. This issue has already been documented by Professor Stoffer, but to make it clear, R reports to us the following information: \begin{aligned} z_{t} &= \mu(1 - \phi_{1}) + \phi_{1} z_{t-1} + a_{t}\\ z_{t} &= 6.0423(1 - 0.6803) + 0.6803 z_{t-1} + a_{t}\\ z_{t} &= 1.931723 + 0.6803 z_{t-1} + a_{t} \end{aligned} We are NOT told, as one may suspect, that: \begin{aligned} z_{t} &= C + \phi_{1} z_{t-1} + a_{t}\\ z_{t} &= 6.0423 + 0.6803 z_{t-1} + a_{t} \end{aligned}
To calculate predicated values, denoted $\hat{z}_{t}$ (that's a "z hat" and not a "z tilde" which I already said represents deviations from the mean), we simply plug values into the following equation. Fortunately, this is done automatically in R, but here are examples using our estimates and the data shown in Table 1 above: \begin{aligned} \hat{z}_{t} &= \hat{C} + \hat{\phi}_{1} z_{t-1}\\ \hat{z}_{2} &= 1.931723 + 0.6803 (4.4) = 4.925043 \\ \hat{z}_{3} &= 1.931723 + 0.6803 (5.8) = 5.877463 \\ \hat{z}_{4} &= 1.931723 + 0.6803 (6.7) = 6.484733 \\ \hat{z}_{60} &= 1.931723 + 0.6803 (13.3) = 10.979713 \end{aligned} Notice that we do not calculate $\hat{z}_{1}$ because in that case we have $\hat{z}_{1} = \hat{C} + \hat{\phi}_{1} z_{0}$ and data is not available for $t=0$.

To calculate the residuals, denoted $\hat{a}_{t}$, for our AR(1) model we subtract the predicted values from the observed values. Again, this is done automatically in R, but here again are some examples: \begin{aligned} \hat{a}_{t} &= z_{t} - \hat{z}_{t}\\ \hat{a}_{2} &= 5.8 - 4.925043 = 0.874957\\ \hat{a}_{3} &= 6.7 - 5.877463 = 0.827537\\ \hat{a}_{4} &= 7.1 - 6.484733 = 0.615267\\ \hat{a}_{60} &= 6.2 - 10.979713 = -4.779713 \end{aligned} It has been shown here that estimating an AR(1) model using a realization with $t=60$ gives rise to 59 residuals. Now here is the problem...

Let's take a look at the residuals that are calculated in R using arima(). To do this, we can assign the model arima(x, c(1,0,0)) using our data to an object called x and then check out the attributes of object x; type in R attributes(x). We find that there is an attribute name called "residuals". These are the residuals of the estimated model. Viewing the residuals using x$residuals we find that 60 residuals have been calculated, not 59 as we should expect given the calculations performed above! R reports to us the residuals shown in Table 2:

Table 2 \begin{array}{|c|c|} \hline t & \hat{a}_{t} \\ \hline 1 & -1.2037423 \\ 2 & 0.8749301 \\ 3 & 0.8225200 \\ 4 & 0.6102564 \\ \vdots & \vdots \\ 60 & -4.7796768 \\ \hline \end{array} Where is this extra residual that R calculates, namely $\hat{a}_{1} = -1.2037423$, coming from? Recall that we do not calculate $\hat{z}_{1}$ and subsequently $\hat{a}_{1}$ because these calculations involve $z_{0}$ which we do not know. See Table 1, it's n.a.!

To decipher what R appears be doing, consider the following: \begin{aligned} \hat{z}_{t} &= \hat{C} + \hat{\phi}_{1} z_{t-1}\\ \hat{z}_{1} &= \hat{C} + \hat{\phi}_{1} z_{0} \end{aligned} If we treat $z_{0}$ as zero we end up with: \[\hat{z}_{1} = \hat{C} = 1.931723\] Now we can calculate a residual for $t=1$ in this way: \begin{aligned} \hat{a}_{t} &= z_{t} - \hat{z}_{t}\\ \hat{a}_{1} &= 4.4 - 1.931723 = 2.468277 \end{aligned} But this value certainly does not correspond to the value of $\hat{a}_{1}$ that R calculates; namely $\hat{a}_{1} = -1.2037423$.

OK, let's go back and again consider the following: \[\hat{z}_{1} = \hat{C} + \hat{\phi}_{1} z_{0}\] but let's rewrite it and work with the equation in this form: \[\hat{z}_{1} = \mu(1 - \phi_{1}) + \phi_{1} z_{0}.\] If we again treat $z_{0}$ as zero we end up with: \begin{aligned} \hat{z}_{1} &= \mu(1 - \phi_{1}) \\ \hat{z}_{1} &= 6.0423(1 - 1.931723) = -5.629749 \end{aligned} Calculating the residual using this predicted value, we get \[\hat{a}_{1} = 4.4 - 5.629749 = -1.229749\] which is close to R's calculation of $\hat{a}_{1}$, namely $\hat{a}_{1} = −1.2037423$.

I'm thus led to believe that this is what's going on when R reports this extra residual for $t=1$. I have tried a few other procedures, but they point in the wrong direction and end up way off the mark.

The question now to be raised is that when using R for time-series analysis, at the diagnostic stage of the Box-Jenkins methodology, should the analyst exclude this extra (non-zero) residual or keep it when performing diagnostic checks? Furthermore, ask yourself what happens to the residuals in R when an AR(2), say, with 60 observations is estimated.

Note: expect some small discrepancies between my calculated figures and R's calculated figures, this is down to rounding and shouldn't be much of an issue here.

How to install the tseries package

I recently installed the tseries package because I wanted to make use of its bds.test() function. This function allows one to carry out the BDS test. It's one of the most popular tests for testing the iid hypothesis and it can also be used to test for non-linearity.

When I tried to install the tseries package in R using install.packages(tseries), I received an error message saying that the package could not be found or something like that there was no such package. Installing packages in R is usually straight forward enough, so I was slightly puzzled on this occasion.

I figured my best bet would be to go directly to the CRAN website and download the tseries source package tarball file then install it using the Terminal in Ubuntu. After trying to install the package from the terminal, I received a second error message. This time I was informed that I'd need to first download and install another package called quadprog.

So, I downloaded and installed the quadprog package and then finally this enabled me to successfully install the tseries package. Some others have encountered the same problem when trying to install the tseries package, so in case you do to,  you may find it useful to know the commands that I used to get it all sorted out.

If you're not used to using the terminal, you may like to know that the wget command is a utility for non-interactive download of files from the Web - you simply type wget followed by the url of the file to be downloaded. The sudo command basically allows you to execute commands as a superuser or administrator on your computer. If you need to, you can read about R CMD INSTALL and how to install R packages here.