Exercise 1:

A committee of researchers is examining the ranking of a set of postgraduate dissertations submitted to a dissertation writing competition. The committee is concerned that Randolph, one of the committee members, tends to be strongly biased in favor of comparative grading. In short, they’re concerned that the grade that Randolph assigns to the current paper is pretty strongly dependent on the grade he assigned to the last paper. Below is the sequence of scores that Randolph has assigned to the dissertations in the prize competition. You’ll be responsible for determining whether the concerns about Randolph are founded.

scores = c(16, 14, 19, 17, 15, 20, 16, 14, 17, 
           15, 20, 17, 18, 19, 15, 17, 16, 18, 
           15, 17, 18, 14, 16, 15, 17, 14, 20)

1. Create a plot of scores. Describe how the score of the papers changes as Randolph grades.²
2. Draw a set of random scores with close to the same average score. Do these exhibit a similar structure? ³

Below, we can identify the number of times Randolph scores a dissertation above the average and then, in the next score, gives it below the average.

low.to.hi = ((scores[1:26] < mean(scores)) * (scores[2:27] > mean(scores)))
hi.to.low = ((scores[1:26] > mean(scores)) * (scores[2:27] < mean(scores)))
crosses.mean = low.to.hi + hi.to.low
total.crossings = sum(crosses.mean)

Thus, Randolph’s series crosses the mean score 17 times in 26 chances to do so.

3. In the same way you made the mock dataset in step 2, generate 1000 mock datasets. Make a histogram/kernel density plot of these crossings. How many times does a random sequence of scores cross the mean more than Randolph’s sequence of scores? ⁴
4. Given your simulations, do you think there is any evidence that Randolph is an unfair grader?

(For more, see Sargent & Stachurski on the basics of Markov chains.)

Recursion & Temporal Autoregressive models

So, let’s say you have this model of an autoregressive process at order \(p\). If we were to express the dependence over time as a graph, we’d be able to do so as follows: \[Y_t \leftarrow Y_{t-1} \leftarrow Y_{t-2} \leftarrow \dots \leftarrow Y_{t-p}\] Thus, when it comes to our model statement, let’s introduce an additional term, \(c\) which models the drift of the series in mean as the time develops. This yields: \[ Y_t =\phi_1 Y_{t-1} + \phi_2 Y_{t-2} + \dots + \phi_p Y_{t-p} + c + \epsilon_t\] But, by our own assumptions, this model applies further back in time, too:

To make the algebra easier, let’s just assume \(p=1\). Then, a full AR(1) system is:

Then, let’s substitute forward. For the first substitution, replacing \(Y_1\) with the equation generating \(Y_1\): \[ Y_2 = \phi \left(\phi Y_{0} + c + \epsilon_1\right) + c + \epsilon_2 = \phi^2 Y_0 + \phi \epsilon_1 + \phi c + c + \epsilon_2 \] which then gets inserted into the equation predicting \(Y_3\): \[ Y_3 = \phi \left( \phi \left(\phi Y_{0} + c + \epsilon_1\right) + c + \epsilon_2 \right) + c + \epsilon_3 = \phi^3Y_0 + \phi^2 c + \phi^2 \epsilon_1 + \phi c + \phi \epsilon_2 + c + \epsilon_3 \] And up on the chain until we reach the current time period: \[ Y_t = \phi \left( \phi \left( \dots \left( \phi \left( \phi Y_{0} + \epsilon_1\right) + \epsilon_2 \right) \dots \right) + \epsilon_t\right) = \phi^t Y_0 + \phi^{t-1}\epsilon_1 + \phi^{t-2} \epsilon_2 + \dots + \epsilon_t = \phi^{t}Y_0 + \sum_{i=0}^{t} \phi^{i}(c + \epsilon_{t-i})\] In the limit, if the signal is infinite, though, we obtain the following representation: \[ Y_t = \frac{c}{1 - \phi} + \sum_{i=0}^{\infty} \phi^{i}\epsilon_{t-i}\] As you can see, this representation involves the infinite sum of errors (sometimes called innovations in time series literature) that accumulate over time, attenuated by \(\phi\) each time period. This makes a few properties clear.⁷

Exercise 2.2

Return to the specifications we worked with in the previous section. Looking only at the ones that were AR(1):

• What happens to the sum of the last term if \(\phi\) is very small/close to zero?
• What happens to the sum of the last term if \(\phi\) is larger than one? what about when \(\phi\) is more negative than -1?

In case the results in Exercise 2.2 are in any doubt, when \(\phi\) is larger than 1 (or smaller than \(-1\)) in an AR(1) model, the error term becomes cumulative. Thinking intuitively, this means that the “total” error in the process (the sum of all “epsilon” terms) becomes larger and larger as time goes on. If \(|\phi| < 1\), the errors after a certain point in the past go to zero.⁸

Estimating AR models & Autocorrelation/covariance

While some elegant theoretical results are available in terms of the Yule Walker equations describing the autocovariance of the process (a good discussion is here) this method has largely been superseded with the estimation through maximum likelihood and information criteria. Essentially, these criteria balance the verbosity of the chosen model (in terms of the number of periods in the past to use for \(p\)) and the fit of the model. Once the order is identified, estimation via maximum likelihood is straightforward.

