Parametric Bootstrap (The Ornstein–Uhlenbeck Model)

Ivan Jacob Agaloos Pesigan

2025-01-18

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} \mathrm{d} \boldsymbol{\eta}_{i, t} = \boldsymbol{\Phi} \left( \boldsymbol{\eta}_{i, t} - \boldsymbol{\mu} \right) \mathrm{d}t + \boldsymbol{\Sigma}^{\frac{1}{2}} \mathrm{d} \mathbf{W}_{i, t} \end{equation}$$ where $\boldsymbol{\mu}$ is the long-term mean or equilibrium level, $\boldsymbol{\Phi}$ is the rate of mean reversion, determining how quickly the variable returns to its mean, $\boldsymbol{\Sigma}$ is the matrix of volatility or randomness in the process, and $\mathrm{d}\boldsymbol{W}$ is a Wiener process or Brownian motion, which represents random fluctuations.

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 long-term mean vector $\boldsymbol{\mu}$ be given by

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

Let the rate of mean reversion matrix $\boldsymbol{\Phi}$ be given by

$$\begin{equation} \boldsymbol{\Phi} = \left( \begin{array}{ccc} -0.357 & 0 & 0 \\ 0.771 & -0.511 & 0 \\ -0.45 & 0.729 & -0.693 \\ \end{array} \right) . \end{equation}$$

Let the dynamic process noise covariance matrix $\boldsymbol{\Sigma}$ be given by

$$\begin{equation} \boldsymbol{\Sigma} = \left( \begin{array}{ccc} 0.2445556 & 0.0220159 & -0.0500476 \\ 0.0220159 & 0.070678 & 0.0153946 \\ -0.0500476 & 0.0153946 & 0.0755306 \\ \end{array} \right) . \end{equation}$$

Let $\Delta t = 0.1$.

R Function Arguments

n
#> [1] 5
time
#> [1] 1000
delta_t
#> [1] 0.1
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
mu
#> [1] 0 0 0
phi
#>        [,1]   [,2]   [,3]
#> [1,] -0.357  0.000  0.000
#> [2,]  0.771 -0.511  0.000
#> [3,] -0.450  0.729 -0.693
sigma
#>             [,1]       [,2]        [,3]
#> [1,]  0.24455556 0.02201587 -0.05004762
#> [2,]  0.02201587 0.07067800  0.01539456
#> [3,] -0.05004762 0.01539456  0.07553061
sigma_l # sigma_l <- t(chol(sigma))
#>             [,1]      [,2]     [,3]
#> [1,]  0.49452559 0.0000000 0.000000
#> [2,]  0.04451917 0.2620993 0.000000
#> [3,] -0.10120330 0.0759256 0.243975
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 <- "ou"

We use the PBSSMOUFixed 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 Ornstein–Uhlenbeck model multiple times is computationally intensive.

