---
output:
 # html_document
  pdf_document
title: "Ch 7: Checking observability of state-space models with TVPs"
author: "Sy-Miin Chow"
date: "June 4, 2021"
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, message = FALSE)
```

```{r, echo = TRUE}
#Set up display options to suppress scientific notations
rm(list=ls())
options(scipen=999)
```

## Background
The goal of this script is to provide several examples of checking the observability of state-space models with time-varying parameters (TVPs). More details can be found in Chapter 7 of Gates, K., Chow, S-M., & Molenaar, P.C.M. Analysis of Intraindividual Variation. Systems Approaches to Human Process Analysis. New York, NY: Taylor & Francis.

## Define a function to check observability
```{r, echo = TRUE}
#Function to get observability matrix - see Equation 7.10 in Chapter 7
#Note that Wikipedia defines O as the transpose of Eq 7.10:
getObservabilityMatrix = function(Jacob_h,Jacob_g){
ne = dim(Jacob_g)[1]
ny = dim(Jacob_h)[1]
J = t(Jacob_h)
O_recursive = matrix(NA,ne,ny*ne)
for (i in 0:(ne-1)){
  O_recursive[,(1+i*(ny)):((i+1)*(ny))] = J
  J = t(Jacob_g)%*%J
}
return(O_recursive)
}

```

## TVP Example 1
This is an example motivated by Molenaar, de Gooijer, & Schmitz, 1992. 
In this examples, we have 3 manifest variables, one common factor, and time-varying factor loadings:
$$y(t) = \Lambda(t)\eta(t) + \epsilon(t); y(t) = [y_1(t), y_2(t), y_3(t)]'$$
$$\Lambda(t) = [l_{11}(t),l_{21}(t),l_{31}(t)]'$$
The TV factor loadings follow the random walk model as:
$$\Lambda(t) = \Lambda(t-1) + \zeta_{\Lambda}(t)$$
The expanded form of this model is obtained by inserting $\Lambda(t)$ into: 
$$\eta^*(t) = [\eta(t), l_{11}(t),  l_{21}(t), l_{31}(t)]'.$$
Let:

$y_1(t)  = l_{11}(t)*\eta(t)$ be defined as $h_1$,

$y_2(t)  = l_{21}(t)*\eta(t)$ be defined as $h_2$,

$y_3(t)  = l_{31}(t)*\eta(t)$ be defined as $h_3$. 

At the dynamic level, let $g = [g_1, g_2, g_3, g_4]'$ be a vector of nonlinear dynamic functions, with:

$\eta(t) = \alpha_1*\eta(t-1) + \zeta_{\eta}(t)$ be defined as $g_1$,

$l_{11}(t) = l_{11}(t-1) + \zeta_{l_{11}}(t)$ be defined as $g_2$,

$l_{21}(t) = l_{21}(t-1) + \zeta_{l_{21}}(t)$ be defined as $g_3$, and

$l_{31}(t) = l_{31}(t-1) + \zeta_{l_{31}}(t)$ be defined as $g_4$.

$Jacob_g = \frac{\partial g}{\partial \eta^*}$ is the Jacobian matrix where the ($j$,$k$)th element of $Jacob_g$ carries partial derivative of the
$j$th dynamic function with respect to the $k$th latent variable
$Jacob_h = \frac{\partial h}{\partial \eta^*}$ is the Jacobian matrix where the ($j$,$k$)th element of $Jacob_h$ carries partial derivative of the
$j$th measurement function with respect to the $k$th latent variable.

To check the observability of this model, we substitute arbitrary numerical values into all the unknown elements, and then check the rank of the resulting observability matrix.
```{r, echo = TRUE}
alpha_1 = .8; eta_t = 2; l_11_t = 1; l_21_t = .9; l_31_t = 1.2
ne = 4 #length of eta*(t)
Jacob_g = matrix(c(alpha_1, 0, 0, 0,
                   0,1,0,0,
                   0,0,1,0,
                   0,0,0,1),ncol=ne,byrow=TRUE)
Jacob_h = matrix(c(l_11_t,eta_t, 0, 0,
                   l_21_t,0,eta_t,0,
                   l_31_t,0,0,eta_t),ncol=ne,byrow=TRUE)

# Call the getObservabilityMatrix function to automatically construct the observability matrix;
# here is one way to check the rank of the matrix.
qr(getObservabilityMatrix(Jacob_h,Jacob_g))$rank 

