library(Ryacas0)
library(Matrix)Set output width:
get_output_width()## [1] 120set_output_width(120)
get_output_width()## [1] 120Consider \(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## 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 1If 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 = "")\[\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 ^{6} + a ^{4} + a ^{2} + 1 & a \left( a ^{4} + a ^{2} + 1\right) & a ^{2} \left( a ^{2} + 1\right) & a ^{3} \\ a \left( a ^{4} + a ^{2} + 1\right) & a ^{4} + a ^{2} + 1 & a \left( a ^{2} + 1\right) & a ^{2} \\ a ^{2} \left( a ^{2} + 1\right) & a \left( a ^{2} + 1\right) & a ^{2} + 1 & a \\ a ^{3} & a ^{2} & a & 1 \end{array} \right) \end{align}\]
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## 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 1K2s <- Simplify(L2s * Transpose(L2s))
V2s <- Simplify(Inverse(K2s))cat(
"\\begin{align} K_2 &= ", TeXForm(K2s), " \\\\
V_2 &= ", TeXForm(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 & a ^{2} b & a ^{3} b \\ 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 & a ^{2} b \\ 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}\]
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.25sparsify(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.09sparsify(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.0000Comparing 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] TRUEall.equal(V2, V2s_eval)## [1] TRUE