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

library(Ryacas0)
library(Matrix)

Set output width:

get_output_width()

## [1] 120

set_output_width(120)
get_output_width()

## [1] 120

Autoregression (AR(1))

Consider AR(1) process: x_i = ax_i − 1 + e_i where i = 1, 2, 3 and where x₀ = e₀. Isolating error terms gives that e = L₁x where e = (e₀, …, e₃) and x = (x₀, …x₃) and where L₁ has the form

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

## Yacas matrix:
##      [,1] [,2] [,3] [,4]
## [1,] 1    0    0    0   
## [2,] -a   1    0    0   
## [3,] 0    -a   1    0   
## [4,] 0    0    -a   1

If error terms have variance 1 then Var(e) = LVar(x)L′ so the covariance matrix V1 = Var(x) = L⁻(L⁻)′ while the concentration matrix is K = LL′

# 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 = bx_i + u_i. Then (e, u) can be expressed in terms of (x, y) as (e, u) = L₂(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

## Yacas matrix:
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1    0    0    0    0    0    0   
## [2,] -a   1    0    0    0    0    0   
## [3,] 0    -a   1    0    0    0    0   
## [4,] 0    0    -a   1    0    0    0   
## [5,] 0    -b   0    0    1    0    0   
## [6,] 0    0    -b   0    0    1    0   
## [7,] 0    0    0    -b   0    0    1

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)

## 4 x 4 sparse Matrix of class "dsCMatrix"
##                            
## [1,]  1.0 -0.50  .     .   
## [2,] -0.5  1.25 -0.50  .   
## [3,]  .   -0.50  1.25 -0.50
## [4,]  .    .    -0.50  1.25

sparsify(V1)

## 4 x 4 sparse Matrix of class "dsCMatrix"
##                                   
## [1,] 1.328125 0.65625 0.3125 0.125
## [2,] 0.656250 1.31250 0.6250 0.250
## [3,] 0.312500 0.62500 1.2500 0.500
## [4,] 0.125000 0.25000 0.5000 1.000

## 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)

## 7 x 7 sparse Matrix of class "dsCMatrix"
##                                             
## [1,]  1.0 -0.50  .     .     .     .    .   
## [2,] -0.5  1.25 -0.50  .     0.30  .    .   
## [3,]  .   -0.50  1.25 -0.50 -0.15  0.30 .   
## [4,]  .    .    -0.50  1.25  .    -0.15 0.30
## [5,]  .    0.30 -0.15  .     1.09  .    .   
## [6,]  .    .     0.30 -0.15  .     1.09 .   
## [7,]  .    .     .     0.30  .     .    1.09

sparsify(V2)

## 7 x 7 sparse Matrix of class "dsCMatrix"
##                                                                   
## [1,]  1.3576563  0.7153125  0.340625  0.13625 -0.15 -0.075 -0.0375
## [2,]  0.7153125  1.4306250  0.681250  0.27250 -0.30 -0.150 -0.0750
## [3,]  0.3406250  0.6812500  1.362500  0.54500  .    -0.300 -0.1500
## [4,]  0.1362500  0.2725000  0.545000  1.09000  .     .     -0.3000
## [5,] -0.1500000 -0.3000000  .         .        1.00  .      .     
## [6,] -0.0750000 -0.1500000 -0.300000  .        .     1.000  .     
## [7,] -0.0375000 -0.0750000 -0.150000 -0.30000  .     .      1.0000

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)

## [1] TRUE

all.equal(V2, V2s_eval)

## [1] TRUE