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.