library(Ryacas0)
library(Matrix)
Set output width:
get_output_width()
## [1] 120
set_output_width(120)
get_output_width()
## [1] 120
Consider AR(1) process: xi=axi−1+ei where i=1,2,3 and where x0=e0. Isolating error terms gives that e=L1x where e=(e0,…,e3) and x=(x0,…x3) and where L1 has the form
<- 3
N <- diag("1", 1 + N)
L1chr cbind(1+(1:N), 1:N)] <- "-a"
L1chr[<- as.Sym(L1chr)
L1s 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 %*%
<- Simplify(L1s * Transpose(L1s))
K1s <- Simplify(Inverse(K1s)) V1s
cat(
"\\begin{align} K_1 &= ", TeXForm(K1s), " \\\\
V_1 &= ", TeXForm(V1s), " \\end{align}", sep = "")
K1=(1−a00−aa2+1−a00−aa2+1−a00−aa2+1)V1=(a6+a4+a2+1a(a4+a2+1)a2(a2+1)a3a(a4+a2+1)a4+a2+1a(a2+1)a2a2(a2+1)a(a2+1)a2+1aa3a2a1)
Augument the AR(1) process above with yi=bxi+ui. Then (e,u) can be expressed in terms of (x,y) as (e,u)=L2(x,y) where
<- 3
N <- diag("1", 1 + 2*N)
L2chr cbind(1+(1:N), 1:N)] <- "-a"
L2chr[cbind(1 + N + (1:N), 1 + 1:N)] <- "-b"
L2chr[<- as.Sym(L2chr)
L2s 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
<- Simplify(L2s * Transpose(L2s))
K2s <- Simplify(Inverse(K2s)) V2s
cat(
"\\begin{align} K_2 &= ", TeXForm(K2s), " \\\\
V_2 &= ", TeXForm(V2s), " \\end{align}", sep = "")
K2=(1−a00000−aa2+1−a0−b000−aa2+1−aab−b000−aa2+10ab−b0−bba0b2+10000−bba0b2+10000−b00b2+1)V2=(a6b2+a6+a4b2+a4+a2b2+a2+1a5b2+a5+a3b2+a3+ab2+aa4b2+a4+a2b2+a2a3b2+a3aba2ba3ba5b2+a5+a3b2+a3+ab2+aa4b2+a4+a2b2+a2+b2+1a3b2+a3+ab2+aa2b2+a2baba2ba4b2+a4+a2b2+a2a3b2+a3+ab2+aa2b2+a2+b2+1ab2+a0baba3b2+a3a2b2+a2ab2+ab2+100bbab00100ba2bab0010ba3ba2bab001)
<- function(x) {
sparsify ::Matrix(x, sparse = TRUE)
Matrix
}
<- 0.5
alpha <- -0.3
beta
## AR(1)
<- 3
N <- diag(1, 1 + N)
L1 cbind(1+(1:N), 1:N)] <- -alpha
L1[<- L1 %*% t(L1)
K1 <- solve(K1)
V1 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
<- 3
N <- diag(1, 1 + 2*N)
L2 cbind(1+(1:N), 1:N)] <- -alpha
L2[cbind(1 + N + (1:N), 1 + 1:N)] <- -beta
L2[<- L2 %*% t(L2)
K2 <- solve(K2)
V2 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:
<- Eval(V1s, list(a = alpha))
V1s_eval <- Eval(V2s, list(a = alpha, b = beta))
V2s_eval
all.equal(V1, V1s_eval)
## [1] TRUE
all.equal(V2, V2s_eval)
## [1] TRUE