Sample from a simulated Ornstein-Uhlenbeck process with time-dependent mean

Description

The data was simulated using a standard Euler-Maruyama method.

The simulated process is governed by the SDE #' dx ~ theta * (mu + u - x) * dt + sigma_x * dw

The parameters used for simulation were theta = 2, mu = 0.5, sigma_x = 1.358, sigma_y = 1e-8

The simulation time-step was 1e-3, and observation time-step 1e-1. The simulation was taken from t = 0..20

Usage

Ornstein

Format

A data frame of 201 rows and 3 columns. The columns represent the variables: t (time), y (observation) and u (input).

Shortcut class constructor

Description

This is just a short-cut class constructor so that users can do

model <- ctsm()

instead of

model <- ctsmTMB$new().

Usage

ctsm()

Examples

library(ctsmTMB)
model <- ctsm()

Methods for the 'ctsmTMB' R6 class

Description

The following public methods are used to construct a stochastic state space model system, consisting of a set of stochastic differential equations (SDEs), and one or more algebraic observation equations (AOEs). The AOEs are used to infer information about the value of the (latent) states governed by the SDEs, and thus must be functions of at least one state.

Value

The function returns an object of class R6 and ctsmTMB, which can be used to define a stochastic state space system.

Methods

Public methods

• ctsmTMB$new()

• ctsmTMB$.private()

• ctsmTMB$addSystem()

• ctsmTMB$addObs()

• ctsmTMB$setVariance()

• ctsmTMB$addInput()

• ctsmTMB$setParameter()

• ctsmTMB$setAlgebraics()

• ctsmTMB$setInitialState()

• ctsmTMB$setInitialVarianceScaling()

• ctsmTMB$setLamperti()

• ctsmTMB$setModelname()

• ctsmTMB$setMAP()

• ctsmTMB$getSystems()

• ctsmTMB$getObservations()

• ctsmTMB$getVariances()

• ctsmTMB$getAlgebraics()

• ctsmTMB$getInitialState()

• ctsmTMB$getParameters()

• ctsmTMB$getEstimate()

• ctsmTMB$getLikelihood()

• ctsmTMB$getPrediction()

• ctsmTMB$getSimulation()

• ctsmTMB$estimate()

• ctsmTMB$likelihood()

• ctsmTMB$predict()

• ctsmTMB$simulate()

• ctsmTMB$print()

• ctsmTMB$clone()

Method new()

Initialize private fields

Usage

ctsmTMB$new()

Method .private()

Extract the private fields of a ctsmTMB model object. Primarily used for debugging.

Usage

ctsmTMB$.private()

Method addSystem()

Define stochastic differential equation(s) on the form

d<state> ~ f(t,<states>, <inputs>) * dt + g(t, <states>, <inputs>) * dw

Usage

ctsmTMB$addSystem(form, ...)

Arguments

Method addObs()

Define algebraic observation equations on the form

<observation> ~ h(t, <states>, <inputs>) + e)

where h is the observation function, and e is normally distributed noise with zero mean.

This function only specifies the observation name, and its mean through h.

Usage

ctsmTMB$addObs(form, ..., obsnames = NULL)

Arguments

Method setVariance()

Specify the variance of an observation equation.

