| Type: | Package |
| Title: | Continuous Time Stochastic Modelling using Template Model Builder |
| Version: | 1.1.1 |
| Date: | 2026-07-12 |
| Description: | Perform state and parameter inference, and forecasting, in stochastic state-space systems using the 'ctsmTMB' R6 class. This class provides a user-friendly interface for working with stochastic state space models. Inference is based on maximum likelihood estimation, with derivatives efficiently computed through automatic differentiation enabled by the 'TMB'/'RTMB' packages (Kristensen et al., 2016) <doi:10.18637/jss.v070.i05>. The available inference methods include Kalman filters, in addition to a Laplace approximation-based smoothing method. For further details of these methods refer to the documentation of the 'CTSMR' package https://ctsm.info/ctsmr-reference.pdf and Thygesen (2025) <doi:10.48550/arXiv.2503.21358>. Forecasting capabilities include moment predictions and stochastic path simulations implemented in 'C++' using 'Rcpp' (Eddelbuettel et al., 2018) <doi:10.1080/00031305.2017.1375990> for computational efficiency. |
| License: | GPL-3 |
| URL: | https://github.com/phillipbvetter/ctsmTMB |
| BugReports: | https://github.com/phillipbvetter/ctsmTMB/issues |
| Depends: | R (≥ 4.0.0) |
| Imports: | Deriv, ggplot2, graphics, grDevices, Matrix, patchwork, R6, Rcpp, rlang, RTMB (≥ 1.7), stats, stringr, TMB, utils |
| Suggests: | knitr, rmarkdown, testthat (≥ 3.0.0) |
| LinkingTo: | Rcpp, RcppArmadillo, RcppEigen, zigg |
| VignetteBuilder: | knitr |
| Config/testthat/edition: | 3 |
| Copyright: | See the file COPYRIGHTS |
| Encoding: | UTF-8 |
| LazyData: | true |
| Config/roxygen2/version: | 8.0.0 |
| NeedsCompilation: | yes |
| Packaged: | 2026-07-14 20:52:38 UTC; pbrve |
| Author: | Phillip Vetter [aut, cre, cph], Jan Møller [ctb], Uffe Thygesen [ctb], Peder Bacher [ctb], Henrik Madsen [ctb] |
| Maintainer: | Phillip Vetter <pbrve@dtu.dk> |
| Repository: | CRAN |
| Date/Publication: | 2026-07-14 21:10:09 UTC |
A list of outputs generated from calling the estimate method with the state reconstruction methods
'ekf', 'lkf', 'ukf', 'laplace' and 'laplace.thygesen'.
Description
A list of outputs generated from calling the estimate method with the state reconstruction methods
'ekf', 'lkf', 'ukf', 'laplace' and 'laplace.thygesen'.
Usage
EstimateReferenceData
Format
An object of class list of length 10.
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 = 5, mu = 3, sigma_x = 1, sigma_y = 0.1
The simulation time-step was 1e-3, and observation time-step 1e-1. The simulation was taken from t = 0..20
The simulated input was given by u.sim = cumsum(rnorm(length(t.sim),sd=0.05))
where t.sim is the simulated time vector.
Usage
Ornstein
Format
A data frame of 201 rows and 3 columns. The columns represent the
variables: t (time), y (observation) and u (input).
Sample from a simulated two-state Ornstein-Uhlenbeck process
Description
The data was simulated using a standard Euler-Maruyama method.
The simulated process is governed by the coupled SDEs
dx1 ~ theta * (mu + u - x1) * dt + sigma_x1 * dw1 dx2 ~ alpha * (x1 - x2) * dt + sigma_x2 * dw2
with observations
y1 ~ x1 y2 ~ x2
x2 mean-reverts toward x1, acting as a lagged/smoothed version of x1. The dataset is intended for testing multi-dimensional (> 1 state) inference.
The parameters used for simulation were theta = 5, mu = 3, alpha = 2, sigma_x1 = 1, sigma_x2 = 0.5, sigma_y1 = 0.1, sigma_y2 = 0.1
The simulation time-step was 1e-3, and observation time-step 1e-1. The simulation was taken from t = 0..20
The simulated input was given by u.sim = cumsum(rnorm(length(t.sim), sd=0.05))
where t.sim is the simulated time vector.
Usage
Ornstein2D
Format
A data frame of 201 rows and 4 columns. The columns represent the
variables: t (time), y1 (observation of x1), y2 (observation of x2),
and u (input).
A list of outputs generated from calling the filter, predict and simulate
method on the Ornstein2D data.
Description
A list of outputs generated from calling the filter, predict and simulate
method on the Ornstein2D data.
Usage
OutputReferenceData
Format
An object of class list of length 3.
Create a 1D (1 state, 1 observation) Ornstein-Uhlenbeck model with input-driven mean value. The state is directly observed.
Description
Create a 1D (1 state, 1 observation) Ornstein-Uhlenbeck model with input-driven mean value. The state is directly observed.
Usage
create.Ornstein1D.model()
Value
a ctsmTMB model object
Create a 2D (2 states, 2 observations) Ornstein-Uhlenbeck model with input-driven mean value in the first state and a lagged second state. The states are directly observed.
Description
Create a 2D (2 states, 2 observations) Ornstein-Uhlenbeck model with input-driven mean value in the first state and a lagged second state. The states are directly observed.
Usage
create.Ornstein2D.model()
Value
a ctsmTMB model object
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. The model object has methods
likelihood and estimate for parameter estimation, filter for state reconstruction,
and predict / simulate for forecasting.
Methods
Public methods
ctsmTMB$new()
Initialize private fields
Usage
ctsmTMB$new()
ctsmTMB$.private()
Extract the private fields of a ctsmTMB model object. Primarily used for debugging.
Usage
ctsmTMB$.private()
ctsmTMB$getPrivateFields()
Extract the private fields of a ctsmTMB model object. Primarily used for debugging.
Usage
ctsmTMB$getPrivateFields()
ctsmTMB$getPrivate()
Extract the private fields of a ctsmTMB model object. Primarily used for debugging.
Usage
ctsmTMB$getPrivate()
ctsmTMB$setTrainingMethod()
Set training method
Usage
ctsmTMB$setTrainingMethod(full.prediction)
Arguments
full.predictionboolean whether or not to train against a full prediction over the data, rather than 1-step predictions
ctsmTMB$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
forma formula specifying a stochastic differential equation
...additional formulas similar to
formfor specifying multiple equations at once.
ctsmTMB$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
forma formula specifying an observation equation
...additional formulas similar to
formfor specifying multiple equations at once.obsnamescharacter vector specifying the name of the observation. This is used when the left-hand side of
formconsists of more than just a single variable (of class 'call').
ctsmTMB$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
formformula class specifying the observation equation to be added to the system.
...additional formulas identical to
formto specify multiple observation equations at a time.
ctsmTMB$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
...variable names that specifies the name of input variables in the defined system.
ctsmTMB$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:
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.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
...a named vector or matrix as described above.
ctsmTMB$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
formalgebraic formula
...additional formulas
ctsmTMB$setInitialState()
Declare the initial state values i.e. mean and covariance for the system states.
Usage
ctsmTMB$setInitialState(initial.state)
Arguments
initial.statea named list of two entries 'x0' and 'p0' containing the initial state and covariance of the state.
ctsmTMB$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
'identity' - when the diffusion is state-independent (default)
'log' - when the diffusion is proportional to to x * dw
'logit' - when the diffusion is proportional to x * (1-x) * dw
'sqrt-logit' - when the diffusion is proportional to sqrt(x * (1-x)) * dw
Usage
ctsmTMB$setLamperti(transforms, states = NULL)
Arguments
transformscharacter vector - one of either "identity, "log", "logit", "sqrt-logit"
statesa vector of the state names for which the specified transformations should be applied to.
ctsmTMB$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
namestring defining the model name.
ctsmTMB$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
meanmean vector of the Gaussian prior parameter distribution
covcovariance matrix of the Gaussian prior parameter distribution
ctsmTMB$setAdvancedSettings()
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$setAdvancedSettings( force.ad = TRUE, rtmb.tapeconfig = NULL, tmb.tapeconfig = NULL, nllreport = FALSE )
Arguments
force.ada boolean indicating whether to use state space functions that take advantage of the RTMB::AD(...,force=TRUE) hack which reduces compilation time call to MakeADFun by 20%. This breaks some functionalities such as REPORT.
rtmb.tapeconfigoptions to be passed to TapeConfig.
tmb.tapeconfigoptions to be passed to config.
nllreporta boolean to disable / enable REPORT quantities out of the likelihood function, such as states and residuals.
ctsmTMB$getSystems()
Retrieve system equations.
Usage
ctsmTMB$getSystems()
ctsmTMB$getObservations()
Retrieve observation equations.
Usage
ctsmTMB$getObservations()
ctsmTMB$getVariances()
Retrieve observation variances
Usage
ctsmTMB$getVariances()
ctsmTMB$getAlgebraics()
Retrieve algebraic relations
Usage
ctsmTMB$getAlgebraics()
ctsmTMB$getInitialState()
Retrieve initially set state and covariance
Usage
ctsmTMB$getInitialState()
ctsmTMB$getParameters()
Get initial (and estimated) parameters.
Usage
ctsmTMB$getParameters(type = "all", value = "all")
Arguments
typeone of "all", "free" or "fixed" parameters.
valueone of "all", "initial", "estimate", "lower" or "upper". When not "all", a named numeric vector is returned.
ctsmTMB$getTimers()
Retrieve initially timers
Usage
ctsmTMB$getTimers()
ctsmTMB$getEstimate()
Retrieve initially set state and covariance
Usage
ctsmTMB$getEstimate()
ctsmTMB$getLikelihood()
Retrieve initially set state and covariance
Usage
ctsmTMB$getLikelihood()
ctsmTMB$getPrediction()
Retrieve initially set state and covariance
Usage
ctsmTMB$getPrediction()
ctsmTMB$getSimulation()
Retrieve initially set state and covariance
Usage
ctsmTMB$getSimulation()
ctsmTMB$filter()
Perform state filtering (or smoothing for the 'laplace' method)
Usage
ctsmTMB$filter( data, pars = NULL, method = "ekf", ode.solver = "rk4", ode.timestep = diff(data$t), first.order.input.hold = FALSE, loss = "quadratic", loss_c = NULL, ukf.hyperpars = c(1, 0, 3), initial.state = self$getInitialState(), laplace.residuals = FALSE, estimate.initial.state = FALSE, use.cpp = TRUE, silent = FALSE, ... )
Arguments
datadata.frame containing time-vector 't', observations and inputs. The observations can take
NA-values.parsfixed parameter vector parsed to the objective function for prediction/filtration. The default parameter values used are the initial parameters provided through
setParameter, unless theestimatemethodcharacter vector specifying the filtering method used for state/likelihood calculations. Must be one of either "lkf", "ekf", "laplace".
ode.solverSets the ODE solver used in the Kalman Filter methods for solving the moment differential equations. The default "euler" is the Forward Euler method, alternatively the classical 4th order Runge Kutta method is available via "rk4".
ode.timestepnumeric value. Sets the time step-size in numerical filtering schemes. The defined step-size is used to calculate the number of steps between observation time-points as defined by the provided
data. If the calculated number of steps is larger than N.01 where N is an integer, then the time-step is reduced such that exactly N+1 steps is taken between observations The step-size is used in the two following ways depending on the chosen method:Kalman filters: The time-step is used as the step-size in the numerical Forward-Euler scheme to compute the prior state mean and covariance estimate as the final time solution to the first and second order moment differential equations.
TMB method: The time-step is used as the step-size in the Euler-Maruyama scheme for simulating a sample path of the stochastic differential equation, which serves to link together the latent (random effects) states.
first.order.input.holdlogical which determines whether the moment equations are solved using first-order (TRUE) or zero-order (FALSE) hold on the inputs when integrating between time-points.
losscharacter vector. Sets the loss function type (only implemented for the kalman filter methods). The loss function is per default quadratic in the one-step residuals as is natural when the Gaussian (negative log) likelihood is evaluated, but if the tails of the distribution is considered too small i.e. outliers are weighted too much, then one can choose loss functions that accounts for this. The three available types available:
Quadratic loss (
quadratic).Quadratic-Linear (
huber)Quadratic-Constant (
tukey)
The cutoff for the Huber and Tukey loss functions are determined from a provided cutoff parameter
loss_c. The implementations of these losses are approximations (pseudo-huber and sigmoid approximation respectively) for smooth derivatives.loss_ccutoff value for huber and tukey loss functions. Defaults to
c=3ukf.hyperparsThe hyperparameters alpha, beta, and kappa used for sigma points and weights construction in the Unscented Kalman Filter.
initial.statea named list of two entries 'x0' and 'p0' containing the initial state and covariance of the state
laplace.residualsboolean - whether or not to calculate one-step ahead residuals using the method of oneStepPredict.
estimate.initial.stateboolean value. When TRUE the initial state and covariance matrices are estimated as the stationary solution of the linearized mean and covariance differential equations. When the system contains time-varying inputs, the first element of these is used.
use.cppa boolean to indicate whether to use C++ to perform calculations
silentlogical value whether or not to suppress printed messages such as 'Checking Data', 'Building Model', etc. Default behaviour (FALSE) is to print the messages.
...additional arguments
ctsmTMB$smooth()
Perform state filtering (or smoothing for the 'laplace' method)
Usage
ctsmTMB$smooth( data, pars = NULL, method = "laplace", ode.solver = "euler", ode.timestep = diff(data$t), first.order.input.hold = FALSE, loss = "quadratic", loss_c = NULL, initial.state = self$getInitialState(), laplace.residuals = FALSE, estimate.initial.state = FALSE, silent = FALSE, ... )
Arguments
datadata.frame containing time-vector 't', observations and inputs. The observations can take
NA-values.parsfixed parameter vector parsed to the objective function for prediction/filtration. The default parameter values used are the initial parameters provided through
setParameter, unless theestimatemethodcharacter vector specifying the filtering method used for state/likelihood calculations. Must be one of either "lkf", "ekf", "laplace".
ode.solverSets the ODE solver used in the Kalman Filter methods for solving the moment differential equations. The default "euler" is the Forward Euler method, alternatively the classical 4th order Runge Kutta method is available via "rk4".
ode.timestepnumeric value. Sets the time step-size in numerical filtering schemes. The defined step-size is used to calculate the number of steps between observation time-points as defined by the provided
data. If the calculated number of steps is larger than N.01 where N is an integer, then the time-step is reduced such that exactly N+1 steps is taken between observations The step-size is used in the two following ways depending on the chosen method:Kalman filters: The time-step is used as the step-size in the numerical Forward-Euler scheme to compute the prior state mean and covariance estimate as the final time solution to the first and second order moment differential equations.
TMB method: The time-step is used as the step-size in the Euler-Maruyama scheme for simulating a sample path of the stochastic differential equation, which serves to link together the latent (random effects) states.
first.order.input.holdlogical which determines whether the moment equations are solved using first-order (TRUE) or zero-order (FALSE) hold on the inputs when integrating between time-points.
losscharacter vector. Sets the loss function type (only implemented for the kalman filter methods). The loss function is per default quadratic in the one-step residuals as is natural when the Gaussian (negative log) likelihood is evaluated, but if the tails of the distribution is considered too small i.e. outliers are weighted too much, then one can choose loss functions that accounts for this. The three available types available:
Quadratic loss (
quadratic).Quadratic-Linear (
huber)Quadratic-Constant (
tukey)
The cutoff for the Huber and Tukey loss functions are determined from a provided cutoff parameter
loss_c. The implementations of these losses are approximations (pseudo-huber and sigmoid approximation respectively) for smooth derivatives.loss_ccutoff value for huber and tukey loss functions. Defaults to
c=3initial.statea named list of two entries 'x0' and 'p0' containing the initial state and covariance of the state
laplace.residualsboolean - whether or not to calculate one-step ahead residuals using the method of oneStepPredict.
estimate.initial.stateboolean value. When TRUE the initial state and covariance matrices are estimated as the stationary solution of the linearized mean and covariance differential equations. When the system contains time-varying inputs, the first element of these is used.
silentlogical value whether or not to suppress printed messages such as 'Checking Data', 'Building Model', etc. Default behaviour (FALSE) is to print the messages.
...additional arguments
ctsmTMB$estimate()
Estimate the fixed effects parameters in the specified model.
Usage
ctsmTMB$estimate( data, method = "ekf", ode.solver = "rk4", ode.timestep = diff(data$t), first.order.input.hold = FALSE, loss = "quadratic", loss_c = NULL, ukf.hyperpars = c(1, 0, 3), initial.state = self$getInitialState(), trace = 10, control = list(trace = trace, iter.max = 1e+05, eval.max = 1e+05), use.hessian = FALSE, report = TRUE, laplace.residuals = FALSE, unconstrained.optim = FALSE, estimate.initial.state = FALSE, silent = FALSE, compile = FALSE, ... )
Arguments
dataa
data.framewith a time-column (must be named"t"), observations and inputs used maximum-likelihood parameter and state estimation. The observations can takeNA-values.methoda
characterstring specifying the likelihood method used for parameter and state estimation.ode.solvera
characterstring to determine the ODE solver scheme for the Kalman methods when solving the moment differential equations. The methods are either Forward Euler ("euler") or 4th order Runge-Kutta ("rk4") (default).ode.timestepa
numericvalue that determines the ODE solver time step-size. The passed value is used to calculate the number of steps between time-points indata$tsuch that an integer number of steps are taken. For details see the "Estimation" vignette on the webpage. The step-size is used for solving the moment differential equations (Kalman methods) or for simulating the stochastic SDE path (Laplace methods).first.order.input.holdlogical which determines whether the moment equations are solved using first-order (TRUE) or zero-order (FALSE) hold on the inputs when integrating between time-points.
losscharacter vector. Sets the loss function type (only implemented for the kalman filter methods). The loss function is per default quadratic in the one-step residuals as is natural when the Gaussian (negative log) likelihood is evaluated, but if the tails of the distribution is considered too small i.e. outliers are weighted too much, then one can choose loss functions that accounts for this. The three available types available:
Quadratic loss (
quadratic).Quadratic-Linear (
huber)Quadratic-Constant (
tukey)
The cutoff for the Huber and Tukey loss functions are determined from a provided cutoff parameter
loss_c. The implementations of these losses are approximations (pseudo-huber and sigmoid approximation respectively) for smooth derivatives.loss_ccutoff value for huber and tukey loss functions. Defaults to
c=3ukf.hyperparsThe hyperparameters alpha, beta, and kappa used for sigma points and weights construction in the Unscented Kalman Filter.
initial.statea named list of two entries 'x0' and 'p0' containing the initial state and covariance of the state
traceinteger passed to
controlwhich determines number of steps between each print-out during optimization (use 0 to disable tracing print-outs).controllist of control parameters parsed to
nlminbas itscontrolargument. See?stats::nlminbfor more informationuse.hessianboolean value. The default (
TRUE) causes the optimization algorithmstats::nlminbto use the fixed effects hessian of the (negative log) likelihood when performing the optimization. This feature is only available for the kalman filter methods without any random effects.reportboolean - whether or not to report filtered states, observations and residuals.
laplace.residualsboolean - whether or not to calculate one-step ahead residuals using the method of oneStepPredict.
unconstrained.optimboolean value. When TRUE then the optimization is carried out unconstrained i.e. without any of the parameter bounds specified during
setParameter.estimate.initial.stateboolean value. When TRUE the initial state and covariance matrices are estimated as the stationary solution of the linearized mean and covariance differential equations. When the system contains time-varying inputs, the first element of these is used.
silentlogical value whether or not to suppress printed messages such as 'Checking Data', 'Building Model', etc. Default behaviour (FALSE) is to print the messages.
compileboolean for (re)compiling the objective C++ file, used for methods ending with
_cpp....additional arguments
ctsmTMB$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 = "rk4", ode.timestep = diff(data$t), first.order.input.hold = FALSE, loss = "quadratic", loss_c = NULL, ukf.hyperpars = c(1, 0, 3), initial.state = self$getInitialState(), estimate.initial.state = FALSE, silent = FALSE, compile = FALSE, ... )
Arguments
dataa data.frame containing time-vector 't', observations and inputs. The observations can take
NA-values.methodcharacter vector specifying the filtering method used for state/likelihood calculations. Must be one of either "lkf", "ekf", "laplace".
ode.solverSets the ODE solver used in the Kalman Filter methods for solving the moment differential equations. The default "euler" is the Forward Euler method, alternatively the classical 4th order Runge Kutta method is available via "rk4".
ode.timestepthe time-step used in the filtering schemes. The time-step has two different uses depending on the chosen method.
Kalman Filters: The time-step is used when numerically solving the moment differential equations.
Laplace Approximation: The time-step is used in the Euler-Maruyama simulation scheme for simulating a sample path of the stochastic differential equation, which serves to link together the latent (random effects) states.
The defined step-size is used to calculate the number of steps between observation time-points as defined by the provided
data. If the calculated number of steps is larger than N.01 where N is an integer, then the time-step is reduced such that exactly N+1 steps is taken between observations The step-size is used in the two following ways depending on the chosen method:Kalman filters: The time-step is used as the step-size in the numerical Forward-Euler scheme to compute the prior state mean and covariance estimate as the final time solution to the first and second order moment differential equations.
TMB method: The time-step is used as the step-size in the Euler-Maruyama scheme for simulating a sample path of the stochastic differential equation, which serves to link together the latent (random effects) states.
first.order.input.holdlogical which determines whether the moment equations are solved using first-order (TRUE) or zero-order (FALSE) hold on the inputs when integrating between time-points.
losscharacter vector. Sets the loss function type (only implemented for the kalman filter methods). The loss function is per default quadratic in the one-step residuals as is natural when the Gaussian (negative log) likelihood is evaluated, but if the tails of the distribution is considered too small i.e. outliers are weighted too much, then one can choose loss functions that accounts for this. The three available types available:
Quadratic loss (
quadratic).Quadratic-Linear (
huber)Quadratic-Constant (
tukey)
The cutoff for the Huber and Tukey loss functions are determined from a provided cutoff parameter
loss_c. The implementations of these losses are approximations (pseudo-huber and sigmoid approximation respectively) for smooth derivatives.loss_ccutoff value for huber and tukey loss functions. Defaults to
c=3ukf.hyperparsThe hyperparameters alpha, beta, and kappa used for sigma points and weights construction in the Unscented Kalman Filter.
initial.statea named list of two entries 'x0' and 'p0' containing the initial state and covariance of the state
estimate.initial.stateboolean value. When TRUE the initial state and covariance matrices are estimated as the stationary solution of the linearized mean and covariance differential equations. When the system contains time-varying inputs, the first element of these is used.
silentlogical value whether or not to suppress printed messages such as 'Checking Data', 'Building Model', etc. Default behaviour (FALSE) is to print the messages.
compileboolean for (re)compiling the objective C++ file, used for methods ending with
_cpp....additional arguments
Returns
The method retunrns a list which includes function handles for the likelihood function.
- par
A vector of initial parameter values. This is the expected input size to
fn,grandhe.- fn
Negative log-likelihood function.
- gr
Negative log-likelihood gradient.
- he
Negative log-likelihood hessian.
ctsmTMB$predict()
Perform prediction/filtration to obtain state mean and covariance estimates. The predictions are
obtained by solving the moment equations k.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 = "rk4",
ode.timestep = diff(data$t),
first.order.input.hold = FALSE,
k.ahead = nrow(data) - 1,
return.k.ahead = 0:min(k.ahead, nrow(data) - 1),
return.variance = c("marginal", "none", "covariance", "correlation"),
ukf.hyperpars = c(1, 0, 3),
initial.state = self$getInitialState(),
estimate.initial.state = private$algo.settings$estimate.initial,
use.cpp = TRUE,
silent = FALSE,
...
)
Arguments
datadata.frame containing time-vector 't', observations and inputs. The observations can take
NA-values.parsfixed parameter vector parsed to the objective function for prediction/filtration. The default parameter values used are the initial parameters provided through
setParameter, unless theestimatefunction has been run, then the default values will be those at the found optimum.methodThe prediction method
ode.solverSets the ODE solver used in the Kalman Filter methods for solving the moment differential equations. The default "euler" is the Forward Euler method, alternatively the classical 4th order Runge Kutta method is available via "rk4".
ode.timestepnumeric value. Sets the time step-size in numerical filtering schemes. The defined step-size is used to calculate the number of steps between observation time-points as defined by the provided
data. If the calculated number of steps is larger than N.01 where N is an integer, then the time-step is reduced such that exactly N+1 steps is taken between observations The step-size is used in the two following ways depending on the chosen method:Kalman filters: The time-step is used as the step-size in the numerical Forward-Euler scheme to compute the prior state mean and covariance estimate as the final time solution to the first and second order moment differential equations.
TMB method: The time-step is used as the step-size in the Euler-Maruyama scheme for simulating a sample path of the stochastic differential equation, which serves to link together the latent (random effects) states.
first.order.input.holdlogical which determines whether the moment equations are solved using first-order (TRUE) or zero-order (FALSE) hold on the inputs when integrating between time-points.
k.aheadinteger specifying the desired number of time-steps (as determined by the provided data time-vector) for which predictions are made (integrating the moment ODEs forward in time without data updates).
return.k.aheadnumeric vector of integers specifying which k.ahead predictions to that should be returned.
return.variancea string to indicate what kind of dispersions should be reported.
ukf.hyperparsThe hyperparameters alpha, beta, and kappa used for sigma points and weights construction in the Unscented Kalman Filter.
initial.statea named list of two entries 'x0' and 'p0' containing the initial state and covariance of the state
estimate.initial.statebool - stationary estimation of initial mean and covariance
use.cppa boolean to indicate whether to use C++ to perform calculations
silentlogical value whether or not to suppress printed messages such as 'Checking Data', 'Building Model', etc. Default behaviour (FALSE) is to print the messages.
...additional 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,...,k.ahead). If return.dispersion = TRUE then the state covariance/correlation
matrix is returned, otherwise only the marginal variances are returned.
ctsmTMB$simulate()
Perform prediction/filtration to obtain state mean and covariance estimates. The predictions are
obtained by solving the moment equations k.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 = TRUE,
cpp.seeds = NULL,
method = c("ekf", "lkf", "ukf", "laplace", "laplace.thygesen"),
ode.solver = "rk4",
ode.timestep = diff(data$t),
first.order.input.hold = FALSE,
simulation.timestep = diff(data$t),
k.ahead = nrow(data) - 1,
return.k.ahead = 0:min(k.ahead, nrow(data) - 1),
n.sims = 100,
ukf.hyperpars = c(1, 0, 3),
initial.state = self$getInitialState(),
estimate.initial.state = private$algo.settings$estimate.initial,
silent = FALSE,
...
)
Arguments
datadata.frame containing time-vector 't', observations and inputs. The observations can take
NA-values.parsfixed parameter vector parsed to the objective function for prediction/filtration. The default parameter values used are the initial parameters provided through
setParameter, unless theestimatefunction has been run, then the default values will be those at the found optimum.use.cppa boolean to indicate whether to use C++ to perform calculations
cpp.seedsan integer seed value to control RNG normal draws on the C++ side.
method-
The natural TMB-style formulation where latent states are considered random effects and are integrated out using the Laplace approximation. This method only yields the gradient of the (negative log) likelihood function with respect to the fixed effects for optimization. The method is slower although probably has some precision advantages, and allows for non-Gaussian observation noise (not yet implemented). One-step / K-step residuals are not yet available in the package.
(Continuous-Discrete) Extended Kalman Filter where the system dynamics are linearized to handle potential non-linearities. This is computationally the fastest method.
(Continuous-Discrete) Unscented Kalman Filter. This is a higher order non-linear Kalman Filter which improves the mean and covariance estimates when the system display high nonlinearity, and circumvents the necessity to compute the Jacobian of the drift and observation functions.
All package features are currently available for the kalman filters, while TMB is limited to parameter estimation. In particular, it is straight-forward to obtain k-step-ahead predictions with these methods (use the
predictS3 method), and stochastic simulation is also available in the cases where long prediction horizons are sought, where the normality assumption will be inaccurate ode.solverSets the ODE solver used in the Kalman Filter methods for solving the moment differential equations. The default "euler" is the Forward Euler method, alternatively the classical 4th order Runge Kutta method is available via "rk4".
ode.timestepnumeric value. Sets the time step-size in numerical filtering schemes. The defined step-size is used to calculate the number of steps between observation time-points as defined by the provided
data. If the calculated number of steps is larger than N.01 where N is an integer, then the time-step is reduced such that exactly N+1 steps is taken between observations The step-size is used in the two following ways depending on the chosen method:Kalman filters: The time-step is used as the step-size in the numerical Forward-Euler scheme to compute the prior state mean and covariance estimate as the final time solution to the first and second order moment differential equations.
TMB method: The time-step is used as the step-size in the Euler-Maruyama scheme for simulating a sample path of the stochastic differential equation, which serves to link together the latent (random effects) states.
first.order.input.holdlogical which determines whether the moment equations are solved using first-order (TRUE) or zero-order (FALSE) hold on the inputs when integrating between time-points.
simulation.timesteptimestep used in the euler-maruyama scheme
k.aheadinteger specifying the desired number of time-steps (as determined by the provided data time-vector) for which predictions are made (integrating the moment ODEs forward in time without data updates).
return.k.aheadnumeric vector of integers specifying which k.ahead predictions to that should be returned.
n.simsnumber of simulations
ukf.hyperparsThe hyperparameters alpha, beta, and kappa used for sigma points and weights construction in the Unscented Kalman Filter.
initial.statea named list of two entries 'x0' and 'p0' containing the initial state and covariance of the state
estimate.initial.statebool - stationary estimation of initial mean and covariance
silentlogical value whether or not to suppress printed messages such as 'Checking Data', 'Building Model', etc. Default behaviour (FALSE) is to print the messages.
...additional arguments
return.dispersionboolean value to indicate whether the covariance (instead of the correlation) should be returned.
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,...,k.ahead). If return.dispersion = TRUE then the state covariance/correlation
matrix is returned, otherwise only the marginal variances are returned.
ctsmTMB$print()
Function to print the model object
Usage
ctsmTMB$print()
ctsmTMB$clone()
The objects of this class are cloneable with this method.
Usage
ctsmTMB$clone(deep = FALSE)
Arguments
deepWhether to make a deep clone.
Create a ctsmTMB model faster avoiding $...
Description
Create a ctsmTMB model faster avoiding $...
Usage
newModel()
Value
Print of ctsmTMB model object
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(),
ylims = c(NA, NA),
standardized.residuals = TRUE,
residual.burnin = 0L,
residual.vs.obs.and.inputs = FALSE,
...
)
Arguments
x |
A R6 ctsmTMB fit object |
print.plot |
a single integer determining which element out of all
states/observations (depending on the argument to |
type |
a character vector either 'residuals' or 'states' determining what to plot. |
state.type |
a character vector either 'prior', 'posterior' or 'smoothed' determining what kind of states to plot. |
against.obs |
name of an observation to plot state predictions against. |
ggtheme |
ggplot2 theme to use for creating the ggplot. |
ylims |
limits on the y-axis for residual time-series plot |
standardized.residuals |
logical whether or not the time-series residual plot uses standardized residuals or not. |
residual.burnin |
integer N to remove the first N residuals |
residual.vs.obs.and.inputs |
the residual plots also include a new window with time-series plots of residuals, associated observations and inputs |
... |
additional arguments |
Value
a (list of) ggplot residual plot(s)
Examples
## Not run:
# Given a fit from estimate i.e.
fit <- model$estimate(data)
# We can generate a standard residual analysis plot by
plot(fit)
# We can remove the first say 10 residuals of the series if they are off due
# to a bad initial prior state
plot(fit, residual.burnin=10)
# We can also plot residuals against inputs and observations instead
plot(fit, residual.vs.obs.and.inputs=TRUE)
# We can also plot the prior or posterior states using the 'type' and
# 'state.type' arguments.
plot(fit, type="states", state.type="prior")
# We can also show observations on that plot
plot(fit, type="states", state.type="prior", against.obs="y")
## End(Not run)
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.obs = NULL,
...
)
Arguments
x |
a ctsmTMB.pred object |
y |
not used |
k.ahead |
an integer indicating which k-ahead predictions to plot |
state.name |
a string indicating which states to plot |
type |
one of 'states' or 'observations', to plot |
against.obs |
name of an observations to plot predictions against |
... |
additional arguments |
Value
A plot of predicted states
Examples
## Not run:
# Given a prediction
pred <- model$predict(data, k.ahead=10)
# We can plot the all k.head predictions in the same plot
plot(pred)
# This is usually messy. We can instead set a specific horizon
plot(pred, k.ahead=10)
# We can select among all states using the 'state.name' argument
plot(pred, k.ahead=10, state.name="x")
# We can also plot against some observations using the 'against' argument
plot(pred, k.ahead=10, state.name="x", against="y")
# We can also show observations on that plot
plot(fit, type="states", state.type="prior", against.obs="y")
## End(Not run)
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
x |
a profile.ctsmTMB object |
y |
not in use |
include.opt |
boolean which indicates whether or not to include the total likelihood optimizer in the plot. |
... |
additional arguments |
Examples
## Not run:
# Given a model fit produced by
fit <- model$estimate(Ornstein)
# We calculate the profile likelihood
prof <- profile(fit, parlist=list(theta=NULL))
# And then simply call plot on that object
plot(prof)
## End(Not run)
Basic print of ctsmTMB objects
Description
Basic print of ctsmTMB objects
Usage
## S3 method for class 'ctsmTMB'
print(x, ...)
Arguments
x |
an object of class 'ctsmTMB' |
... |
additional arguments (not in use) |
Value
Print of ctsmTMB model object
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
x |
a ctsmTMB fit object |
... |
additional arguments |
Value
Print of ctsmTMB fit object
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
fitted |
a ctsmTMB fit object |
parlist |
a named-list of parameters to profile over. The user can either supply grid-values in the list or leave it empty. If the any one list is empty then grid-values will be calculated using the estimated parameter mean value and standard deviation. |
grid.size |
a vector of |
grid.qnt |
a vector of |
hessian |
a boolean indicating whether to use the hessian or not during the profile optimization. |
silent |
boolean whether or not to mute current iteration number
the |
control |
a list of optimization output controls (see nlminb) |
... |
various arguments (not in use) |
Note
The implementation was modified from that of https://github.com/kaskr/adcomp/blob/master/TMB/R/tmbprofile.R
Examples
## Not run:
# Given a model fit produced by
fit <- model$estimate(Ornstein)
# We call profile with theta=NULL to use the standard range theta +- 3 std. devs.
prof <- profile(fit, parlist=list(theta=NULL))
## End(Not run)
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
object |
a ctsmTMB fit object |
correlation |
boolean indicating whether or not to display the parameter correlation structure |
... |
additional arguments |
Value
a summary of the estimated ctsmTMB model fit