{'fraction_non_alphanumeric': 0.041910092705964665, 'fraction_numerical': 0.011509532971838377, 'mean_word_length': 4.929475212611491, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
bssm: Bayesian Inference of Non-linear and Non-Gaussian State Space Models in R

We present an R package bssm for Bayesian non-linear/non-Gaussian state space modelling. Unlike the existing packages, bssm allows for easy-to-use approximate inference based on Gaussian approximations such as the Laplace approximation and the extended Kalman filter. The package accommodates also discretely observed latent diffusion processes. The inference is based on fully automatic, adaptive Markov chain Monte Carlo (MCMC) on the hyperparameters, with optional importance sampling post-correction to eliminate any approximation bias. The package implements also a direct pseudo-marginal MCMC and a delayed acceptance pseudo-marginal MCMC using intermediate approximations. The package offers an easy-to-use interface to define models with linear-Gaussian state dynamics with non-Gaussian observation models, and has an Rcpp interface for specifying custom non-linear and diffusion models.


Introduction
State space models (SSM) are a flexible class of latent variable models commonly used in analysing time series data (cf. Durbin and Koopman, 2012). There are a number of packages available for state space modelling for R, especially for two special cases: a linear-Gaussian SSM (LGSSM) where both the observation and state densities are Gaussian with linear relationships with the states, and an SSM with discrete state space, which is sometimes called a hidden Markov model (HMM). These classes admit analytically tractable marginal likelihood functions and conditional state distributions (conditioned on the observations), making inference relatively straightforward. See for example (Petris and Petrone, 2011;Tusell, 2011;Helske, 2017;Helske and Helske, 2019) for review of some of the R packages dealing with these type of models. The present R package bssm is designed for Bayesian inference of general state space models with non-Gaussian and/or non-linear observational and state equations. The package primary aim is to provide easy-to-use and fast functions for fully Bayesian inference with common time series models such as basic structural time series model (Harvey, 1989) with exogenous covariates and simple stochastic volatility models. The package accomodates also custom non-linear models and discretised diffusion models.
When extending the state space modelling to non-linear or non-Gaussian models, some difficulties arise. As the likelihood is no longer analytically tractable, computing the latent state distributions, as well as hyperparameter estimation of the model becomes more challenging. One general option is to use Markov chain Monte Carlo (MCMC) methods targeting the full joint posterior of hyperparameters and the latent states, for example by Gibbs sampling or Hamiltonian Monte Carlo. Unfortunately, the joint posterior is typically very high dimensional and due to the strong autocorrelation structures of the state densities, the efficiency of such methods can be relatively poor. Another asymptotically exact approach is based on the pseudo-marginal particle MCMC approach (Andrieu et al., 2010), where the likelihood function and the state distributions are estimated using sequential Monte Carlo (SMC) i.e. the particle filter (PF) algorithm. Instead of computationally demanding Monte Carlo methods, approximation-based methods such extended and unscented Kalman filters may be used, as well as Laplace approximations, which are provided for example by the INLA (Lindgren and Rue, 2015) R package. The latter are computationally appealing, but may lead to hard-to-quantify biases of the posterior.
Some of the R packages suitable for Bayesian state space modelling include pomp (King et al., 2016), rbi (Jacob and Funk, 2020), nimbleSMC (Michaud et al., 2020;NIMBLE Development Team, 2020), and rstan (Stan Development Team, 2020). With the package pomp, user defines the model using R or C snippets for simulation from and evaluation of the latent state and observation level densities, allowing flexible model construction. The rbi package is an interface to LibBi (Murray, 2015), a standalone software with a focus on Bayesian state space modelling on high-performance computers.
The pomp package provides several simulation-based inference methods mainly based on iterated filtering and maximum likelihood, whereas rbi is typically used for Bayesian inference via particle MCMC. For a more detailed comparison of differences of rbi/LibBi and pomp with examples, see (Funk and King, 2020). The nimbleSMC package contains some particle filtering algorithms which can be used in the general Nimble modelling system (de Valpine et al., 2017), whereas the rstan package provides an R interface to the Stan C++ package, a general statistical modelling platform (Carpenter et al., 2017).
The key difference to the aforementioned packages and motivation behind the present bssm package is to combine the use of fast approximation-based methods with Monte Carlo correction step, leading to computationally efficient and unbiased (approximation error free) inference of the joint posterior of hyperparameters and latent states, as suggested in . In a nutshell, the method uses MCMC which targets an approximate marginal posterior of the hyperparameters, and an importance sampling type weighting which provides asymptotically exact inference on the joint posterior of hyperparameters and the latent states. In addition to this two-stage procedure, the bssm supports also delayed acceptance pseudo-marginal MCMC (Christen and Fox, 2005) using the approximations, and direct pseudo-marginal MCMC. To our knowledge, importance sampling and delayed acceptance in this form are not available in other Bayesian state space modelling packages in R.

Supported models
We denote the sequence of observations (y 1 , . . . , y T ) as y, and the sequence of latent state variables (α 1 , . . . , α T ) as α. The latent states α t ∈ R d are typically vector-valued, whereas we focus mainly on scalar observations y t ∈ R (vector-valued observations are also supported, assuming conditional independence (given α t ) in case of non-Gaussian observations).
A general state space model consists of two parts: observation level densities g (θ) t (y t |α t ) and latent state transition densities µ (θ) t (α t+1 |α t ). Typically both g (θ) t and µ (θ) t depend on unknown parameter vector θ for which we can define arbitrary prior p(θ).
In a linear-Gaussian SSM, both g (θ) t and µ (θ) t are Gaussian densities and they depend linearly on the current and previous state vectors, respectively. Section Models with linear-Gaussian state dynamics describes a common extension to these models supported by bssm, which relaxes the assumptions on observational density g (θ) t , by allowing exponential family links, and stochastic volatility models. While the main focus of bssm is in state space models with linear-Gaussian state dynamics, there is also support for more general non-linear models, discussed briefly in Section Other state space models. Section Using the bssm package describes how arbitrary models based on these definitions are constructed in bssm.

Models with linear-Gaussian state dynamics
The primary class of models supported by bssm consists of SSMs with linear-Gaussian state dynamics of form where c t ∈ R d , T t ∈ R d×d , and R t ∈ R d×k can depend on the unknown parameters θ and covariates. The noise terms η t ∼ N(0, I k ) and α 1 ∼ N(a 1 , P 1 ) are independent. These state dynamics can be combined with the observational level density g t of form where parameters φ and the known vector u t are distribution specific and can be omitted in some cases. Currently, following observational level distributions are supported: , where u t is the known exposure at time t.
• Binomial distribution: g t (y t |d t + Z t α t , u t ) = B(u t , logit −1 (d t + Z t α t )), where u t is the number of trials and logit −1 (d t + Z t α t ) is the probability of the success. • Negative binomial distribution: is the expected value, φ is the dispersion parameter, and u t is a known offset term.
For multivariate models, these distributions can be combined arbitrarily, except the stochastic volatility model case which is currently handled separately. Also for fully Gaussian model, the observational level errors t can be correlated across time series.

Other state space models
The general non-linear Gaussian model in the bssm has following form: The bssm package also supports models where the state equation is defined as a continuous-time diffusion model of the form where B t is a Brownian motion and where µ and σ are scalar-valued functions, with the univariate observation density p(y k |α k ) defined at integer times k = 1 . . . , n.

Inference methods
The main goal of bssm is to facilitate easy-to-use full Bayesian inference of the joint posterior p(α, θ|y) for models discussed in Section Supported models. The inference methods implemented in bssm are based on a factorised approach where the joint posterior of hyperparameters θ and latent states α is given as where p(y|θ) is the parameter marginal likelihood and p(α|y, θ) is the smoothing distribution.
All the inference algorithms are based on a Markov chain Monte Carlo on the parameters θ, whose single iteration may be summarised as follows: Algorithm 1 One iteration of MCMC algorithm for sampling p(θ|y).

Direct inference: marginal algorithm and particle MCMC
The simplest case is with a linear-Gaussian SSM, where we can use the exact marginal likelihood p(y|θ) = p(y|θ), in which case Algorithm 1 reduces to (an adaptive) random-walk Metropolis algorithm targeting the posterior marginal of the parameters θ. Inference from the full posterior may be done using the simulation smoothing algorithm (Durbin and Koopman, 2002) conditional to the sampled hyperparameters.
The other 'direct' option, which can be used with any model, is using the bootstrap particle filter (BSF) (Gordon et al., 1993), which leads to a randomp(y|θ) which is an unbiased estimator of p(y|θ). In this case, Algorithm 1 reduces to (an adaptive) particle marginal Metropolis-Hastings (Andrieu et al., 2010). Full posterior inference is achieved simultaneously, by picking particle trajectories based on their ancestries as in the filter-smoother algorithm (Kitagawa, 1996). Note that with BSF, the desired acceptance rate needs to be lower, depending on the number of particles used (Doucet et al., 2015).

Approximate inference: Laplace approximation and the extended Kalman filter
The direct BSF discussed above may be used with any non-linear and/or non-Gaussian model, but may be slow and/or poor mixing. To alleviate this, bssm provides computationally efficient (intermediate) approximate inference in case of non-Gaussian observation models in Section Models with linear-Gaussian state dynamics, and in case of non-linear dynamics in Section Other state space models.
With non-Gaussian models of Section Models with linear-Gaussian state dynamics, we use an approximating Gaussian modelp(y, α|θ) which is a Laplace approximation of p(α, y|θ) following (Durbin and Koopman, 2000). We write the likelihood as follows wherep(y|θ) is the likelihood of the Laplace approximation and the expectation is taken with respect to its conditionalp(α|y, θ) (Durbin and Koopman, 2012). Indeed, denotingα as the mode ofp(α|θ, y), we may write Ifp resembles p with typical values of α, the latter logarithm of expectation is zero. We takep(y|θ) as the expression on the right, dropping the expectation.
Whenp is approximate, the MCMC algorithm targets an approximate posterior marginal. Approximate full inference may be done analogously as in Section Direct inference: marginal algorithm and particle MCMC, by simulating trajectories conditional to the sampled parameter configurations θ i . We believe that approximate inference is often good enough for model development, but strongly recommend using post-correction as discussed in Section Post-processing by importance weighting to check the validity of the final inference.
In addition to these algorithms, bssm also supportsp(y|θ) based on the extended KF (EKF) or iterated EKF (IEKF) (Jazwinski, 1970) which can be used for models with non-linear dynamics (Section Other state space models). Approximate smoothing based on (iterated) EKF is also supported. It is also possible to perform direct inference as in Section Direct inference: marginal algorithm and particle MCMC, but instead of the BSF, employ particle filter based on EKF (Van Der Merwe et al., 2001).

Post-processing by importance weighting
The inference methods in Section Approximate inference: Laplace approximation and the extended Kalman filter are computationally efficient, but come with a bias. The bssm implements importancesampling type post-correction as discussed in . Indeed, having MCMC samples (θ i ) from the approximate posterior constructed as in Section Approximate inference: Laplace approximation and the extended Kalman filter, we may produce (random) weights and latent states, such that the weighted samples form estimators which are consistent with respect to the true posterior p(α, θ|y).
The primary approach which we recommend for post-correction is based on a "ψ-APF' ' -a particle filter using the Gaussian approximations of Section Approximate inference: Laplace approximation and the extended Kalman filter. In essence, this particle filter employs the dynamics and a look-ahead strategy coming from the approximation, which leads to low-variance estimators; see  and package vignettes 1 for more detailed description. Naturally ψ-APF can also be used in place of BSF in direct inference of Section Direct inference: marginal algorithm and particle MCMC.

Direct inference using approximation-based delayed acceptance
An alternative to approximate MCMC and post-correction, bssm also supports an analogous delayed acceptance method (Christen and Fox, 2005;Banterle et al., 2019) (here denoted by DA-MCMC). This algorithm is similar to 1, but in case of 'acceptance', leads to second-stage acceptance using the same weights as the post-correction would; see  for details. Note that as in direct approach for non-Gaussian/non-linear models, the desired acceptance rate with DA-MCMC should be lower than the default 0.234.
The DA-MCMC also leads to consistent posterior estimators, and often outperforms the direct particle marginal Metropolis-Hastings. However, empirical findings  and theoretical considerations  suggest that approximate inference with post-correction may often be preferable. The bssm supports parallelisation with post-correction using OpenMP, which may further promote the latter.

Inference with diffusion state dynamics
For general continuous-time diffusion models, the transition densities are intractable. The bssm uses Millstein time-discretisation scheme for approximate simulation, and inference is based on the corresponding BSF. Fine time-discretisation mesh gives less bias than the coarser one, with increased computational complexity. The DA and IS approaches can be used to speed up the inference by using coarse discretisation in the first stage and then using more fine mesh in the second stage. For comparison of DA and IS approaches in case of geometric Brownian motion model, see .

Using the bssm package
Main functions of bssm related to the MCMC sampling, approximations, and particle filtering are written in C++, with help of Rcpp (Eddelbuettel and François, 2011) and RcppArmadillo (Eddelbuettel and Sanderson, 2014) packages. On the R side, the package uses S3 methods to provide a relatively unified workflow independent of the type of the model one is working with. The model building functions such as bsm_ng and svm are used to construct the model objects of same name which can be then passed to other methods, such as logLik and run_mcmc which compute the log-likelihood value and run MCMC algorithm respectively. We will now briefly describe the main functionality of bssm. For more detailed descriptions of different functions and their arguments, see the corresponding documentation in R and the package vignettes.

Constructing the model
For models with linear-Gaussian state dynamics, bssm includes some predefined models such as bsm_lg and bsm_ng for univariate Gaussian and non-Gaussian structural time series models with external covariates, for which user only needs to supply the data and priors for unknown model parameters. In addition, bssm supports general model building functions ssm_ulg, ssm_mlg for custom univariate and multivariate Gaussian models and ssm_ung, and ssm_mng for their non-Gaussian counterparts. For these models, users need to supply their own R functions for the evaluation of the log prior density and for updating the model matrices given the current value of the parameter vector θ. It is also possible to avoid defining the matrices manually by leveraging the formula interface of the KFAS package (Helske, 2017) together with as_bssm function which converts KFAS model to a bssm equivalent model object. This is especially useful in case of complex multivariate models with covariates.
As an example, consider a Gaussian local linear trend model of the form with zero-mean Gaussian noise terms t , η t , ξ t with unknown standard deviations. Using the time series of the mean annual temperature (in Fahrenheit) in New Haven, Connecticut, from 1912 to 1971 (available in the datasets package) as an example, this model can be built with bsm function as library("bssm") data("nhtemp", package = "datasets") prior <-halfnormal(1, 10) bsm_model <-bsm_lg(y = nhtemp, sd_y = prior, sd_level = prior, sd_slope = prior) Here we use helper function halfnormal which defines half-Normal prior distribution for the standard deviation parameters, with the first argument defining the initial value of the parameter, and second defines the scale parameter of the half-Normal distribution. Other prior options are normal, tnormal (truncated normal), gamma, and uniform.
For models where the state equation is no longer linear-Gaussian, we use pointer-based interface by defining all model components as well as functions defining the Jacobians of Z(·) and T(·) needed by the extended Kalman filter as C++ snippets. General non-linear Gaussian model can be defined with the function ssm_nlg. Discretely observed diffusion models where the state process is assumed to be continuous stochastic process can be constructed using the ssm_sde function, which takes pointers to C++ functions defining the drift, diffusion, the derivative of the diffusion function, and the log-densities of the observations and the prior. As an example of the latter, let us consider an Ornstein-Uhlenbeck process with parameters θ = (φ, ν, σ) = (0.5, 2, 1) and the initial condition α 0 = 1. For observation density, we use Poisson distribution with parameter exp(α k ). We first simulate a trajectory x 0 , . . . , x n using the sde.sim function from the sde package (Iacus, 2016) and use that for the simulation of observations y: library("sde") x <-sde.sim(t0 = 0, T = 100, X0 = 1, N = 100, drift = expression(0.5 * (2 -x)), sigma = expression(1), sigma.x = expression(0)) y <-rpois(100, exp(x[-1])) We then compile and build the model as Rcpp::sourceCpp("ssm_sde_template.cpp") pntrs <-create_xptrs() sde_model <-ssm_sde(y, pntrs$drift, pntrs$diffusion, pntrs$ddiffusion, pntrs$obs_density, pntrs$prior, c(0.5, 2, 1), 1, FALSE) The templates for the C++ functions for SDE and non-linear Gaussian models can be found from the package vignettes on the CRAN 2 .

Markov chain Monte Carlo in bssm
The main purpose of the bssm is to allow computationally efficient MCMC-based inference for various state space models. For this task, a method run_mcmc can be used. The function takes a number of arguments, depending on the model class, but for many of these, default values are provided. For linear-Gaussian models, we only need to supply the number of iterations. Using the previously created local linear trend model for the New Haven temperature data of Section Constructing the model, we run an MCMC with 100,000 iterations where first 10,000 is discarded as a burn-in (burn-in phase is also used for the adaptation of the proposal distribution): mcmc_bsm <-run_mcmc(bsm_model, iter = 1e5, burnin = 1e4) The print method for the output of the MCMC algorithms gives a summary of the results, and detailed summaries for θ and α can be obtained using summary function. For all MCMC algorithms, bssm uses so-called jump chain representation of the Markov chain X 1 , . . . , X n , where we only store each accepted X k and the number of steps we stayed on the same state. So for example if X 1:n = (1, 2, 2, 1, 1, 1), we present such chain asX = (1, 2, 1), N = (1, 2, 3). This approach reduces the storage space and makes it more computationally efficient to use importance sampling type correction algorithms. One drawback of this approach is that the results from the MCMC run correspond to weighted samples from the target posterior, so some of the commonly used postprocessing tools need to be adjusted. Of course, in case of other methods than IS-weighting, the simplest option is to just expand the samples to typical Markov chain using the stored counts N. This can be done using the function expand_sample which returns an object of class "mcmc" of the coda package (Plummer et al., 2006) (thus the plotting and diagnostic methods of coda can also be used). We can also directly transform the posterior samples to a "data.frame" object by using as.data.frame method for the MCMC output (for IS-weighting, the returned data frame contains additional column weights). This is useful for example for visualization purposes with the ggplot2 (Wickham, 2016) package: library("ggplot2") d <-as.data.frame(mcmc_bsm, variable = "theta") ggplot(d, aes(x = value)) + geom_density(bw = 0.1, fill = "#9ebcda") + facet_wrap(~variable, scales = "free") + theme_bw() Figure 1 shows the estimated posterior densities of the three standard deviation parameters of the model. The relatively large observational level standard deviation σ y suggests that the underlying latent temperature series is much smoother than the observed series, which can also be seen from d <-as.data.frame(mcmc_bsm, variable = "states") summary_y <-d %>% filter(variable == "level") %>% group_by(time) %>% summarise(mean = mean(value), lwr = quantile(value, 0.025), upr = quantile(value, 0.975)) ggplot(summary_y, aes(x = time, y = mean)) + geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.25) + geom_line() + geom_point(data = data.frame(mean = nhtemp, time = time(nhtemp))) + theme_bw() + xlab("Year") + ylab("Mean annual temperature in New Haven") For non-Gaussian models the default MCMC algorithm is approximate inference based on Laplace approximation combined with importance sampling post-correction (Section Post-processing by importance weighting). It is also possible to perform first approximate MCMC using the argument mcmc_type = "approx", and then perform the post-correction step using the results from the approximate MCMC. In doing so, we can also use the function suggest_N to find a suitable number of particles N for ψ-APF in the spirit of Doucet et al. (2015): out_approx <-run_mcmc(mng_model, mcmc_type = "approx", iter = 50000) est_N <-suggest_N(mng_model, out_approx) out_exact <-post_correct(mng_model, out_approx, particles = est_N$N) The function suggest_N computes the standard deviation of the logarithm of the post-correction weights (i.e. the random part of log-likelihood of ψ-APF) at the approximate MAP estimator of θ using a range of N and returns a list with component N which is the smallest number of particles where the standard deviation was less than one. For small and moderate problems typically 10-20 particles is enough.

