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

#### 2020-04-14

library(Ryacas)
library(Matrix)

## 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 <- ysym(L1chr)
L1s
## {{ 1,  0,  0,  0},
##  {-a,  1,  0,  0},
##  { 0, -a,  1,  0},
##  { 0,  0, -a,  1}}

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'$$

K1s <- L1s %*% t(L1s)
V1s <- solve(K1s)
cat(
"\\begin{align} K_1 &= ", tex(K1s), " \\\\
V_1 &= ", tex(V1s), " \\end{align}", sep = "")

\begin{align} K_1 &= \left( \begin{array}{cccc} 1 & - a & 0 & 0 \\ - a & a ^{2} + 1 & - a & 0 \\ 0 & - a & a ^{2} + 1 & - a \\ 0 & 0 & - a & a ^{2} + 1 \end{array} \right) \\ V_1 &= \left( \begin{array}{cccc} a ^{2} + a ^{4} + a ^{6} + 1 & a + a ^{3} + a ^{5} & a ^{2} + a ^{4} & a ^{3} \\ a + a ^{3} + a ^{5} & a ^{2} + a ^{4} + 1 & a + a ^{3} & a ^{2} \\ a ^{2} + a ^{4} & a + a ^{3} & a ^{2} + 1 & a \\ a ^{3} & a ^{2} & a & 1 \end{array} \right) \end{align}

## 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 <- ysym(L2chr)
L2s
## {{ 1,  0,  0,  0,  0,  0,  0},
##  {-a,  1,  0,  0,  0,  0,  0},
##  { 0, -a,  1,  0,  0,  0,  0},
##  { 0,  0, -a,  1,  0,  0,  0},
##  { 0, -b,  0,  0,  1,  0,  0},
##  { 0,  0, -b,  0,  0,  1,  0},
##  { 0,  0,  0, -b,  0,  0,  1}}
K2s <- L2s %*% t(L2s)
V2s <- simplify(solve(K2s))
cat(
"\\begin{align} K_2 &= ", tex(K2s), " \\\\
V_2 &= ", tex(V2s), " \\end{align}", sep = "")

\begin{align} K_2 &= \left( \begin{array}{ccccccc} 1 & - a & 0 & 0 & 0 & 0 & 0 \\ - a & a ^{2} + 1 & - a & 0 & - b & 0 & 0 \\ 0 & - a & a ^{2} + 1 & - a & a b & - b & 0 \\ 0 & 0 & - a & a ^{2} + 1 & 0 & a b & - b \\ 0 & - b & b a & 0 & b ^{2} + 1 & 0 & 0 \\ 0 & 0 & - b & b a & 0 & b ^{2} + 1 & 0 \\ 0 & 0 & 0 & - b & 0 & 0 & b ^{2} + 1 \end{array} \right) \\ V_2 &= \left( \begin{array}{ccccccc} a ^{6} b ^{2} + a ^{6} + a ^{4} b ^{2} + a ^{4} + a ^{2} b ^{2} + a ^{2} + 1 & a ^{5} b ^{2} + a ^{5} + a ^{3} b ^{2} + a ^{3} + a b ^{2} + a & a ^{4} b ^{2} + a ^{4} + a ^{2} b ^{2} + a ^{2} & a ^{3} b ^{2} + a ^{3} & a b & b a ^{2} & b a ^{3} \\ a ^{5} b ^{2} + a ^{5} + a ^{3} b ^{2} + a ^{3} + a b ^{2} + a & a ^{4} b ^{2} + a ^{4} + a ^{2} b ^{2} + a ^{2} + b ^{2} + 1 & a ^{3} b ^{2} + a ^{3} + a b ^{2} + a & a ^{2} b ^{2} + a ^{2} & b & a b & b a ^{2} \\ a ^{4} b ^{2} + a ^{4} + a ^{2} b ^{2} + a ^{2} & a ^{3} b ^{2} + a ^{3} + a b ^{2} + a & a ^{2} b ^{2} + a ^{2} + b ^{2} + 1 & a b ^{2} + a & 0 & b & a b \\ a ^{3} b ^{2} + a ^{3} & a ^{2} b ^{2} + a ^{2} & a b ^{2} + a & b ^{2} + 1 & 0 & 0 & b \\ b a & b & 0 & 0 & 1 & 0 & 0 \\ b a ^{2} & b a & b & 0 & 0 & 1 & 0 \\ b a ^{3} & b a ^{2} & b a & b & 0 & 0 & 1 \end{array} \right) \end{align}

## Numerical evalation in R

sparsify <- function(x) {
if (requireNamespace("Matrix", quietly = TRUE)) {
library(Matrix)

return(Matrix::Matrix(x, sparse = TRUE))
}

return(x)
}

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(yac_expr(V1s), list(a = alpha))
V2s_eval <- eval(yac_expr(V2s), list(a = alpha, b = beta))
all.equal(V1, V1s_eval)
## [1] TRUE
all.equal(V2, V2s_eval)
## [1] TRUE