The State Space Model

Ivan Jacob Agaloos Pesigan

2025-02-14

Model

The measurement model is given by $$\begin{equation} \mathbf{y}_{i, t} = \boldsymbol{\nu} + \boldsymbol{\Lambda} \boldsymbol{\eta}_{i, t} + \boldsymbol{\varepsilon}_{i, t}, \quad \mathrm{with} \quad \boldsymbol{\varepsilon}_{i, t} \sim \mathcal{N} \left( \mathbf{0}, \boldsymbol{\Theta} \right) \end{equation}$$ where $\mathbf{y}_{i, t}$, $\boldsymbol{\eta}_{i, t}$, and $\boldsymbol{\varepsilon}_{i, t}$ are random variables and $\boldsymbol{\nu}$, $\boldsymbol{\Lambda}$, and $\boldsymbol{\Theta}$ are model parameters. $\mathbf{y}_{i, t}$ represents a vector of observed random variables, $\boldsymbol{\eta}_{i, t}$ a vector of latent random variables, and $\boldsymbol{\varepsilon}_{i, t}$ a vector of random measurement errors, at time $t$ and individual $i$. $\boldsymbol{\nu}$ denotes a vector of intercepts, $\boldsymbol{\Lambda}$ a matrix of factor loadings, and $\boldsymbol{\Theta}$ the covariance matrix of $\boldsymbol{\varepsilon}$.

An alternative representation of the measurement error is given by $$\begin{equation} \boldsymbol{\varepsilon}_{i, t} = \boldsymbol{\Theta}^{\frac{1}{2}} \mathbf{z}_{i, t}, \quad \mathrm{with} \quad \mathbf{z}_{i, t} \sim \mathcal{N} \left( \mathbf{0}, \mathbf{I} \right) \end{equation}$$ where $\mathbf{z}_{i, t}$ is a vector of independent standard normal random variables and $\left( \boldsymbol{\Theta}^{\frac{1}{2}} \right) \left( \boldsymbol{\Theta}^{\frac{1}{2}} \right)^{\prime} = \boldsymbol{\Theta}$ .

The dynamic structure is given by $$\begin{equation} \boldsymbol{\eta}_{i, t} = \boldsymbol{\alpha} + \boldsymbol{\beta} \boldsymbol{\eta}_{i, t - 1} + \boldsymbol{\zeta}_{i, t}, \quad \mathrm{with} \quad \boldsymbol{\zeta}_{i, t} \sim \mathcal{N} \left( \mathbf{0}, \boldsymbol{\Psi} \right) \end{equation}$$ where $\boldsymbol{\eta}_{i, t}$, $\boldsymbol{\eta}_{i, t - 1}$, and $\boldsymbol{\zeta}_{i, t}$ are random variables, and $\boldsymbol{\alpha}$, $\boldsymbol{\beta}$, and $\boldsymbol{\Psi}$ are model parameters. Here, $\boldsymbol{\eta}_{i, t}$ is a vector of latent variables at time $t$ and individual $i$, $\boldsymbol{\eta}_{i, t - 1}$ represents a vector of latent variables at time $t - 1$ and individual $i$, and $\boldsymbol{\zeta}_{i, t}$ represents a vector of dynamic noise at time $t$ and individual $i$. $\boldsymbol{\alpha}$ denotes a vector of intercepts, $\boldsymbol{\beta}$ a matrix of autoregression and cross regression coefficients, and $\boldsymbol{\Psi}$ the covariance matrix of $\boldsymbol{\zeta}_{i, t}$.

An alternative representation of the dynamic noise is given by $$\begin{equation} \boldsymbol{\zeta}_{i, t} = \boldsymbol{\Psi}^{\frac{1}{2}} \mathbf{z}_{i, t}, \quad \mathrm{with} \quad \mathbf{z}_{i, t} \sim \mathcal{N} \left( \mathbf{0}, \mathbf{I} \right) \end{equation}$$ where $\left( \boldsymbol{\Psi}^{\frac{1}{2}} \right) \left( \boldsymbol{\Psi}^{\frac{1}{2}} \right)^{\prime} = \boldsymbol{\Psi}$ .