A defined observation variable y in e.g. addObs(y ~ h(t,<states>,<inputs>) is perturbed by Gaussian noise with zero mean and variance to-be specified using setVariance(y ~ p(t,<states>,<inputs>). We can for instance declare setVariance(y ~ sigma_x^2 where sigma_x is a fixed effect parameter to be declared through setParameter.

Usage

ctsmTMB$setVariance(form, ...)

Arguments

Method addInput()

Declare variables as data inputs

Declare whether a variable contained in system, observation or observation variance equations is an input variable. If e.g. the system equation contains an input variable u then it is declared using addInput(u). The input u must be contained in the data.frame .data provided when calling the estimate or predict methods.

Usage

ctsmTMB$addInput(...)

Arguments

Method setParameter()

Declare which variables that are (fixed effects) parameters in the specified model, and specify the initial optimizer guess, as well as lower / upper bounds during optimization. There are two ways to declare parameters:

1. You can declare parameters using formulas i.e. setParameter( theta = c(1,0,10), mu = c(0,-10,10) ). The first value is the initial value for the optimizer, the second value is the lower optimization bound and the third value is the upper optimization bound.

2. You can provide a 3-column matrix where rows corresponds to different parameters, and the parameter names are provided as rownames of the matrix. The columns values corresponds to the description in the vector format above.

Usage

ctsmTMB$setParameter(...)

Arguments

Method setAlgebraics()

Add algebraic relations.

Algebraic relations is a convenient way to transform parameters in your equations. In the Ornstein-Uhlenbeck process the rate parameter theta is always positive, so estimation in the log-domain is a good idea. Instead of writing exp(theta) directly in the system equation one can transform into the log domain using the algebraic relation setAlgebraics(theta ~ exp(logtheta)). All instances of theta is replaced by exp(logtheta) when compiling the C++ function. Note that you must provide values for logtheta now instead of theta when declaring parameters through setParameter

Usage

ctsmTMB$setAlgebraics(form, ...)

Arguments

Method setInitialState()

Declare the initial state values i.e. mean and covariance for the system states.

Usage

ctsmTMB$setInitialState(initial.state)

Arguments

Method setInitialVarianceScaling()

A scalar value that is multiplied onto the estimated initial state covariance matrix. The scaling is only applied when the initial state/cov is estimated, not when it is set by the user.

Usage

ctsmTMB$setInitialVarianceScaling(scaling)

Arguments

Method setLamperti()

Set a Lamperti Transformation

If the provided system equations have state dependent diffusion in a few available ways then it is advantageous to perform a transformation to remove the state dependence. This comes at the cost of a more complicated drift function. The following types of state-dependence is currently supported

1. 'identity' - when the diffusion is state-independent (default)

2. 'log' - when the diffusion is proportional to to x * dw

3. 'logit' - when the diffusion is proportional to x * (1-x) * dw

4. 'sqrt-logit' - when the diffusion is proportional to sqrt(x * (1-x)) * dw

Usage

ctsmTMB$setLamperti(transforms, states = NULL)

Arguments

Method setModelname()

Set modelname used to create the C++ file for TMB

When calling TMB::MakeADFun the (negative log) likelihood function is created in the directory specified by the setCppfilesDirectory method with name <modelname>.cpp

Usage

ctsmTMB$setModelname(name)

Arguments

Method setMAP()

Enable maximum a posterior (MAP) estimation.

Adds a maximum a posterior contribution to the (negative log) likelihood function by evaluating the fixed effects parameters in a multivariate Gaussian with mean and covariance as provided.

Usage

ctsmTMB$setMAP(mean, cov)

Arguments

Method getSystems()

Retrieve system equations.

Usage

ctsmTMB$getSystems()

Method getObservations()

Retrieve observation equations.

Usage

ctsmTMB$getObservations()

Method getVariances()

Retrieve observation variances

Usage

ctsmTMB$getVariances()

Method getAlgebraics()

Retrieve algebraic relations

Usage

ctsmTMB$getAlgebraics()

Method getInitialState()

Retrieve initially set state and covariance

Usage

ctsmTMB$getInitialState()

Method getParameters()

Get initial (and estimated) parameters.

Usage

ctsmTMB$getParameters(type = "all", value = "all")

Arguments

Method getEstimate()

Retrieve initially set state and covariance

Usage

ctsmTMB$getEstimate()

Method getLikelihood()

Retrieve initially set state and covariance

Usage

ctsmTMB$getLikelihood()

Method getPrediction()

Retrieve initially set state and covariance

Usage

ctsmTMB$getPrediction()

Method getSimulation()

Retrieve initially set state and covariance

Usage

ctsmTMB$getSimulation()

Method estimate()

Estimate the fixed effects parameters in the specified model.

Usage

ctsmTMB$estimate(
  data,
  method = "ekf",
  ode.solver = "euler",
  ode.timestep = diff(data$t),
  loss = "quadratic",
  loss_c = NULL,
  control = list(trace = 1, iter.max = 1e+05, eval.max = 1e+05),
  use.hessian = FALSE,
  laplace.residuals = FALSE,
  unconstrained.optim = FALSE,
  estimate.initial.state = FALSE,
  silent = FALSE
)

Arguments

Method likelihood()

Construct and extract function handlers for the negative log likelihood function.

The handlers from TMB's MakeADFun are constructed and returned. This enables the user to e.g. choose their own optimization algorithm, or just have more control of the optimization workflow.

Usage

ctsmTMB$likelihood(
  data,
  method = "ekf",
  ode.solver = "euler",
  ode.timestep = diff(data$t),
  loss = "quadratic",
  loss_c = NULL,
  estimate.initial.state = FALSE,
  silent = FALSE
)

Arguments

Method predict()

Perform prediction/filtration to obtain state mean and covariance estimates. The predictions are obtained by solving the moment equations n.ahead steps forward in time when using the current step posterior state estimate as the initial condition.

Usage

ctsmTMB$predict(
  data,
  pars = NULL,
  method = "ekf",
  ode.solver = "euler",
  ode.timestep = diff(data$t),
  k.ahead = nrow(data) - 1,
  return.k.ahead = 0:k.ahead,
  return.covariance = TRUE,
  initial.state = self$getInitialState(),
  estimate.initial.state = private$estimate.initial,
  silent = FALSE,
  use.cpp = FALSE
)

Arguments

Returns

A data.frame that contains for each time step the posterior state estimate at that time.step (k = 0), and the prior state predictions (k = 1,...,n.ahead). If return.covariance = TRUE then the state covariance/correlation matrix is returned, otherwise only the marginal variances are returned.

Method simulate()

Perform prediction/filtration to obtain state mean and covariance estimates. The predictions are obtained by solving the moment equations n.ahead steps forward in time when using the current step posterior state estimate as the initial condition.

Usage

ctsmTMB$simulate(
  data,
  pars = NULL,
  use.cpp = FALSE,
  method = "ekf",
  ode.solver = "rk4",
  ode.timestep = diff(data$t),
  simulation.timestep = diff(data$t),
  k.ahead = nrow(data) - 1,
  return.k.ahead = 0:k.ahead,
  n.sims = 100,
  initial.state = self$getInitialState(),
  estimate.initial.state = private$estimate.initial,
  silent = FALSE
)

Arguments

Returns

A data.frame that contains for each time step the posterior state estimate at that time.step (k = 0), and the prior state predictions (k = 1,...,n.ahead). If return.covariance = TRUE then the state covariance/correlation matrix is returned, otherwise only the marginal variances are returned.

Method print()

Function to print the model object

Usage

ctsmTMB$print()

Method clone()

The objects of this class are cloneable with this method.

Usage

ctsmTMB$clone(deep = FALSE)

Arguments

Examples

library(ctsmTMB)
model <- ctsmTMB$new()

# adding a single system equations
model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw)

# adding an observation equation and setting variance
model$addObs(y ~ x)
model$setVariance(y ~ sigma_y^2)

# add model input
model$addInput(u)

# add parameters
model$setParameter(
  theta   = c(initial = 1, lower=1e-5, upper=50),
  mu      = c(initial=1.5, lower=0, upper=5),
  sigma_x = c(initial=1, lower=1e-10, upper=30),
  sigma_y = 1e-2
)

# set the model initial state
model$setInitialState(list(1,1e-1))

# extract the likelihood handlers
nll <- model$likelihood(data=Ornstein)

# calculate likelihood, gradient and hessian w.r.t parameters
nll$fn(nll$par)
nll$gr(nll$par)
nll$he(nll$par)

# estimate the parameters using an extended kalman filter
fit <- model$estimate(data=Ornstein)

# perform moment predictions
pred <- model$predict(data=Ornstein)

# perform stochatic simulations
sim <- model$simulate(data=Ornstein, n.sims=10)

This function creates residual plots for an estimated ctsmTMB object

Description

This function creates residual plots for an estimated ctsmTMB object

Usage

## S3 method for class 'ctsmTMB.fit'
plot(
  x,
  print.plot = 1,
  type = "residuals",
  state.type = "prior",
  against.obs = NULL,
  ggtheme = getggplot2theme(),
  ...
)

Arguments

Value

a (list of) ggplot residual plot(s)

Examples

library(ctsmTMB)
model <- ctsmTMB$new()

# create model
model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw)
model$addObs(y ~ x)
model$setVariance(y ~ sigma_y^2)
model$addInput(u)
model$setParameter(
  theta   = c(initial = 1, lower=1e-5, upper=50),
  mu      = c(initial=1.5, lower=0, upper=5),
  sigma_x = c(initial=1, lower=1e-10, upper=30),
  sigma_y = 1e-2
)
model$setInitialState(list(1,1e-1))

# fit model to data
fit <- model$estimate(Ornstein)

# plot residuals
plot(fit)

# plot filtered states
plot(fit, type="states", against="y")

Plot of k-step predictions from a ctsmTMB prediction object

Description

Plot of k-step predictions from a ctsmTMB prediction object

Usage

## S3 method for class 'ctsmTMB.pred'
plot(
  x,
  y,
  k.ahead = unique(x[["states"]][["k.ahead"]]),
  state.name = NULL,
  type = "states",
  against = NULL,
  ...
)

Arguments

Value

A plot of predicted states

Examples

library(ctsmTMB)
model <- ctsmTMB$new()

# create model
model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw)
model$addObs(y ~ x)
model$setVariance(y ~ sigma_y^2)
model$addInput(u)
model$setParameter(
  theta   = c(initial = 1, lower=1e-5, upper=50),
  mu      = c(initial=1.5, lower=0, upper=5),
  sigma_x = c(initial=1, lower=1e-10, upper=30),
  sigma_y = 1e-2
)
model$setInitialState(list(1,1e-1))

