ctsmTMB

Welcome to the GitHub repository of the ctsmTMB (Continuous Time Stochastic Modelling using Template Model Builder) R package, the intended successor of the ctsmr (Continuous Time Stochastic Modelling for R) package.

Overview

The purpose of ctsmTMB is to offer a user-friendly tool for inference and forecasting in linear and mildly non-linear continuous-discrete stochastic state space systems, i.e. systems on the form

$$ dx_{t} = f\left( t, x_t, u_t, \theta \right) , dt + g\left( t, x_t, u_t, \theta \right) , d\omega_{t} $$ $$ y_{k} = h\left( t_k, x_{t_k}, u_{t_k}, \theta \right) + \varepsilon_{k} $$

In this formulation the state $x_t$ is a stochastic continous-time process whose evolution is governed by the specified stochastic differential equation characterized by its drift $f$ and diffusion $g$. The associated algebraic equation relates measurements $y_{k}$ obtained at discrete times $t_k$ to the state values $x_{t_k}$ through the link function $h$. The measurement process is subject to zero-mean Gaussian noise via $\varepsilon_{k}$ i.e.

$$ \varepsilon_{k} \sim \mathcal{N} \left(0, , \Sigma(t_k, x_{t_k}, u_{t_k}, \theta) \right) $$

State and parameter inference is based on maximum-likelihood estimation, using Gaussian approximations (Kalman filtering and methods based on the Laplace approximation). The likelihood optimization procedure is improved by the use of automatic differentiation made available via the RTMB package due to Kristensen et. al (2016).

ctsmTMB aims to provide a very user-friendly modelling interface, where users interact with the package via certain exposed methods of the custom R6 ctsmTMB class. A state space model is defined by passing symbolic expressions (formulas, e.g. dx ~ f * dt + g * dw) to the addSystem, addObs and setVariance methods, specifying the functions $f$, $g$, $h$ and $\Sigma$. This allows for relatively easy specification, but limits the supported operations to (compositions of) elementary functions on scalars.

Once a model is created inference and forecasting is performed using the following available methods:

  1. estimate:

    Performs state and parameter inference by minimizing the (negative log) likelihood using the stats::nlminb (quasi-Newton) optimizer.

  2. filter:

    Performs state filtration (only for Kalman filter methods).

  3. smooth:

    Performs state smoothing (only for Laplace approximation methods).

  4. predict:

    Performs k-step forecasting of the system mean and variance (only for Kalman filter methods).

  5. simulate:

    Performs k-step forecasting through stochastic path simulations (only for Kalman filter methods).

  6. likelihood:

    Construct the negative log-likelihood function handles for the function value, gradient and hessian.

The likelihood method returns the object passed directly from RTMB::MakeADFun. Having the raw functon handlers enables the user to choose their own optimizer, modify the likelihood calculations, or combine likelihoods from multiple time-series.

Estimation Methods

The following state reconstruction algorithms are currently available:

  1. The Linear Kalman Filter, lkf.

  2. The Extended Kalman Filter, ekf.

  3. The Unscented Kalman Filter, ukf.

  4. The Laplace Smoothers laplace and laplace.thygesen.

Kalman Filters

The package is currently mostly tailored towards the Kalman filter methods. The advantages of the methods are:

  1. The hessian of the likelihood function (w.r.t the fixed parameters $\theta$) is available.

  2. The model residuals are easier to compute for e.g. model validation.

  3. Multi-step predictions / simulations with state updates are easier to compute.

The Unscented Kalman Filter implementation is based on Algorithm 4.7 in S. Särkkä, 2007.

Laplace Smoother

The state-reconstructions based on the laplace method are smoothed estimates, meaning that states are optimized jointly conditioned on all observations. The Laplace approximation is natively built-into and completely handled by RTMB. The additional method laplace.thygesen is an implementation of the stability-improved Laplace approximation for state-dependent diffusion due to Thygesen, 2025. This method prevents biases in the state estimates towards roots in the diffusion function $g(x_t)$ at intermediate time points where no measurements are available.

The Laplace methods have the advantage that they allow for (unimodal) non-Gaussian observation densities, to accommodate the need for e.g. heavier tails in the underlying distribution, but this feature is not yet supported in the package.

Compiler Installation

NOTE: YOU NEED TO MAKE SURE YOU HAVE A WORKING C++ COMPILER AVAILABLE IN R TO INSTALL THE PACKAGE!!

Windows

C++ compilation in R requires Rtools:

  1. Go to https://cran.r-project.org/bin/windows/Rtools/ and find the latest version.

  2. Go to “Control Panel -> System ->”Advanced” (tab) -> “Environment Variables” -> Highlight “Path” -> “Edit” -> Add to the character string in “Variable Value” the path to your Rtools folder C:\Rtools\bin;C:\Rtools\MinGW\bin.

Mac / Unix

Mac users should install Command Line Tools. Run the following command in the Terminal

xcode-select --install

Package Installation

The package can be installed from source from CRAN:

install.packages("ctsmTMB", type="source")

The development version is available on GitHub and R-universe:

remotes::install_github(repo="phillipbvetter/ctsmTMB")

install.packages('ctsmTMB', repos = c('https://phillipbvetter.r-universe.dev'), type="source")

If you encounter problems with a locked folder 00LOCK-ctsmTMB due to a previously failed installation remove it before reinstalling

# Enter the path on your system
enter.your.path <- "/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/00LOCK-ctsmTMB"
unlink(enter.your.path, recursive = TRUE)

Documentation

You can access the documentation for all methods with

?ctsmTMB