Data Generation

Notation

Let $t = 100$ be the number of time points and $n = 5$ be the number of individuals.

Let the measurement model intecept vector $\boldsymbol{\nu}$ be given by

$$\begin{equation} \boldsymbol{\nu} = \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \end{array} \right) . \end{equation}$$

Let the factor loadings matrix $\boldsymbol{\Lambda}$ be given by

$$\begin{equation} \boldsymbol{\Lambda} = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) . \end{equation}$$

Let the measurement error covariance matrix $\boldsymbol{\Theta}$ be given by

$$\begin{equation} \boldsymbol{\Theta} = \left( \begin{array}{ccc} 0.2 & 0 & 0 \\ 0 & 0.2 & 0 \\ 0 & 0 & 0.2 \\ \end{array} \right) . \end{equation}$$

Let the initial condition $\boldsymbol{\eta}_{0}$ be given by

$$\begin{equation} \boldsymbol{\eta}_{0} \sim \mathcal{N} \left( \boldsymbol{\mu}_{\boldsymbol{\eta} \mid 0}, \boldsymbol{\Sigma}_{\boldsymbol{\eta} \mid 0} \right) \end{equation}$$

$$\begin{equation} \boldsymbol{\mu}_{\boldsymbol{\eta} \mid 0} = \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \end{array} \right) \end{equation}$$

$$\begin{equation} \boldsymbol{\Sigma}_{\boldsymbol{\eta} \mid 0} = \left( \begin{array}{ccc} 1 & 0.2 & 0.2 \\ 0.2 & 1 & 0.2 \\ 0.2 & 0.2 & 1 \\ \end{array} \right) . \end{equation}$$

Let the constant vector $\boldsymbol{\alpha}$ be given by

$$\begin{equation} \boldsymbol{\alpha} = \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \end{array} \right) . \end{equation}$$

Let the transition matrix $\boldsymbol{\beta}$ be given by

$$\begin{equation} \boldsymbol{\beta} = \left( \begin{array}{ccc} 0.7 & 0 & 0 \\ 0.5 & 0.6 & 0 \\ -0.1 & 0.4 & 0.5 \\ \end{array} \right) . \end{equation}$$

Let the dynamic process noise $\boldsymbol{\Psi}$ be given by

$$\begin{equation} \boldsymbol{\Psi} = \left( \begin{array}{ccc} 0.1 & 0 & 0 \\ 0 & 0.1 & 0 \\ 0 & 0 & 0.1 \\ \end{array} \right) . \end{equation}$$

R Function Arguments

n
#> [1] 5
time
#> [1] 100
mu0
#> [1] 0 0 0
sigma0
#>      [,1] [,2] [,3]
#> [1,]  1.0  0.2  0.2
#> [2,]  0.2  1.0  0.2
#> [3,]  0.2  0.2  1.0
sigma0_l # sigma0_l <- t(chol(sigma0))
#>      [,1]      [,2]      [,3]
#> [1,]  1.0 0.0000000 0.0000000
#> [2,]  0.2 0.9797959 0.0000000
#> [3,]  0.2 0.1632993 0.9660918
alpha
#> [1] 0 0 0
beta
#>      [,1] [,2] [,3]
#> [1,]  0.7  0.0  0.0
#> [2,]  0.5  0.6  0.0
#> [3,] -0.1  0.4  0.5
psi
#>      [,1] [,2] [,3]
#> [1,]  0.1  0.0  0.0
#> [2,]  0.0  0.1  0.0
#> [3,]  0.0  0.0  0.1
psi_l # psi_l <- t(chol(psi))
#>           [,1]      [,2]      [,3]
#> [1,] 0.3162278 0.0000000 0.0000000
#> [2,] 0.0000000 0.3162278 0.0000000
#> [3,] 0.0000000 0.0000000 0.3162278
nu
#> [1] 0 0 0
lambda
#>      [,1] [,2] [,3]
#> [1,]    1    0    0
#> [2,]    0    1    0
#> [3,]    0    0    1
theta
#>      [,1] [,2] [,3]
#> [1,]  0.2  0.0  0.0
#> [2,]  0.0  0.2  0.0
#> [3,]  0.0  0.0  0.2
theta_l # theta_l <- t(chol(theta))
#>           [,1]      [,2]      [,3]
#> [1,] 0.4472136 0.0000000 0.0000000
#> [2,] 0.0000000 0.4472136 0.0000000
#> [3,] 0.0000000 0.0000000 0.4472136

