This example illustrates how to fit a linear Stochastic differential
equation (SDE) model in the dynr
package.
A damped linear oscillator model is estimated. Details of the model set up is as follows:
d2xdt2=−kx−cdxdt+ζ
This model can also be written in the form of a vector of first-order systems:
The dynamic model is :
(dxdtd2xdt2)=(01−k−c)(xdxdt)+(0ζ)
The measurement model is:
y=(10)(xdxdt)+ϵ
We first load dynr package.
The dynr.data command is used to create a dynr data object. Here the observed variable we are modeling is “y1” . (The simulated data are loaded as part of the package.)
## Loading required package: dynr
## Loading required package: ggplot2
The prep.measurement command
is used to specify the factor loadings matrix. The
values.load
argument specifies that the factor loadings are
fixed at 1 and 0. In this model, all parameters are fixed, which is
indicated by the params.load
argument. The
*.names
argument gives names to the latent and observed
variables.
The prep.noise command is
used to specify dynamic (values.latent
) and observation
(values.observed
) noise components. The latent noise is the
dynamic noise. The observed noise is the measurement noise.
The prep.matrixDynamics command is used to define the differental equation. In this model, the first row of the parameter matrix is fixed at 0 and 1 and the starting values for the second row of the matrix is set at -0.1 and -0.2.
ecov <- prep.noise(
values.latent=diag(c(0, 1), 2), params.latent=diag(c('fixed', 'dnoise'), 2),
values.observed=diag(1.5, 1), params.observed=diag('mnoise', 1))
dynamics <- prep.matrixDynamics(
values.dyn=matrix(c(0, -0.1, 1, -0.2), 2, 2),
params.dyn=matrix(c('fixed', 'spring', 'fixed', 'friction'), 2, 2),
isContinuousTime=TRUE)
This step specifies the covariances and latent state values at t=0. These initialize the recursive algorithm (extended Kalman filter) that dynr uses.
Now we put together everything we’ve previously specified in dynr.model. This code connects the recipes we’ve written up with our data and writes a c file in our working directory. We can inspect c functions that go with each recipe in the c file.
We can check our model specifications in a neatly printed pdf file using the following code.
The printex command is used
to write the model into a Latex file, with a name given by the
outFile
argument. Then, the tools::texi2pdf command generates a
pdf file from the latex file we just created. The system command prints out the pdf
file:
We can also print out the model in R, instead of generating a Latex file, using the command plotFormula.