```
The rank of the returned observability matrix is equal to $ne$, which suggests that the observability matrix is of full rank. This common factor model with TV factor loadings is observable!

## TVP Example 2 - Expanding Example 1 to a two-factor model
In this example, we expanded Example 1 to 6 manifest variables and two latent factors with time-varying factor loadings:
$$y(t) = \Lambda(t)\eta(t) + \epsilon(t); y(t) = [y_1(t), y_2(t), y_3(t), \ldots, y_6(t)]'$$
\begin{eqnarray}
\Lambda(t) = \begin{bmatrix}
l_{11}(t) & 0\\
l_{21}(t)& 0\\
l_{31}(t)& 0\\
0  &     l_{42}(t)\\
0  &     l_{52}(t)\\
0  &     l_{62}(t)
\end{bmatrix}
\end{eqnarray}

The expanded form of this model is obtained by inserting $\Lambda(t)$ into: 
$$\eta^*(t) = [\eta_1(t), \eta_2(t), l_{11}(t),  l_{21}(t), l_{31}(t),\ldots,l_{62}(t)]',$$
and letting the TV factor loadings follow the random walk model as: 
$$y_1(t)  = l_{11}(t)*\eta_1(t) \text{ be defined as }h_1$$,
$$y_2(t)  = l_{21}(t)*\eta_1(t) \text{ be defined as }h_2$$,
$$y_3(t)  = l_{31}(t)*\eta_1(t) \text{ be defined as }h_3$$, 
$$y_4(t)  = l_{42}(t)*\eta_2(t) \text{ be defined as }h_4$$. 
$$y_5(t)  = l_{52}(t)*\eta_2(t) \text{ be defined as }h_5$$. 
$$y_6(t)  = l_{62}(t)*\eta_2(t) \text{ be defined as }h_6$$. 

At the dynamic level, let $g = [g_1, g_2, g_3, g_4, g_5, g_6, g_7, g_8]'$ be a vector of nonlinear dynamic functions, with:

$$\eta_1(t) = \alpha_{11}*\eta_1(t-1) + \alpha_{12}*\eta_2(t-1) + \zeta_{\eta_1}(t), \text{ defined as }g_1$$
$$\eta_2(t) = \alpha_{22}*\eta_2(t-1) + \alpha_{21}*\eta_1(t-1) + \zeta_{\eta_2}(t) \text{, defined as } g_2$$

$$l_{11}(t) = l_{11}(t-1) + \zeta_{l_{11}}(t) \text{ defined as }g_3$$,

$$l_{21}(t) = l_{21}(t-1) + \zeta_{l_{21}}(t) \text{ defined as }g_4$$,

$$l_{31}(t) = l_{31}(t-1) + \zeta_{l_{31}}(t)\text{ defined as }g_5$$,

$$\vdots$$

$$l_{62}(t) = l_{62}(t-1) + \zeta_{l_{62}}(t)\text{ be defined as }g_8$$.

$Jacob_g = \frac{\partial g}{\partial \eta^*}$ is the Jacobian matrix where the ($j$,$k$)th element of $Jacob_g$ carries partial derivative of the
$j$th dynamic function with respect to the $k$th latent variable
$Jacob_h = \frac{\partial h}{\partial \eta^*}$ is the Jacobian matrix where the ($j$,$k$)th element of $Jacob_h$ carries partial derivative of the
$j$th measurement function with respect to the $k$th latent variable.

To check the observability of this model, we substitute arbitrary numerical values into all the unknown elements, and then check the rank of the resulting observability matrix.
```{r, echo = TRUE}
alpha_11 = .8; alpha_12 = .1; alpha_21 = -.3; alpha_22 = .7

eta1_t = 2; eta2_t = -3
l_11_t = 1; l_21_t = .9; l_31_t = 1.2
l_42_t = 1; l_52_t = .8; l_62_t = 0.7

ne = 8 # length of expanded eta*(t)
Jacob_g = matrix(c(alpha_11, alpha_12, rep(0,6),
                   alpha_21, alpha_22, rep(0,6),
                   0,0,1,0,0,0,0,0,
                   0,0,0,1,0,0,0,0,
                   0,0,0,0,1,0,0,0,
                   0,0,0,0,0,1,0,0,
                   0,0,0,0,0,0,1,0,
                   0,0,0,0,0,0,0,1),ncol=ne,byrow=TRUE)
Jacob_h = matrix(c(l_11_t,0,eta1_t,0,0,0,0,0,
                   l_21_t,0,0,eta1_t,0,0,0,0,
                   l_31_t,0,0,0,eta1_t,0,0,0,
                   0,l_42_t,0,0,0,eta2_t,0,0,
                   0,l_52_t,0,0,0,0,eta2_t,0,
                   0,l_62_t,0,0,0,0,0,eta2_t                   
                   ),ncol=ne,byrow=TRUE)

qr(getObservabilityMatrix(Jacob_h,Jacob_g))$rank

