library(Ryacas)
Consider this model: xi=ax0+ei,i=1,…,4 and x0=e0. All terms e0,…,e3 are independent and N(0,1) distributed. Let e=(e0,…,e3) and x=(x0,…x3). Isolating error terms gives that e=L1x where L1 has the form
<- diag(4)
L1chr 2:4, 1] <- "-a"
L1chr[<- ysym(L1chr)
L1 L1
## {{ 1, 0, 0, 0},
## {-a, 1, 0, 0},
## {-a, 0, 1, 0},
## {-a, 0, 0, 1}}
If error terms have variance 1 then Var(e)=LVar(x)L′ so the covariance matrix is V1=Var(x)=L−(L−)′ while the concentration matrix (the inverse covariances matrix) is K=L′L.
<- solve(L1)
L1inv <- t(L1) %*% L1
K1 <- L1inv %*% t(L1inv) V1
cat(
"\\begin{align}
K_1 &= ", tex(K1), " \\\\
V_1 &= ", tex(V1), "
\\end{align}", sep = "")
K1=(3a2+1−a−a−a−a100−a010−a001)V1=(1aaaaa2+1a2a2aa2a2+1a2aa2a2a2+1)
Slightly more elaborate:
<- diag(4)
L2chr 2:4, 1] <- c("-a1", "-a2", "-a3")
L2chr[<- ysym(L2chr)
L2 L2
## {{ 1, 0, 0, 0},
## {-a1, 1, 0, 0},
## {-a2, 0, 1, 0},
## {-a3, 0, 0, 1}}
<- diag(4)
Vechr cbind(1:4, 1:4)] <- c("w1", "w2", "w2", "w2")
Vechr[<- ysym(Vechr)
Ve Ve
## {{w1, 0, 0, 0},
## { 0, w2, 0, 0},
## { 0, 0, w2, 0},
## { 0, 0, 0, w2}}
<- solve(L2)
L2inv <- t(L2) %*% solve(Ve) %*% L2
K2 <- L2inv %*% Ve %*% t(L2inv) V2
cat(
"\\begin{align}
K_2 &= ", tex(K2), " \\\\
V_2 &= ", tex(V2), "
\\end{align}", sep = "")
K2=(1w1+a21w2+a22w2+a23w2−a1w2−a2w2−a3w2−a1w21w200−a2w201w20−a3w2001w2)V2=(w1w1a1w1a2w1a3a1w1w1a21+w2a1w1a2a1w1a3a2w1a2w1a1w1a22+w2a2w1a3a3w1a3w1a1a3w1a2w1a23+w2)