$~$

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.

ctsmTMB offers a user-friendly tool for inference and forecasting in linear and mildly non-linear (multi-dimensional) continuous-discrete stochastic state space systems, i.e. systems on the form

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

That is, 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$. We may estimate the latent state distribution and the fixed effects parameters $p$ of the system through likelihood inference using the discrete-time measurements $\mathcal{Y}{k} = {y{t_0}, y_{t_1},...,y_{t_k}}$. These measurements are related to the latent state through the link function $h$ , but they are also contaminated by zero-mean Gaussian measurement noise $\varepsilon_{t} \sim \mathcal{N}\left(0, \Sigma(t_k, x_{t_k}, u_{t_k}, p) \right)$.

Users interact with the ctsmTMB package via the available methods of the exported R6 ctsmTMB class. The most important of these are addSystem and addObs/setVariance used to specify the functions $f$, $g$, and $h$, $\Sigma$ respectively. The functions are specified indirectly by writing symbolic equation expressions, making it relatively easy to do, but limiting the allowed operations to compositions of simple mathematical functions on scalars.

State and parameter inference and forecasting is carried out using the following available methods:

1. likelihood:

   This method is used to extract the likelihood function handles for the function value, its gradient and hessian, direcly passed back from RTMB::MakeADFun. The is useful for e.g. using other optimizers, modifying the likelihood calculations, or adding together derivatives from several data series.

2. estimate:

   This is the primary work-horse method, which carries out inference by minimizing the (negative log) likelihood using the stats::nlminb quasi-Newton optimizer. The resulting object contains the maximum likelihood parameter and state estimates, and associated marginal uncertainties. See the following section for a list of available inference methods.

3. filter:

   Performs state filtration (Kalman filters only).

4. smoother:

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

5. predict:

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

6. simulate:

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

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.

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

1. The hessian of the likelihood function (w.r.t the fixed parameters) 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.

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 the stability-improved Laplace approximation for state-dependent diffusion due to Thygesen, 2025.

A distinct advantage of the laplace methods over the Kalman filters is the possibility for unimodal non-Gaussian observation densities to accommodate the need for e.g. heavier distribution tails (not yet supported in the package though).

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

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 users should install Command Line Tools. Run the following command in the Terminal

xcode-select --install

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)

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.

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

$$ dx_{t} = \theta \left( \mu + u_t - x_t \right) dt + \sigma_x dB_{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))
df.obs <- data.frame(
  t = t.sim,
  u = u.sim,
  y = NA
)

# set seed for simulate c++. 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

# Now we are ready to see if we can estimate the parameters

############################################################
# 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)

• U. H. Thygesen and K. Kristensen (2025), “Inference in stochastic differential equations using the Laplace approximation: Demonstration and examples”. In: arXiv:2503.21358v2.

• S. Särkkä, “On Unscented Kalman Filtering for State Estimation of Continuous-Time Nonlinear Systems”. In: IEEE Transactions on Automatic Control, 52(9), 1631-1641.