Gaussian dynamic linear model

The local level model

1 The local level model

The local level model, also known as first order dynamic linear model (aka DLM) where \[\begin{eqnarray*} y_t &=& \theta_t + \epsilon_t\\ \theta_t &=& \theta_{t-1} + \omega_t, \end{eqnarray*}\] with observational errors \(\epsilon_t\) iid \(N(0,\sigma^2)\) and systems errors \(\omega_t\) iid \(N(0,\tau^2)\) for \(t=1,\ldots,n\). Let \(y=(y_1,\ldots,y_n)'\) and \(\theta=(\theta_1,\ldots,\theta_n)'\), then we can show that

\[ p(\theta|\sigma^2, \tau^2) \propto \exp\left\{ -\frac{1}{2\tau^2} \sum_{t=1}^n (\theta_t - \theta_{t-1})^2 \right\} \]

An alternative representation based on the precision matrix \(K\): \[ \theta \sim N(0, \tau^2 K^{-1}), \] where the covariance matrix is written as \[ K^{-1} = \begin{pmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & 2 & 2 & \cdots & 2 \\ 1 & 2 & 3 & \cdots & 3 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & 2 & 3 & \cdots & n \end{pmatrix} \]

This is a very common way to present the random walk prior in the closed-form versus MC-based inference, say, in Gaussian Markov Random Field (GMRF) models.

It can be easily shown that the precision matrix is given by \[ K = \begin{pmatrix} 2 & -1 & 0 & \cdots & 0 \\ -1 & 2 & -1 & \cdots & 0 \\ 0 & -1 & 2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & -1 \\ 0 & 0 & 0 & -1 & 1 \end{pmatrix} \]

To sample from the posterior \(p(\theta|\sigma^2,\tau^2,y)\), we leverage the fact that for a Gaussian state-space model, the posterior is also a multivariate Gaussian distribution. In fact, the posterior precision matrix is given by \[ P = \frac{1}{\sigma^2}I + \frac{1}{\tau^2}K, \] while the posterior mean \(\widehat{\theta}\) solves \[ P\widehat{\theta} = \frac{1}{\sigma^2}y. \] To sample efficiently, we use the Cholesky decomposition \(P = LL'\).

2 Simulation exercise

library(Matrix)

# 1. Setup parameters and synthetic data
set.seed(42)
n      = 100
sigma2 = 1.0
tau2   = 0.5
theta_true = cumsum(rnorm(n, 0, sqrt(tau2)))
y          = theta_true + rnorm(n, 0, sqrt(sigma2))

# 2. Construct the Sparse Precision Matrix K
# Diagonals: 2 for all except the last element which is 1
diag_elements = c(rep(2, n-1), 1)
off_diag = rep(-1, n-1)

K = sparseMatrix(
  i = c(1:n, 1:(n-1), 2:n),
  j = c(1:n, 2:n, 1:(n-1)),
  x = c(diag_elements, off_diag, off_diag)
)

# 3. Compute Posterior Precision P and Mean
# P = (1/sigma^2)*I + (1/tau^2)*K
P = (1/sigma2) * Diagonal(n) + (1/tau2) * K

# 4. Cholesky Decomposition
# P = L %*% t(L)
L = chol(P) 

# 5. Sample from the Posterior
# Step A: Solve for the posterior mean (P * m = (1/sigma2) * y)
m = solve(P, y / sigma2)

# Step B: Generate Z ~ N(0, I) and solve L' * v = Z
z = rnorm(n)
v = solve(t(L), z)

# Step C: theta_sample = m + v
theta_sample = as.vector(m + v)

# 6. Quick Plot
plot(y, col = "grey", pch = 16, main = "Posterior Sample via Cholesky")
lines(theta_true, col = "black", lwd = 2, lty = 2)
lines(theta_sample, col = "blue", lwd = 2)
legend("topleft", legend=c("True", "Sample"), col=c("black", "blue"), lty=c(2,1))

3 Sampling from the posterior

To generate \(M=1000\) draws efficiently, we only perform the Cholesky decomposition once. Since the precision matrix \(P\) depends only on the hyperparameters (\(\sigma^2, \tau^2\)) and the model structure, it remains constant for all draws of \(\theta\).

library(Matrix)

M      = 1000

# Construct Sparse Precision Matrix K
# K is tridiagonal with 2 on diag (1 at the end) and -1 on off-diagonals
K = bandSparse(n,k=0:1,diagonals=list(c(rep(2,n-1),1),rep(-1, n-1)),symmetric=TRUE)

# Posterior Precision P and its Cholesky Decomposition
# P = (1/sigma^2)*I + (1/tau^2)*K
# L is upper triangular such that t(L) %*% L = P
P = (1/sigma2) * Diagonal(n) + (1/tau2) * K
L = chol(P) 

# Posterior Mean (m)
# Solve P*m = (1/sigma2)*y
m = solve(P, y / sigma2)

# Efficient Sampling: theta = m + inv(L) * Z
# We generate an n x M matrix of standard normal draws
Z = matrix(rnorm(n * M), nrow = n, ncol = M)

# Solve L*v = Z for all columns (since L is upper triangular)
# Note: backsolve/solve is very fast for sparse triangular matrices
v = solve(L, Z)

# Final draws: Add the mean vector to each column of v
# theta_draws: (n x M) matrix where each column is a draw
theta_draws = as.matrix(m + v)

We use the property that if \(Z \sim N(0, I)\), then \(\theta = m + L^{-1}Z\) results in \(\theta \sim N(m, (L'L)^{-1})\), which is exactly \(N(m, P^{-1})\).

plot(y, type = "n", main = "1,000 Posterior Draws of Theta")
matlines(theta_draws, col = rgb(0.1, 0.1, 0.9, 0.05), lty = 1) # Transparent blue lines
lines(as.numeric(m), col = "red", lwd = 2)                   # Posterior mean

4 Learning \(\theta\), \(\sigma^2\) and \(\tau^2\)

library(Matrix)

# 1. Priors and Initial Values
a_sig = b_sig = 0.01
a_tau = b_tau = 0.01

burnin = 5000
M = burnin + 5000
n = length(y)
theta_samples = matrix(0, n, M)
sig2_samples  = numeric(M)
tau2_samples  = numeric(M)

# Initial values
curr_sig2 = var(y)
curr_tau2 = var(diff(y))
K = bandSparse(n, k = 0:1, diagonals = list(c(rep(2, n-1), 1), rep(-1, n-1)), symmetric = TRUE)

# 2. Gibbs Loop
for (m in 1:M) {
  # --- Step A: Sample Theta ---
  P = (1/curr_sig2) * Diagonal(n) + (1/curr_tau2) * K
  L = chol(P)
  post_mean = solve(P, y / curr_sig2)
  theta_samples[, m] = as.vector(post_mean + solve(L, rnorm(n)))
  
  curr_theta = theta_samples[, m]
  
  # --- Step B: Sample sigma^2 ---
  sse_y = sum((y - curr_theta)^2)
  curr_sig2 = 1 / rgamma(1, a_sig + n/2, b_sig + sse_y/2)
  sig2_samples[m] = curr_sig2
  
  # --- Step C: Sample tau^2 ---
  # Using the quadratic form theta' * K * theta
  sse_theta = as.numeric(t(curr_theta) %*% K %*% curr_theta)
  curr_tau2 = 1 / rgamma(1, a_tau + n/2, b_tau + sse_theta/2)
  tau2_samples[m] = curr_tau2
}


theta_samples = theta_samples[,(burnin+1):M]
sig2_samples = sig2_samples[(burnin+1):M]
tau2_samples = tau2_samples[(burnin+1):M]

# Setup a 3x2 plotting area
# Layout: Rows = (Trace, ACF, Histogram), Cols = (sigma2, tau2)
par(mfrow = c(3, 2), mar = c(4, 4, 2, 1), oma = c(0, 0, 2, 0))

# --- Row 1: Trace Plots ---
plot(sig2_samples, type = "l", col = "steelblue", 
     main = expression(paste("Trace of ", sigma^2)), ylab = "Value")
plot(tau2_samples, type = "l", col = "darkorange", 
     main = expression(paste("Trace of ", tau^2)), ylab = "Value")

# --- Row 2: Autocorrelation Functions (ACF) ---
acf(sig2_samples, main = expression(paste("ACF of ", sigma^2)))
acf(tau2_samples, main = expression(paste("ACF of ", tau^2)))

# --- Row 3: Histograms ---
hist(sig2_samples, breaks = 30, col = "steelblue", border = "white",
     main = expression(paste("Posterior of ", sigma^2)), xlab = "Value")
abline(v = mean(sig2_samples), col = "red", lwd = 2, lty = 2)

hist(tau2_samples, breaks = 30, col = "darkorange", border = "white",
     main = expression(paste("Posterior of ", tau^2)), xlab = "Value")
abline(v = mean(tau2_samples), col = "red", lwd = 2, lty = 2)

# Add a main title for the article
mtext("MCMC Diagnostics for Variance Parameters", outer = TRUE, cex = 1.2, font = 2)

plot(y,type="n",xlab="Time",main=paste(M-burnin," posterior draws of theta",sep=""))
matlines(theta_draws, col = rgb(0.1, 0.1, 0.9, 0.05), lty = 1) # Transparent blue lines
points(y,pch=16)
lines(apply(theta_draws,1,mean),col="red",lwd=2)