Filtering and smoothing
The bssm also offers separate methods for performing (approximate) state filtering and smoothing which may be useful in some custom settings.
For LGSSM, methods kfilter and smoother perform Kalman filtering and smoothing. For non-Gaussian models with linear-Gaussian dynamics, approximate filtering and smoothing estimates can be obtained by calls to kfilter and smoother, in which case these functions first construct an approximating Gaussian model for which the Kalman filter/smoother is then applied. For non-linear models defined by nlg_ssm we can run approximate filtering using the extended Kalman filter with the function ekf, the unscented Kalman filter with the function ukf, or the iterated EKF (IEKF) by changing the argument iekf_iter of the ekf function. Function ekf_smoother can be used for smoothing based on EKF/IEKF.
For particle filtering the bssm package supports general bootstrap particle filter for all model classes of the bssm (function bootstrap_filter). For nlg_ssm, extended Kalman particle filtering (Van Der Merwe et al., 2001) is also supported (function ekpf_filter).
For particle smoothing, function particle_smoother with the smoothing based on BSF is available for all models. In addition, ψ-APF (using argument method = "psi") is available for all models except of ssm_sde class. Currently, only filter-smoother approach (Kitagawa, 1996) for particle smoothing is supported.  compared the computational efficiency of delayed acceptance MCMC and importance sampling type MCMC approaches in various settings. Here we make a small experiment comparing the generic Hamiltonian Monte Carlo using the NUTS sampler (Hoffman and Gelman, 2014) with rstan, and IS-MCMC with bssm. Given that the bssm is specialized for state space models whereas Stan is a general purpose tool suitable for wider range of problems, it is to be expected that bssm performs better in terms of computational efficiency. The purpose of this experiment is to illustrate this fact, i.e., that there is still demand for specialized algorithms for various types of statistical models. For complete code of the experiment, see Appendix: Code for section Comparison of IS-MCMC and HMC.

