Chapter 06: Dependent Data - Time Series and Spatial Statistics

This note is a high-fidelity Markdown migration of the Dependent Data: Time Series and Spatial Statistics chapter from the LaTeX source.

Parent map: index Prerequisites: probability-and-mathstats, linear-regression, maximum-likelihood-and-machine-learning

Concept map

flowchart TD
  A[Dependent Data] --> B[Time Series]
  B --> C[Stationarity]
  C --> D[Ergodicity]
  B --> E[AR MA ARMA]
  B --> F[Unit Root]
  F --> G[Cointegration]
  B --> H[HAC Inference]
  A --> I[Spatial Statistics]
  I --> J[Kriging]
  I --> K[Spatial Autocorrelation]
  K --> L[Variogram]
  I --> M[Spatial Regression]
  I --> N[GMRF GP CAR]

Dependent Data: Time series and spatial statistics

Time Series

A time series is a sequence of data points ${\left\{w_t\right\}}_{t=1}^T$ observed over time. In a random sample, points are iid, so the joint distribution $f(w_1,\dots ,w_T)={\prod }_{t=1}^{T}f(w_t)$. In time series, this is clearly violated, since observations that are temporally close to each other tend to be more similar.

A stochastic process is a sequence of random variables $\left\{\dots ,Y_{-1},Y_0,Y_1,\dots \right\}$ indexed by elements in a set of indices $\left\{Y_t:t\in \mathcal{T}\right\}$. Hypothetical repeated realisations of a stochastic process look like

$${\left\{w_t^{(1)},w_t^{(2)},\dots ,w_t^{(n)}\right\}}_{t=-\infty }^{\infty }$$

The index set $\mathcal{T}$ may be either countable, in which case we get a discrete time process or uncountable, in which case we get a continuous time process.

State Space

We assume $\exists$ a set $\mathcal{Y}\in \mathbb{R}\text{ s.t. }\forall t\in \mathcal{T},Y_t\in \mathcal{Y}$. Then, $\mathcal{Y}$ is called the State Space of the stochastic process.

Consider a random process ${\left\{Y_t\right\}}_{t=1}^{\infty }$ and an increasing sequence of information sets ${\left\{{\mathcal{F}}_t\right\}}_{t=1}^{\infty }$ i.e. collection of $\sigma -$fields s.t. ${\mathcal{F}}_0\subset {\mathcal{F}}_1\dots {\mathcal{F}}_{\infty }\subset \mathcal{F}$. If $Y_t$ belongs to the information set ${\mathcal{F}}_t$ and is absolutely integrable [i.e. $y_t\in L_0({\mathcal{F}}_t)\cap L_1(\mathcal{F})$], and $\mathbb{E}\left[Y_{t+1}|{\mathcal{F}}_t\right]=Y_t\forall t<\infty$ then ${\left\{Y_t\right\}}_{t=0}^{\infty }$ is called a martingale. In words, the conditional expected value of the next observation, given all the past observations, is equal to the most recent observation.

The autocovariance of $Y_t$ is the covariance between $Y_t$ and its $j^{th}$ lagged value

$${\gamma }_{jt}:=\mathbb{E}\left[Y_t-{\mu }_t\right][Y_{t-j}-{\mu }_{t-j}]$$

the variance-covariance matrix of $y=\left\{y_t\right\}$ has Toeplitz form:

$$\mathbb{V}\left[y\right]={\left[{\gamma }_{|i-j|}\right]}_{i,j=1}^T.$$

the $j^{th}$ order correlation coefficient ${\rho }_j:={\gamma }_j/{\gamma }_0$.

A random process is said to be stationary if the distribution functions of $(X_{t_1},X_{t_2}\dots )$ and $(X_{t_1+j},X_{t_2+j}\dots )$ are the same $\forall t_1,\dots ,t_k,h\in \mathbb{Z}$.

A process is said to be covariance (or weakly) stationary if

$\mathbb{E}\left[Y_t\right]=\mu \forall t\in T$

${\gamma }_{jt}=\mathbb{E}\left[Y_t-{\mu }_t\right][Y_{t-j}-{\mu }_{t-j}]={\gamma }_j\forall t\in T$

i.e. neither the mean nor the autocovariances depend on the date $t$; stationary expectation, variance, and covariance. Most relevant variables aren’t stationary, but their detrended or first-differenced versions may be.

If $X_0,X_1,\dots$ is a Markov Process,

$$\mathbf{Pr}\left(X_{n+1}\le x|X_1,\dots ,X_n\right)=\mathbf{Pr}\left(X_{n+1}\le x|X_n\right)$$

that is, the conditional distribution of $X_{n+1}$ given $X_1,\dots ,X_n$ does not depend on $X_1,\dots ,X_{n-1}$.

Markov Chain

A Markov chain is simply a Markov process in which the state-space is a countable set. Since a Markov chain is a markov process, the conditional distribution of $X_{t+1}|X_1,\dots ,X_t$ depends only on $X_t$. The conditional distribution is often represented by a Transition matrix where

$${\mathbf{P}}_{ij}^{t,t+1}=\mathbf{Pr}\left(X_{t+1}=j|X_t=i\right);i,j=1,\dots J$$

If $\mathbf{P}$ is the same $\forall t$, we say the Markov chain has stationary transition probabilities.

A stationary process is ergodic if any two variables positioned far apart in the sequence are almost independently distributed.

$\left\{x_t\right\}$ is ergodic if, for any two bounded functions $f(.)$ in $k+1$ variables and $g(.)$ in $l+1$ variables,

$$\begin{array}{rl}\underset{N\to \infty }{\lim }& \left|\mathbb{E}\left[f(x_t,\dots ,x_{t+k})g(x_{t+N},\dots ,x_{t+l+N})\right]\right|-\\ & \left|\mathbb{E}\left[f(x_t,\dots ,x_{t+k})\right]\right|\left|\mathbb{E}\left[g(x_{t+N},\dots ,x_{t+l+N})\right]\right|=0\end{array}$$

