CRAN Task Views

In the previous post, I provided a list of R packages that I find useful for time-series analysis, economic and financial research. There are actually expert summaries of R packages by subject on CRAN. This is a really useful resource and I'd highly recommend checking it out.

Visit: http://cran.r-project.org/web/views/

R Packages for Econometrics & Time-Series Analysis

During my masters I learned how to use STATA. The STATA documentation is very good, particularly the [TS] manual. I like using STATA, but over time, I came across some tasks that I could not readily perform in STATA. Often, I wanted to make a unusual calculation or estimate a particular model, but I found that these tasks would be difficult or simply impossible to perform in STATA.

Essentially, STATA is more or less a "canned packaged" statistical tool and it has its limitations. There is, of course, the programming language version of STATA called MATA, but at the time I was fairly uninformed about its capabilities. Having since programmed using MATA, I'd still say that the data analysis and modelling tasks that interest me would be difficult to perform in MATA - even for an advanced user.

For the reasons above (and more), I started using the R programming language. R is not a "canned packaged" statistical software; it is a programming language. It is used for statistical analysis, but it has many other uses and different people use it for different reasons. Some people use it for econometrics, others for biometrics, mathematical modelling, statistics, and so on. On the whole, the limitations of R basically come down to the capabilities of the user; if you can program then the sky is the limit.

My general approach to using R and modelling data is that I like to write the code myself. That is, I try not to rely on add-in's or "packages" that have been written by other users of the language. In this way, I know exactly all of the formulas being used in my analysis since I program them myself. This gives me great confidence to be able to stand over the figures and results of my analysis. It would be difficult to ask me a question about where a particular figure is coming from or how it has been made that I could not answer with confidence.

Having said that, there are some really smart people out there writing R packages; I'm talking about leaders of the field in time-series analysis, statistics, etc. So, given that such people's work ought to be deemed reliable, an R user should probably use their packages if they are well documented, tested and trusted.

In other words, writing your own code some of the time is great, but writing it ALL of the time is probably an inefficient way to allocate your time. For instance, if an algorithm has been translated into R code, say, for estimating ARIMA or VAR models then you may as well take advantage of it - not without having read the source documentation, though, of course! Using a black-box for data analysis or forecasting could spell trouble for numerous reasons!

The problem is then, since there are bucket-loads of R packages which ones should you use? Well, here I'd like to suggest a few R packages that I have come across and found useful for time-series analysis and econometrics. In the next few posts I will provide some examples of what they can do.

In no particular order, they are as follows:

• Ecdat: Data sets for econometrics
• pear: Package for Periodic Autoregression Analysis
• vars: VAR Modelling
• strucchange: Testing, Monitoring, and Dating Structural Changes
• dse: Dynamic Systems Estimation (time series package)
• KFAS: Kalman Filter and Smoother for Exponential Family State Space Models
• dlm: Bayesian and Likelihood Analysis of Dynamic Linear Models
• tsDyn: Nonlinear time series models with regime switching
• ccgarch: Conditional Correlation GARCH models
• rugarch: Univariate GARCH models
• dlnm: Distributed Lag Non-linear Models
• urca: Unit root and cointegration tests for time series data
• apt: Asymmetric Price Transmission (apt)
• quantmod: Quantitative Financial Modelling Framework
• FinancialInstrument: Financial Instrument Model Infrastructure for R
• TTR: Technical Trading Rules
• RandomFields: Simulation and Analysis of Random Fields

Non-linearity in the Irish Stock Market

In this post, I'm going to go through an application of the BDS test using the R programming language. To be specific, I will test for non-linearity in the Irish Stock Market.

The application will involve daily data for the ISEQ index. The data spans from 1999-10-04 to 2013-05-03. There are a total of 3444 observations. Recall from my previous post that for the BDS test to have a reliable performance, the sample size should be greater than 500 observations. The data was obtained from Yahoo! Finance and you can download a .txt file of the data here.

The method that I am going to use involves taking differences of the logs of the data and then applying the BDS test to this transformed series. Note, however, that there is another method of using the BDS test to test for non-linearity which is to fit an ARIMA model (or say a GARCH model) and then to apply the BDS test to the residuals of the model. For further examples of the method that I am going to use in this post, see Chu, P.K.K. (2001) or these course notes by Eric Zivot from the University of Washington.

Let's begin by plotting the data to see what we're working with. The figure below shows a plot of the raw data.

    

The data is shown in logged form in the next figure. Generally speaking, taking logs has the effect of stabilizing the variance.

And here is the data in differences of the logs (which can be interpreted as daily percent change):

