Forecasting Seasonal Data in R

A few posts ago, I talked about visualizing seasonal data in R. In that post I plotted data for the retail sector in Ireland and showed that December was an important month in the year for retailers. To recap, below is a plot of the data, obtained from the CSO, once again.

 

Economic statistics can tell us about the past and present, but if we are to plan for the future ahead of time, say, to control inventory, we're going to need forecasts of some sort. Forecasts can, of course, be qualitative or quantitative in nature, but here I will focus on producing point and interval forecasts using the univariate Box-Jenkins method. It's worth noting that properly constructed ARIMA models produce optimal univariate forecasts; the mean-squared forecast error is smaller than for any other univariate forecast (produced by univariate, linear, fixed coefficient models).

If we're going to construct a model we better be able to write it out, so let's have a quick review of some notation used for expressing seasonal non-seasonal multiplicative models. There are also seasonal non-seasonal additive models, but let's leave them aside for the purposes of this blog post.

The general form of seasonal non-seasonal multiplicative models can be written compactly as: \[ \phi_{p}(B) \Phi_{p} (B^{s}) \nabla^{d} \nabla_{s}^{D} \tilde{z} = \Theta_{Q}(B^{s})\theta_{q}(B)a_{t} \] where $\phi_{p}(B)$ is the nonseasonal AR operator, $\theta_{q}(B)$ is the nonseasonal MA operator, $\Phi_{p} (B^{s})$ is the seasonal AR operator, $\Theta_{Q}(B^{s})$ is the seasonal MA operator, $\nabla^{d}$ is the nonseasonal differencing operator, $\nabla_{s}^{D}$ is the seasonal differencing operator, and $\tilde{z} = z_{t} - \mu$, is the data expressed in deviations from the mean.

This sort of process is often referred to as an ARIMA(p,d,q)(P,D,Q)$_{s}$ process. Lower case letters p, d, and q denote the nonseasonal orders while the upper case letter P, D, Q denote the seasonal orders. The letter s denotes the type of periodicity; 4 for quarterly; 12 for monthly.

Using the Box-Jenkins method, I identified and estimated the following tentative model for the aforementioned CSO data: an ARIMA(1,1,0)(0,1,0)$_{12}$ model. The model can be expressed in the following form: \[(1 - \phi_{1}B)(1 - B)(1 - B^{12}) \tilde{z} = a_{t}\] or equivalently as: \[ \phi_{1}(B) \nabla \nabla_{12} \tilde{z} = a_{t} \] It's easy to see from either of these equations that the model consists of one non-seasonal autoregressive term, it is in first differences, and it has been seasonally differenced once with monthly periodicity. Furthermore, there are no non-seasonal moving average terms, nor are there any seasonal autoregressive terms or seasonal moving average terms.

I should say that this compact notation may look a bit daunting at first, but in actual fact it's far clearer than the common algebraic form normally used to express difference equations. For argument's sake, let's write the above model in common algebraic form. Using a bit of algebra, we get the following series of equations: \begin{aligned} (1 - \phi_{1}B) (1 - B) z_{t} - z_{t-12} &= a_{t}\\ (1 - \phi_{1}B) z_{t} - z_{t-12} - z_{t-1} + z_{t-13} &= a_{t}\\ z_{t} - z_{t-12} - z_{t-1} + z_{t-13} - \phi_{1}z_{t-1} + \phi_{1}z_{t-13} + \phi_{1}z_{t-2} - \phi_{1}z_{t-14} &= a_{t} \end{aligned} \begin{aligned} (z_{t} - z_{t-12}) - (z_{t-1} - z_{t-13}) &= \phi_{1}z_{t-1} - \phi_{1}z_{t-13} - \phi_{1}z_{t-2} + \phi_{1}z_{t-14} + a_{t}\\ (z_{t} - z_{t-12}) - (z_{t-1} - z_{t-13}) &= \phi_{1}(z_{t-1} - z_{t-13} - z_{t-2} + z_{t-14}) + a_{t} \end{aligned} Expressed in common algebraic form, this model is less easy to identify hence the preference for using the compact notation (it's not only about saving paper!).

After performing diagnostic checks on the ARIMA(1,1,0)(0,1,0)$_{12}$, I decided that, for the purposes of this blog post, the model could be used to produce forecasts. I have plotted a fit of the ARIMA(1,1,0)(0,1,0)$_{12}$ model below and the figure also shows both point and interval out-of-sample forecasts. Since I omitted a legend from the plot, the dashed line is the model fit, the solid black line is the actual data, and you should be able to make sense of the line (point forecasts) and shading (interval forecasts) in the out-of-sample period.

Rather than go through all of the steps of the analysis here, I've instead written an R script to show how I carried out the analysis using R. The script also contains the code that I used to colour in between the forecast intervals in gray. This wasn't so easy, but I managed to do it without using any packages from the library. In fact, to carry out the analysis I only resorted to using my own R function (presented in a previous post). The plotting scheme that I used to present the fit of the model and the out-of-sample forecasts above is one that I created using the base package graphics only.