i.e. ${\lim }_{j\to \infty }{\gamma }_j=0$

Sufficient condition for ergodicity is $x_t$ be covariance stationary and ${\sum }_{j=0}^{\infty }\left|{\gamma }_j\right|<\infty$

Ergodic processes have the following property

$$\mathbb{V}\left[\sum _{t=1}^T{x}_t\right]=\sum _{j=1-T}^{T-1}(T-\left|j\right|){\gamma }_j$$

this result implies that

$$\underset{T\to \infty }{\lim }\mathbb{V}\left[\frac{1}{\sqrt{T}}\sum _{t=1}^T{x}_t\right]=\sum _{j=-\infty }^{\infty }{\gamma }_j<\infty$$

This permits us to swap $i$s for $t$s and derive Asymptotic theory with dependent observations, such as LLN and CLT.

A family of r.v.s $\left\{X_t\right\}$ indexed by a continuous variable $t$ over $[0,\infty )$ is a Brownian Motion iff

$X(0)=0$

$\left\{X(s_i+t_i)-X_{(}{S}_i)\right\}$ over an arbitrary collection of disjoint intervals $(s_i,s_i+t_i$ are independent r.v.s

$\forall s,t\ge 0,X(s+t)-X(s)\sim \mathcal{N}\left(0,t\right)$

White noise is a sequence $\{{\epsilon }_t\}$ whose elements have mean zero and variance ${\sigma }^2$, and for which ${\epsilon }_t$‘s are uncorrelated over time

$\mathbb{E}\left[{\epsilon }_t\right]=0$

$\mathbb{E}\left[{\epsilon }_t^2\right]={\sigma }^2$

$\mathbb{E}\left[{\epsilon }_t{\epsilon }_{t-j}\right]=0\forall j\ne 0$

A moving average of order $q$, $MA(q)$ is a weighted average of the $q$ most recent values of a white noise defined as

$$Y_t=\mu +{\epsilon }_t+{\theta }_1{\epsilon }_{t-1}+\cdots +{\theta }_q+{\epsilon }_{t-q}$$

An autoregressive process of order $p$, $AR(p)$ is given by $Y_t$ as a linear combination of $p$ lags of itself and one white noise

$$Y_t=\mu +{\varphi }_1{Y}_{t-1}+\cdots +{\varphi }_p{Y}_{t-p}+{\epsilon }_t$$

ARMA(p, q) combines AR(p) and MA(q)

$$Y_t=\mu +\underset{\text{MA(q)}}{\underbrace{{\epsilon }_t+{\epsilon }_t+{\theta }_1{\epsilon }_{t-1}+\cdots +{\theta }_q+{\epsilon }_{t-q}}}+\underset{\text{AR(p)}}{\underbrace{{\varphi }_1{Y}_{t-1}+\cdots +{\varphi }_p{Y}_{t-p}}}$$

Consider AR(1): $Y_t=\rho Y_{t-1}+{\epsilon }_t$. Since this holds at $t$, it holds at $t-1⟹Y_{t-1}=\rho Y_{t-2}+{\epsilon }_{t-1}$. Substitute into original to get $Y_t=\rho (\rho Y_{t-2}+{\epsilon }_{t-1})+{\epsilon }_t$. Repeat ad infinitum to obtain, as long as $\rho <1$

$$Y_t=\sum _{s=0}^{\infty }{\rho }^s{\epsilon }_{t-s}$$

In other words, AR(1) $\equiv$ MA($\infty$) ; they are different representations of the same underlying stochastic process.

Wold Representation: All covariance-stationary time series processes can be represented by / decomposed into a deterministic component and a $MA(\infty )$

In a stationary process, $\mathbb{E}\left[x_t\right]=\mu$, which is seldom true. A less restrictive assumption that allows for nonstationarity is to specify the mean as a function of time.

$$\mathbb{E}\left[x_t\right]=\alpha +\beta t\text{specify}{x}_t=\alpha +\beta t+{\epsilon }_t;\epsilon \text{stationary}$$

A random walk ($\equiv I(1)$) is a process such that $\mathbb{E}\left[x_t|x_{t-1},x_{t-2},\dots \right]=x_{t-1}$.

$x_t=x_{t-1}+{\epsilon }_t,{\epsilon }_t\sim \mathcal{N}\left(0,{\sigma }^2\right)$ = AR(1) process with $\varphi =1=:$ Unit Root. Rewrite as

$$x_t=x_{t-1}+{\epsilon }_t=x_0+\sum _j^t{\epsilon }_j$$

Random walk with drift

$$x_t=x_{t-1}+\delta +{\epsilon }_t=\delta t+{\epsilon }_t+{\epsilon }_{t-1}\dots {\epsilon }_1+x_0=x_0+\delta t+\sum _j^t{\epsilon }_j$$

For the following model $x_t=\rho x_{t-1}+{\epsilon }_t,{\epsilon }_t\sim \mathcal{N}\left(0,{\sigma }^2\right)$ = AR(1)

test $\rho =1$. Distribution of $\hat{\rho }$ under the null $\rho =1$ is non-standard: CLT not valid. test to use: Dickey Fuller, Augmented Dickey Fuller, Phillips-Perron.

Let $y_t\sim I(1)\wedge x_t\sim I(1)$. $y_t$ and $x_t$ are said to be cointegrated if $\exists \psi \text{ s.t. }{y}_t-\psi x_t\sim I(0)$. For example, let

$$\begin{array}{rl}{y}_t& =\beta x_t+{\epsilon }_t\\ x_t& =x_{t-1}+{\upsilon }_t\end{array}$$

where $({\epsilon }_t,{\upsilon }_t)$ is white noise. Then, $y_t,x_t\sim I(1)$, but $y_t-\beta x_t\sim I(0)$, with cointegration vector $a=(1,-\beta )$.

decomposes an observed time series $X_t,t=1,2,\dots ,n$ into a trend ${\overline{X}}_t$ and a stationary component ${\tilde{X}}_t=X_t-{\overline{X}}_t$ so that the trend ${\left\{{\overline{X}}_t\right\}}_{t=1}^n$ minimises

$$\frac{1}{n}\sum _{t=1}^n(X_t-{\overline{X}}_t{)}^2+\underset{\text{Penalty for incorporating fluctuations}}{\underbrace{w\frac{1}{n}\sum _{t=2}^{n-1}{\left(({\overline{X}}_{t+1}-{\overline{X}}_t)-({\overline{X}}_t-{\overline{X}}_{t-1})\right)}^2}}$$

$w$ is a tuning parameter. In quarterly data, $w=1600$.

Regression with time series

Basic assumption in conventional OLS with time series is $\mathbb{E}\left[y_t|x_1,\dots ,x_T\right]=\mathbb{E}\left[y_t|x_t\right]=x_t^{\prime }\beta$. Equivalently, $y_t=x_t^{\prime }\beta +{\epsilon }_t\mathbb{E}\left[{\epsilon }_t|X\right]=0$ where $X=(x_1,\dots ,x_T{)}^{\prime }$. The second classical assumption is $\mathbb{E}\left[{\epsilon }_t^2|x\right]={\sigma }^2\forall t;\mathbb{E}\left[{\epsilon }_t{\epsilon }_{t-j}\right]=0\forall t,j$.

$\mathbb{E}\left[u_t,u_{t-j}|X\right]\ne 0$ is called autocorrelation. Fix: Newey-West HAC consistent variance estimator ‘meat’

$$\hat{V}={\hat{\Gamma }}_0+\sum _{i=1}^m\left(1-\frac{j}{m+1}\right)({\hat{\Gamma }}_j+{\hat{\Gamma }}_j^{\prime })\text{where }{\hat{\Gamma }}_j=\frac{1}{T-j}\sum _{t=j+1}^T{\hat{\epsilon }}_j{\hat{\epsilon }}_{t-j}{x}_t{x}_{t-j}^{\prime }$$

with variance estimated the normal way

$$\hat{W}={\left(\frac{1}{T}\sum _{t=1}^T{x}_t{x}_t^{\prime }\right)}^{-1}\hat{V}{\left(\frac{1}{T}\sum _{t=1}^T{x}_t{x}_t^{\prime }\right)}^{-1}$$

Consider the model

$$y_t=\delta +\alpha y_{t-1}{\beta }_0{x}_t+{\beta }_1{x}_{t-1}+{\epsilon }_t$$

Subtract $y_{t-1}$ and add $-{\beta }_0{x}_{t-1}+{\beta }_0{x}_{t-1}$ to the l.h.s. we get

$$y_t-y_{t-1}=\delta -(1-\alpha )y_{t-1}+{\beta }_0(x_t-x_{t-1})+({\beta }_0+{\beta }_1)x_{t-1}+{\epsilon }_t$$

and

$$\Delta y_t=\delta +{\beta }_0\Delta x_t-(1-\alpha )(y_{t-1}-\gamma x_{t-1})+{\epsilon }_t$$

where $\gamma$ is the long run effect

$$\gamma =\frac{{\beta }_0+{\beta }_1}{1-\alpha }$$

A Quandt Likelihood ratio test begins with no knowledge of when the trend break occurs [although researchers typically know of the timing for substantive reasons], and sequentially estimates the following model

$$\Delta Y_t=\log Y_t-\log Y_{t-1}=\alpha +{\delta }_0{D}_t(\tau )+{\epsilon }_t$$

where $\Delta Y_t$ is the first difference of the outcome, and $D_t(\tau )$ is an indicator variable equal to zero for all years before $\tau$ and one for all subsequent years. The researcher varies $\tau$ and tests the null that ${\delta }_0=0$, and the largest F-statistic is used to determine the best possible break point. Use Andrews (2003) critical values to account for multiple-testing.

Spatial Statistics

A spatial stochastic process is a collection of random variables $y(u)$ indexed by location $u$: $\left\{y_i,i\in \mathcal{D}\subset {\mathbb{R}}^d\right\}$, where $\mathcal{D}$ is either a continuous surface of a finite set of discrete locations.

For each location $u$, $y(u)$ is a random variable, and thus needs to be modeled. Basic approach is to assume $\mathbb{E}\left[y(u)\right],\mathbb{V}\left[y(u)\right]$ exist, and decompose

$$y(u)=\underset{\text{mean function}}{\underbrace{m(u)}}+\underset{\text{error}}{\underbrace{e(u)}}$$

mean function $m(u)=\mathbb{E}\left[y(u)\right]$ and stochastic error process $e(u)\text{ s.t. }\mathbb{E}\left[e(u)\right]=0$.

Kriging - modeling $m(u)$

Main reference: Christensen (2019, ch. 8).

Assume linear structure for $m(u)$. $p$ known functions of $u$, $x_1(u),\dots ,x_p(u)$ s.t.

$$m(u)=\sum _{j=1}^p{\beta }_j{x}_j(u)$$

A special case of this is the Ordinary Kriging model where

$$m(u)=\mu$$

for unknown $\mu$. The most basic model is Simple Kriging where

$$m(u)={\mu }_0$$

with known ${\mu }_0$.

Assume the universal kriging model $m(u)={\sum }_j^p{\beta }_j{x}_j(u)$ holds, we have data on locations $u_1,\dots ,u_n$, and that we wish to predict the value of $y(u_0)$. The model can be written

$$\begin{array}{rl}\mathbf{Y}& =\mathbf{X}\mathbf{\beta }+\mathbf{e}\\ \mathbb{E}\left[e\right]& =0\\ \mathrm{Cov}\left(e\right)& =\mathbf{\Sigma }=[{\sigma }_{ij}]=\sigma (u_i,u_j)i,j=1,\dots ,n\end{array}$$

Let ${\Sigma }_{Y0}:=\left[\begin{array}{c}{\sigma }_{10}\\ ⋮\\ {\sigma }_{n0}\end{array}\right]$

The best linear unbiased predictor of $y_0$ is

$${\hat{y}}_0=x_0^{\prime }\hat{\beta }+{\delta }^{\prime }(Y-X\hat{\beta })$$

where $\hat{\beta }={\left(X^{\prime }{\Sigma }^{-1}X\right)}^{-1}{X}^{\prime }{\Sigma }^{-1}Y$ and $\delta ={\Sigma }^{-1}{\Sigma }_{Y0}$.

Spatial Autocorrelation: Modelling $e(u)$

Spatial Autocorrelation is expressed as

$$\begin{array}{rl}\sigma (u,w)& :=\mathrm{Cov}\left(y(u),y(w)\right)=\mathrm{Cov}\left(e(u),e(w)\right)=\sigma (w,u)\\ & =\mathbb{E}\left[y(u)y(w)\right]-\mathbb{E}\left[y(u)\right]\mathbb{E}\left[y(w)\right]\ne 0\forall i\ne j\end{array}$$

Covariance is often modelled in terms of an unknown parameter $\mathbf{\theta }$, in which case we write $\sigma (u,w;\mathbf{\theta })$. Assumptions made about $e(u)$ include:

• second-order stationarity,
• strict stationarity,
• intrinsic stationarity,
• increment stationarity,
• isotropy.

Covariance functions can be modelled in three basic ways:

• specify a functional form for the stochastic process generating $\left\{y_i,i\in \mathcal{D}\right\}$, and derive covariance from that process,
• model covariance directly as a function of a small number of parameters,
• leave covariance unspecified and estimate nonparametrically.

A process $y(u)$ is strictly stationary if for all $k$, locations $u_1,\dots ,u_k$, Borel sets $C_1,\dots ,C_k$, and shifts $h\in {\mathbb{R}}^d$,

$$\mathbf{Pr}\left(y(u_1)\in C_1,\dots ,y(u_k)\in C_k\right)=\mathbf{Pr}\left(y(u_1+h)\in C_1,\dots ,y(u_k+h)\in C_k\right).$$

This implies translation invariance of the joint law. In particular:

$$m(u)=\mu ,\sigma (u,w)=\sigma (u-w)=\sigma (h).$$

If, in addition, the finite-dimensional distributions are multivariate Gaussian, the process is a Gaussian Process.

Second-order (weak) stationarity imposes the same constant mean and covariance depending only on distance, but does not require full strict stationarity.

Increment-stationarity requires invariant increment laws:

$$\mathbf{Pr}\left(y(u_2)-y(u_1)\in C_1,\dots ,y(u_k)-y(u_{k-1})\in C_k\right)=\mathbf{Pr}\left(y(u_2+h)-y(u_1+h)\in C_1,\dots ,y(u_k+h)-y(u_{k-1}+h)\in C_k\right).$$

Brownian motion is increment-stationary but not strictly stationary.

For increment-stationary processes, the semivariogram is

$$\gamma (u,w)=\frac{1}{2}\mathbb{E}\left[(y(u)-y(w){)}^2\right]=\frac{1}{2}\mathbb{V}\left[y(u)-y(w)\right]=\frac{1}{2}\left(\sigma (u,u)+\sigma (w,w)-2\sigma (w,u)\right).$$

The variogram is $2\gamma (u,w)$. Under increment-stationarity:

$$\gamma (u,w)=\gamma (u+h,w+h)=\gamma (u-w,0)=\gamma (u-w).$$

An intrinsically-stationary process satisfies the constant-mean restriction and this semivariogram invariance condition. All second-order stationary processes are intrinsically stationary, but not vice versa.

For a linear model, stipulate a nonnegative definate weighting matrix, and fit

$$Y=X\beta +e,\mathbb{E}\left[e\right]=0,\mathrm{Cov}\left(e\right)={\Sigma }_0.$$

to obtain residuals ${\hat{e}}_0=Y-X\hat{\beta }$. For any vector $h$, there is a finite number $N_h$ of pairs of observations $y_i,y_j$ for which $u_i-u_j=h$. For each of these pairs, list the corresponding residual pairs, $({\hat{e}}_{0i},{\hat{e}}_{0i(h)}),i=1,\dots ,N_h$. If $N_h\ge 1$, the traditional empirical covariance estimator is

$$\hat{\sigma }(h)=\hat{\sigma }(-h)=\frac{1}{N_h}\sum _{i=1}^{N_h}{\hat{e}}_i{\hat{e}}_{i(h)}$$

The traditional empirical semivariogram estimator in ordinary kriging (no covariates) is

$$\hat{\gamma }(h)=\frac{1}{2N_h}\sum _{i=1}^{N_h}(y_i-y_{i(h)}{)}^2$$

A second-order stationary process is said to be isotropic if

$$\sigma (u-w)=\sigma (\left|\left|u-w\right|\right|)$$

An intrinsically stationary process is isotropic if

$$\gamma (u-w)=\gamma (\left|\left|u-w\right|\right|)$$ $$\mathbf{y}-\mu \mathbf{\iota }=\rho \mathbf{W}(\mathbf{y}-\mu \mathbf{\iota })+\mathbf{\epsilon }={\left(\mathbf{I}-\rho \mathbf{W}\right)}^{-1}+\mathbf{\epsilon }$$

where $\mathbf{W}$ is a $N\times N$ weight matrix. a spatial lag for $y_i$

$$\mathbf{W}{y}_i=\sum _j{w}_{ij}{y}_j$$

A parsimonious specification of a small number of parameters for the covariance matrix is typically presumed.

$$\mathrm{Cov}\left({\epsilon }_i,{\epsilon }_j\right)={\sigma }^{2}f(d_{ij},\phi )$$

where ${\epsilon }_i,{\epsilon }_j$ are residuals, ${\sigma }^2$ is the error variance, $d_{ij}$ is the distance between $i,j$, and $f$ is a distance decay function such that $\frac{\partial f}{\partial d}<0$ and $\left|f(d_{ij},\phi )\right|\le 1$, with $\phi \in \Phi$ being a $p\times 1$ vector.

The generalised Moran’s I is a weighted, scaled cross-product

$$\mathcal{I}:=\frac{n{\sum }_{i=1}^n{\sum }_{j\ne i}{w}_{ij}(y_i-\overline{y})(y_j-\overline{y})}{{\sum }_{i=1}^n{\sum }_{j\ne i}{w}_{ij}{\sum }_i(y_i-\overline{y}{)}^2}$$

Its expected value is $\frac{-1}{n-1}$.

A test for Moran’s I involves shuffling the locations of points and computing $\mathcal{I}$ $S$ times. This produces a randomization distribution under $H_0$.

A Monte-carlo P-value is

$$\hat{p}=\frac{1+{\sum }_{s=1}^S{1}_{I_s^{\ast }\ge I_{obs}}}{S+1}$$

Spatial Linear Regression

A simple spatial regression is

$$\mathbf{y}=\rho \mathbf{Wy}+\mathbf{X}\mathbf{\beta }+\mathbf{\epsilon }$$

the solution is

$$\hat{\mathbf{\beta }}={\left(\mathbf{X}{\mathbf{X}}^{\top }\right)}^{-1}\mathbf{X}{(}^{\top }\mathbf{I}-\rho \mathbf{W})\mathbf{y}$$

Its reduced form is

$$\mathbf{y}={\left(\mathbf{I}-\rho \mathbf{W}\right)}^{-1}\mathbf{X}\mathbf{\beta }+{\left(\mathbf{I}-\rho \mathbf{W}\right)}^{-1}\mathbf{\epsilon }$$

The spatial lag term induces correlation between the error and explanatory variables, and thus must be treated as an endogenous variable.

A spatial error model is simply an linear model with a non-spherical but typically parametric structure in the error covariance matrix.

$$\mathbf{y}=\mathbf{X}\mathbf{\beta }+\underset{\text{Composite error ε}}{\underbrace{\lambda \mathbf{W}\mathbf{\xi }+\mathbf{\eta }}},\eta \sim \mathcal{N}\left(0,{\sigma }^2\mathbf{I}\right)$$ $$\mathbb{E}\left[\epsilon {\epsilon }^{\prime }\right]=\mathbf{\Omega }(\mathbf{\theta })$$

A covariance function decomposes into a systematic part and idiosyncratic noise as follows

$${\mathbf{\Sigma }}_{ij}={\sigma }^{2}C(\Vert i-j\Vert ,\mathbf{\pi })+{\tau }^2{1}_{ij}\equiv {\sigma }^2{\mathbf{C}}_{UU}+{\tau }^2\mathbf{I}$$

where $C$ is a correlation function, $\Vert i-j\Vert$ is the distance between points $i,j$.

Kelly recommends using a Whittle-Matern function defined next. These parameters can be fitted on the error distribution to estimate the covariance matrix.

A covariance function $\mathrm{Cov}\left(Y(\mathbf{x}),Y\left({\mathbf{x}}^{\prime }\right)\right)$ describes the joint variability between a stochastic process $Y(\cdot )$ at two locations $\mathbf{x}$ and ${\mathbf{x}}^{\prime }$. This covariance function is vital in spatial prediction. The fields package includes common parametric covariance families (e.g. exponential and Matern) as well as nonparametric models (e.g. radial and tensor basis functions).

When modeling $\mathrm{Cov}\left(Y(\mathbf{x}),Y\left({\mathbf{x}}^{\prime }\right)\right),$ we are often forced to make simplifying assumptions.

Stationarity assumes we can represent the covariance function as

$$\mathrm{Cov}(Y(\mathbf{x}+\mathbf{h}),Y(\mathbf{x}))=C(\mathbf{h})$$

for some function $C:{\mathbb{R}}^d\to \mathbb{R}$ where $\dim (\mathbf{x})=d$.

• Isotropy assumes we can represent the covariance function as

$$\mathrm{Cov}(Y(\mathbf{x}+\mathbf{h}),Y(\mathbf{x}))=C(\Vert \mathbf{h}\Vert )$$

for some function $C:\mathbb{R}\to \mathbb{R},$ where $\Vert \cdot \Vert$ is a vector norm.

Exponential :

$$\mathrm{Cov}\left(Y(\mathbf{x}),Y\left({\mathbf{x}}^{\prime }\right)\right)=C(r)=\rho e^{-r/\theta }+{\sigma }^2{1}_{\mathbf{x}={\mathbf{x}}^{\prime }}$$

Matern:

$$\mathrm{Cov}\left(Y(\mathbf{x}),Y\left({\mathbf{x}}^{\prime }\right)\right)=C(r)=\rho \left(\frac{2^{1-\nu }}{\Gamma (\nu )}{\left(\frac{r}{\theta }\right)}^{\nu }{K}_{\nu }\left(\frac{r}{\theta }\right)\right)+{\sigma }^2{1}_{\mathbf{x}={\mathbf{x}}^{\prime }}$$

where $K_{\nu }$ is a modified Bessel function of the second kind, of order $\nu$

Matern covariance depends on $(\rho ,\theta ,\nu ,{\sigma }^2)$, while exponential depends on $\rho ,\theta ,{\sigma }^2)$, where

