---
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.