Now let's test the data expressed in differences of the logs for non-linearity. We can do this in R by using the bds.test() function from the tseries package. The hypothesis of the test is as follows:

• H$_{0}$ : the series is i.i.d. 
• H$_{1}$ : the series is not i.i.d. and this implies that non-linearity is present. 

The BDS test statistics are shown in the first table below and the second table contains the corresponding p-values. Note that we reject the null hypothesis at the 5% significance level when the BDS test statistic is greater in absolute value than 1.96.

BDS Test Statistics 

\begin{array}{rrrrr} \hline m & \epsilon = 0.5 & \epsilon = 1.0 & \epsilon = 1.5 & \epsilon = 2.0 \\ \hline 2 & 13.81 & 15.47 & 15.59 & 14.15 \\ 3 & 19.99 & 21.41 & 21.49 & 20.08 \\ 4 & 25.16 & 25.43 & 24.62 & 22.88 \\ 5 & 30.00 & 28.94 & 27.12 & 25.07 \\ \hline \end{array}

p-values

\begin{array}{rrrrr} \hline m & \epsilon = 0.5 & \epsilon = 1.0 & \epsilon = 1.5 & \epsilon = 2.0 \\ \hline 2 & 0.00 & 0.00 & 0.00 & 0.00 \\ 3 & 0.00 & 0.00 & 0.00 & 0.00 \\ 4 & 0.00 & 0.00 & 0.00 & 0.00 \\ 5 & 0.00 & 0.00 & 0.00 & 0.00 \\ \hline \end{array} The tables above show that all of the test statistics are greater than 1.96, so we should reject the null hypothesis of i.i.d. The p-values provide further evidence for rejecting the null.

These results suggest that the time-series data for the ISEQ index is non-linearly dependent. This evidence could be used to justify exploring the possibility of building a non-linear model for the series. Keep in mind that we'd normally like to build simpler linear models and not necessarily complicate matters by building non-linear models, but based on the evidence presented here, perhaps a non-linear approach to modelling is appropriate for this particular series.

Here is my R script for this analysis. You can use it to replicate what I've done in this post. I've also included in the script the code used for creating the LaTeX tables.

BDS Test for Non-Linearity

A widely used test for non-linearity is the BDS test due to Brock et al. (1987), but later published as Brock et al. (1996). Strictly speaking, the BDS test is a test of independence (a test of the i.i.d hypothesis), but it can be used to test for non-linearity and also as a test of model misspecification.

In this post I'd like to give a brief overview of the BDS test and at the end I'll provide an application of carrying out the BDS test in R. To introduce the BDS test, let me begin by saying that the method of the BDS test has its roots in the concept of spatial correlation from chaos theory.

Let's assume that $\{y_{t}\}$ is a strictly stationary process. Now consider the following vector called the m-history of the process: \[{\bf y}_{t}^{m} = (y_{t}, y_{t-1}, \ldots, y_{t+m-1})\] The construction of the m-history means that the univariate series $y_{t}$ has been embedded into a higher dimensional space. The purpose of this operation is to be able to examine the spatial correlation of the time-series. Notice that the vector, ${\bf y}_{t}^{m}$, consists of overlapping series. This becomes clear if we write out the following equations: \begin{aligned} y_{1}^{m} &= (y_{1}, y_{2}, \ldots, y_{m})\\ y_{2}^{m} &= (y_{2}, y_{3}, \ldots, y_{m+1})\\ & ~\vdots\\ y_{T-m+1}^{m} &= (y_{T-m}, y_{T-m+1}, \ldots, y_{T}) \end{aligned} To demonstrate the concept of embedding (which is due to Dutch mathematician Floris Takens), here is an example from Cromwell et al. (1994). Suppose we have a time-series, $y_{t}$, where $t = 1, 2, \ldots, 5$ and $m = 3$. Embedding this series into higher dimensions results in the following sequence of overlapping series: \begin{aligned} y_{1}^{3} &= (y_{1}, y_{2}, y_{3})\\ y_{2}^{3} &= (y_{2}, y_{3}, y_{4})\\ y_{3}^{3} &= (y_{3}, y_{4}, y_{5}) \end{aligned} In this example, the embedding operation creates three new three dimensional vectors: $y_{1}^{3}, y_{2}^{3}, y_{3}^{3}$. A requirement is that the vectors have the same length. Once the series has been embedded, we need an instrument to allows us to measure (a metric).