$\theta$: is the range of the process at which observations become uncorrelated

$\rho$ : marginal variance / ‘sil’

${\sigma }^2$ : small scale variation such as measurement error

$\nu$ : smoothness

$$\mathbf{y}=\mathbf{X}\mathbf{\gamma }+\stackrel{\text{AC in Y}}{\overbrace{\beta \mathbf{Wy}}}+\stackrel{\text{AC in X}}{\overbrace{\mathbf{WX}\mathbf{\theta }}}+\stackrel{\text{AC in errors}}{\overbrace{\mathbf{W}\mathbf{\nu }\mathbf{\lambda }}}+\mathbf{\epsilon }$$

Here, $\mathbf{W}$ is a weight matrix (typically row-standardised), so $\mathbf{WM}$ is a spatial lag. In spatial econometrics, the above form nests many popular regressions

Spatially Autoregressive (SAR) Model : $\mathbf{\lambda }=\mathbf{\theta }=0$

Spatially lagged $x$ : $\mathbf{\beta }=\mathbf{\lambda }=0$

Spatial Durbin Model : $\mathbf{\lambda }=0$

Spatial Error model : $\mathbf{\beta }=\mathbf{\lambda }=0$

In the Social Interactions literature (e.g., Manski 1993), the above expression is written in the form of conditional expectations

