ctsmTMB

CRAN status R-CMD-check

Overview

Welcome to this GitHub repository which hosts the R package ctsmTMB (Continuous Time Stochastic Modelling using Template Model Builder), the intended successor of (and heavily inspired by) the CTSM (Continuous Time Stochastic Modelling) package.

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 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_{t_k} = h\left( t_k, x_{t_k}, u_{t_k}, \theta \right) + \varepsilon_{t}\]

The latent state \(x_t\) is a continuous-time process whose evolution is governed by an Itô stochastic differential equation with drift \(f\) and diffusion \(g\). The state is measured at discrete times \(y_{t_k}\) through the link function \(h\), either directly (in which case \(h\) is the identity map) or indirectly. These measurements are however also contaminated by zero-mean Gaussian noise i.e. \[\varepsilon_{t} \sim \mathcal{N} \left(0, \, \Sigma(t_k, x_{t_k}, u_{t_k}, \theta) \right)\]

We estimate the latent state distribution and the fixed effects parameters \(p\) using maximum likelihood methods based on the collection of discrete-time measurements \(\mathcal{Y}_{k} = \left\{ y_{t_0}, y_{t_1},...,y_{t_k} \right\}\). The package implements various approximate Gaussian methods to perform this inference, specifically Kalman filtering schemes and via Laplace approximations. The inference optimization is carried out using automatic differentiation of the likelihood function for its gradient and hessian, enabled by the use of the TMB/RTMB packages due to Kristensen et. al (2016).

Users interact with the ctsmTMB package via the available methods of the exported R6 ctsmTMB class. The primary methods for defining a state space model are addSystem, addObs and setVariance. These methods takes as input symbolic expressions (R formulas) specifying \(f\), \(g\), \(h\) and \(\Sigma\) respectively. This makes specification relatively easy, but limits the allowed 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 (Kalman filters only).

  3. smooth:

    Performs state smoothing (only available for the Laplace approximation methods).

  4. predict:

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

  5. simulate:

    Performs k-step forecasting through stochastic path simulations (Kalman filters only).

  6. likelihood:

    Construct the (negative log) likelihood function handles for the function value, gradient and hessian. This is directly passed back from RTMB::MakeADFun. This method allows for the use of other optimizers, modifying the likelihood calculations, or combining several data 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 TMB. The additional method laplace.thygesen is an implementation of the stability-improved Laplace approximation for state-dependent diffusion due to Thygesen, 2025.

While not yet supported in the package, a distinct advantage of the Laplace methods is the possibility for (unimodal) non-Gaussian observation densities to accommodate the need for e.g. heavier tails in the underlying distribution.

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:;C:*.

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}\). The observations remain linear and Gaussian:

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

The code chunk below simulates data from this process using an Euler-Maruyama scheme, generates an appropriate ctsmTMB model object, performs parameter estimation using an Extended Kalman Filter (the Linear Kalman Filter method='lkf' could also be used) and inspects the resulting residuals, moment predictions and stochastic simulations.

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)
                      simulation.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.est <- cbind(fit$states$mean$prior, x.sd=fit$states$sd$prior[,"x"])
t <- df.est[,"t"]
x <- df.est[,"x"]
x.sd <- df.est[,"x.sd"]
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(df.obs$t, df.obs$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)

# Plot 10 step predictions vs data
dfp <- pred.10step$states[,c("t.j","x","var.x")]
t <- dfp[,"t.j"]
x <- dfp[,"x"]
x.sd <- sqrt(dfp[,"var.x"])
plot(x=t, y=x, type="n", main="10 Step Predictions 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(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
dfp <- model$predict(df.obs, method="ekf")
xpred <- dfp$states[,"x"]
sdf <- model$simulate(df.obs, method="ekf", n.sims=10)
t <- sdf$times$i0[,"t.j"]
xsim <- sdf$states$x$i0

matplot(t, xsim, type="l", lty="solid", col="grey70", main="Full Prediction/Simulations vs Observations", xlab="Time", ylim=c(0,3))
lines(t, xpred, col="steelblue", lwd=2)
points(df.obs$t, df.obs$y, col="tomato", pch=16, cex=0.7)

# Perform residual analysis
p2 <- plot(fit)

Bibliography