Note of caution: admittedly, I identified other possible ARIMA(p,d,q)(P,D,Q)$_{s}$ models that I did not subsequently estimate or diagnostically check. If this were a thorough analysis of the data, these alternative models would also have to be entertained. In any event, this post was just meant to be a quick taster. As a larger exercise, I might compare this model to some of the other tentative models in a later post. In other words, the model presented here is only a candidate for the final model.

Univariate Box-Jenkins Analysis in R

Two graphical devices used in UBJ-analysis are the autocorrelation function and the partial autocorrelation. These tools are used at the identification stage of the Box-Jenkins method in order to identify tentative ARIMA models. This process is done by comparing the estimated acf and estimated pacf to a set of theoretical acf's and pacf's. The models chosen at this stage are only candidates for the final model - the Box-Jenkins method is an iterative process.

Inspection of the autocorrelation function at the identification stage can also help detect whether or not the time-series is stationary. If the autocorrelation function decays slowly then some order of differencing may be required.

At the diagnostic checking stage of the Box-Jenkins method, the residual autocorrelation function can be used as a guide to determining whether or not the residuals are independent. If they are not independent, it may be necessary to return to the identification stage.

To determine the statistical significance of the autocorrelation coefficients, we can test them individually (under the null that they equal zero) using t-tests. To carry out the t-tests we need to first calculate standard errors of the autocorrelation coefficient sampling distribution. Bartlett (1946) derived an approximate expression for doing this. The same can be done to the partial autocorrelation coefficients, although using slightly different standard errors.

The autocorrelation coefficients can also be tested for joint significance using the Q-statistic developed by Box and Pierce. The null this time is that all of the coefficients up to a certain number of lags are all equal to zero.

So, to perform a Box-Jenkins analysis of a time-series, we are going to need (at the least) to program the computer to do quite a few things for us. We'll need: autocorrelation coefficients, partial autocorrelation coefficients, standard errors, t-values, Q-statistics and p-values. We'll also need to have plots of the acf and pacf.

I'm not aware of any ready-made way of gathering all of this information in R, so to make Box-Jenkins analysis easier in R, I've designed a function that will do a variety of things for us automatically. The R function that I've created does the following:

• Calculates autocorrelation coefficients
• Calculates partial autocorrelation coefficients
• Calculates t-values based on Bartlett's approximations
• Calculates the Box-Pierce Q statistic
• Plots the acf and pacf in a single frame

The details of how the autocorrelation and partial autocrrelation coefficients are calculated along with the formulas that I used to calculate the standard errors, t-values and Q-stats will be documented in a separate post.

For now, a sample output of the program looks similar to the table below. In practice, the results appear as a dataframe in R, but the output can easily be exported to $LaTeX$, as shown here.

\begin{array}{c|cccc|cc} \hline Lag & ACF & T-VAL & Q-Stat & Prob & PACF & T-VAL \\ \hline 1 & -0.51 & -3.94 & 3.88 & 0.05 & -0.51 & -3.94 \\ 2 & ~~~0.22 & ~~~1.41 & 4.63 & 0.10 & -0.05 & -0.36 \\ 3 & -0.25 & -1.50 & 5.54 & 0.14 & -0.20 & -1.59 \\ 4 & ~~~0.16 & ~~~0.92 & 5.91 & 0.21 & -0.06 & -0.50 \\ 5 & -0.09 & -0.50 & 6.02 & 0.30 & -0.03 & -0.21 \\ 6 & ~~~0.09 & ~~~0.53 & 6.15 & 0.41 & ~~~0.02 & ~~~0.17 \\ 7 & -0.02 & -0.13 & 6.16 & 0.52 & ~~~0.07 & ~~~0.51 \\ 8 & ~~~0.05 & ~~~0.31 & 6.20 & 0.62 & ~~~0.09 & ~~~0.73 \\ 9 & ~~~0.05 & ~~~0.26 & 6.23 & 0.72 & ~~~0.18 & ~~~1.41 \\ 10 & -0.12 & -0.71 & 6.46 & 0.78 & -0.03 & -0.19 \\ 11 & ~~~0.10 & ~~~0.57 & 6.61 & 0.83 & ~~~0.05 & ~~~0.37 \\ 12 & -0.11 & -0.64 & 6.80 & 0.87 & -0.04 & -0.31 \\ 13 & ~~~0.16 & ~~~0.88 & 7.17 & 0.89 & ~~~0.05 & ~~~0.41 \\ 14 & -0.21 & -1.14 & 7.80 & 0.90 & -0.14 & -1.07 \\ 15 & ~~~0.22 & ~~~1.19 & 8.53 & 0.90 & ~~~0.03 & ~~~0.24 \\ \hline \end{array}
There were some issues with the plotting scheme associated with acf() and pacf(), so I've revamped them to produce figures like the one shown below. This plotting scheme has the advantage of showing both the acf and the pacf in the same figure; handy for comparing them to theoretical acf's and pacf's. I have also removed the zeroth lag from the acf because that was confusing for some people. The y-axis on both plots have been set to have the same scale (from -1 to +1) and the labeling is a slight improvement on the original too.

Note: since I prefer minimalism when presenting data for analysis, I was thinking about removing the blue color from the dashed horizontal lines. I may or may not make this change.

