# The structure of the concentration and covariance matrix in a simple state-space model In Ryacas: R Interface to the 'Yacas' Computer Algebra System

knitr::opts_chunk$set(echo = TRUE)  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  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 = "")  ## Dynamic linear model Augument theAR(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

K2s <- L2s %*% t(L2s)
V2s <- solve(K2s)
# Try simplify; causes timeout on CRAN Fedora, hence in try() call.
# Else, just use
# V2s <- simplify(solve(K2s))
try(V2s <- simplify(V2s), silent = TRUE)

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


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


## Try the Ryacas package in your browser

Any scripts or data that you put into this service are public.

Ryacas documentation built on Jan. 17, 2023, 1:11 a.m.