knitr::opts_chunk$set(echo = TRUE)
library(Ryacas) library(Matrix)
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 = "")
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
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 = "")
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.