I should also stress that the t-values deserve far more attention than the dashed blue line. In fact, Box and Jenkins suggested a variety of practical warning levels to be wary of. For example, t-values of 1.6 at short lags during the identification stage and t-values of 1.25 for short lags in the residual acf!

These details ought not to be overlooked! It's the t-values that should be inspected closely and compared to the warning levels suggested by Box and Jenkins. Just having a quick gander at the plots to see if there are any spikes sticking out beyond the dashed blue lines simply does not cut it for a proper application of the Box-Jenkins methodology.

R Time Series Issue: Autocorrelation Function

In two earlier blog posts, I talked about using R for time series analysis and in particular I pointed out that the arima() function does some slightly strange things when calculating residuals.

In this post, I'd like to draw attention to the acf() function which plots the autocorrelation function for a particular time series. In actual fact, the default setting is that acf() will plot the autocorrelation function, but the plotting can be suppressed by specifying the argument plot=FALSE. If the plotting is suppressed then R will instead report the values of autocorrelation coefficients for up to 30 lags. If you want to analyse more than 30 lags, say 50, then you need to specify max.lags=50.

The standard formula for calculating autocorrelation coefficients is \[r_{k} = \frac{\sum_{t=1}^{n-k}(z_{t} - \bar{z})(z_{t+k} - \bar{z})}{\sum_{t=1}^{n} (z_{t} - \bar{z})^{2}}\] where $r_{k}$ is the autocorrelation coefficient at lag $k$, $z_{t}$ is a time-series, and $\bar{z}$ is the mean of the time-series.

The formula above is not the forumla for calculating autocorrelation coefficients because there are other ways of calculating $r_{k}$. Pankratz (1983) (p.35) refers me to Jenkins and Watts (1968) for a discussion and evaluation of the alternatives.

A problem is that the exact equation used by acf() itself is given in Venables and Ripley's Modern Applied Statistics with S and since I do not have access to this book, how do I know which formula is being computed when I call acf() in R? Reviewing the source code in RStudio by pressing F2 does not readily reveal the formula used by the function. Not knowing the exact details behind the acf() function is kind of like working with a black-box!

To bring the point further, if I were to present a piece of research that involved estimating autocorrelation coefficients using the acf() function and then someone came along and quizzed me about the formulas used to carry out the research, how could I respond? Er... I just used a function and I have no clue how it was programmed, sorry about that!? That wouldn't be a nice situation to find oneself in.

The basic point I'm trying to make is that it's good practice to be aware of the formulas being used in the analysis that you're carrying out. Having at least a fair idea about how the formulas could be used "by hand" is a welcomed bonus.

I'm not alone in trying to decipher what acf() actually does. Professor of Statistics Galit Shmuéli asked a question on Cross Validated also with the intention of finding out exactly what it does.

To get around the problem I decided to write my own function in R that will calculate autocorrelation coefficients for a series of data using the aforementioned standard formula. I've tested the function in order to confirm that it works correctly by replicating various results for estimated acf's reported in Pankratz (1983) - as well as some other books and papers.

Here is an R script which contains the function that I wrote to calculate autocorrelation coefficients. I've written the code in such a way that the code and the actual formula look similar to each other. To test the function, I've used two datasets from Box and Jenkins (1976). After you run the script, you can compare the results to those found in the book if you have a copy of it.

I've also devised a basic plotting scheme that overcomes a weird plotting convention associated with acf(). Professor Stoffer at the Department of Statistics at the University of Pittsburgh has already shone light on this issue on his website (see Issue #6).

How to simulate a Random Walk in R

I've put together this R script to demonstrate how to simulate a random walk in R. Actually, the script shows how to do a lot more than simulate a random walk. To give an overview of the content of the R script, I wrote code to simulate the following processes:

• Gaussian white noise
• Random walk
• AR(1) with positive phi
• AR(1) with negative phi
• AR(2) with two positive phi's
• MA(1) with positive theta
• MA(1) with negative theta
• MA(2) with two positive theta's
• ARMA(1,1) with positive phi and positive theta

I suggest that you edit the script and play around with the values and signs on the coefficients. For example, try to find out what happens when some of the models do not satisfy the stationarity conditions. What would the autocorrelation function look like in that case?

Keep in mind that the pure MA processes are always stationary even if the coefficients on theta are greater than one in absolute value. Despite the fact that pure MA processes are always stationary because they contain no AR terms, they may, however, not be invertible.

In case you need some quick revision... the conditions for stationarity and invertibility are...

As I already said, tinker around with the coefficients and don't forget to inspect the autocorrelation functions and partial autocorrelation functions to see how each process behaves.

Note that there is a function from the stats package for simulating ARIMA models. The function is called arima.sim(). I did not use this function in the script that I compiled for this blog post. The reason being I wanted to write code that would make sense to those who have not used R previously and to those who understand the theory behind the models, but who are perhaps new to programming languages.

Visualizing Seasonal Data in R

The end of the year is a very important time for retailers. A quick look at some data from the CSO database shows why. It's clear in the figure below that there is a significant spike in retail sales every December.

Now let's imagine that you want to analyse this time series in R, but rather than looking at the whole series, you want to focus on the data for the month of December. In that case, you might like to make a plot similar to the one I've created below.