$$y_i={\mathbf{x}}_i^{\prime }\mathbf{\gamma }+\mathbb{E}{\left[y|{\mathbf{w}}_i\right]}^{\prime }\stackrel{\text{Endogenous}}{\overbrace{\beta }}+\mathbb{E}{\left[{\mathbf{x}}_i|{\mathbf{w}}_i\right]}^{\prime }\stackrel{\text{Contextual}}{\overbrace{\mathbf{\theta }}}+\mathbb{E}{\left[\nu |{\mathbf{w}}_i\right]}^{\prime }\underset{\text{Correlated}}{\underbrace{\lambda }}+{\epsilon }_i$$

in practice, the expectations are replaced with empirical counterparts $\hat{\mathbb{E}}(y|{\mathbf{w}}_i)=\mathbf{W}y$ and so on, so the estimation steps are isomorphic.

Define unobservables as $\mathbf{\upsilon }=\mathbf{W}\nu +\mathbf{\epsilon }$, and assume they are uncorrelated with observables $\mathbf{x}$; that is, there is no sorting and no omitted spatial variables. Then, we can write

$$\mathbf{y}=\mathbf{X}\gamma +\mathbf{Wy}\beta +\mathbf{WX}\mathbf{\theta }+\mathbf{\upsilon }$$

