25. Linear State Space Models#

25.2.1. Primitives#

The primitives of the model are

1. the matrices $A,C,G$

2. shock distribution, which we have specialized to $N(0,I)$

3. the distribution of the initial condition $x_0$, which we have set to $N({\mu }_0,{\mathrm{\Sigma }}_0)$

Given $A,C,G$ and draws of $x_0$ and $w_1,w_2,\dots$, the model (25.1) pins down the values of the sequences $\{x_t\}$ and $\{y_t\}$.

Even without these draws, the primitives 1–3 pin down the probability distributions of $\{x_t\}$ and $\{y_t\}$.

Later we’ll see how to compute these distributions and their moments.

25.2.2. Examples#

By appropriate choice of the primitives, a variety of dynamics can be represented in terms of the linear state space model.

The following examples help to highlight this point.

They also illustrate the wise dictum finding the state is an art.

25.2.2.1. Second-order difference equation#

Let $\{y_t\}$ be a deterministic sequence that satisfies

(25.2)#$$y_{t+1}={\varphi }_0+{\varphi }_1{y}_t+{\varphi }_2{y}_{t-1}\text{s.t.}{y}_0,y_{-1}\text{ given}$$

To map (25.2) into our state space system (25.1), we set

$$\begin{array}{r}{x}_t=\left[\begin{array}{c}1\\ y_t\\ y_{t-1}\end{array}\right]A=\left[\begin{array}{ccc}1& 0& 0\\ {\varphi }_0& {\varphi }_1& {\varphi }_2\\ 0& 1& 0\end{array}\right]C=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]G=\left[\begin{array}{ccc}0& 1& 0\end{array}\right]\end{array}$$

You can confirm that under these definitions, (25.1) and (25.2) agree.

The next figure shows dynamics of this process when ${\varphi }_0=1.1,{\varphi }_1=0.8,{\varphi }_2=-0.8,y_0=y_{-1}=1$

Later you’ll be asked to recreate this figure.

25.2.2.2. Univariate Autoregressive Processes#

We can use (25.1) to represent the model

(25.3)#$$y_{t+1}={\varphi }_1{y}_t+{\varphi }_2{y}_{t-1}+{\varphi }_3{y}_{t-2}+{\varphi }_4{y}_{t-3}+\sigma w_{t+1}$$

where $\{w_t\}$ is iid and standard normal.

To put this in the linear state space format we take $x_t={\left[\begin{array}{cccc}{y}_t& y_{t-1}& y_{t-2}& y_{t-3}\end{array}\right]}^{\prime }$ and

$$\begin{array}{r}A=\left[\begin{array}{cccc}{\varphi }_1& {\varphi }_2& {\varphi }_3& {\varphi }_4\\ 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right]C=\left[\begin{array}{c}\sigma \\ 0\\ 0\\ 0\end{array}\right]G=\left[\begin{array}{cccc}1& 0& 0& 0\end{array}\right]\end{array}$$

The matrix $A$ has the form of the companion matrix to the vector $\left[\begin{array}{cccc}{\varphi }_1& {\varphi }_2& {\varphi }_3& {\varphi }_4\end{array}\right]$.

The next figure shows dynamics of this process when

$${\varphi }_1=0.5,{\varphi }_2=-0.2,{\varphi }_3=0,{\varphi }_4=0.5,\sigma =0.2,y_0=y_{-1}=y_{-2}=y_{-3}=1$$

25.2.3. Moving Average Representations#

A nonrecursive expression for $x_t$ as a function of $x_0,w_1,w_2,\dots ,w_t$ can be found by using (25.1) repeatedly to obtain

Representation (25.5) is a moving average representation.

It expresses $\{x_t\}$ as a linear function of

1. current and past values of the process $\{w_t\}$ and

2. the initial condition $x_0$

As an example of a moving average representation, let the model be

You will be able to show that $A^t=\left[\begin{array}{cc}1& t\\ 0& 1\end{array}\right]$ and $A^{j}C={\left[\begin{array}{cc}1& 0\end{array}\right]}^{\prime }$.

Substituting into the moving average representation (25.5), we obtain

where $x_{1t}$ is the first entry of $x_t$.

The first term on the right is a cumulated sum of martingale differences, and is therefore a martingale.

The second term is a translated linear function of time.

For this reason, $x_{1t}$ is called a martingale with drift.

25.3.1. Unconditional Moments#

