Movatterモバイル変換

Welcome to this GitHub repository which hostsctsmTMB(ContinuousTime Stochastic Modelling using Template Model Builder), theintended successor of, and heavily inspired by, theCTSM package(ContinuousTime Stochastic Modelling). The purpose of the package is tofacilitate a user-friendly tool for (state and parameter) inference, andforecasting, in (multi-dimensional) continuous-discrete stochastic statespace systems, i.e. systems on the form

\[dx_{t} = f(t, x_t, u_t, \theta) \, dt + g(t, x_t, u_t, \theta) \, dB_{t}\]\[y_{t_k} = h(t_k, x_{t_k}, u_{t_k}, \theta)\]

Here the latent state\(x_t\)evolves continuously in time, governed by a set of stochasticdifferential equations, and information about the system is available atdiscrete times through the observations\(y_{t_k}\).

ThectsmTMB package is essentially wrapper aroundtheTMB/RTMB packages(Template Model Builder) thatprovide automatic differentiation of the likelihood function, and accessto other computational tools such as the Laplace approximation. Thelikelihood function is constructed based on the (symbolic) user-providedstate space equations, which may be specified using the implementedOOP-styleR6ctsmTMB class, with methodssuch asaddSystem (for defining system equations), andaddObs (for defining observation equations).

The primary work-horse ofctsmTMB is theestimate method, which carries out inference by minimizingthe (negative log) likelihood using thestats::nlminbquasi-Newton optimizer. The resulting object contains the maximumlikelihood parameter and state estimates, and associated marginaluncertainties. The available inference methods are the Linear andExtended Kalman filters in addition to filtration (actually smoothing)using a Laplace approximation approach.

The package facilities forecasting through thepredict(moment forecasts) andsimulate (stochastic pathsimulations) methods. The calculations for these may be carried out ineitherR (default) or for additional speed inC++ usingRcpp.

Estimation Methods

Kalman Filters

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

The Unscented Kalman Filter implementation is based onAlgorithm4.7 inS.Särkkä, 2007.

Laplace Smoother

The state-reconstructions based on thelaplace(approximation) method aresmoothed estimates, meaning thatstates are optimized jointly conditioned on all observations. TheLaplace approximation is natively built-into and completely handled byTMB. The additional methodlaplace.thygesen is an implementation of the thestability-improved laplace method for systems with state-dependentdiffusion and is due toThygesen, 2025.

Installation

You may experience an issue with a PAT token in which case you cantry r-universe instead

If you encounter problems with a locked folder00LOCK-ctsmTMB due to a previously failed installationremove it with the command below before installating

It is important to note that users must have working C++ compilers toinstall and usectsmTMB.

Windows

Mac / Unix

Mac users should installCommand Line Tools. Run thefollowing command in the Terminal

Test the Installation

Once you have installed the package is it a good idea to test whetherTMB and C++ compilation works. You should be able torun all examples without compilation errors:

Getting Started

You can visit the packagewebpageand browse the vignettes for example uses, in particular seeGettingStarted.

Documentation

CodeExample - Inference in 1D Ornstein-Uhlenbeck Process

library(ggplot2)library(patchwork)library(dplyr)library(reshape2)library(ctsmTMB)############################################################# Data simulation############################################################# Simulate data using Euler Maruyamaset.seed(20)pars=c(theta=10,mu=1,sigma_x=1,sigma_y=0.1)#dt.sim=1e-3t.sim=seq(0,5,by=dt.sim)dw=rnorm(length(t.sim)-1,sd=sqrt(dt.sim))u.sim=cumsum(rnorm(length(t.sim),sd=0.05))x=3for(iin1:(length(t.sim)-1)) {  x[i+1]= x[i]+ pars[1]*(pars[2]-x[i]+u.sim[i])*dt.sim+ pars[3]*dw[i]}# Extract observations and add noisedt.obs=1e-2ids=seq(1,length(t.sim),by=round(dt.obs/ dt.sim))t.obs= t.sim[ids]y= x[ids]+ pars[4]*rnorm(length(t.obs))# forcing inputu= u.sim[ids]# Create data.data=data.frame(t = t.obs,y = y,u = u)############################################################# Model creation and estimation############################################################# Create model objectmodel= ctsmTMB$new()# Add system equationsmodel$addSystem(  dx~ theta* (mu-x+u)* dt+ sigma_x*dw)# Add observation equationsmodel$addObs(  y~ x)# Set observation equation variancesmodel$setVariance(  y~ sigma_y^2)# Add vector inputmodel$addInput(u)# Specify parameter initial values and lower/upper bounds in estimationmodel$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 =c(initial=1e-1,lower=1e-10,upper=30))# Set initial state mean and covariancemodel$setInitialState(list(x[1],1e-1*diag(1)))# Carry out estimation with default settings (extended kalman filter)fit<- model$estimate(data=.data,method="ekf")# Check parameter estimates against truthp0= fit$par.fixedcbind(p0,pars)# Create plot of one-step predictions, simulated states and observationst.est= fit$states$mean$prior$tx.mean= fit$states$mean$prior$xx.sd= fit$states$sd$prior$xplot1=ggplot()+geom_ribbon(aes(x=t.est,ymin=x.mean-2*x.sd,ymax=x.mean+2*x.sd),fill="grey",alpha=0.9)+geom_line(aes(x=t.est, x.mean),col="steelblue",lwd=1)+geom_line(aes(x=t.sim,y=x))+geom_point(aes(x=t.obs,y=y),col="tomato",size=0.5)+labs(title="1-Step State Estimates vs Observations",x="Time",y="")+theme_minimal()# Predict to obtain k-step-ahead predictions to see model forecasting abilitypred.list= model$predict(data=.data,k.ahead=10,method="ekf",)# Create plot all 10-step predictions against datapred= pred.list$statespred10step= pred%>% dplyr::filter(k.ahead==10)plot2=ggplot()+geom_ribbon(aes(x=pred10step$t.j,ymin=pred10step$x-2*sqrt(pred10step$var.x),ymax=pred10step$x+2*sqrt(pred10step$var.x)),fill="grey",alpha=0.9)+geom_line(aes(x=pred10step$t.j,pred10step$x),color="steelblue",lwd=1)+geom_point(aes(x=t.obs,y=y),color="tomato",size=0.5)+labs(title="10 Step Predictions vs Observations",x="Time",y="")+theme_minimal()# Perform prediction ignoring all datapred.list= model$predict(data=.data,method="ekf")# Perform simulation ignoring all datasim.list= model$simulate(data=.data,method="ekf",n.sims=10)# Collapse simulation data for easy use with ggplotsim.df= sim.list$states$x$i0%>%select(!c("i","j","t.i","k.ahead"))%>%  reshape2::melt(.,id.var="t.j")# Plot all simulations and the prediction against observationsplot3=ggplot()+geom_line(data=sim.df,aes(x=t.j,y=value,group=variable),color="grey")+geom_line(aes(x=pred.list$states$t.j,y=pred.list$states$x),color="steelblue")+geom_point(aes(x=t.obs,y=y),color="tomato",size=0.5)+labs(title="No Update Prediction and Simulations vs Observations",x="Time",y="")+theme_minimal()+theme(legend.position ="none")# Create plotp1<- patchwork::wrap_plots(plot1, plot2, plot3,ncol=1)# Create one-step-ahead residual analysis plotp2<-plot(fit)# Wrap both plotspatchwork::wrap_plots(p1,p2[[1]],ncol=2)