Visualizing the Dynamics Without Process Noise (n = 5 with Different Initial Condition)

Using the SimSSMFixed Function from the simStateSpace Package to Simulate Data

library(simStateSpace)
sim <- SimSSMFixed(
  n = n,
  time = time,
  mu0 = mu0,
  sigma0_l = sigma0_l,
  alpha = alpha,
  beta = beta,
  psi_l = psi_l,
  nu = nu,
  lambda = lambda,
  theta_l = theta_l,
  type = 0
)
data <- as.data.frame(sim)
head(data)
#>   id time          y1          y2         y3
#> 1  1    0 -0.68686529 -0.23269186  0.5864176
#> 2  1    1 -0.24266216 -0.01684359 -0.2646096
#> 3  1    2  0.92818305 -0.05859969 -1.0290302
#> 4  1    3  0.03836036  0.57871750 -0.2909122
#> 5  1    4  0.14986876  0.76497073  0.8038175
#> 6  1    5 -0.10242480  0.63283996  0.2144810
summary(data)
#>        id         time             y1                 y2           
#>  Min.   :1   Min.   : 0.00   Min.   :-2.20038   Min.   :-2.343454  
#>  1st Qu.:2   1st Qu.:24.75   1st Qu.:-0.45945   1st Qu.:-0.487994  
#>  Median :3   Median :49.50   Median :-0.03995   Median :-0.003589  
#>  Mean   :3   Mean   :49.50   Mean   :-0.04093   Mean   :-0.016301  
#>  3rd Qu.:4   3rd Qu.:74.25   3rd Qu.: 0.38653   3rd Qu.: 0.496862  
#>  Max.   :5   Max.   :99.00   Max.   : 1.44938   Max.   : 2.361441  
#>        y3          
#>  Min.   :-2.03584  
#>  1st Qu.:-0.39631  
#>  Median : 0.02100  
#>  Mean   : 0.02513  
#>  3rd Qu.: 0.44660  
#>  Max.   : 2.81589
plot(sim)

Model Fitting

Prepare Data

dynr_data <- dynr::dynr.data(
  data = data,
  id = "id",
  time = "time",
  observed = c("y1", "y2", "y3")
)

Prepare Initial Condition

dynr_initial <- dynr::prep.initial(
  values.inistate = mu0,
  params.inistate = c("mu0_1_1", "mu0_2_1", "mu0_3_1"),
  values.inicov = sigma0,
  params.inicov = matrix(
    data = c(
      "sigma0_1_1", "sigma0_2_1", "sigma0_3_1",
      "sigma0_2_1", "sigma0_2_2", "sigma0_3_2",
      "sigma0_3_1", "sigma0_3_2", "sigma0_3_3"
    ),
    nrow = 3
  )
)

Prepare Measurement Model

dynr_measurement <- dynr::prep.measurement(
  values.load = diag(3),
  params.load = matrix(data = "fixed", nrow = 3, ncol = 3),
  state.names = c("eta_1", "eta_2", "eta_3"),
  obs.names = c("y1", "y2", "y3")
)

Prepare Dynamic Process

dynr_dynamics <- dynr::prep.formulaDynamics(
  formula = list(
    eta_1 ~ alpha_1_1 * 1 + beta_1_1 * eta_1 + beta_1_2 * eta_2 + beta_1_3 * eta_3,
    eta_2 ~ alpha_2_1 * 1 + beta_2_1 * eta_1 + beta_2_2 * eta_2 + beta_2_3 * eta_3,
    eta_3 ~ alpha_3_1 * 1 + beta_3_1 * eta_1 + beta_3_2 * eta_2 + beta_3_3 * eta_3
  ),
  startval = c(
    alpha_1_1 = alpha[1], alpha_2_1 = alpha[2], alpha_3_1 = alpha[3],
    beta_1_1 = beta[1, 1], beta_1_2 = beta[1, 2], beta_1_3 = beta[1, 3],
    beta_2_1 = beta[2, 1], beta_2_2 = beta[2, 2], beta_2_3 = beta[2, 3],
    beta_3_1 = beta[3, 1], beta_3_2 = beta[3, 2], beta_3_3 = beta[3, 3]
  ),
  isContinuousTime = FALSE
)