```
The rank of the returned observability matrix is equal to $ne$. So this model is observable too.


## TVP Example 3 - Two-factor model, TV auto- and cross-regression weights
In this example, we have a two-factor model with time-invariant factor loadings, but TV auto- and cross-regression weights at the factor level as:
$$y(t) = \Lambda\eta(t) + \epsilon(t); y(t) = [y_1(t), y_2(t), y_3(t), \ldots, y_6(t)]'$$
\begin{eqnarray}
\Lambda(t) = \begin{bmatrix}
l_{11} & 0\\
l_{21}& 0\\
l_{31}& 0\\
0  &     l_{42}\\
0  &     l_{52}\\
0  &     l_{62}
\end{bmatrix}
\end{eqnarray}

At the dynamic level, let $g = [g_1, g_2, g_3, g_4, g_5, g_6]'$ be a vector of nonlinear dynamic functions, with:

$\eta_1(t) = \alpha_{11}(t-1)\eta_1(t-1) + \alpha_{12}(t-1)\eta_2(t-1) + \zeta_{\eta_1}(t)$, defined as $g_1$

$\eta_2(t) = \alpha_{22}(t-1)\eta_2(t-1) + \alpha_{21}(t-1)\eta_1(t-1) + \zeta_{\eta_2}(t)$, defined as $g_2$

$\alpha_{11}(t) = \alpha_{11}(t-1) + \zeta_{\alpha_{11}}(t)$, defined as $g_3$,

$\alpha_{12}(t) = \alpha_{12}(t-1) + \zeta_{\alpha_{12}}(t)$, defined as $g_4$,

$\alpha_{21}(t) = \alpha_{21}(t-1) + \zeta_{\alpha_{21}}(t)$, defined as $g_5$,

$\alpha_{22}(t) = \alpha_{22}(t-1) + \zeta_{\alpha_{22}}(t)$, defined as $g_6$,

As before, $Jacob_g = \frac{\partial g}{\partial \eta^*}$ is the Jacobian matrix where the ($j$,$k$)th element of $Jacob_g$ carries partial derivative of the
$j$th dynamic function with respect to the $k$th latent variable; 
$Jacob_h = \frac{\partial h}{\partial \eta^*}$ is the Jacobian matrix where the ($j$,$k$)th element of $Jacob_h$ carries partial derivative of the
$j$th measurement function with respect to the $k$th latent variable.

To check the observability of this model, we substitute arbitrary numerical values into all the unknown elements, and then check the rank of the resulting observability matrix.
```{r, echo = TRUE}
a_11_t = .8; a_12_t = .1; a_21_t = -.3; a_22_t = .7

eta1_t = 2; eta2_t = -3
l_11 = 1; l_21 = .9; l_31 = 1.2
l_42 = 1; l_52 = .8; l_62 = 0.7

ne = 6 # length of expanded eta*(t)
Jacob_g = matrix(c(a_11_t, a_12_t, rep(0,4),
                   a_21_t, a_22_t, rep(0,4),
                   0,0,1,0,0,0,
                   0,0,0,1,0,0,
                   0,0,0,0,1,0,
                   0,0,0,0,0,1),ncol=ne,byrow=TRUE)
Jacob_h = matrix(c(l_11,0,rep(0,4),
                   l_21,0,rep(0,4),
                   l_31,0,rep(0,4),
                   0,l_42,rep(0,4),
                   0,l_52,rep(0,4),                   
                   0,l_62,rep(0,4)
                   ),ncol=ne,byrow=TRUE)

qr(getObservabilityMatrix(Jacob_h,Jacob_g))$rank

```
This model is not observable.


## TVP Example 4 - Two-factor model, TV cross-regression weights
Seeing that the model in Example 3 is not observable, we modify it so only the two cross-regression parameters were allowed to be TVPs:
$$\eta_1(t) = \alpha_{11}\eta_1(t-1) + \alpha_{12}(t-1)\eta_2(t-1) + \zeta_{\eta_1}(t), \text{ defined as }g_1$$

$$\eta_2(t) = \alpha_{22}\eta_2(t-1) + \alpha_{21}(t-1)\eta_1(t-1) + \zeta_{\eta_2}(t), \text{ defined as } g_2$$

$$\alpha_{12}(t) = \alpha_{12}(t-1) + \zeta_{\alpha_{12}}(t), \text{ defined as }g_3$$,

$$\alpha_{21}(t) = \alpha_{21}(t-1) + \zeta_{\alpha_{21}}(t), \text{ defined as }g_4$$.

Now let's repeat the same process.
```{r, echo = TRUE}
a_11 = .8; a_12_t = .1; a_21_t = -.3; a_22 = .7

eta1_t = 2; eta2_t = -3
l_11 = 1; l_21 = .9; l_31 = 1.2
l_42 = 1; l_52 = .8; l_62 = 0.7

ne = 4 # length of expanded eta*(t)
Jacob_g = matrix(c(a_11, a_12_t, eta1_t, 0,
                   a_21_t, a_22, 0, eta2_t,
                   0,0,1,0,
                   0,0,0,1),ncol=ne,byrow=TRUE)
Jacob_h = matrix(c(l_11,0,rep(0,2),
                   l_21,0,rep(0,2),
                   l_31,0,rep(0,2),
                   0,l_42,rep(0,2),
                   0,l_52,rep(0,2),                   
                   0,l_62,rep(0,2)
),ncol=ne,byrow=TRUE)

qr(getObservabilityMatrix(Jacob_h,Jacob_g))$rank

```
This streamlined model is now observable.