Comparison of IS-MCMC and HMC
We consider the case of a random walk with drift model with negative binomial observations and some known covariate x t , defined as with zero-mean Gaussian noise term η t with unknown standard deviation σ µ . Based on this we simulate one realization of y and x with n = 200, φ = 5, β = −0.9, ν = 0.01, σ µ = 0.1.
For the IS approach we use ng_bsm function for model building, with prior variances 100 and 0.01 for the initial states µ 1 and ν 1 . For hyperparameters, we used fairly uninformative half-Normal distribution with standard deviation 0.5 for σ µ and 0.1 for σ ν . We then ran the IS-MCMC algorithm with run_mcmc using a burn-in phase of length 10,000 and run 50,000 iterations after the burn-in, with 10 particles per SMC.
Using the same set up, we ran the MCMC with rstan using 15,000 iterations (with first 5000 used for warm-up). Note that in order to avoid sampling problems, it was necessary to tweak the default control parameters of the sampler (see Appendix). Table 1 shows the results. We see both methods produce identical results (within the Monte Carlo error), but while rstan produces similar Monte Carlo standard errors with smaller amount of total iterations than bssm, the total computation time of rstan is almost 80 times higher than with bssm (58 minutes versus 45 seconds), which suggests that for these type of problems it is highly beneficial to take advantage of the known model structure and available approximations versus general Bayesian software such as Stan which makes no distinction between latent states α and hyperparameters θ.

Conclusions
State space models are a flexible tool for analysing a variety of time series data. Here we introduced the R package bssm for fully Bayesian state space modelling for a large class of models with several alternative MCMC sampling strategies. All computationally intensive parts of the package are
Compared to early versions of the bssm package, the option to define R functions for model updating and prior evaluation have lowered the bar for analysing custom models. The package is also written in a way that it is relatively easy to extend to new model types similar to current bsm_lg in future. The bssm package could be expanded to allow other proposal adaptation schemes such as adaptive Metropolis algorithm by Haario et al. (2001), as well as support for multivariate SDE models and automatic differentiation for EKF-type algorithms.