--- title: "Example: A Linear Stochastic Differential Equation Model" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Example: A Linear Stochastic Differential Equation Model} %\VignetteEngine{knitr::rmarkdown} \usepackage[utf8]{inputenc} --- Introduction --------------- 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: $\frac{d^2x}{dt^2} = -kx - c\frac{dx}{dt} + \zeta$ This model can also be written in the form of a vector of first-order systems: The dynamic model is : $$\begin{pmatrix} \frac{dx}{dt}\\ \frac{d^2x}{dt^2} \end{pmatrix} = \begin{pmatrix} 0 & 1\\ -k & -c \end{pmatrix}\begin{pmatrix} x\\ \frac{dx}{dt} \end{pmatrix} + \begin{pmatrix} 0\\ \zeta \end{pmatrix}$$ The measurement model is: $$y = \begin{pmatrix} 1 & 0 \end{pmatrix}\begin{pmatrix} x\\ \frac{dx}{dt} \end{pmatrix} + \epsilon$$ Data --------------- 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.) ```{r data, results="hide"} require(dynr) # Data data(Oscillator) data <- dynr.data(Oscillator, id="id", time="times", observed="y1") ``` Measurement Model --------------- 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. ```{r measurement, results="hide"} meas <- prep.measurement( values.load=matrix(c(1, 0), 1, 2), params.load=matrix(c('fixed', 'fixed'), 1, 2), state.names=c("Position","Velocity"), obs.names=c("y1")) ``` Dynamic Model --------------- 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. ```{r noise cov, results="hide"} 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) ``` Initial Values --------------- This step specifies the covariances and latent state values at t=0. These initialize the recursive algorithm (extended Kalman filter) that dynr uses. ```{r initials, results="hide"} initial <- prep.initial( values.inistate=c(0, 1), params.inistate=c('inipos', 'fixed'), values.inicov=diag(1, 2), params.inicov=diag('fixed', 2)) ``` Dynr Model --------------- 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. ```{r model, results="hide"} model <- dynr.model(dynamics=dynamics, measurement=meas, noise=ecov, initial=initial, data=data, outfile="LinearSDE.c") ``` Tex Options --------------- 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: ![](LinearSDE.png) We can also print out the model in R, instead of generating a Latex file, using the command _plotFormula_. ```{r tex, results="hide",eval=FALSE} printex(model,ParameterAs=model$param.names,show=FALSE,printInit=TRUE, outFile="LinearSDE.tex") tools::texi2pdf("LinearSDE.tex") system(paste(getOption("pdfviewer"), "LinearSDE.pdf")) ``` Optimization Step and Results --------------- Now the model is specified and we are ready to cook dynr! The _summary_ command gives the model fitting results. ```{r cook, results="hide", eval=FALSE} res <- dynr.cook(model, verbose=FALSE) ``` ```{r serve, eval=FALSE} summary(res) ```