First, note that the model for all of \(Y_t, Y_{t-1},\dots,Y_{0}\) can be decomposed using the markov property. For an AR(1), this looks like: \[ f(Y_t, Y_{t-1},\dots,Y_{0}) = f(Y_0) * f(Y_1|Y_0) * f(Y_2 | Y_1) * \dots *f(Y_t | Y_{t-1})\] Since there is a “root” of this dependence chain at \(Y_0\), estimation can simply proceed through the sum like before. Assuming that innovations are normally-distributed with a time-constant variance, we can obtain the following likelihood:

This likelihood can be maximized using the conventional approaches.

Thinking about the autocovariance functions, there are two relevant ideas that are also important to understand. The autocovariance function (usually denoted \(\gamma(k)\)) is the covariance between \(Y_t\) and \(Y_{t-k}\) for some arbitrary value, \(k\). The autocovariance function standardizes this by the inherent variance of \(Y_0\), as is typical for correlation measures, to bring the domain between \((-1,1)\). Further, a partial autocovariance/autocorrelation provides the covariance between \(Y\) and \(Y_{t-k}\) for an arbitrary lag \(k\) while removing the contribution of all lags. Thus, if an autoregressive process is of order \(p\), its partial autocorrelation at lag \(p\) should be approximately zero, and this should be the first lag at which this becomes zero.

Exercise 2.3: fitting an AR

Some basic time-series functionality is built directly into R’s default stats module. For one of the series you simulated above, estimate an AR(p) model using the ar function in R:

model = ar(outcome)

Note that, by default, this will fit an autoregressive process where \(p\) is estimated using the Akaike information criterion. Since you know the “right” order of each of the outcomes we worked with in the previous qustions, fit an AR to them and see if the parameters are close to correct given their standard errors and point estimates. ⁹

Further, since you do know the “right” order, force the ar function to fit an autoregressive model at the “wrong order.” For instance, if I wanted to force the AR function to fit a model with order k, I’d have to specify that I do not want to fit using AIC maximization for the order and specify that order k should be fit:

k=3
model.bad.specification = ar(outcome, aic=FALSE, order.max=k)

• Fit an AR on the AR(5) where you force the model estimated to be an AR(20). Then, read about what the model’s $aic attribute means in the docstring of the ar function by running ?ar. Plot this attribute. Do you see a tradeoff between model fit and model complexity?

• There is an extensive set of tools for examining the autoregressive structure of these kinds of sequences. Try plotting the first 20 periods of the autocorrelation function and partial autocorrelation function for the stable AR(5) series discussed before using the acf and pacf functions in R.
  • Does the partial autocorrelation function go to zero at the 5th lag?
  • Do the signs of the autocorrelation function match the signs of the terms in \(\phi\)?

(There are many more models in this family. If you’re very interested in these kinds of econometric autoregressive models, a strong overview of the concepts + applications in R is provided by the author of the forecast package in R, Robert Hyndman, for free. This would be an excellent resource to learn to fit these models directly in R. Also, see Sargent & Stachurski on univariate and multivariate linear state space models, a generalization of these processes)

The spatial lag polynomial

Indeed, the “final” structure of the spatial autoregressive lag model contains higher-order contributions regardless of the assumed order or structure of \(\mathbf{W}\). To show this, though, we must reduce this model. Recall how we changed the AR(p) statement from being about \(Y_{t-k}\) to being a potentially-inifinite accumulation of errors.

A similar structure can be recognized in these models, and must be done in order to estimate these models using maximum likelihood. Here, since \(Y\) is a term on both sides of the equation, we must attempt to factor it into one side alone in order to express the rest as a “generating function” for \(Y\). So, a simple application of matrix algebra yields a similar style of “accumulation of errors” form:

Practically speaking, this can be understood by recognizing that the inverted matrix component \((I - \rho \mathbf{W})^{-1}\), sometimes called a “spatial filter matrix,” is actually a “Leontief inverse” that arises from a similarly-structured model of the recursive, mutually-implicative dynamics governing how goods and services are produced in an economy. For our purposes, this means we can restate our expression arrived at above into the sum of an infinite series of accumulating errors, just like in the time series case:

\[ Y = (I - \rho \mathbf{W})^{-1} \epsilon = (I + \rho \mathbf{W} + \rho^2 \mathbf{W}^2 + \dots)\epsilon = \sum_{i=0}^\infty\rho^i\mathbf{W}^i\epsilon\] Thus, again, there is an equivalence in the statement of \(Y\) as accumulating error and \(Y\) as directly autoregressive. Both are equivalent spellings of the same mathematical sentence. Further, we might think that \(\rho\) again should again be bound between \(-1\) and \(1\) in order the sum to vanish as \(i\) increases. However, since \(\rho\) is mediated by \(\mathbf{W}^i\), it is necessary both that \(\mathbf{W}\) and \(\rho\) exhibit special properties in order for this to be stable. That is, the eigenspectrum of the connectivity graph provided by \(\mathbf{W}\) will govern the stable ranges of \(\rho\). This can get quite complex at higher orders, too, as discussed in Elhorst (2014), sections 2.3-2.5.