The correlation integral, introduced by Grassberger and Procaccia (1983), is a measure of spatial correlation. It measures the distance between the points in the m-dimensional space constructed by the embedding procedure; in the example above, it measures the distance between points in three dimensions. The correlation integral can also be interpreted as a measure of the frequency with which temporal patterns are repeated in the data.

The correlation integral at embedding dimension m is defined as: \[ C_{m, T} (\epsilon) = T_{m} (T_{m} - 1)^{-1} \sum_{t \neq s} I[y_{t}^{m}, y_{s}^{m}, \epsilon ] \] where $T_{m} = T - m + 1$, $t$ and $s$ both range from $1$ to $T - m + 1$ in the summation, m is the embedding dimension, and $\epsilon$ is the threshold value. $I(\cdot)$ is called the indicator function and is defined as: \begin{aligned} I[y_{t}^{m}, y_{s}^{m}, \epsilon ] &= 1 ~~~~~\text{if}~ || y_{t}^{m} - y_{s}^{m} || \leq \epsilon \\ &= 0 ~~~~~\text{otherwise} \end{aligned} The term, $|| y_{t}^{m} - y_{s}^{m} ||$, is a metric for the distance between two pairs. An example using the aforementioned hypothetical data set would be that the distance between $y_{1}^{3} = [y_{1}, y_{2}, y_{3}]$ and $y_{3}^{3} = [y_{3}, y_{4}, y_{5}]$ is given by: \[y_{1}^{3} - y_{3}^{3} = [ y_{1} - y_{3}, y_{2} - y_{4}, y_{3} - y_{5} ] \] The next step of the test involves choosing the maximum element from this resulting vector and then comparing it to the threshold value, denoted, $\epsilon$. If the maximum value is greater than $\epsilon$ then the pair is counted. Otherwise, the pair is not counted. The counting continues as the values for $t$ and $s$ are varied. Then, the number of counted pairs is divided by the total number of possible pairs. This fraction of pairs that are within a distance of $\epsilon$ of each other is the computed value for the correlation integral. The extreme cases are that all pairs are chosen, $C_{m, T} (\epsilon) = 1$, and that no pairs are chosen, $C_{m, T} (\epsilon) = 0$.

The BDS statistic is: \[ BDS = \sqrt{T} \frac{C_{m,T}(\epsilon) - C_{1,T}^{m}(\epsilon)}{\sigma_{m}(\epsilon)}\] where the variance, $\sigma_{m}^{2}(\epsilon)$, has the following complicated definition: \[ \sigma_{m}^{2}(\epsilon) = 4\{K(\epsilon)^{m} + 2 [\sum_{i=1}^{m-1} K(\epsilon)^{m-i} C(\epsilon)^{2i} + (m-1)^{2} C(\epsilon)^{2m} - m^{2} K (\epsilon) C(\epsilon)^{2m-2}]\} \] The indicator functions $C(\epsilon) = E\{I[y_{i}, y_{j}, \epsilon ]\}$ and $K(\epsilon) = E\{I[y_{i}, y_{j}, \epsilon ] I[y_{j}, y_{k}, \epsilon ]\}$ are defined in a similar way. For a complete definition, see Brock et al. (1996).

The null hypothesis of the BDS test is $H_{0}: y_{t}$ is i.i.d. The test is a two tailed test, so we should reject the null hypothesis of i.i.d at the 5% significance level when $|BDS| > 1.96$.

Notice that the statistic is a function of two parameters; the embedding dimension, m, and the threshold value, $\epsilon$. Values for these two parameters must be chosen, so this raises the question of how to choose values for m and $\epsilon$.

Various rules of thumb have been used to choose a value or values for $\epsilon$ and a short list of references can be found be referring to the concluding section of Brock (1995). The default values using the bds.test() function in R are $\epsilon$ = (0.5, 1.0, 1.5, 2.0) and this seems like the norm based on what I've read about a couple of statistics packages. These values for $\epsilon$ are usually converted back into units of the original data, so the may appear differently in the R output.

Granger et al. (2010) points out that some authors have reported the value of the test statistic for a number of $\epsilon$ values and that this is not particularly helpful. My reading of the literature tells me that it is, indeed, the norm to report results using a number of values for $\epsilon$.

To test for non-linearity in practice, we must first build a linear model for a time-series $y_{t}$. For instance, we could follow the Box-Jenkins methodology and construct an ARIMA model. The BDS test is then applied to the residuals of this model. If the null hypothesis is rejected, one can conclude that there is non-linearity in the data, but the form of non-linearity cannot be determined by the test. Therefore, if non-linearity is present, we can't say whether it is chaos or a nonlinear stochastic process.

