The Ornstein–Uhlenbeck 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} \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.

Data Generation

Notation

Let $t = 1000$ be the number of time points and $n = 1000$ 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] 1000
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

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

Using the `SimSSMOUFixed` Function from the `simStateSpace` Package to Simulate Data

library(simStateSpace)
sim <- SimSSMOUFixed(
  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,
  type = 0
)
data <- as.data.frame(sim)
head(data)
#>   id time          y1          y2        y3
#> 1  1  0.0  0.29937539 -1.37581548 1.3779071
#> 2  1  0.1 -0.98770381 -0.03632195 0.8363080
#> 3  1  0.2  0.33221051 -0.40321664 1.2054318
#> 4  1  0.3 -0.09485392 -0.82030556 1.0272653
#> 5  1  0.4 -1.50322069 -0.36841853 0.1821731
#> 6  1  0.5 -0.75049839  0.35752476 0.2862544
summary(data)
#>        id              time             y1                  y2           
#>  Min.   :   1.0   Min.   : 0.00   Min.   :-3.981741   Min.   :-4.091222  
#>  1st Qu.: 250.8   1st Qu.:24.98   1st Qu.:-0.490953   1st Qu.:-0.584801  
#>  Median : 500.5   Median :49.95   Median : 0.005798   Median : 0.009320  
#>  Mean   : 500.5   Mean   :49.95   Mean   : 0.005499   Mean   : 0.007007  
#>  3rd Qu.: 750.2   3rd Qu.:74.92   3rd Qu.: 0.504090   3rd Qu.: 0.600191  
#>  Max.   :1000.0   Max.   :99.90   Max.   : 4.319747   Max.   : 4.187849  
#>        y3           
#>  Min.   :-3.673010  
#>  1st Qu.:-0.479035  
#>  Median : 0.002762  
#>  Mean   : 0.001891  
#>  3rd Qu.: 0.482969  
#>  Max.   : 4.363632
plot(sim)

Model Fitting

Prepare Data

dynr_data <- dynr::dynr.data(
  dataframe = 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 ~ (phi_1_1 * (eta_1 - mu_1_1)) + (phi_1_2 * (eta_2 - mu_2_1)) + (phi_1_3 * (eta_3 - mu_3_1)),
    eta_2 ~ (phi_2_1 * (eta_1 - mu_1_1)) + (phi_2_2 * (eta_2 - mu_2_1)) + (phi_2_3 * (eta_3 - mu_3_1)),
    eta_3 ~ (phi_3_1 * (eta_1 - mu_1_1)) + (phi_3_2 * (eta_2 - mu_2_1)) + (phi_3_3 * (eta_3 - mu_3_1))
  ),
  startval = c(
    mu_1_1 = mu[1], mu_2_1 = mu[2], mu_3_1 = mu[3],
    phi_1_1 = phi[1, 1], phi_1_2 = phi[1, 2], phi_1_3 = phi[1, 3],
    phi_2_1 = phi[2, 1], phi_2_2 = phi[2, 2], phi_2_3 = phi[2, 3],
    phi_3_1 = phi[3, 1], phi_3_2 = phi[3, 2], phi_3_3 = phi[3, 3]
  ),
  isContinuousTime = TRUE
)

Prepare Process Noise