Using (25.1), it’s easy to obtain expressions for the (unconditional) means of $x_t$ and $y_t$.

We’ll explain what unconditional and conditional mean soon.

Letting ${\mu }_t:=\mathbb{E}[x_t]$ and using linearity of expectations, we find that

Here ${\mu }_0$ is a primitive given in (25.1).

The variance-covariance matrix of $x_t$ is ${\mathrm{\Sigma }}_t:=\mathbb{E}[(x_t-{\mu }_t)(x_t-{\mu }_t{)}^{\prime }]$.

Using $x_{t+1}-{\mu }_{t+1}=A(x_t-{\mu }_t)+C{w}_{t+1}$, we can determine this matrix recursively via

As with ${\mu }_0$, the matrix ${\mathrm{\Sigma }}_0$ is a primitive given in (25.1).

As a matter of terminology, we will sometimes call

• ${\mu }_t$ the unconditional mean of $x_t$

• ${\mathrm{\Sigma }}_t$ the unconditional variance-convariance matrix of $x_t$

This is to distinguish ${\mu }_t$ and ${\mathrm{\Sigma }}_t$ from related objects that use conditioning information, to be defined below.

However, you should be aware that these “unconditional” moments do depend on the initial distribution $N({\mu }_0,{\mathrm{\Sigma }}_0)$.

25.3.2. Distributions#

In general, knowing the mean and variance-covariance matrix of a random vector is not quite as good as knowing the full distribution.

However, there are some situations where these moments alone tell us all we need to know.

These are situations in which the mean vector and covariance matrix are sufficient statistics for the population distribution.

(Sufficient statistics form a list of objects that characterize a population distribution)

One such situation is when the vector in question is Gaussian (i.e., normally distributed).

This is the case here, given

1. our Gaussian assumptions on the primitives

2. the fact that normality is preserved under linear operations

In fact, it’s well-known that

In particular, given our Gaussian assumptions on the primitives and the linearity of (25.1) we can see immediately that both $x_t$ and $y_t$ are Gaussian for all $t\ge 0$ 2.

Since $x_t$ is Gaussian, to find the distribution, all we need to do is find its mean and variance-covariance matrix.

But in fact we’ve already done this, in (25.6) and (25.7).

Letting ${\mu }_t$ and ${\mathrm{\Sigma }}_t$ be as defined by these equations, we have

By similar reasoning combined with (25.8) and (25.9),

25.3.3. Ensemble Interpretations#

How should we interpret the distributions defined by (25.11)–(25.12)?

Intuitively, the probabilities in a distribution correspond to relative frequencies in a large population drawn from that distribution.

Let’s apply this idea to our setting, focusing on the distribution of $y_T$ for fixed $T$.

We can generate independent draws of $y_T$ by repeatedly simulating the evolution of the system up to time $T$, using an independent set of shocks each time.

The next figure shows 20 simulations, producing 20 time series for $\{y_t\}$, and hence 20 draws of $y_T$.

The system in question is the univariate autoregressive model (25.3).

The values of $y_T$ are represented by black dots in the left-hand figure

In the right-hand figure, these values are converted into a rotated histogram that shows relative frequencies from our sample of 20 $y_T$’s.

(The parameters and source code for the figures can be found in file linear_models/paths_and_hist.jl)

Here is another figure, this time with 100 observations

Let’s now try with 500,000 observations, showing only the histogram (without rotation)

The black line is the population density of $y_T$ calculated from (25.12).

The histogram and population distribution are close, as expected.

By looking at the figures and experimenting with parameters, you will gain a feel for how the population distribution depends on the model primitives listed above, as intermediated by the distribution’s sufficient statistics.

25.3.3.1. Ensemble means#

In the preceding figure we approximated the population distribution of $y_T$ by

1. generating $I$ sample paths (i.e., time series) where $I$ is a large number

2. recording each observation $y_T^i$

3. histogramming this sample

Just as the histogram approximates the population distribution, the ensemble or cross-sectional average

$${\overline{y}}_T:=\frac{1}{I}\sum _{i=1}^I{y}_T^i$$

approximates the expectation $\mathbb{E}[y_T]=G{\mu }_T$ (as implied by the law of large numbers).

Here’s a simulation comparing the ensemble averages and population means at time points $t=0,\dots ,50$.

The parameters are the same as for the preceding figures, and the sample size is relatively small ($I=20$).

The ensemble mean for $x_t$ is