You can also visit the package webpage and browse the vignettes for further details, see in particular Getting Started.

The methods documentation is also available on the package homepage.

Code Example - Inference in 1D Ornstein-Uhlenbeck Process

We consider estimating the parameters of the modified Ornstein-Uhlenbeck process

$$ dx_{t} = \theta \left( \mu + u_t - x_t \right) dt + \sigma_x d\omega_{t} $$

where the stationary mean $\mu$ has been augmented with the addition of a time-varying input $u_{t}$. We observe the state directly with noise i.e.:

$$ y_{k} = x_{t_k} + \varepsilon_{k} \qquad \varepsilon_{k} \sim \mathcal{N}\left(0, \sigma_{y}^2 \right) $$

In the code chunk below we demonstate how to:

  1. create a ctsmTMB model object for the specified model.

  2. generate data from the stochastic process by sampling paths from joint distribution of the SDE using the simulate method via an euler-maruyama scheme.

  3. estimate model parameters and prior/posterior states using an extended kalman filter.

  4. perform post-analysis of the model residuals for model verification.

  5. forecast the fitted model 10-steps ahead, and also generate full moment and stochastic forecasts.

library(ctsmTMB)

############################################################
# Model creation
############################################################

# Create model object
model <- newModel()

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

# Add observation equations
model$addObs(
  y ~ x
)

# Set observation equation variances
model$setVariance(
  y ~ sigma_y^2
)

# Add vector input
model$addInput(u)

# Specify parameter initial values and lower/upper bounds in estimation
model$setParameter(
  theta   = c(initial=2, lower=1e-5, upper=50),
  mu      = c(initial=2, lower=0, upper=5),
  sigma_x = c(initial=1e-2, lower=1e-10, upper=30),
  sigma_y = c(initial=1e-2, lower=1e-10, upper=30)
)

# Set initial state mean and covariance
initial.state <- list(x0=3, p0=0.01*diag(1))
model$setInitialState(initial.state)

############################################################
# Simulate observations
############################################################

# First we create data (no observations)
r.seed <- 20 # this seed for for u.sim rnorm draws
set.seed(r.seed)
true.pars = c(theta=10, mu=1, sigma_x=1, sigma_y=0.1)
dt.sim = 1e-2
t.sim = seq(0, 5, by=dt.sim)
u.sim = cumsum(rnorm(length(t.sim),sd=0.05)) # Create arbitrary input signal
df.obs <- data.frame(
  t = t.sim,
  u = u.sim,
  y = NA
)

# set seeds for simulation. The first seed is brownian motion rng,
# the second is for observations noise rng
cpp.seeds <- c(20,20)
sim <- model$simulate(data=df.obs,
                      # use the true parameters 
                      pars = true.pars,
                      # ekf method is default
                      method="ekf", 
                      # use a smaller stepsize than diff(df.obs$t)
                      sim.timestep = 1e-3,
                      # we only need 1 simulation
                      n.sims = 1,
                      # apply rng seeds
                      cpp.seeds = cpp.seeds)

# store simulated observations in the data.frame
df.obs$y <- sim$observations$y$i0

############################################################
# Model estimation
############################################################

# Carry out estimation with default settings (extended kalman filter)
fit <- model$estimate(df.obs, method="ekf")

# Check parameter estimates against truth
fitted.pars <- fit$par.fixed
cbind(true.pars, fitted.pars, difference=true.pars-fitted.pars)

par(mfrow=c(3,1))
# Plot prior predictions (1-step predictions) against simulation (truth) and observations (data)
df.prior <- fit$states$prior$mean_and_sd
with(cbind(df.prior,y=df.obs$y),
     {
       plot(x=t, y=x, type="n", main="1-Step State Estimates vs Observations", xlab="Time", ylab="",  ylim=c(0,3))
       polygon(c(t, rev(t)), c(x+1.96*x.sd, rev(x-1.96*x.sd)), col="grey70", border=NA)
       lines(t, x, col="steelblue", lwd=2)
       points(t, y, col="tomato", pch=16, cex=0.7)
     }
)

# Predict to obtain 10-step-ahead predictions to see model forecasting ability
pred.10step <- model$predict(df.obs, k.ahead=10, method="ekf", return.k.ahead = 10)
pred.10step$states$mean_and_sd %>% head

# Plot 10 step predictions vs data
with(pred.10step$states$mean_and_sd,{
  plot(x=t.j, y=x, type="n", main="10 Step Predictions vs Observations", xlab="Time", ylab="", ylim=c(0,3))
  polygon(c(t.j, rev(t.j)), c(x+1.96*x.sd, rev(x-1.96*x.sd)), col="grey70", border=NA)
  lines(t.j, x, col="steelblue", lwd=2)
}
)
points(df.obs$t, df.obs$y, col="tomato", pch=16, cex=0.7)

# Compare full (i.e. without updating to data along the way) prediction and simulation
pred <- model$predict(df.obs, method="ekf")
sim <- model$simulate(df.obs, method="ekf", n.sims=10, sim.timestep=1e-3)
xsim <- sim$states$x$i0
x <- pred$states$mean_and_sd$x
# All use the same time vector since all series are full predictions
matplot(df.obs$t, xsim, 
        type="l", lty="solid", col="grey70",  ylim=c(0,3),
        xlab="Time", main="Full Prediction/Simulations vs Observations")
lines(df.obs$t, x, col="steelblue", lwd=2)
points(df.obs$t, df.obs$y, col="tomato", pch=16, cex=0.7)

# Perform residual analysis
plot(fit)

Bibliography