dynr_noise <- dynr::prep.noise(
  values.latent = sigma,
  params.latent = matrix(
    data = c(
      "sigma_1_1", "sigma_2_1", "sigma_3_1",
      "sigma_2_1", "sigma_2_2", "sigma_3_2",
      "sigma_3_1", "sigma_3_2", "sigma_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 = "ou.c"
)

Add lower and upper bounds to aid in the optimization.

model$lb[
  c(
    "phi_1_1",
    "phi_1_2",
    "phi_1_3",
    "phi_2_1",
    "phi_2_2",
    "phi_2_3",
    "phi_3_1",
    "phi_3_2",
    "phi_3_3"
  )
] <- -1.5
model$ub[
  c(
    "phi_1_1",
    "phi_1_2",
    "phi_1_3",
    "phi_2_1",
    "phi_2_2",
    "phi_2_3",
    "phi_3_1",
    "phi_3_2",
    "phi_3_3"
  )
] <- +1.5

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.003618424 0.003958073 -0.0001418471 -0.3492072 
#> -0.01051393 0.006002923 0.7704453 -0.5100404 -0.0006556109 -0.4511334 0.7283834 
#> -0.6934193 -1.407162 0.0922747 -0.2015102 -2.660961 0.2845036 -2.818866 
#> -1.611389 -1.611329 -1.609704 -0.02470139 0.007251663 0.03819751 -0.006521297 
#> 0.08963148 0.1948516 0.04919685 0.1787325 -0.03387155 
#> 
#> Transformed fitted parameters:  0.003618424 0.003958073 -0.0001418471 
#> -0.3492072 -0.01051393 0.006002923 0.7704453 -0.5100404 -0.0006556109 
#> -0.4511334 0.7283834 -0.6934193 0.2448372 0.02259228 -0.0493372 0.07196571 
#> 0.01532883 0.07527187 0.1996101 0.1996222 0.1999469 -0.02470139 0.007251663 
#> 0.03819751 0.9934999 0.08904887 0.1935851 1.058409 0.2050968 1.037972 
#> 
#> Doing end processing
#> Warning in sqrt(diag(iHess)): NaNs produced
#> Warning in sqrt(diag(x$inv.hessian)): NaNs produced
#> Warning: These parameters may have untrustworthy standard errors: mu0_2_1,
#> sigma0_1_1, sigma0_2_2, sigma0_3_3.
#> Total Time: 1.089781 
#> Backend Time: 1.089773

Summary

summary(results)
#> Coefficients:
#>              Estimate Std. Error  t value   ci.lower   ci.upper Pr(>|t|)    
#> mu_1_1      0.0036184  0.0042587    0.850 -0.0047285  0.0119653   0.1978    
#> mu_2_1      0.0039581  0.0068552    0.577 -0.0094778  0.0173940   0.2818    
#> mu_3_1     -0.0001418  0.0044739   -0.032 -0.0089106  0.0086269   0.4874    
#> phi_1_1    -0.3492072  0.0122732  -28.453 -0.3732622 -0.3251522   <2e-16 ***
#> phi_1_2    -0.0105139  0.0107170   -0.981 -0.0315188  0.0104909   0.1633    
#> phi_1_3     0.0060029  0.0080773    0.743 -0.0098283  0.0218341   0.2287    
#> phi_2_1     0.7704453  0.0061490  125.297  0.7583936  0.7824970   <2e-16 ***
#> phi_2_2    -0.5100404  0.0054032  -94.396 -0.5206305 -0.4994503   <2e-16 ***
#> phi_2_3    -0.0006556  0.0042162   -0.155 -0.0089192  0.0076079   0.4382    
#> phi_3_1    -0.4511334  0.0083562  -53.988 -0.4675113 -0.4347555   <2e-16 ***
#> phi_3_2     0.7283834  0.0074170   98.204  0.7138463  0.7429205   <2e-16 ***
#> phi_3_3    -0.6934193  0.0056349 -123.057 -0.7044636 -0.6823750   <2e-16 ***
#> sigma_1_1   0.2448372  0.0022716  107.783  0.2403850  0.2492894   <2e-16 ***
#> sigma_2_1   0.0225923  0.0007997   28.252  0.0210249  0.0241596   <2e-16 ***
#> sigma_3_1  -0.0493372  0.0009188  -53.695 -0.0511381 -0.0475363   <2e-16 ***
#> sigma_2_2   0.0719657  0.0006040  119.151  0.0707819  0.0731495   <2e-16 ***
#> sigma_3_2   0.0153288  0.0004449   34.453  0.0144568  0.0162009   <2e-16 ***
#> sigma_3_3   0.0752719  0.0006862  109.689  0.0739269  0.0766169   <2e-16 ***
#> theta_1_1   0.1996101  0.0003745  532.954  0.1988760  0.2003442   <2e-16 ***
#> theta_2_2   0.1996222  0.0003167  630.224  0.1990014  0.2002430   <2e-16 ***
#> theta_3_3   0.1999469  0.0003210  622.902  0.1993177  0.2005760   <2e-16 ***
#> mu0_1_1    -0.0247014  0.0219017   -1.128 -0.0676278  0.0182251   0.1297    
#> mu0_2_1     0.0072517        NaN       NA        NaN        NaN       NA    
#> mu0_3_1     0.0381975  0.0178003    2.146  0.0033096  0.0730855   0.0159 *  
#> sigma0_1_1  0.9934999        NaN       NA        NaN        NaN       NA    
#> sigma0_2_1  0.0890489  0.0325965    2.732  0.0251609  0.1529368   0.0031 ** 
#> sigma0_3_1  0.1935851  0.0331723    5.836  0.1285685  0.2586016   <2e-16 ***
#> sigma0_2_2  1.0584087        NaN       NA        NaN        NaN       NA    
#> sigma0_3_2  0.2050968  0.0309378    6.629  0.1444598  0.2657338   <2e-16 ***
#> sigma0_3_3  1.0379722        NaN       NA        NaN        NaN       NA    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log-likelihood value at convergence = 4297047.59
#> AIC = 4297107.59
#> BIC = 4297462.06

#> [1] -0.024701391  0.007251663  0.038197509

Parameter Estimates

mu_hat
#> [1]  0.0036184237  0.0039580727 -0.0001418471
phi_hat
#>            [,1]        [,2]          [,3]
#> [1,] -0.3492072 -0.01051393  0.0060029228
#> [2,]  0.7704453 -0.51004038 -0.0006556109
#> [3,] -0.4511334  0.72838338 -0.6934192985
sigma_hat
#>             [,1]       [,2]        [,3]
#> [1,]  0.24483722 0.02259228 -0.04933720
#> [2,]  0.02259228 0.07196571  0.01532883
#> [3,] -0.04933720 0.01532883  0.07527187
theta_hat
#>           [,1]      [,2]      [,3]
#> [1,] 0.1996101 0.0000000 0.0000000
#> [2,] 0.0000000 0.1996222 0.0000000
#> [3,] 0.0000000 0.0000000 0.1999469
mu0_hat
#> [1] -0.024701391  0.007251663  0.038197509
sigma0_hat
#>            [,1]       [,2]      [,3]
#> [1,] 0.99349992 0.08904887 0.1935851
#> [2,] 0.08904887 1.05840869 0.2050968
#> [3,] 0.19358507 0.20509676 1.0379722
beta_var1_hat <- expm::expm(phi_hat)
beta_var1_hat
#>            [,1]         [,2]        [,3]
#> [1,]  0.7020358 -0.005532608 0.003577427
#> [2,]  0.5011358  0.598106144 0.001021309
#> [3,] -0.1012113  0.400281020 0.499296528

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

Ivan Jacob Agaloos Pesigan

2025-03-25

Model

Data Generation

Notation

R Function Arguments

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

Using the `SimSSMOUFixed` Function from the `simStateSpace` Package to Simulate Data

Model Fitting

Prepare Data

Prepare Initial Condition

Prepare Measurement Model

Prepare Dynamic Process

Prepare Process Noise

Prepare the Model

Fit the Model

Summary

Parameter Estimates

References

The Ornstein–Uhlenbeck Model

Ivan Jacob Agaloos Pesigan

2025-03-25

Model

Data Generation

Notation

R Function Arguments

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

Using the SimSSMOUFixed Function from the simStateSpace Package to Simulate Data

Model Fitting

Prepare Data

Prepare Initial Condition

Prepare Measurement Model

Prepare Dynamic Process

Prepare Process Noise

Prepare the Model

Fit the Model

Summary

Parameter Estimates

References

Using the `SimSSMOUFixed` Function from the `simStateSpace` Package to Simulate Data