bssm (version 2.0.2)

ssm_sde: Univariate state space model with continuous SDE dynamics

Description

Constructs an object of class ssm_sde by defining the functions for the drift, diffusion and derivative of diffusion terms of univariate SDE, as well as the log-density of observation equation. We assume that the observations are measured at integer times (missing values are allowed).

Usage

ssm_sde(
  y,
  drift,
  diffusion,
  ddiffusion,
  obs_pdf,
  prior_pdf,
  theta,
  x0,
  positive
)

Value

An object of class ssm_sde.

Arguments

y

Observations as univariate time series (or vector) of length \(n\).

drift, diffusion, ddiffusion

An external pointers for the C++ functions which define the drift, diffusion and derivative of diffusion functions of SDE.

obs_pdf

An external pointer for the C++ function which computes the observational log-density given the the states and parameter vector theta.

prior_pdf

An external pointer for the C++ function which computes the prior log-density given the parameter vector theta.

theta

Parameter vector passed to all model functions.

x0

Fixed initial value for SDE at time 0.

positive

If TRUE, positivity constraint is forced by abs in Milstein scheme.

Details

As in case of ssm_nlg models, these general models need a bit more effort from the user, as you must provide the several small C++ snippets which define the model structure. See vignettes for an example and cpp_example_model.

Examples

 # Takes a while on CRAN
library("sde")
set.seed(1)
# theta_0 = rho = 0.5
# theta_1 = nu = 2
# theta_2 = sigma = 0.3
x <- sde.sim(t0 = 0, T = 50, X0 = 1, N = 50,
       drift = expression(0.5 * (2 - x)),
       sigma = expression(0.3),
       sigma.x = expression(0))
y <- rpois(50, exp(x[-1]))

# source c++ snippets
pntrs <- cpp_example_model("sde_poisson_OU")

sde_model <- ssm_sde(y, pntrs$drift, pntrs$diffusion,
 pntrs$ddiffusion, pntrs$obs_density, pntrs$prior,
 c(rho = 0.5, nu = 2, sigma = 0.3), 1, positive = FALSE)

est <- particle_smoother(sde_model, L = 12, particles = 500)

ts.plot(cbind(x, est$alphahat, 
  est$alphahat - 2*sqrt(c(est$Vt)), 
  est$alphahat + 2*sqrt(c(est$Vt))), 
  col = c(2, 1, 1, 1), lty = c(1, 1, 2, 2))


# Takes time with finer mesh, parallelization with IS-MCMC helps a lot
out <- run_mcmc(sde_model, L_c = 4, L_f = 8, 
  particles = 50, iter = 2e4,
  threads = 4L)