Premultiplying by $\mathbf{Wy}$ gives

$$\mathbf{Wy}=\mathbf{WX}\gamma +\mathbf{WWy}\beta +\mathbf{W}\mathbf{WX}\mathbf{\theta }+\mathbf{W}\mathbf{\upsilon }$$

This shows that $\mathbf{Gy}$ is correlated with $\mathbf{\upsilon }$, i.e. $\mathbb{E}\left[\mathbf{\upsilon }|\mathbf{Wy}\right]\ne 0$, and least square estimates of the above regression are biased.

If we assume $\mathbf{W}$ is idempotent (by constructing a block-diagonal, transitive matrix), we can simplify the above expression to

$$\begin{array}{rl}\mathbf{W}y& =\mathbf{WX}\frac{\mathbf{\gamma }+\mathbf{\theta }}{1-\mathbf{\beta }}+\mathbf{W}\upsilon /(1-\beta )\text{Plugging in definition for Wy}\\ y& =\mathbf{X}\underset{\tilde{\gamma }}{\underbrace{\gamma /(1-\beta )}}+\mathbf{WX}\underset{\tilde{\theta }}{\underbrace{(\gamma \beta +\theta )/(1-\beta )}}+\underset{\tilde{\upsilon }}{\underbrace{\upsilon +\mathbf{W}\upsilon \beta /(1-\beta )}}\end{array}$$