To be able to create a graph like the one above you need to extract the data points for each December from the whole series. There are probably many ways of doing it, but I've written a short R script that demonstrates just two possible methods. The first method is a bit hackish and the second method is a bit more systematic. The script also shows how I labelled the x-axis.

Note that the basic idea behind the first method is to extract every nth element of a vector in R while the second method involves sub-setting a ts object.

How to run a Hausman Specification Test in R

Consider the following structural equations: \begin{aligned} M_{t}^{d} &= \beta_{0} + \beta_{1} Y_{t} + \beta_{2} R_{t} + \beta_{3} P_{t} + u_{1t}\\ M_{t}^{s} &= \alpha_{0} + \alpha_{1} Y_{t} + u_{2t} \end{aligned} where $M$ = money, $Y$ = income, $R$ = rate of interest, $P$ = price, and the $u$'s are error terms. Assume that $R$ and $P$ are exogenous and $M$ and $Y$ are endogenous.

Let's begin by deriving the reduced form equations for the model.

By the equilibrium condition, $M_{t}^{d} = M_{t}^{s}$, we have: \begin{aligned} \beta_{0} + \beta_{1} Y_{t} + \beta_{2} R_{t} + \beta_{3} P_{t} + u_{1t} &= \alpha_{0} + \alpha_{1} Y_{t} + u_{2t}\\ \beta_{1} Y_{t} - \alpha_{1} Y_{t} & = \alpha_{0} - \beta_{0} - \beta_{2} R_{t} - \beta_{3} P_{t} + u_{2t} - u_{1t}\\ Y_{t} (\beta_{1} - \alpha_{1}) &= \alpha_{0} - \beta_{0} - \beta_{2} R_{t} - \beta_{3} P_{t} + u_{2t} - u_{1t}\\ Y_{t} &= \frac{\alpha_{0} - \beta_{0}}{\beta_{1} - \alpha_{1}} - \frac{ \beta_{2} R_{t}}{\beta_{1} - \alpha_{1}} - \frac{\beta_{3} P_{t}}{\beta_{1} - \alpha_{1}} + \frac{u_{2t} - u_{1t}}{\beta_{1} - \alpha_{1}}\\ Y_{t} &= \frac{\alpha_{0} - \beta_{0}}{\beta_{1} - \alpha_{1}} - \frac{ \beta_{2}}{\beta_{1} - \alpha_{1}}R_{t} - \frac{\beta_{3} }{\beta_{1} - \alpha_{1}}P_{t} + \frac{u_{2t} - u_{1t}}{\beta_{1} - \alpha_{1}}\\ Y_{t} &= \Pi_{0} + \Pi_{1}R_{t} + \Pi_{2}P_{t} + w_{t} \end{aligned} where $\Pi_{0} = \frac{\alpha_{0} - \beta_{0}}{\beta_{1} - \alpha_{1}}$ , $\Pi_{1} = - \frac{ \beta_{2}}{\beta_{1} - \alpha_{1}}$, $\Pi_{2} = - \frac{\beta_{3} }{\beta_{1} - \alpha_{1}}$ , and $w_{t} = \frac{u_{2t} - u_{1t}}{\beta_{1} - \alpha_{1}}$.

Substituting $Y_{t} = \Pi_{0} + \Pi_{1}R_{t} + \Pi_{2}P_{t} + w_{t}$ into either the supply function or the demand function, we can solve for the second reduced form equation which is: \[M_{t} = \Pi_{3} + \Pi_{4}R_{t} + \Pi_{5}P_{t} + v_{t}.\] Now let's find out if the structural equations are identified or not identified.

