demo/SDETVP.R

In dynr: Dynamic Models with Regime-Switching

#--------------------------------------------------------
# Filename: LogisticSetpoint.R
# Purpose: An illustrative example of using dynr to fit
#   a stochastic damped oscillator model with time-
#   varying (logistic growth) setpoint parameter
#   (Chen, Chow & Hunter, expected 2018)
#--------------------------------------------------------

#---- (1) Load packages ----
library(dynr)
#---- (2) Read in data ----
data(LogisticSetPointSDE)
dataLog <- dynr.data(LogisticSetPointSDE, id="id", time="times", observed="obsy")
#---- (3) Specify recipes for all model pieces ----
#---- (3a) Measurement ----
measLog <- prep.measurement(
  values.load=matrix(c(1, 0, 0), 1, 3), # starting values and fixed values
  params.load=matrix("fixed", 1, 3),
  state.names=c("x","y","z"),
  obs.names=c("obsy")) # parameter numbers or indication that parameter is fixed

#---- (3b) Dynamic and measurement noise cov structures----
ecovLog <- prep.noise(
  values.latent=diag(c(0,0.1,0), 3), params.latent=diag(c('fixed', 'pnoise','fixed'), 3), 
  values.observed=diag(0.1, 1), params.observed=diag('mnoise', 1)) 

#---- (3c) Initial condition specification ----
initialLog <- prep.initial(
  values.inistate=c(-5, 0 ,.1),
  params.inistate=c('fixed', 'fixed','fixed'), 
  values.inicov=diag(0, 3),
  params.inicov=diag('fixed', 3)) 


#---- (3d) Dynamic model ----
formulaLogistic<-list(
  x~y,
  y~eta*(x-z)+damp*y,
  z~r*z*(1-b*z)
)

dynamicsLog<-prep.formulaDynamics(formula=formulaLogistic,startval=c(eta=-.1,damp=-.05,b=0.5,r=0.3),
                                  isContinuousTime=TRUE)

#---- (4) Create model and cook it all up ----

modelLog <- dynr.model(dynamics=dynamicsLog, measurement=measLog, 
                       noise=ecovLog, initial=initialLog, #transform=trans,
                       data=dataLog, outfile="setpointLog.c")


modelLog$ub[c("eta",'r', 'pnoise')] <- c(0,10,10)


resLog <- dynr.cook(modelLog, verbose=FALSE)
#---- (5) Serve it! ----

summary(resLog)

round(coef(resLog),3)
# True parameters: freq=-1; damp=-.1; b=0.1; r=.5; sigma1=0.1, mnoise=1