$${\overline{x}}_T:=\frac{1}{I}\sum _{i=1}^I{x}_T^i\to {\mu }_T(I\to \mathrm{\infty })$$

The limit ${\mu }_T$ is a “long-run average”.

(By long-run average we mean the average for an infinite ($I=\mathrm{\infty }$) number of sample $x_T$’s)

Another application of the law of large numbers assures us that

$$\frac{1}{I}\sum _{i=1}^I(x_T^i-{\overline{x}}_T)(x_T^i-{\overline{x}}_T{)}^{\prime }\to {\mathrm{\Sigma }}_T(I\to \mathrm{\infty })$$

25.3.4. Joint Distributions#

In the preceding discussion we looked at the distributions of $x_t$ and $y_t$ in isolation.

This gives us useful information, but doesn’t allow us to answer questions like

• what’s the probability that $x_t\ge 0$ for all $t$?

• what’s the probability that the process $\{y_t\}$ exceeds some value $a$ before falling below $b$?

• etc., etc.

Such questions concern the joint distributions of these sequences.

To compute the joint distribution of $x_0,x_1,\dots ,x_T$, recall that joint and conditional densities are linked by the rule

From this rule we get $p(x_0,x_1)=p(x_1|x_0)p(x_0)$.

The Markov property $p(x_t|x_{t-1},\dots ,x_0)=p(x_t|x_{t-1})$ and repeated applications of the preceding rule lead us to

The marginal $p(x_0)$ is just the primitive $N({\mu }_0,{\mathrm{\Sigma }}_0)$.

In view of (25.1), the conditional densities are

25.4.1. Visualizing Stability#

Let’s look at some more time series from the same model that we analyzed above.

This picture shows cross-sectional distributions for $y$ at times $T,T^{\prime },T^{″}$

Note how the time series “settle down” in the sense that the distributions at $T^{\prime }$ and $T^{″}$ are relatively similar to each other — but unlike the distribution at $T$.

Apparently, the distributions of $y_t$ converge to a fixed long-run distribution as $t\to \mathrm{\infty }$.

When such a distribution exists it is called a stationary distribution.

25.4.2. Stationary Distributions#

In our setting, a distribution ${\psi }_{\mathrm{\infty }}$ is said to be stationary for $x_t$ if

Since

1. in the present case all distributions are Gaussian

2. a Gaussian distribution is pinned down by its mean and variance-covariance matrix

we can restate the definition as follows: ${\psi }_{\mathrm{\infty }}$ is stationary for $x_t$ if

where ${\mu }_{\mathrm{\infty }}$ and ${\mathrm{\Sigma }}_{\mathrm{\infty }}$ are fixed points of (25.6) and (25.7) respectively.

25.4.3. Covariance Stationary Processes#

Let’s see what happens to the preceding figure if we start $x_0$ at the stationary distribution.

Now the differences in the observed distributions at $T,T^{\prime }$ and $T^{″}$ come entirely from random fluctuations due to the finite sample size.

By

• our choosing $x_0\sim N({\mu }_{\mathrm{\infty }},{\mathrm{\Sigma }}_{\mathrm{\infty }})$

• the definitions of ${\mu }_{\mathrm{\infty }}$ and ${\mathrm{\Sigma }}_{\mathrm{\infty }}$ as fixed points of (25.6) and (25.7) respectively

we’ve ensured that

Moreover, in view of (25.14), the autocovariance function takes the form ${\mathrm{\Sigma }}_{t+j,t}=A^j{\mathrm{\Sigma }}_{\mathrm{\infty }}$, which depends on $j$ but not on $t$.

This motivates the following definition.

A process $\{x_t\}$ is said to be covariance stationary if

• both ${\mu }_t$ and ${\mathrm{\Sigma }}_t$ are constant in $t$

• ${\mathrm{\Sigma }}_{t+j,t}$ depends on the time gap $j$ but not on time $t$

In our setting, $\{x_t\}$ will be covariance stationary if ${\mu }_0,{\mathrm{\Sigma }}_0,A,C$ assume values that imply that none of ${\mu }_t,{\mathrm{\Sigma }}_t,{\mathrm{\Sigma }}_{t+j,t}$ depends on $t$.

25.4.4. Conditions for Stationarity#

25.4.5. Ergodicity#

Let’s suppose that we’re working with a covariance stationary process.

In this case we know that the ensemble mean will converge to ${\mu }_{\mathrm{\infty }}$ as the sample size $I$ approaches infinity.