A structural equation can be said to be identified or unidentified. Furthermore, identified equations can be classified as just-identified, overidentified, or, underidentified. To determine the identification of structural equations we can use the order and rank conditions of identification. The order condition of identifiability is defined (in Gujarati's Basic Econometrics) as follows:

In a model of $M$ simultaneous equations, in order for an equation to be identified, it must exclude at least $M - 1$ variables (endogenous as well as predetermined) appearing in the model. If it excludes exactly $M - 1$ variables, the equation is just identified. If it excludes more than $M - 1$ variables, it is overidentified.

There are two endogenous variables in the model, namely, $Y_{t}$ and $M_{t}$, so in this case $M = 2$. This means that the demand function should omit at least $2 - 1 = 1$ variable in order to be identified. Since the demand function excludes $0$ variables from the model, it is unidentified. The supply function, however, omits $2$ variables, and is thus classified as overidentified.

With the preliminaries out of the way, we can now focus on the Hausman specification test.

To carry out the test on the supply function, there are two steps.

The first step is to regress $Y_{t}$ on $R_{t}$ and $P_{t}$. So, estimate by OLS the reduced form equation: $Y_{t} = \Pi_{0} + \Pi_{1}R_{t} + \Pi_{2}P_{t} + w_{t}$. Then, obtain the residuals, denoted $\hat{w_{t}}$, from this regression.

The second step is to regress $M_{t}$ on $Y_{t}$ and $\hat{w_{t}}$, and then perform a t-test on the coefficient of $\hat{w_{t}}$. The null hypothesis is that there is no simultaneity.

OK. Let's do this in R using some actual data.

Here is an R script that I've written to show how it can be done. It happens to be rather easy.

The dataset that I've used comes from an empirical exercise at the end of Chapter 19 in Gujarati's Basic Econometrics. The dataset is downloadable as a zip file from the book's website and the specific .txt file for the data is named Table 19.4.txt. You can go fetch it yourself or to save you the hassle, here is the data.

Before you run my R script you will need to change "/home/name/Desktop/hausmanData.txt" in the scan() function to the location of the dataset saved on your computer.

The term "pure economics"

In this earlier blog post, I used the term pure economics. Admittedly, I did not make it explicit as to what I was referring to when I used the term. Nevertheless, I trust the reader could follow my argument.

Associated with the term, pure economics, is the French economist, Léon Walras. In fact, arguably his most important work is titled: Elements of Pure Economics. He also used the terms applied economics and social economics and made distinctions between each of them.

According to Wikipedia (forgive me on this one, I do not have access to the actual source):

Léon Walras planned to organize his main work into volumes on "pure," "applied," and "social" economics. Jaffé (1983) describes Walras's plan as involving making a distinction between that which is true, is useful, and is just. In using the term true, Walras referred to propositions that necessarily followed from the nature of things. Pure economics then involves pure logic. Applied economics involves examining ways to achieve practical goals and requires the making judgments about whether or not the logic of pure economics was relevant to the real world. Social economics also presumed pure economics, but dealt with a different range of questions than did applied economics.

Now, when I used the term, I did so much in the same way as did Léon Walras. Interestingly, Vilfredo Pareto, an Italian economist, used the term in a similar vain.

To go over what I said earlier, pure economics involves developing an axiomatic system and as long as that system has a logical development, it holds true. It holds true in itself, but may, however, not be true about the real world. What we observe may be on the contrary to or unexplained by the implications asserted by the axiomatic system. In this way, pure economics - an axiomatic system - only tells us that an economy with such and such properties, subsequently has such and such other properties. All propositions about the real world, belong to applied economics, or, econometrics. Again, pure economics illustrates while econometrics discovers.

Having mentioned Léon Walras, it's probably a good time to say that his name is pronoucned with the letter 's'. Some have claimed that his name is pronounced 'Val-rah' or 'Wal-rah', but both instances are incorrect.

The source that tells me that this is a footnote to the translator's forward of Elements of Pure Economics! The footnote says:

I was told by Walras's daughter, the late Mlle Aline Walras, that the final s is sounded in the correct pronunciation of the name.

You simply never know what you'll find out when reading a book!

Transforming imaginary and complex numbers

Difference equations are the topic of interest in Chapter 1 of Hamilton's Time Series Analysis. The chapter includes a section that outlines the general solution of pth-order difference equations with distinct eigenvalues. On page 14, Hamilton says that:

interesting possibilities arise if some of the eigenvalues are complex

and then a series of equations follow this statement. To guide the reader through the equations, footnotes are provided at the bottom of the page that refer the reader to the mathematical appendix at the back of the book.

It can be a little bit difficult to figure out what's going on in the equations on pages 14 to 16 of the book, so I thought it would be a good idea to write up some notes in order to make the material more understandable.

The equations involve transforming imaginary and complex numbers to trigonometric functions. It's straightforward stuff really, but since no reference is made to the names of the formulas being used, it may leave you scratching your head as to what's going on.

So, our aim here is to understand how to transform imaginary and complex numbers to trigonometric functions. To do this, we're going to need to familiarize ourselves with the following:

1. Expressing Cartesian coordinates in polar coordinate form; 
2. Euler's formula; 
3. De Moivre's Theorem.

Let's begin with expressing Cartesian coordinates in polar coordinate form.

$a$ and $b$ in the diagram below are called Cartesian coordinates. These can be expressed in terms of $\theta$ and $R = \sqrt{a^{2}+b^{2}}$ which are called polar coordinates. To do this, we use the following formulas: \[a = R~cos~\theta,~ ~~~~~ b = R~sin~\theta,~ ~~~~~ R>0.\] Applying these formulas to the conjugate complex number $(a \pm b i)$, we have: \[a \pm b i = R~cos~\theta \pm i~R~sin~\theta = R(cos~\theta \pm i~sin~\theta) \] and this explains Equation [1.2.36] in Hamilton.

Second on the list is Euler's formula. This tells us that: \[ e^{\pm i \theta} = cos~\theta \pm i~sin~\theta.\] Thus we can express $(a \pm bi) = R(cos~\theta \pm i~sin~\theta)$ as \[a \pm bi = R(cos~\theta \pm i~sin~\theta) = R e^{\pm i \theta}.\] The last item on the list is De Moivre's theorem, but first let's see how a conjugate complex number is raised to the jth power. It is done as follows: \[(a \pm bi)^{j} = (Re^{\pm i \theta})^{j} = R^{j} e^{\pm i j \theta}.\] Or, using De Moivre's theorem we have: \[(a \pm bi)^{j} = R^{j}(cos~j~\theta \pm i~sin~j~\theta).\] and now we've arrived at Equation [1.2.37] in Hamilton.

Although not explicitly mentioned in the textbook, it's Euler's formula and De Moivre's theorem that lead to both Equation [1.2.37] and Equation [1.2.38] in Chapter 1 of Hamilton.

The rest of the section on the general solution of pth-order difference equations with distinct eigenvalues involves no more than basic substitution and a little bit of algebra. Despite the fact that the notation looks a bit "heavy", it's easy to follow once you know to use Euler's formula and De Moivre's theorem.

Vector Algebra (complete proofs)

Open up a decent undergraduate Linear Algebra textbook (my favourites are Fraleigh & Beauregard and Strang) and you'll quickly come across definitions of vector addition, vector subtraction, and scalar multiplication. A list of the properties of these vector operations usually follow and while one or two proofs are given, the rest are left up to the reader as exercises. That shouldn't be too much of a problem to any mathematics student.

You will, however, encounter such proofs - again left as exercises - in graduate textbooks for economics and finance; one such example being Harrison & Waldron. They also appear in this Advanced Engineering Mathematics textbook that I happened to pick up at a bargain price in Oxfam!

Familiarity with vector algebra is pretty essential to anyone studying a numerical subject and up to a certain level students can advance in their chosen field without being able to write proofs.

There comes a time, though, when it is necessary for the student to be able to understand and write proofs. That's got to be a nightmare for anyone who has not yet grasped how to do it. Even worse, there's probably no course designed to teach it to someone for the first time. Gone are the days of the trivium.

Mathematics students at Virginia Commonwealth University are quite the lucky ones because they have a whole course designed to teach them how to write mathematical proofs. The course is called Introduction to Mathematical Reasoning and the book that's used for the course can be found here (it's an excellent textbook!).