# fit model to data
fit <- model$estimate(Ornstein)

# perform moment predictions
pred <- model$predict(Ornstein)

# plot the k.ahead=10 predictions
plot(pred, against="y.data")


# plot filtered states
plot(fit, type="states", against="y")

#' Plot a profile likelihood ctsmTMB object

Description

#' Plot a profile likelihood ctsmTMB object

Usage

## S3 method for class 'ctsmTMB.profile'
plot(x, y, include.opt = TRUE, ...)

Arguments

Examples

library(ctsmTMB)
model <- ctsmTMB$new()

# create model
model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw)
model$addObs(y ~ x)
model$setVariance(y ~ sigma_y^2)
model$addInput(u)
model$setParameter(
  theta   = c(initial = 1, lower=1e-5, upper=50),
  mu      = c(initial=1.5, lower=0, upper=5),
  sigma_x = c(initial=1, lower=1e-10, upper=30),
  sigma_y = 1e-2
)
model$setInitialState(list(1,1e-1))

# fit model to data
fit <- model$estimate(Ornstein)

# calculate profile likelihood
out <- profile(fit,parlist=list(theta=NULL))

# plot profile
plot(out)

Basic print of ctsmTMB objects

Description

Basic print of ctsmTMB objects

Usage

## S3 method for class 'ctsmTMB'
print(x, ...)