25.6.1. Forecasting Formulas – Conditional Means#

The natural way to predict variables is to use conditional distributions.

For example, the optimal forecast of $x_{t+1}$ given information known at time $t$ is

The right-hand side follows from $x_{t+1}=A{x}_t+C{w}_{t+1}$ and the fact that $w_{t+1}$ is zero mean and independent of $x_t,x_{t-1},\dots ,x_0$.

That ${\mathbb{E}}_t[x_{t+1}]=\mathbb{E}[x_{t+1}\mid x_t]$ is an implication of $\{x_t\}$ having the Markov property.

The one-step-ahead forecast error is

The covariance matrix of the forecast error is

More generally, we’d like to compute the $j$-step ahead forecasts ${\mathbb{E}}_t[x_{t+j}]$ and ${\mathbb{E}}_t[y_{t+j}]$.

With a bit of algebra we obtain

In view of the iid property, current and past state values provide no information about future values of the shock.

Hence ${\mathbb{E}}_t[w_{t+k}]=\mathbb{E}[w_{t+k}]=0$.

It now follows from linearity of expectations that the $j$-step ahead forecast of $x$ is

The $j$-step ahead forecast of $y$ is therefore

25.6.2. Covariance of Prediction Errors#

It is useful to obtain the covariance matrix of the vector of $j$-step-ahead prediction errors

Evidently,

$V_j$ defined in (25.21) can be calculated recursively via $V_1=C{C}^{\prime }$ and

$V_j$ is the conditional covariance matrix of the errors in forecasting $x_{t+j}$, conditioned on time $t$ information $x_t$.

Under particular conditions, $V_j$ converges to

Equation (25.23) is an example of a discrete Lyapunov equation in the covariance matrix $V_{\mathrm{\infty }}$.

A sufficient condition for $V_j$ to converge is that the eigenvalues of $A$ be strictly less than one in modulus.

Weaker sufficient conditions for convergence associate eigenvalues equaling or exceeding one in modulus with elements of $C$ that equal $0$.

25.6.3. Forecasts of Geometric Sums#

In several contexts, we want to compute forecasts of geometric sums of future random variables governed by the linear state-space system (25.1).

We want the following objects

• Forecast of a geometric sum of future $x$’s, or ${\mathbb{E}}_t\left[\sum _{j=0}^{\mathrm{\infty }}{\beta }^j{x}_{t+j}\right]$.

• Forecast of a geometric sum of future $y$’s, or ${\mathbb{E}}_t\left[\sum _{j=0}^{\mathrm{\infty }}{\beta }^j{y}_{t+j}\right]$.

These objects are important components of some famous and interesting dynamic models.

For example,

• if $\{y_t\}$ is a stream of dividends, then $\mathbb{E}\left[\sum _{j=0}^{\mathrm{\infty }}{\beta }^j{y}_{t+j}|x_t\right]$ is a model of a stock price

• if $\{y_t\}$ is the money supply, then $\mathbb{E}\left[\sum _{j=0}^{\mathrm{\infty }}{\beta }^j{y}_{t+j}|x_t\right]$ is a model of the price level

25.8.1. Exercise 1#

Replicate this figure using the LSS type from lss.jl.

25.8.2. Exercise 2#

Replicate this figure modulo randomness using the same type.

25.8.3. Exercise 3#

Replicate this figure modulo randomness using the same type.

The state space model and parameters are the same as for the preceding exercise.

25.8.4. Exercise 4#

Replicate this figure modulo randomness using the same type.

The state space model and parameters are the same as for the preceding exercise, except that the initial condition is the stationary distribution.

Hint: You can use the stationary_distributions method to get the initial conditions.

The number of sample paths is 80, and the time horizon in the figure is 100.

Producing the vertical bars and dots is optional, but if you wish to try, the bars are at dates 10, 50 and 75.

25.9.1. Exercise 1#

phi0, phi1, phi2 = 1.1, 0.8, -0.8

A = [1.0 0.0 0
     phi0 phi1 phi2
     0.0 1.0 0.0]
C = zeros(3, 1)
G = [0.0 1.0 0.0]
mu_0 = ones(3)

lss = LSS(A, C, G; mu_0)

x, y = simulate(lss, 50)
plot(dropdims(y, dims = 1), color = :blue, linewidth = 2, alpha = 0.7)
plot!(xlabel = "time", ylabel = L"y_t", legend = :none)

25.9.2. Exercise 2#

