Parametric Bootstrap (The State Space Model)

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

Parameters

Notation

Let $t = 1000$ 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] 1000
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

Parametric Bootstrap

R <- 1000L # use at least 1000 in actual research
path <- getwd()
prefix <- "ssm"

We use the PBSSMFixed function from the bootStateSpace package to perform parametric bootstraping using the parameters described above. The argument R specifies the number of bootstrap replications. The generated data and model estimates are stored in path using the specified prefix for the file names. The ncores = parallel::detectCores() argument instructs the function to use all available CPU cores in the system.

NOTE: Fitting the state space model multiple times is computationally intensive.

library(bootStateSpace)
pb <- PBSSMFixed(
  R = R,
  path = path,
  prefix = prefix,
  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,
  ncores = parallel::detectCores(),
  seed = 42
)
summary(pb)
#> Call:
#> PBSSMFixed(R = R, path = path, prefix = prefix, 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, 
#>     ncores = parallel::detectCores(), seed = 42)
#>             est     se    R    2.5%   97.5%
#> beta_1_1    0.7 0.0430 1000  0.6116  0.7813
#> beta_2_1    0.5 0.0378 1000  0.4310  0.5792
#> beta_3_1   -0.1 0.0269 1000 -0.1551 -0.0512
#> beta_1_2    0.0 0.0272 1000 -0.0507  0.0549
#> beta_2_2    0.6 0.0317 1000  0.5371  0.6575
#> beta_3_2    0.4 0.0284 1000  0.3493  0.4595
#> beta_1_3    0.0 0.0185 1000 -0.0374  0.0346
#> beta_2_3    0.0 0.0230 1000 -0.0436  0.0450
#> beta_3_3    0.5 0.0286 1000  0.4418  0.5523
#> psi_1_1     0.1 0.0136 1000  0.0758  0.1299
#> psi_2_2     0.1 0.0113 1000  0.0798  0.1232
#> psi_3_3     0.1 0.0130 1000  0.0778  0.1286
#> theta_1_1   0.2 0.0115 1000  0.1740  0.2206
#> theta_2_2   0.2 0.0106 1000  0.1781  0.2183
#> theta_3_3   0.2 0.0130 1000  0.1714  0.2226
#> mu0_1_1     0.0 0.4736 1000 -0.9850  0.9076
#> mu0_2_1     0.0 0.4644 1000 -0.8640  1.0137
#> mu0_3_1     0.0 0.4684 1000 -0.9131  0.9377
#> sigma0_1_1  1.0 0.6458 1000  0.0349  2.3562
#> sigma0_2_1  0.2 0.4791 1000 -0.6613  1.2473
#> sigma0_3_1  0.2 0.4459 1000 -0.6209  1.1667
#> sigma0_2_2  1.0 0.6609 1000  0.0351  2.5214
#> sigma0_3_2  0.2 0.4653 1000 -0.6812  1.1830
#> sigma0_3_3  1.0 0.6390 1000  0.0421  2.4159
summary(pb, type = "bc")
#> Call:
#> PBSSMFixed(R = R, path = path, prefix = prefix, 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, 
#>     ncores = parallel::detectCores(), seed = 42)
#>             est     se    R    2.5%   97.5%
#> beta_1_1    0.7 0.0430 1000  0.6144  0.7831
#> beta_2_1    0.5 0.0378 1000  0.4313  0.5808
#> beta_3_1   -0.1 0.0269 1000 -0.1551 -0.0512
#> beta_1_2    0.0 0.0272 1000 -0.0534  0.0530
#> beta_2_2    0.6 0.0317 1000  0.5355  0.6563
#> beta_3_2    0.4 0.0284 1000  0.3452  0.4561
#> beta_1_3    0.0 0.0185 1000 -0.0388  0.0327
#> beta_2_3    0.0 0.0230 1000 -0.0429  0.0457
#> beta_3_3    0.5 0.0286 1000  0.4455  0.5557
#> psi_1_1     0.1 0.0136 1000  0.0761  0.1308
#> psi_2_2     0.1 0.0113 1000  0.0808  0.1244
#> psi_3_3     0.1 0.0130 1000  0.0768  0.1283
#> theta_1_1   0.2 0.0115 1000  0.1738  0.2206
#> theta_2_2   0.2 0.0106 1000  0.1773  0.2182
#> theta_3_3   0.2 0.0130 1000  0.1735  0.2234
#> mu0_1_1     0.0 0.4736 1000 -0.9757  0.9435
#> mu0_2_1     0.0 0.4644 1000 -0.8997  0.9845
#> mu0_3_1     0.0 0.4684 1000 -0.9102  0.9393
#> sigma0_1_1  1.0 0.6458 1000  0.2545  4.1916
#> sigma0_2_1  0.2 0.4791 1000 -0.4479  1.5940
#> sigma0_3_1  0.2 0.4459 1000 -0.4620  1.4333
#> sigma0_2_2  1.0 0.6609 1000  0.2410  3.7139
#> sigma0_3_2  0.2 0.4653 1000 -0.4630  1.6344
#> sigma0_3_3  1.0 0.6390 1000  0.2078  3.3093

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/