In summary, $\beta ,\theta$ cannot be separately identified from the composite parameters $\tilde{\beta },\tilde{\theta }$. This is the reflection problem discussed by Manski (1993).

Spatial Modelling

Based on Rue and Held (2005) and lecture notes.

$x_1,x_2$ are conditionally independent given $x_3$ if, for a given value of $x_3$, learning $x_2$ gives one no additional information about $x_1$. The density representation is therefore

$$\mathsf{f}\left(\mathbf{x}\right)=\mathsf{f}\left(x_1|x_3\right)\mathsf{f}\left(x_2|x_3\right)\mathsf{f}\left(x_3\right)$$

which is a simplification of the general representation.

$$\mathsf{f}\left(\mathbf{x}\right)=\mathsf{f}\left(x_1|x_2,x_3\right)\mathsf{f}\left(x_2|x_3\right)\mathsf{f}\left(x_3\right)$$ $$x\perp \perp y|z⟺\mathsf{f}\left(x,y,z\right)=g(x,z)h(y,z)$$

for some functions $f,g$ , and $\forall z\text{ with }\mathsf{f}\left(z\right)>0$

$$x_t=\varphi x_{t-1}+{\epsilon }_t;\epsilon \stackrel{\mathrm{iid}}{\sim }\mathcal{N}\left(0,1\right),\left|\varphi \right|<1$$

This can be re-expressed as

$$x_t|x_1,\dots ,x_{t-1}\sim \mathcal{N}\left(\varphi x_{t-1},1\right)\forall t=2,\dots ,n$$

So, for $x_s,x_t,1\le s<t\le n$,

$$x_s\perp \perp x_t|\left\{x_{s+1},\dots ,x_{t-1}\right\}\text{ if}t-s>1$$

In addition to the conditional distribution, also assume the marginal distribution of $x_1\stackrel{\mathrm{iid}}{\sim }\mathcal{N}\left(0,1/(1-\varphi {)}^2\right)$, which is the stationary distribution of this process. Then, the join distribution of $\mathbf{x}$ is

$$\begin{array}{rl}\mathsf{f}\left(\mathbf{x}\right)& =\mathsf{f}\left(x_1\right)\mathsf{f}\left(x_2|x_1\right)\dots ,\mathsf{f}\left(x_n|x_{n-1}\right)\\ & =\frac{1}{(2\pi {)}^{n/2}}{\left|\mathbf{Q}\right|}^{1/2}\exp \left(-\frac{1}{2}{\mathbf{x}}^{\prime }\mathbf{Qx}\right)\end{array}$$

where $\mathbf{Q}$ is a precision matrix of the form

