-
Notifications
You must be signed in to change notification settings - Fork 7
/
ssm-matrix.Rmd
executable file
·126 lines (104 loc) · 2.64 KB
/
ssm-matrix.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
---
title: "The structure of the concentration and covariance matrix in a simple state-space model"
author: "Mikkel Meyer Andersen and Søren Højsgaard"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{The structure of the concentration and covariance matrix in a simple state-space model}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r, message=FALSE}
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
```{r}
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'$
```{r}
K1s <- L1s %*% t(L1s)
V1s <- solve(K1s)
```
```{r, results="asis"}
cat(
"\\begin{align} K_1 &= ", tex(K1s), " \\\\
V_1 &= ", tex(V1s), " \\end{align}", sep = "")
```
## Dynamic linear model
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
```{r}
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
```
```{r}
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)
```
```{r, results="asis"}
cat(
"\\begin{align} K_2 &= ", tex(K2s), " \\\\
V_2 &= ", tex(V2s), " \\end{align}", sep = "")
```
## Numerical evalation in R
```{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:
```{r}
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)
```