Prepare Process Noise

dynr_noise <- dynr::prep.noise(
  values.latent = psi,
  params.latent = matrix(
    data = c(
      "psi_1_1", "psi_2_1", "psi_3_1",
      "psi_2_1", "psi_2_2", "psi_3_2",
      "psi_3_1", "psi_3_2", "psi_3_3"
    ),
    nrow = 3
  ),
  values.observed = theta,
  params.observed = matrix(
    data = c(
      "theta_1_1", "fixed", "fixed",
      "fixed", "theta_2_2", "fixed",
      "fixed", "fixed", "theta_3_3"
    ),
    nrow = 3
  )
)

Prepare the Model

model <- dynr::dynr.model(
  data = dynr_data,
  initial = dynr_initial,
  measurement = dynr_measurement,
  dynamics = dynr_dynamics,
  noise = dynr_noise,
  outfile = "ssm.c"
)

Fit the Model

results <- dynr::dynr.cook(
  model,
  debug_flag = TRUE,
  verbose = FALSE
)
#> [1] "Get ready!!!!"
#> using C compiler: ‘gcc (Ubuntu 13.3.0-6ubuntu2~24.04) 13.3.0’
#> Optimization function called.
#> Starting Hessian calculation ...
#> Finished Hessian calculation.
#> Original exit flag:  3 
#> Modified exit flag:  3 
#> Optimization terminated successfully: ftol_rel or ftol_abs was reached. 
#> Original fitted parameters:  -0.01111479 0.006754815 0.01695684 0.570367 
#> 0.1114704 -0.05800102 0.5171025 0.4670397 0.0101843 0.04070013 0.2466306 
#> 0.5618936 -2.089745 -0.1198123 -0.1608378 -1.462141 0.03533652 -2.585696 
#> -1.712114 -2.493556 -1.567645 -0.31268 0.5428907 0.08120753 -2.256296 0.6458685 
#> 1.699779 -1.547901 -2.477312 -0.1063313 
#> 
#> Transformed fitted parameters:  -0.01111479 0.006754815 0.01695684 0.570367 
#> 0.1114704 -0.05800102 0.5171025 0.4670397 0.0101843 0.04070013 0.2466306 
#> 0.5618936 0.1237187 -0.01482302 -0.01989864 0.2335157 0.01057298 0.07883343 
#> 0.1804839 0.08261565 0.2085357 -0.31268 0.5428907 0.08120753 0.1047377 
#> 0.06764678 0.178031 0.256385 -0.4119249 2.507059 
#> 
#> Doing end processing
#> Successful trial
#> Total Time: 27.33181 
#> Backend Time: 26.70152

Summary