Students of economics are, perhaps, not so lucky. Some may already lack math skills (or rather, the confidence in their math ability), and I've not really come across any material that aims to teach them how to write proofs. It's kind of like they're left to their own devices to try and work it out for themselves. A tad unfair - remember, some 3rd year math students take a whole course in this art!

Now, to write proofs for the properties of vector algebra is very easy and as basic as you can get, but since someone out there (maybe a student of economics) might be having trouble with writing them, I've written up complete proofs and added comments to them that explain how each step in the proof was made. I think they're quite instructive and definitely more detailed than any I've seen in a variety of textbooks. Feedback welcome.

Note: proofs can be written for the properties of matrix algebra in a similar fashion. Following the tradition, I'll leave them as an exercise. Follow the technique I used for vector algebra and you'll have no problem.

The notion of "encompassing"

Christopher L. Gilbert in his paper Professor Hendry's Econometric Methodology gives an outline of six criteria for model selection. These six criteria originate from Hendry and Richard's paper The Econometric Analysis of Economic Time Series. The six criteria are:

1. Models must be data admissible.
2. Models must be consistent with theory.
3. Satisfactory models should have regressors that are (at least) weakly exogenous. 
4. Satisfactory models must  exhibit parameter constancy.
5. A satisfactory model must be data coherent.
6. A satisfactory model should encompass a wide range of rival models.

A model that satisfies all six of these criteria is said to be congruent with the data.

The sixth criterion about encompassing is the one which I want to focus on in this post. Gilbert explains encompassing in the following paragraph:

One model may be said to encompass a second model if it can explain the second model's results. A good model should not only explain the data, but should also explain the successes and the failures of rival models in accounting for the same data. The Special Theory of Relativity encompasses Newtonian mechanics (which could not explain the Michelson-Morley results), but is itself encompassed by the General Theory of Relativity (as it cannot explain the visibility of the planet Mercury behind the Sun during a solar eclipse).

To bring the reference to physics up to date, the search is now on for theoretical physicists to find a Theory of Everything that encompasses both General Relativity and Quantum Mechanics. One contender for the job is a theory fathered by Professor Leonard Susskind at Stanford University called string theory.

A basic example of encompassing would be taking a model with a variable whose coefficient is set equal to one and then developing a second model in which the coefficient is allowed to take on a range of values, say, between zero and one. In this case, the second model encompasses the first model because when the coefficient in the second model happens to equal to one, it reduces to the first model. The basic idea is that it's possible to introduce parameters into existing equations to give them a little bit more of a kick.

University of Yale Professor Charles Bailyn has spent some of his lecture time discussing scientific revolutions and the notion of encompassing - see this Astrophysics lecture. He does an excellent job at explaining it, so it's worth repeating. In the lecture, he said:

But Einstein didn't say that Newton was wrong. That wasn't what happened. That wasn't how a revolution works. You revolutionize science only if the new theory encompasses what was previously known--the previous theories. And this is this business of the parameter, gamma, where under appropriate circumstances, gamma is one. And under those circumstances, you reproduce the theory that was already in existence. And so, you know, there's been hundreds and hundreds of years of Newtonian physics, all of which is basically right, one way or another. So whatever your new idea is has to encompass all the things that were right about the old idea. And so, you don't overturn Newton. What you do is, you demonstrate that Newton's theory is only a part of a bigger theory. And then, the bigger theory is the new theory and it encompasses new kinds of data and explains new kinds of things, while simultaneously not overturning all the things that you know are right about the old theory. 