using Random
Random.seed!(42) # For deterministic results.

phi1, phi2, phi3, phi4 = 0.5, -0.2, 0, 0.5
sigma = 0.2

A = [phi1 phi2 phi3 phi4
     1.0 0.0 0.0 0.0
     0.0 1.0 0.0 0.0
     0.0 0.0 1.0 0.0]
C = [sigma
     0.0
     0.0
     0.0]
G = [1.0 0.0 0.0 0.0]

ar = LSS(A, C, G; mu_0 = ones(4))
x, y = simulate(ar, 200)

plot(dropdims(y, dims = 1), color = :blue, linewidth = 2, alpha = 0.7)
plot!(xlabel = "time", ylabel = L"y_t", legend = :none)

25.9.3. Exercise 3#

phi1, phi2, phi3, phi4 = 0.5, -0.2, 0, 0.5
sigma = 0.1

A = [phi1 phi2 phi3 phi4
     1.0 0.0 0.0 0.0
     0.0 1.0 0.0 0.0
     0.0 0.0 1.0 0.0]
C = [sigma
     0.0
     0.0
     0.0]
G = [1.0 0.0 0.0 0.0]
I = 20
T = 50
ar = LSS(A, C, G; mu_0 = ones(4))
ymin, ymax = -0.5, 1.15

ensemble_mean = zeros(T)
ys = []
for i in 1:I
    x, y = simulate(ar, T)
    y = dropdims(y, dims = 1)
    push!(ys, y)
    ensemble_mean .+= y
end

ensemble_mean = ensemble_mean ./ I
plot(ys, color = :blue, alpha = 0.2, linewidth = 0.8, label = "")
plot!(ensemble_mean, color = :blue, linewidth = 2, label = L"\overline{y_t}")
m = moment_sequence(ar)
pop_means = zeros(0)
for (i, t) in enumerate(m)
    (mu_x, mu_y, Sigma_x, Sigma_y) = t
    push!(pop_means, mu_y[1])
    i == 50 && break
end
plot!(pop_means, color = :green, linewidth = 2, label = L"G \mu_t")
plot!(ylims = (ymin, ymax), xlabel = "time", ylabel = L"y_t", legendfont = font(12))

25.9.4. Exercise 4#

phi1, phi2, phi3, phi4 = 0.5, -0.2, 0, 0.5
sigma = 0.1

A = [phi1 phi2 phi3 phi4
     1.0 0.0 0.0 0.0
     0.0 1.0 0.0 0.0
     0.0 0.0 1.0 0.0]
C = [sigma
     0.0
     0.0
     0.0]
G = [1.0 0.0 0.0 0.0]

T0 = 10
T1 = 50
T2 = 75
T4 = 100

ar = LSS(A, C, G; mu_0 = ones(4))
ymin, ymax = -0.6, 0.6

mu_x, mu_y, Sigma_x, Sigma_y = stationary_distributions(ar)
ar = LSS(A, C, G; mu_0 = mu_x, Sigma_0 = Sigma_x)
colors = ["c", "g", "b"]

ys = []
x_scatter = []
y_scatter = []
for i in 1:80
    rcolor = colors[rand(1:3)]
    x, y = simulate(ar, T4)
    y = dropdims(y, dims = 1)
    push!(ys, y)
    x_scatter = [x_scatter; T0; T1; T2]
    y_scatter = [y_scatter; y[T0]; y[T1]; y[T2]]
end

plot(ys, linewidth = 0.8, alpha = 0.5)
plot!([T0 T1 T2; T0 T1 T2], [-1 -1 -1; 1 1 1], color = :black, legend = :none)
scatter!(x_scatter, y_scatter, color = :black, alpha = 0.5)
plot!(ylims = (ymin, ymax), ylabel = L"y_t", xticks = [], yticks = ymin:0.2:ymax)
plot!(annotations = [(T0 + 1, -0.55, L"T"); (T1 + 1, -0.55, L"T^\prime");
                     (T2 + 1, -0.55, L"T^{\prime\prime}")])

1

The eigenvalues of $A$ are $(1,-1,i,-i)$.

2

The correct way to argue this is by induction. Suppose that $x_t$ is Gaussian. Then (25.1) and (25.10) imply that $x_{t+1}$ is Gaussian. Since $x_0$ is assumed to be Gaussian, it follows that every $x_t$ is Gaussian. Evidently this implies that each $y_t$ is Gaussian.