These models are usually fit using maximum likelihood methods, but can also be fit using Bayesian methods, or multi-stage least squares algorithms, depending on the assumed variance structure of \(\nu\).

Autoregressive error specifications

Many other specifications focus on specifying an autoregressive process for \(\epsilon\), the error term, rather than for the \(Y\) term directly. Broadly speaking, these spatial error models come in an incredibly variety and many exhibit distinct properties. I’ll run through a few of them below.

Exercise 3.1

The STAR model

These kinds of spatial- and temporal-autoregressive models can be combined. In these space-time autoregressive models, the direct autoregressive and temporal autoregressive models are combined. These were suggested first by Cliff & Ord in early work, extended by Pfeifer & Deutsch (1980), see a very rigorous treatment in RJ Bennett’s book and are also discussed in Cressie & Wikle (2013). While theoretically useful, they’re not often used in practice for a few reasons.

First, just from a practical perspective, it’s usually the case that either space or time is a significant factor, but usually their joint autoregressive effects are unnecessary.

Second, there are significant stability issues, since space and time reverberate. To notice this, consider a STAR(1,1) while assuming \(\mathbf{W}\) is constant over time. In practice, this is spelled like a typical spatial lag model but with a single term from the past introduced as also being relevant: \[Y_t = \rho \mathbf{W} Y_t + \phi Y_{t-1} + \epsilon_t\] We do the same factor-and-reduce strategy here focusing on \(Y_t\): \[Y_t = (I - \rho \mathbf{W})^{-1}(\phi Y_{t-1} + \epsilon_t)\] but, the same equation applies in the previous period: \[ Y_{t-1} = (I - \rho \mathbf{W})^{-1}(\phi Y_{t-2} + \epsilon_{t-1})\] And the same recursive substitution strategy makes clear that this becomes a rapidly-increasing tower of spatial filter applications. Letting \(\mathbf{F} = (I - \rho \mathbf{W})^{-1}\), this can be stated more concisely as: \[ Y_{t} = \mathbf{F}\left[\phi\left(\mathbf{F}\left[\phi \left(\mathbf{F}\left[\phi \dots + \epsilon_{t-2}\right]\right) + \epsilon_{t-1}\right]\right) + \epsilon_t \right] = \sum_{j=0}^p(I - \rho \mathbf{W})^{-(j+1)}\phi^j\epsilon_{t-j}\] where a negative matrix power indicates the matrix power of the inverse of that matrix. The stability of this statement now depends both on \(\rho\) and on \(\phi\), so stability conditions must be assessed by considering the whole lag polynomial including \(\phi\) and \(\mathbf{F}\). The Leontief-expanded version also still applies, which allows the solutions for any single spatial filter term to be expanded using the infinite series representation, raised to the appropriate power: \[ Y_{t} = \sum_{j=0}^p(I - \rho \mathbf{W})^{-(j+1)}\phi^j\epsilon_{t-j} = \sum_{j=0}^p\left(I + \sum_k^\infty\rho^k\mathbf{W}^k\right)^{j+1}\phi^j\epsilon_{t-j}\]

Further, let’s focus in on a single observation \(y_{it}\): \[ y_{it} = \rho \sum_{j \in \eta(i)} w_{ij}y_{j,t} + \phi y_{i,t-1} + \epsilon_{i,t}\] where now the summation involves all \(j\) observations in the neighborhood of site \(i\) (denoted by \(\eta(i)\)).

Let’s pick a specific one of these \(y_{j,t}\) who is in the neighbourhood of \(i\) and call it \(y_{k,t}\). Now, consider its generating equation: \[ y_{kt} = \rho \sum_{\ell \in \eta(k), \ell \neq i} w_{k\ell}y_{\ell,t} + \rho w_{ik} y_{it} + \phi y_{k,t-1} + \epsilon_{k,t}\] substitute: \[ y_{kt} = \rho \sum_{\ell \in \eta(k), \ell \neq i} w_{k\ell}y_{\ell,t} + \rho^2 w_{ik}\sum_{j \in \eta(i)} w_{ij}y_{j,t} + \rho\phi w_{ik} y_{i,t-1} + \phi y_{k,t-1} + \rho w_{ik} \epsilon_{i,t} + \epsilon_{k,t}\] Looking closely, you’ll note that two \(t-1\) terms are now present; \(\rho \phi w_{ik} y_{i,t-1}\) and \(\phi y_{k,t-1}\). This means that despite the fact that \(y_{kt}\) only includes direct autoregressive structures for of its current neighbor’s values and its own past value, these combine to make its neighbours past values relevant as well. Thus, there is “bleed” between the spatial and temporal autoregressive structures. For the model, the time- and space-autoregressive parameter \(\rho\) and \(\phi\) will interact, complicating the determination of process stability.

Usually, these models are fit with an iterative backfitting process discussed by Pfeffer & Deutsch (1980), and recent work on the problem tends to focus on Kalman filtering¹³.