And this is the thing that the people who wish they were Einstein, but aren't, don't understand. And so, anytime you see something where you're--where they're going to revolutionize all of science, it's certainly wrong. Because there's a large fraction of science in current existence that is correct. And you have to understand that. And you have to encompass it and go beyond that--but it's not going to get overturned. The planets aren't going to stop going in Keplerian orbits because somebody has a new idea. And so, it has to encompass the old idea.

Note the reference to the parameter gamma which equals one in a Newtonian framework, whereas in General Relativity, gamma can be some other number, a number other than one. The same idea of introducing new parameters into old equations can, of course, be applied to economics.

History shows us that encompassing has played an important role (and will continue to do so) in advancing our understanding of physics and of the world in general. Furthermore, there's little doubt that encompassing is important for the future of economics and it's reassuring to me that it's recognized, voiced, and supported by those at the top of the field.

Unfortunately, however, it's easy to be critical of economics and econometrics and there will always be cranks who'll stand up to say that everything we've done up to now has been wrong!

Average Economic Regression

There's a sort of textbook or stylized view of econometrics that is often taught to students. It's quite simplified and is referred to as average economic regression (AER) by Christopher L. Gilbert in his paper Professor Hendry's Econometric Methodoogy. It goes a little like this (in my own words):

To build an econometric model we first consider what economic theory has to say. Having done that, we throw some variables into a single equation, add an error term, and then all we have to do is obtain estimates for the beta coefficients. Job done. Well, nearly. Before we finish, we must solve some potential problems, such as low Durbin-Watson statistics, wrong signs, and insignificant coefficients. To deal with these problems, we chuck some variables out of the model and possibly apply the Cochrane-Orcutt transformation to help deal with serial correlation in the error term. One eventually emerges from the econometrics basement serving on a golden platter the final model that has all of the correct signs, with t-values greater than 2, and an R-squared that is acceptably high for the model builder.

Gilbert explains what's wrong with this. He says:

The problem with the AER approach is that we are using econometrics to illustrate the theories which we believe independently. The alternative might be to use econometrics to discover which views of the economy (or market) are tenable and to test, scientifically, the rival views. Our normal concept of science tends to be based on the natural, and in particular, the physical sciences, and does not obviously translate to the social sciences. This alternative programme for econometrics raises the question of what is required for economics to be scientific. 

Professor Edward Leamer also views textbook econometrics with some skepticism. He says:

Theoretical econometricians have interpreted scientific objectivity to mean that an economist must identify exactly the variables in the model, the functional form, and the distribution of the errors. Given these assumptions, and given a data set, the econometric method produces an objective inference from a data set, unencumbered by the subjective opinions of the researcher. 

This advice could be treated as ludicrous, except that it fills all the econometric textbooks. Fortunately, it is ignored by applied econometricians. The econometric art as it is practiced at the computer terminal involves fitting many, perhaps thousands, of statistical models. One or several that the researcher finds pleasing are selected for reporting purposes. This searching for a model is often well intentioned, but there can be no doubt that such a specification search invalidates the traditional theories of inference. The concepts of unbiasedness, consistency, efficiency, maximum likelihood estimation, in fact, all the concepts of traditional theory, utterly lose their meaning by the time an applied researcher pulls from the bramble of computer output the one thorn of a model he likes best, the one he chooses to portray as a rose. The consuming public is hardly fooled by this chicanery. The econometrician's shabby art is humorously and disparagingly labelled 'data mining', 'fishing', 'grubbing', 'number crunching'. A joke evokes the Inquisition: 'If you torture the data long enough, Nature will confess' (Coase). Another suggests methodological fickleness: 'Econometricians, like artists, tend to fall in love with their models' (wag unknown). Or how about: 'There are two things you are better off not watching in the making: sausages and econometric estimates'. 

This is a sad and decidedly unscientific state of affairs we find ourselves in. Hardly anyone takes data analyses seriously. Or perhaps more accurately, hardly anyone takes anyone else's data analyses seriously. Like elaborately plumed birds who have long since lost the ability to procreate but not the desire, we preen and strut and display our t-values.

My interpretation of the relationship between economic theory and econometrics tells me that a distinction should be drawn between what economic theory tells us and what econometrics can tell us. To lay the foundations of my argument, let me introduce the notion of a formal implication via a passage from Bertrand Russell's Principles of Mathematics:

[But] since the growth of non-Euclidean Geometry, it has appeared that pure mathematics has no concern with the question whether the axioms and propositions of Euclid hold of actual space or not: this is a question for applied mathematics, to be decided, so far as any decision is possible, by experiment and observation. What pure mathematics asserts is merely that the Euclidean propositions follow from the Euclidean axioms - i.e. it asserts an implication: any space which has such and such properties has also such and such other properties. Thus, as dealt with in pure mathematics, the Euclidean and non-Euclidean Geometries are equally true: in each nothing is affirmed except implications. All propositions as to what actually exists, like the space we live in, belong to experimental or empirical science, not to mathematics; when they belong to applied mathematics, they arise from giving to one or more of the variables in a proposition of pure mathematics some constant value satisfying the hypothesis, and thus enabling us, for that value of the variable, actually to assert both hypothesis and consequent instead of asserting merely the implication. We assert always in mathematics that if a certain assertion p is true of any entity x, or of any set of entities x, y, z, . . ., then some other assertion q is true of those entities; but we do not assert either p or q separately of our entities. We assert a relation between the assertions p and q, which I shall call formal implication. (Chapter 1, p.5)

