enhance data generation script - Oscillator

[ares7823]

Hi @mhunter1,
I used the data simulation codes from 'Oscillator' in dataDoc.R and I found that it is only suitable for equal time interval t = 1.
The main relevant part is

tx[,1] <- x0
for(i in 2:(tdim+1)){
q <- t(rmvnorm(1, rep(0, xdim), tQ))
tdx[,i] <- tA %*% tx[,i-1] + tB %*% tu[,i-1] + q
expA <- as.matrix(expm(tA * (tT[,i]-tT[,i-1])))
intA <- solve(tA) %*% (expA - tI)
tx[,i] <- expA %*% tx[, i-1] + intA %*% tB %*% tu[,i-1] + intA %*% q
ty[,i-1] <- tC %*% tx[,i] + tD %*% tu[,i-1] + t(rmvnorm(1, rep(0, ydim), tR))
}

Here, I think tQ matrix and pulled error q should reflect on the time step in the continuous-time modeling.
Based on Voelkle & Oud (2013, p109, Eq 4), I tried to have the code be flexible to time step, such as:

for(i in 2:(tdim+1)){
  deltat <- tT[,i]-tT[,i-1]
  tA_pnd <- kronecker(tA, diag(1,xdim))+ kronecker(diag(1,xdim),tA)
  expA_pnd <- as.matrix(expm(tA_pnd * deltat))
  intA_pnd <- solve(tA_pnd) %*% (expA_pnd - diag(1,dim(expA_pnd)[1])) 
  tQ_d <- matrix(intA_pnd %*% c(t(tQ)), nrow = xdim, byrow = T)
  
  q <- t(rmvnorm(1, rep(0, xdim), tQ_d)) 
  expA <- as.matrix(expm(tA * deltat)) 
  intA <- solve(tA) %*% (expA - tI) 
  tx[,i] <- expA %*% tx[, i-1] + intA %*% tB %*% tu[,i-1] + q 
  ty[,i-1] <- tC %*% tx[,i] + tD %*% tu[,i-1] + t(rmvnorm(1, rep(0, ydim), tR)) # measurement model
}

I checked that the codes worked but I would appreciate it if you can verify this.

[mhunter1]

@ares7823 You are correct that the continuous-time dynamic noise matrix (Q) needs to be discretized (1) to allow for unequal observations intervals, and (2) to actually be correct.

I'm less sure about the way you're discretizing. The fundamental problem is computing the integral

int_0^{deltaT} expm(A tau) %*% Q %*% t(expm(A tau))  dtau

I'm assuming the Kronecker products and the vectorizations from Voelkle & Oud are intended for that purpose. However, I haven't reviewed them closely and the math looks like that of the asymptotic latent variable covariance. So, I'm suspicious. The dynamic noise discretization we use in dynr uses RK4 on an extended differential equation. The dynamic noise discretization we use in OpenMx uses an analytic solution based on block matrices of a matrix exponential from Van Loan (1978). I have an R implementation of the latter at https://github.com/OpenMx/OpenMx/blob/master/R/MxExpectationStateSpace.R#L769:L795

[ares7823]

@mhunter1 Thanks for sharing your thought. The way I used is for the case that deltat is any fixed number besides 1 but I haven't thought about whether it would work as well for the unequal space. The problem I had is getting small q when I directly use the code from dynr with deltat = .1