A second method of testing for non-linearity using the BDS test involves taking differences of the logs of the time-series and then applying the BDS test to this transformed series. For an example of this method, see Chu, P.K.K. (2001).

I've decided that this post is long enough, so I will pick up the application in R in the next post. Before wrapping up, though, a couple of final comments about the BDS test.

Granger et al. (2010) lend the following advice: at least in small samples up to 200 observations, the BDS test may not be the ideal choice. Interestingly, these course notes by Eric Zivot say that 500 observations is usually thought to be the minimal sample size for the BDS test to have reliable performance. Finally, one advantage of the BDS test is that it does not require any distributional assumption on the data being tested.

(Unlinked) Reference:
Brock, W., & Dechert, W. D. J. Scheinkman, 1987,“A Test for Independence Based on the Correlation Dimension.”. Economics Working Paper.

Exponentially Weighted Moving Average models

Open up any finance textbook and you won't have to look far before you come across exponentially weighted average (EWMA) models. The EWMA is typically used to measure historical volatility, but it can, of course, be used in other contexts. For example, see Muth (1960) for an early discussion of the statistical properties of EWMA models and its use in macroeconomics, in particular, in relation to the adaptive expectations mechanism associated with Milton Friedman.

Based on the amount of online video tutorials and other online content providing explanations of the EWMA, I can only infer that the EWMA is a popular tool for those working in the finance industry. And why not, eh? It's a simple model that can be understood and interpreted by mangers and the like.

Indeed, for some time-series, the EWMA can provide an adequate representation of the data generating mechanism and produce optimal forecasts, however, the EWMA can also produce sub-optimal forecasts when (mis)applied to other series of data. For this reason, using the EWMA simply out of habit may be a bad habit worth kicking.

Before using the EWMA for forecasting, one ought to grant some time to considering whether or not the EWMA is appropriate for the data in question. That is, of course, under the assumption that time is not of the essence and that the practitioners' ultimate goal is to produce optimal forecasts. In other words, let's assume that the practitioner is not under a severe time constraint and faced with the task of updating or developing hundreds of forecasting models for a company.

The claim thus far is that strict use of the EWMA is not the highest standard of forecasting. A justification of this claim now follows.

Consider the following fact:

• The exponentially weighted moving average is algebraically equivalent to an ARIMA(0,1,1) model. 

To put it another way, the EWMA is a particular model within the class of ARIMA models. In fact, there are various types of EWMA models and these happen to be included in the class of ARIMA(0,d,q) models - see Cogger (1974).

Now, recall the fact that properly constructed ARIMA models produce optimal univariate forecasts - meaning that the mean-squared forecast error is smaller than for any other univariate forecast (produced by univariate, linear, fixed coefficient models).

Suppose now that rather than using the EWMA by habit, we instead use the Box-Jenkins method and let the data talk to us. In this situation, there are two possible outcomes; the appropriate model may, or may not, be an ARIMA(0,1,1), or, equivalently, an EWMA.

If it turns out that the appropriate model is an ARIMA(0,1,1) then our forecasts will be optimal. On the other hand, if the appropriate model is not an ARIMA(0,1,1) then our model has not been properly constructed and thus it must be admitted that our forecasts could be improved by specifying a different model from the suite of ARIMA models.

R Time Series Issue: Seasonal ARIMA Models

In my previous post, I mentioned that there are seasonal non-seasonal (S-NS) multiplicative models and seasonal non-seasonal additive models.

About these two types of models, Pankratz (1983) tells us that:

Multiplicative models are an excellent starting place for modeling S-NS realizations, and they often give good results. Most computer programs for estimating UBJ-ARIMA models provide at least an option for multiplicative seasonality; some simply assume that seasonal models are multiplicative.

It is then said that:

Occasionally, the assumption that the seasonal and nonseasonal elements are multiplicative is inappropriate.

This led me to question whether or not S-NS models in R using the arima() function are multiplicative by assumption or if there is an option to choose between multiplicative and additive models.

To fit a S-NS model in R, we can specify the function using code in the following form:

arima(x, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = s)

Other options can be specified, but there doesn't appear to be any reference to say whether or not the models are multiplicative by assumption.

I'm not the first to raise this question about S-NS models in R. The question was asked on the R mailing list here, but no answer was given. A follow up question was asked and in the email it is revealed that the arima() function assumes the models to be multiplicative.

The question now is: how to model S-NS additive models in R?

For the moment, one way around the problem of estimating S-NS additive models would be to use the lm() function. A proper solution would be to write code to slightly revamp arima(). I'm not sure if anyone is working on this problem, so I might try and write my own function that will do the trick.

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.