Arguments

Value

Print of ctsmTMB model object

Examples

library(ctsmTMB)
model <- ctsmTMB$new()

# print empty model
print(model)

# add elements to model and see new print
model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw)
model$addObs(y ~ x)
model$setVariance(y ~ sigma_y^2)
model$addInput(u)
model$setParameter(
  theta   = c(initial = 1, lower=1e-5, upper=50),
  mu      = c(initial=1.5, lower=0, upper=5),
  sigma_x = c(initial=1, lower=1e-10, upper=30),
  sigma_y = 1e-2
)
print(model)

Basic print of objects ctsmTMB fit objects

Description

Basic print of objects ctsmTMB fit objects

Usage

## S3 method for class 'ctsmTMB.fit'
print(x, ...)

Arguments

Value

Print of ctsmTMB fit object

Examples

library(ctsmTMB)
model <- ctsmTMB$new()

# create model
model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw)
model$addObs(y ~ x)
model$setVariance(y ~ sigma_y^2)
model$addInput(u)
model$setParameter(
  theta   = c(initial = 1, lower=1e-5, upper=50),
  mu      = c(initial=1.5, lower=0, upper=5),
  sigma_x = c(initial=1, lower=1e-10, upper=30),
  sigma_y = 1e-2
)
model$setInitialState(list(1,1e-1))

# fit model to data
fit <- model$estimate(Ornstein)

# print fit
print(fit)

#' Performs full multi-dimensional profile likelihood calculations

Description

#' Performs full multi-dimensional profile likelihood calculations

Usage

## S3 method for class 'ctsmTMB.fit'
profile(
  fitted,
  parlist,
  grid.size = rep(10, length(parlist)),
  grid.qnt = rep(3, length(parlist)),
  hessian = FALSE,
  silent = FALSE,
  control = list(trace = 0, iter.max = 1000, eval.max = 1000),
  ...
)

Arguments

Note

The implementation was modified from that of https://github.com/kaskr/adcomp/blob/master/TMB/R/tmbprofile.R

Examples

library(ctsmTMB)
model <- ctsmTMB$new()

# create model
model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw)
model$addObs(y ~ x)
model$setVariance(y ~ sigma_y^2)
model$addInput(u)
model$setParameter(
  theta   = c(initial = 1, lower=1e-5, upper=50),
  mu      = c(initial=1.5, lower=0, upper=5),
  sigma_x = c(initial=1, lower=1e-10, upper=30),
  sigma_y = 1e-2
)
model$setInitialState(list(1,1e-1))

# fit model to data
fit <- model$estimate(Ornstein)

# calculate profile likelihood
out <- profile(fit,parlist=list(theta=NULL))

Basic summary of ctsmTMB fit object

Description

Basic summary of ctsmTMB fit object

Usage

## S3 method for class 'ctsmTMB.fit'
summary(object, correlation = FALSE, ...)

Arguments

Value

a summary of the estimated ctsmTMB model fit

Examples

library(ctsmTMB)
model <- ctsmTMB$new()

# create model
model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw)
model$addObs(y ~ x)
model$setVariance(y ~ sigma_y^2)
model$addInput(u)
model$setParameter(
  theta   = c(initial = 1, lower=1e-5, upper=50),
  mu      = c(initial=1.5, lower=0, upper=5),
  sigma_x = c(initial=1, lower=1e-10, upper=30),
  sigma_y = 1e-2
)
model$setInitialState(list(1,1e-1))

# fit model to data
fit <- model$estimate(Ornstein)

# print model summary
summary(fit, correlation=TRUE)