$$\mathbf{Q}=\left(\begin{array}{ccccccc}1& -\varphi & & & & & \\ -\varphi & 1+{\varphi }^2& -\varphi & & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & & -\varphi & 1+{\varphi }^2& -\varphi \\ & & & & & -\varphi & 1\end{array}\right)$$

This tridiagonal form is due to the fact that $x_i\perp \perp x_j$ if $\left|i-j\right|>1$ given the rest of the sequence. This is generally true for any GMRF: $Q_{ij}=0,i\ne j⟹x_i\perp \perp x_j|\left\{x_k:k\ne i,j\right\}$.

While the conditional independence structure is readily apparent from the precision matrix, it isn’t evident in the covariance matrix $\mathbf{\Sigma }={\mathbf{Q}}^{-1}$, which is completely dense with entries

$${\sigma }_{ij}=\frac{1}{1-{\varphi }^2}{\varphi }^{|i-j|}$$

Entries of the covariance matrix $\mathbf{\Sigma }$ only give direct information about the marginal dependence structure, not the conditional one.

A spatial process $Y(\mathbf{s})s\in \mathcal{D}\subset {\mathbb{R}}^2$ is said to follow a Gaussian Process if any realisation $\mathbf{Y}=(Y(\mathbf{s}{)}_1,\dots ,Y(\mathbf{s}{)}_n{)}^{\prime }$ at the finite number of locations ${\mathbf{s}}_1,\dots ,{\mathbf{s}}_n$ follows an $N-$ variate Gaussian. More precisely, let $\mu (\mathbf{s}:\mathcal{D}\to \mathbb{R}$ denote a mean function returning a mean at location $\mathbf{s}$ (typically assumed to be linear in covariates $\mathbf{X}(\mathbf{s})=(1,X_1(\mathbf{s}),\dots ,X_p(\mathbf{s}){)}^{\prime }$) and $\mathbb{C}({\mathbf{s}}_1,{\mathbf{s}}_2):{\mathcal{D}}^2\to {\mathbb{R}}^{+}$ denote a covariance function. Then, $Y(\mathbf{s})$ follows a spatial Gaussian process, and $\mathbf{Y}$ has a density

$${\mathsf{f}}_{\mathbf{Y}}\left(\mathbf{y}\right)=\left(\frac{1}{\sqrt{2\pi }}\right){\left|\mathbf{\Sigma }\right|}^{-1/2}\exp \left\{-\frac{1}{2}(\mathbf{y}-\mathbf{\mu }{)}^{\prime }{\left(\mathbf{\Sigma }\right)}^{-1}(\mathbf{y}-\mathbf{\mu })\right\}$$

Where $\mathbf{\mu }=(\mu ({\mathbf{s}}_1),\dots ,\mu ({\mathbf{s}}_N){)}^{\prime }$ is the mean vector and $\mathbf{\Sigma }={\left\{\mathbb{C}({\mathbf{s}}_i,{\mathbf{s}}_j)\right\}}_{ij}$ is the $N\times N$ covariance matrix. Evaluating this density requires $\mathcal{O}(N^3)$ operations and $\mathcal{O}(N^2)$ memory, which means it does not scale well with large datasets. See Heaton et al. (2019) for an overview of alternatives.

Let $\mathbf{x}$ be associated with some property of points (typically location), with no natural ordering of the indices. The joint density of a zero-mean GMRF is specified by each of the $n$ full-conditionals

$$x_i|{\mathbf{x}}_{-i}\sim \mathcal{N}\left(\sum _{j:j\ne i}{\beta }_{ij}{x}_j,{\left(\kappa \right)}_i^{-1}\right)$$

these are called CAR models. The associated precision matrix is

$$\mathbf{Q}=Q_{ij}=\{\begin{array}{ll}{\kappa }_i& i=j\\ -{\kappa }_i{\beta }_{ij}& i\ne j\end{array}$$

which is symmetric and positive-definite.

A random vector $\mathbf{x}=(x_1,\dots ,x_n{)}^{\prime }\in {\mathbb{R}}^n$ is called a GMRF wrt a labelled graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with mean $\mathbf{\mu }$ and precision matrix $\mathbf{Q}>0$ iff its density has the form

$$\mathsf{f}\left(\mathbf{x}\right)=(2\pi {)}^{-n/2}{\left|\mathbf{Q}\right|}^{1/2}\exp \left(-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^{\prime }\mathbf{Q}(\mathbf{x}-\mathbf{\mu })\right)$$

and $Q_{ij}\ne 0⟺\{i,j\}\in \mathcal{E}\forall i\ne j$. If $\mathbf{Q}$ is completely dense, $\mathcal{G}$ is completely connected. In spatial settings, $\mathbf{Q}$ is typically sparse [depending on how neighbours are defined.]

Key summary quantities

$$\mathbb{E}\left[x_i|{\mathbf{x}}_{-i}\right]={\mu }_i-\frac{1}{Q_{ii}}\sum _{j:j\sim i}{Q}_{ij}(x_j-{\mu }_j)$$

$\text{Prec}(x_i,{\mathbf{x}}_{-i})=Q_{ii}$ and

$$\text{Corr}(x_i,x_j|{\mathbf{x}}_{-ij})=\frac{-Q_{ij}}{\sqrt{Q_{ii}{Q}_{jj}}},i\ne j$$

Let $\mathbf{x}$ be a GMRF wrt $\mathcal{G}=(\mathcal{V},\mathcal{E})$. The following are equivalent

Pairwise Markov Property: $x_i\perp x_j|{\mathbf{x}}_{-ij}\text{ if}\left\{i,j\right\}\notin \mathcal{E}\wedge i\ne j$

Local Markov Property; $x_i\perp {\mathbf{x}}_{\left\{i,ne(i)\right\}}|{\mathbf{x}}_{ne(i)}\forall i\in \mathcal{V}$

Global Markov: ${\mathbf{x}}_A\perp {\mathbf{x}}_B|{\mathbf{x}}_C$ for disjoint sets $A,B,C$ where $C$ separates $A,B$ and $A\text{ and }B$ are nonempty.

let the spatial process at location $\mathbf{s}\in \mathcal{D}$ be

$$Z(\mathbf{s})=\mathbf{X}(\mathbf{s})\mathbf{\beta }+w(\mathbf{s}),\forall \mathbf{s}\in \mathcal{D}$$

where $\mathbf{X}(\mathbf{s})$ collects a $p-$ vectors of covariates for site $\mathbf{s}$, and $\mathbf{\beta }$ is a p-vector of coefficients. Spatial dependence can be imposed by modelling $\left\{w(\mathbf{s}):\mathbf{s}\in \mathcal{D}\right\}$ as a zero-mean stationary Gaussian Process. Distributionally, this implies that for any ${\mathbf{s}}_1,\dots ,{\mathbf{s}}_n\in D$, if we let $\mathbf{w}=(w({\mathbf{s}}_1),\dots ,w({\mathbf{s}}_n){)}^{\prime }$, and $\mathbf{\Theta }$ be the parameters of the model

$$\mathbf{w}|\mathbf{\Theta }\sim \mathcal{N}\left(0,\mathbf{\Sigma }(\mathbf{\Theta })\right)$$

where $\mathbf{\Sigma }(\mathbf{\Theta })$ is the covariance matrix of a n-dimensional normal density. We need $\mathbf{\Sigma }(\mathbf{\Theta })$ to be Symmetric, PD for this distribution to be proper.

Special cases:

Exponential Covariance Matrix: $\mathbf{\Theta }=(\psi ,\varphi ,\kappa )$ $\mathbf{\Sigma }(\mathbf{\Theta })=\psi \mathbf{I}+\kappa \mathbf{H}(\varphi )$, where the $i,j$th element of $\mathbf{H}(\varphi )=\exp (-\Vert {\mathbf{s}}_i-{\mathbf{s}}_j\Vert /\varphi )$. The ‘nugget’ $\psi$ is the variance of the non-spatial error, $\kappa$ dictates the scale, and $\varphi$ dictates the range of the spatial dependence.

Matern Covariance: $\mathbf{\Theta }=(\psi ,\kappa ,\varphi ,\nu )>0$ for distance $x:=\Vert {\mathbf{s}}_i-{\mathbf{s}}_j\Vert$.

$$\mathrm{Cov}\left(x;\varphi ,\psi ,\kappa ,\nu \right)=\{\begin{array}{ll}\frac{\kappa }{2^{\nu -1}\Gamma (\nu )}(2\sqrt{\nu }x/\varphi {)}^{\nu }{K}_{\nu }(2\sqrt{\nu }x/\varphi )& \text{ if}x>0\\ \psi +\kappa & \text{ if}x=0\end{array}$$

where $K_{\nu }(x)$ is a modified Bessel function of order $\nu$.

Specifying $\Sigma$ directly can be awkward when dealing with irregular spatial data [i.e. every real use case].

So, random effects $\mathbf{w}$ are modelled conditionally. Let ${\mathbf{w}}_{-i}$ denote the vector of $\mathbf{w}$ excluding $\mathbf{w}({\mathbf{s}}_i)$. Model $w({\mathbf{s}}_i)$ in terms of its full-conditional.

$$w({\mathbf{s}}_i)|{\mathbf{w}}_{-i},\mathbf{\Theta }\sim \mathcal{N}\left(\sum _{j=1}^n{c}_{ij}w({\mathbf{s}}_j),{\kappa }_i^{-1}\right),i=1,\dots ,n$$

where $c_{ij}$ describes the neighbourhood structure.

Besag (1974) proved that if $\mathbf{Q}$ is symmetric PD, with ${\kappa }_i$ in the diagonals and $-{\kappa }_i{c}_{ij}$ in the off-diagonals. $\mathbf{w}|\mathbf{\Theta }\sim \mathcal{N}\left(0,{\mathbf{Q}}^{-1}\right)$. Simplest version assumes common precision parameter ${\kappa }_i=\tau$.

Intrinsic GMRF: $f(\mathbf{w}|\mathbf{\Theta })\sim {\tau }^{(N-1)/2}\exp (-{\mathbf{w}}^{\prime }\mathbf{Q}(\tau )\mathbf{w})$. When $c_{ij}=1$ for neighbours (i.e. adjacency matrix instead of distances), it simplifies further to

$$(\mathbf{w}|\mathbf{\Theta })\sim {\tau }^{(N-1)/2}\exp \left(-\frac{1}{2}\sum _{i\sim j}(w({\mathbf{s}}_i)-w({\mathbf{s}}_j){)}^2\right)$$

Let $\left\{Z(\mathbf{s}):\mathbf{s}\in D\right\}$ and $\left\{w(\mathbf{s}:\mathbf{s}\in D)\right\}$ be two spatial processes on $D\subset {\mathbb{R}}^d(d\in {\mathbb{Z}}^{+})$. Assume $Z({\mathbf{s}}_i)$s are conditional independent given random effects $w({\mathbf{s}}_1),\dots ,w({\mathbf{s}}_n)$, and that $Z({\mathbf{s}}_i)$ follow some common distributional form, and

$$\mathbb{E}\left[Z({\mathbf{s}}_i)|\mathbf{w}\right]=\mu ({\mathbf{s}}_i)\forall i=1,\dots n$$

Let $\eta (\mathbf{s})=h(\mu (\mathbf{s}))$ for some known link function $h(\cdot )$ e.g. $h(x)=\log \left(\frac{x}{1-x}\right)$ for logit. Assume linear form for projection

$\eta (\mathbf{s})=\mathbf{X}(\mathbf{s})\mathbf{\beta }+w(\mathbf{s})$. Spatial dependence via $\mathbf{w}|\mathbf{\Theta }\sim \mathcal{N}\left(0,\mathbf{\Sigma }(\mathbf{\Theta })\right)$, where $\Sigma$ is often Matern.