The structure of the concentration and covariance matrix in a simple state-space model

In Ryacas0: Legacy 'Ryacas' (Interface to 'Yacas' Computer Algebra System)

knitr::opts_chunk$set(echo = TRUE)

library(Ryacas0)
library(Matrix)

Set output width:

get_output_width()
set_output_width(120)
get_output_width()

Autoregression ($AR(1)$)

Consider $AR(1)$ process: $x_i = a x_{i-1} + e_i$ where $i=1,2,3$ and where $x_0=e_0$. Isolating error terms gives that $$ e = L_1 x $$ where $e=(e_0, \dots, e_3)$ and $x=(x_0, \dots x_3)$ and where $L_1$ has the form

N <- 3
L1chr <- diag("1", 1 + N)
L1chr[cbind(1+(1:N), 1:N)] <- "-a"
L1s <- as.Sym(L1chr)
L1s

If error terms have variance $1$ then $\mathbf{Var}(e)=L \mathbf{Var}(x) L'$ so the covariance matrix $V1=\mathbf{Var}(x) = L^- (L^-)'$ while the concentration matrix is $K=L L'$

# FIXME: * vs %*%
K1s <- Simplify(L1s * Transpose(L1s))
V1s <- Simplify(Inverse(K1s))

cat(
  "\\begin{align} K_1 &= ", TeXForm(K1s), " \\\\ 
                  V_1 &= ", TeXForm(V1s), " \\end{align}", sep = "")

Dynamic linear model

Augument the $AR(1)$ process above with $y_i=b x_i + u_i$. Then $(e,u)$ can be expressed in terms of $(x,y)$ as $$ (e,u) = L_2(x,y) $$ where

N <- 3
L2chr <- diag("1", 1 + 2*N)
L2chr[cbind(1+(1:N), 1:N)] <- "-a"
L2chr[cbind(1 + N + (1:N), 1 + 1:N)] <- "-b"
L2s <- as.Sym(L2chr)
L2s

K2s <- Simplify(L2s * Transpose(L2s))
V2s <- Simplify(Inverse(K2s))

cat(
  "\\begin{align} K_2 &= ", TeXForm(K2s), " \\\\ 
                  V_2 &= ", TeXForm(V2s), " \\end{align}", sep = "")

Numerical evalation in R

sparsify <- function(x) {
  Matrix::Matrix(x, sparse = TRUE)
}

alpha <- 0.5
beta <- -0.3

## AR(1)
N <- 3
L1 <- diag(1, 1 + N)
L1[cbind(1+(1:N), 1:N)] <- -alpha
K1 <- L1 %*% t(L1)
V1 <- solve(K1)
sparsify(K1)
sparsify(V1)

## Dynamic linear models
N <- 3
L2 <- diag(1, 1 + 2*N)
L2[cbind(1+(1:N), 1:N)] <- -alpha
L2[cbind(1 + N + (1:N), 1 + 1:N)] <- -beta
K2 <- L2 %*% t(L2)
V2 <- solve(K2)
sparsify(K2)
sparsify(V2)

Comparing with results calculated by yacas:

V1s_eval <- Eval(V1s, list(a = alpha))
V2s_eval <- Eval(V2s, list(a = alpha, b = beta))

all.equal(V1, V1s_eval)
all.equal(V2, V2s_eval)

Ryacas0 documentation built on Jan. 12, 2023, 5:11 p.m.