library(bootStateSpace)
pb <- PBSSMOUFixed(
  R = R,
  path = path,
  prefix = prefix,
  n = n,
  time = time,
  delta_t = delta_t,
  mu0 = mu0,
  sigma0_l = sigma0_l,
  mu = mu,
  phi = phi,
  sigma_l = sigma_l,
  nu = nu,
  lambda = lambda,
  theta_l = theta_l,
  ncores = parallel::detectCores(),
  seed = 42
)
summary(pb)
#> Call:
#> PBSSMOUFixed(R = R, path = path, prefix = prefix, n = n, time = time, 
#>     delta_t = delta_t, mu0 = mu0, sigma0_l = sigma0_l, mu = mu, 
#>     phi = phi, sigma_l = sigma_l, nu = nu, lambda = lambda, theta_l = theta_l, 
#>     ncores = parallel::detectCores(), seed = 42)
#>                est     se    R    2.5%   97.5%
#> phi_1_1    -0.3570 0.1675 1000 -0.7096 -0.0697
#> phi_2_1     0.7710 0.1069 1000  0.5697  0.9811
#> phi_3_1    -0.4500 0.1009 1000 -0.6561 -0.2732
#> phi_1_2     0.0000 0.1468 1000 -0.2557  0.3007
#> phi_2_2    -0.5110 0.0943 1000 -0.7053 -0.3323
#> phi_3_2     0.7290 0.0894 1000  0.5624  0.9165
#> phi_1_3     0.0000 0.1128 1000 -0.2537  0.1974
#> phi_2_3     0.0000 0.0708 1000 -0.1377  0.1327
#> phi_3_3    -0.6930 0.0704 1000 -0.8356 -0.5639
#> sigma_1_1   0.2446 0.0310 1000  0.1909  0.3117
#> sigma_2_1   0.0220 0.0118 1000 -0.0017  0.0431
#> sigma_3_1  -0.0500 0.0124 1000 -0.0741 -0.0255
#> sigma_2_2   0.0707 0.0091 1000  0.0536  0.0892
#> sigma_3_2   0.0154 0.0064 1000  0.0017  0.0267
#> sigma_3_3   0.0755 0.0098 1000  0.0570  0.0950
#> theta_1_1   0.2000 0.0051 1000  0.1900  0.2103
#> theta_2_2   0.2000 0.0046 1000  0.1915  0.2095
#> theta_3_3   0.2000 0.0046 1000  0.1912  0.2085
#> mu0_1_1     0.0000 0.4713 1000 -0.9073  0.8767
#> mu0_2_1     0.0000 0.4589 1000 -0.9450  0.8877
#> mu0_3_1     0.0000 0.4455 1000 -0.8783  0.8600
#> sigma0_1_1  1.0000 0.6262 1000  0.0491  2.3055
#> sigma0_2_1  0.2000 0.4099 1000 -0.6189  1.0723
#> sigma0_3_1  0.2000 0.4364 1000 -0.6363  1.1324
#> sigma0_2_2  1.0000 0.5931 1000  0.0567  2.1617
#> sigma0_3_2  0.2000 0.4286 1000 -0.5908  1.1979
#> sigma0_3_3  1.0000 0.6133 1000  0.0455  2.3477
summary(pb, type = "bc")
#> Call:
#> PBSSMOUFixed(R = R, path = path, prefix = prefix, n = n, time = time, 
#>     delta_t = delta_t, mu0 = mu0, sigma0_l = sigma0_l, mu = mu, 
#>     phi = phi, sigma_l = sigma_l, nu = nu, lambda = lambda, theta_l = theta_l, 
#>     ncores = parallel::detectCores(), seed = 42)
#>                est     se    R    2.5%   97.5%
#> phi_1_1    -0.3570 0.1675 1000 -0.6991 -0.0571
#> phi_2_1     0.7710 0.1069 1000  0.5650  0.9795
#> phi_3_1    -0.4500 0.1009 1000 -0.6402 -0.2653
#> phi_1_2     0.0000 0.1468 1000 -0.2514  0.3026
#> phi_2_2    -0.5110 0.0943 1000 -0.7051 -0.3323
#> phi_3_2     0.7290 0.0894 1000  0.5530  0.8976
#> phi_1_3     0.0000 0.1128 1000 -0.2496  0.2018
#> phi_2_3     0.0000 0.0708 1000 -0.1370  0.1340
#> phi_3_3    -0.6930 0.0704 1000 -0.8275 -0.5506
#> sigma_1_1   0.2446 0.0310 1000  0.1934  0.3157
#> sigma_2_1   0.0220 0.0118 1000 -0.0005  0.0443
#> sigma_3_1  -0.0500 0.0124 1000 -0.0741 -0.0257
#> sigma_2_2   0.0707 0.0091 1000  0.0564  0.0922
#> sigma_3_2   0.0154 0.0064 1000  0.0015  0.0266
#> sigma_3_3   0.0755 0.0098 1000  0.0587  0.0976
#> theta_1_1   0.2000 0.0051 1000  0.1900  0.2102
#> theta_2_2   0.2000 0.0046 1000  0.1912  0.2093
#> theta_3_3   0.2000 0.0046 1000  0.1912  0.2086
#> mu0_1_1     0.0000 0.4713 1000 -0.8934  0.8838
#> mu0_2_1     0.0000 0.4589 1000 -0.8479  0.9346
#> mu0_3_1     0.0000 0.4455 1000 -0.8713  0.8708
#> sigma0_1_1  1.0000 0.6262 1000  0.2581  3.8576
#> sigma0_2_1  0.2000 0.4099 1000 -0.4145  1.3483
#> sigma0_3_1  0.2000 0.4364 1000 -0.4470  1.3556
#> sigma0_2_2  1.0000 0.5931 1000  0.2967  3.9061
#> sigma0_3_2  0.2000 0.4286 1000 -0.4373  1.3701
#> sigma0_3_3  1.0000 0.6133 1000  0.2395  3.3381

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

Chow, S.-M., Losardo, D., Park, J., & Molenaar, P. C. M. (2023). Continuous-time dynamic models: Connections to structural equation models and other discrete-time models. In R. H. Hoyle (Ed.), Handbook of structural equation modeling (2nd ed.). The Guilford Press.

Deboeck, P. R., & Preacher, K. J. (2015). No need to be discrete: A method for continuous time mediation analysis. Structural Equation Modeling: A Multidisciplinary Journal, 23(1), 61–75. https://doi.org/10.1080/10705511.2014.973960

Oravecz, Z., Tuerlinckx, F., & Vandekerckhove, J. (2011). A hierarchical latent stochastic differential equation model for affective dynamics. Psychological Methods, 16(4), 468–490. https://doi.org/10.1037/a0024375

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/

Uhlenbeck, G. E., & Ornstein, L. S. (1930). On the theory of the brownian motion. Physical Review, 36(5), 823–841. https://doi.org/10.1103/physrev.36.823