---
title: "21 - Concentration and covariance matrix in an autoregressive model and in a dynamic linear model"
author: Mikkel Meyer Andersen and Søren Højsgaard
date: "`r date()`"
output: 
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
    number_sections: true
vignette: >
  %\VignetteIndexEntry{21 - Concentration and covariance matrix in an autoregressive model and in a dynamic linear model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r, message=FALSE}
library(caracas)
```

```{r, include = FALSE}
inline_code <- function(x) {
  x
}

if (!has_sympy()) {
  # SymPy not available, so the chunks shall not be evaluated
  knitr::opts_chunk$set(eval = FALSE)
  
  inline_code <- function(x) {
    deparse(substitute(x))
  }
}
```


## Autoregressive model ($AR(1)$)

```{r ar1, echo=F}
N <- 3
L1 <- diag(4)
L1[cbind(1 + (1:N), 1:N)] <- "-a"
L1 <- as_sym(L1)
```

```{r, echo=F}
e <- as_sym(paste0("e", 0:3))
x <- as_sym(paste0("x", 0:3))
u <- as_sym(paste0("u", 1:3))
y <- as_sym(paste0("y", 1:3))
eu <- c(e, u)
xy <- c(x, y)
```


Consider this model: 
$$
x_i = a x_{i-1} + e_i, \quad i=1, \dots, 3
$$ 
and $x_0=e_0$. All terms $e_0, \dots, e_3$ are independent and $N(0,v^2)$ distributed. 
Let $e=(e_0, \dots, e_3)$ and $x=(x_0, \dots x_3)$. Hence $e \sim N(0, v^2 I)$.
Isolating error terms gives
$$
e= `r inline_code(tex(e))` = `r inline_code(tex(L1))` `r inline_code(tex(x))` = L_1 x
$$

```{r}
<<ar1>>
```

Since $\mathbf{Var}(e)=v^2 I$ we have 
$\mathbf{Var}(e)=v^2 I=L \mathbf{Var}(x) L'$ so the covariance matrix of $x$ is $V_1=\mathbf{Var}(x) = v^2 L^- (L^-)'$ 
 while the concentration matrix (the inverse covariances matrix) is $K=v^{-2}L' L$.

```{r}
def_sym(v2)
L1inv <- inv(L1)
V1 <- v2 * L1inv %*% t(L1inv)
K1 <- (t(L1) %*% L1) / v2
```

```{r, echo=F, results="asis"}
cat(
  "\\begin{align} 
    K_1 &= ", tex(K1), " \\\\ 
   V_1 &= ", tex(V1), " 
  \\end{align}", sep = "")
```


## Dynamic linear model

```{r L2, echo=F}
N <- 3
L2 <- diag("1", 1 + 2*N)
L2[cbind(1 + (1:N), 1:N)] <- "-a"
L2[cbind(1 + N + (1:N), 1 + 1:N)] <- "-b"
L2 <- as_sym(L2)
```

Augment the $AR(1)$ process above with $y_i=b x_i + u_i$ for
$i=1,2,3$.  Suppose $u_i\sim N(0, w^2)$ and all $u_i$ are independent
and inpendent of $e$.
Then
$(e,u)$ can be expressed in terms of $(x,y)$ as
$$
(e,u) = `r inline_code(tex(eu))` = `r inline_code(tex(L2))` `r inline_code(tex(xy))` = L_2 (x,y)
$$
where
```{r}
<<L2>>
```


```{r}
Veu <- diag(1, 7)
diag(Veu)[1:4] <- "v2"
diag(Veu)[5:7] <- "w2"
Veu
Veu <- as_sym(Veu)
Veu
L2inv <- inv(L2) 
V2 <- L2inv %*% Veu %*% t(L2inv) 
K2 <- t(L2) %*% inv(Veu) %*% L2
```

```{r, results="asis", echo=F}
cat(
  "\\begin{align} K_2 &= ", tex(K2), " \\\\ 
                  V_2 &= ", tex(V2), " \\end{align}", sep = "")
```