summary(results)
#> Coefficients:
#>             Estimate Std. Error t value  ci.lower  ci.upper Pr(>|t|)    
#> alpha_1_1  -0.011115   0.018420  -0.603 -0.047218  0.024988   0.2733    
#> alpha_2_1   0.006755   0.025266   0.267 -0.042765  0.056275   0.3947    
#> alpha_3_1   0.016957   0.016181   1.048 -0.014758  0.048672   0.1476    
#> beta_1_1    0.570367   0.114848   4.966  0.345269  0.795465   <2e-16 ***
#> beta_1_2    0.111470   0.059513   1.873 -0.005174  0.228115   0.0308 *  
#> beta_1_3   -0.058001   0.059092  -0.982 -0.173819  0.057817   0.1634    
#> beta_2_1    0.517103   0.148811   3.475  0.225437  0.808768   0.0003 ***
#> beta_2_2    0.467040   0.123224   3.790  0.225525  0.708555   0.0001 ***
#> beta_2_3    0.010184   0.087063   0.117 -0.160457  0.180825   0.4535    
#> beta_3_1    0.040700   0.087530   0.465 -0.130855  0.212256   0.3211    
#> beta_3_2    0.246631   0.075807   3.253  0.098052  0.395209   0.0006 ***
#> beta_3_3    0.561894   0.074405   7.552  0.416062  0.707725   <2e-16 ***
#> psi_1_1     0.123719   0.047152   2.624  0.031303  0.216134   0.0045 ** 
#> psi_2_1    -0.014823   0.015467  -0.958 -0.045138  0.015492   0.1692    
#> psi_3_1    -0.019899   0.013849  -1.437 -0.047042  0.007245   0.0757 .  
#> psi_2_2     0.233516   0.086871   2.688  0.063251  0.403780   0.0037 ** 
#> psi_3_2     0.010573   0.015306   0.691 -0.019426  0.040572   0.2450    
#> psi_3_3     0.078833   0.032477   2.427  0.015179  0.142488   0.0078 ** 
#> theta_1_1   0.180484   0.041291   4.371  0.099555  0.261413   <2e-16 ***
#> theta_2_2   0.082616   0.075990   1.087 -0.066322  0.231554   0.1388    
#> theta_3_3   0.208536   0.035254   5.915  0.139439  0.277632   <2e-16 ***
#> mu0_1_1    -0.312680   0.219527  -1.424 -0.742946  0.117586   0.0775 .  
#> mu0_2_1     0.542891   0.259385   2.093  0.034505  1.051276   0.0184 *  
#> mu0_3_1     0.081208   0.733100   0.111 -1.355642  1.518057   0.4559    
#> sigma0_1_1  0.104738   0.166334   0.630 -0.221272  0.430747   0.2646    
#> sigma0_2_1  0.067647   0.129939   0.521 -0.187029  0.322322   0.3014    
#> sigma0_3_1  0.178031   0.384164   0.463 -0.574917  0.930979   0.3216    
#> sigma0_2_2  0.256385   0.235130   1.090 -0.204462  0.717232   0.1380    
#> sigma0_3_2 -0.411925   0.468374  -0.879 -1.329922  0.506072   0.1898    
#> sigma0_3_3  2.507059   1.695344   1.479 -0.815755  5.829873   0.0699 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log-likelihood value at convergence = 2609.40
#> AIC = 2669.40
#> BIC = 2795.84

Parameter Estimates

alpha_hat
#> [1] -0.011114788  0.006754815  0.016956841
beta_hat
#>            [,1]      [,2]        [,3]
#> [1,] 0.57036703 0.1114704 -0.05800102
#> [2,] 0.51710252 0.4670397  0.01018430
#> [3,] 0.04070013 0.2466306  0.56189364
psi_hat
#>             [,1]        [,2]        [,3]
#> [1,]  0.12371870 -0.01482302 -0.01989864
#> [2,] -0.01482302  0.23351568  0.01057298
#> [3,] -0.01989864  0.01057298  0.07883343
theta_hat
#>           [,1]       [,2]      [,3]
#> [1,] 0.1804839 0.00000000 0.0000000
#> [2,] 0.0000000 0.08261565 0.0000000
#> [3,] 0.0000000 0.00000000 0.2085357
mu0_hat
#> [1] -0.31267997  0.54289065  0.08120753
sigma0_hat
#>            [,1]        [,2]       [,3]
#> [1,] 0.10473770  0.06764678  0.1780310
#> [2,] 0.06764678  0.25638497 -0.4119249
#> [3,] 0.17803098 -0.41192489  2.5070593

References

Chow, S.-M., Ho, M. R., Hamaker, E. L., & Dolan, C. V. (2010). Equivalence and differences between structural equation modeling and state-space modeling techniques. Structural Equation Modeling: A Multidisciplinary Journal, 17(2), 303–332. https://doi.org/10.1080/10705511003661553

Ou, L., Hunter, M. D., & Chow, S.-M. (2019). What’s for dynr: A package for linear and nonlinear dynamic modeling in R. The R Journal, 11(1), 91. https://doi.org/10.32614/rj-2019-012

R Core Team. (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/