Although Bertrand Russell is addressing the nature and bounds of pure mathematics and applied mathematics, an analogous story could be told of pure economics and econometrics.

Like pure mathematics, what pure economics deals with are sets of axioms and from those axioms follow propositions. Those familiar with Andreu Mas-Colell et. al's Microeconomic Theory (or Euclid's Elements for that matter!) will no doubt be comfortable with this sort of axiomatic system and the use of deductive logic.

The pure economists' concern is to come up with axioms and ensure that the progression from axioms to propositions will be carried out in a logical manner - in the parlance of formal logic, the inferential steps should be "watertight".

So, with one set of axioms we have a set of accompanying propositions and formal implications. As long as the axiomatic system is developed in a logical manner, it may be said to be true, albeit within its own scope - or space. Another set of axioms will generate another set of accompanying propositions and formal implications. Some of these may intersect the sets associated with the first axiomatic system. Indeed, one axiomatic system may encompass another one. Again, this second axiomatic system may be said to be true within its space.

That's what axiomatic systems in pure economics tell us. They assert implications that hold within a given space, or rather, in a particular model of an economy. Whether or not these implications conform to actual space - or actual economies - is another thing. The question to be asked then is why should an econometric model which aims to help us discover facts about actual economies be constrained to corroborate with an axiomatic system which may not necessarily hold true or provide a satisfactory description of actual economic phenomena? It probably shouldn't be that way because then we're using econometrics to illustrate rather than to discover - a reiteration of Gilbert.

Lawrence R. Klein, puts it nicely in Chapter 1 of his Economic Fluctuations in the United States (1921–1941). He says:

The purpose in building econometric models is to describe the way in which the system actually operates. [...] We want also to discovered the best possible theory or theories which explain the fluctuations that we observe. If we know the quantitative characteristics of the economic system, we shall be able to forecast with a specified level of probability the course of certain economic magnitudes such as employment, output, or income; and we shall also be able to forecast with a specified level of probability the effect upon the system of various economic policies. In the course of the search for models which are suitable for our purposes of forecasting and making policy recommendations, we shall inevitably have to consider several alternative theories as admissible hypotheses. The acceptance or  rejection of these hypotheses in the course of our search for truth will be our contribution to the problem of testing business-cycle theories. This is an important problem for which contributions are needed, because too often writers have made bold statements about the operation of the system without examining the factual data to determine whether or not their statements are true.* As a result we have had  continuing economic fluctuations and continuing publication of business cycle theories.

[*the emphasis is my own and does not appear in the original source]

It seems that an econometrican who builds models that are forced to conform to theories derived from an axiomatic system of pure economics is going astray from finding out about how actual economies work. Think for moment that Platonic solids and polygons from Euclid's Geometry can be described as having such and such properties and that's a useful exercise in itself, but believing that such solids and shapes occur in reality, is perhaps, not so wise. Go find me a perfect triangle!

Non-Euclidean Geometries, which permit shapes such as triangles whose angles sum to less than 180 degrees (hyperbolic geometry) or to greater than 180 degrees (elliptic geometry), may conform more closely to actual space than what Euclid's system tells us. As an example, an implication of Einstein's theory of general relativity is that physical space itself is not Euclidean, and Euclidean space is a good approximation for it only where the gravitational field is weak.* Importantly, each geometry, Euclidean or non-Euclidean, is true within its own framework, but their implications may hold differently in actual space. This is just a restatement of what Bertrand Russell said.

Coming back to the average economic regression and a situation in which a coefficient turns out to have the wrong sign. Well, if that's what the data tells us, maybe that's just the way it is for this particular economy - data and observation assumed to be accurate, i.e. assuming no identification problem. If economic theory suggests the sign of the coefficient should be otherwise, then the implications from that axiomatic system may not hold for the actual economy under scrutiny. If this is the case, it should follow that a new axiomatic system is built (which encompasses the old one) to take into account the new information we've discovered from econometrics. This shouldn't be thought of as a step back for economic theory, but an opportunity to advance the subject! Again, Newton didn't overturn Kepler, and Einstein didn't overturn Newton; each encompassed the other.

Similar to how applied mathematics deals with informing us about actual space, econometrics ought to aim at discovering information about economies that actually exist - freed from any bounds that pure economics may set upon itself. A man of considerable stature in economics once said: when the facts change, I change my mind, what do you do, sir?

Well, it's a bit like that. If revisions (or extensions) need to be made to economic theory, be open to making them. Now, it's not a case of scrap what we've already got. It's more so a case of expanding on what some great minds have already developed and kindly handed down to us. Nor is it a case that we're permitted to run any old nonsensical regressions. Sensible regressions that can be explained in terms of some (possibly hitherto unknown) axiomatic system are what we should be seeking to build.

I'll conclude the argument with four points.

• Pure mathematics illustrates. 
• Pure economics illustrates. 
• Applied mathematics discovers. 
• Econometrics discovers.

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.