## Introduction

Dynamic regression models are extension to basic linear regression models where instead of constant but unknown regression coefficients, the underlying coefficients are assumed to vary over “time” according to random walk. In their simplest form these models can be used to model regression models with additional time series component, and they also allow robust modelling of phenomenas where the effect size of the predictor variables can vary during the period of the study. The R (R Core Team 2017) package walker provides an efficient method for fully Bayesian inference of such models, where the main computations are performed using the Markov chain Monte Carlo (MCMC) algorithms provided by Stan (Stan Development Team 2016b, Stan Development Team (2016a)). This also allows the use of many general diagnostic and graphical tools provided by several Stan related R packages such as ShinyStan (Stan Development Team 2017).

More specifically, the dynamic linear regression model is defined as \begin{aligned} y_t &= x'_t \beta_t + \epsilon_t, \quad t = 1,\ldots, n\\ \beta_{t+1} &= \beta_t + \eta_t, \end{aligned} where $$y_t$$ is the observation at time $$t$$, $$x_t$$ contains the corresponding predictor variables, $$\beta_t$$ is a $$k$$ dimensional vector of regression coefficients at time $$t$$, $$\epsilon_t \sim N(0, \sigma^2_{\epsilon})$$, and $$\eta_t \sim N(0, D)$$, with $$D$$ being $$k \times k$$ diagonal matrix with diagonal elements $$\sigma^2_{i,\eta}$$, $$i=1,\ldots,k$$. Denote the unknown parameters of the model by $$\beta = (\beta'_1, \ldots, \beta'_n)'$$ and $$\sigma = (\sigma_{\epsilon}, \sigma_{1, \eta}, \ldots, \sigma_{k, \eta})$$. We define priors for first $$\beta_1$$ as $$N(\mu_{\beta_1}, \sigma_{\beta_1})$$, and for $$\sigma_i \sim N(\mu_{\sigma_i}, \sigma_{\sigma_i})$$, $$i=1,\ldots,k+1$$, truncated to positive values, with slighly awful notation.

Although in principle writing dynamic linear regression model as above in Stan or BUGS (Lunn et al. 2000) language is straighforward, these standard implementations are computationally inefficient and prone to severe problems related to convergence of the underlying MCMC algorithm. For example in (block) Gibbs sampling approach we target the joint posterior $$p(\beta, \sigma | y)$$ by sampling from $$p(\beta | \sigma, y)$$ and $$p(\sigma | \beta, y)$$. But because of strong autocorrelations between the coeffients $$\beta$$ at different time points, as well as with the associated standard error parameters, this MCMC scheme can be very inefficient. Also, the total number of parameters to sample increase with number of data points $$n$$. Although the Hamiltonian Monte Carlo algoritms offered by Stan are typically more efficient exploring the posterior than Gibbs-type algorithms, we still encounter similar problems (see the illustration in next Section).

An alternative solution used by walker is based on the marginalization of the regression coefficients $$\beta$$ during the MCMC sampling by using the Kalman filter. This provides fast and accurate inference of marginal posterior $$p(\sigma | y)$$. Then, the corresponding joint posterior $$p(\sigma, \beta | y) = p(\beta | \sigma, y)p(\sigma | y)$$ can be obtained by simulating the regression coefficients given marginal posterior of standard deviations. This sampling can be performed for example by simulation smoothing algorithm by Durbin and Koopman (2002). Note that we have opted to sample the $$\beta$$ parameters given $$\sigma$$’s, but it is also possible to obtain somewhat more accurate summary statistics such as mean and variance of these parameters by using the standard Kalman smoother for compution of $$\textrm{E}(\beta| \sigma, y)$$ and $$\textrm{Var}(\beta| \sigma, y)$$, and using the law of total expectation.

In this vignette we first introduce the basic use using simple linear regression model with first order random walk, and then discuss extensions to second order random walk, as well as how to deal with non-Gaussian data.

## Illustration

Let us consider a observations $$y$$ of length $$n=250$$, generated by random walk (i.e. time varying intercept) and two predictors. This is rather small problem, but it was chosen in order to make possible comparisons with the “naive” Stan implementation. For larger problems (in terms of number of observations and especially number of predictors) it is very difficult to get naive implementation to work properly, as even after tweaking several parameters of the underlying MCMC sampler, one typically ends up with divergent transitions or low BMFI index, meaning that the results are not to be trusted.

First we simulate the coefficients and the predictors:

set.seed(1)
n <- 100
beta1 <- cumsum(c(0.5, rnorm(n - 1, 0, sd = 0.05)))
beta2 <- cumsum(c(-1, rnorm(n - 1, 0, sd = 0.15)))
x1 <- rnorm(n, mean = 2)
x2 <- cos(1:n)
rw <- cumsum(rnorm(n, 0, 0.5))
ts.plot(cbind(rw, beta1 * x1, beta2 * x2), col = 1:3)

signal <- rw + beta1 * x1 + beta2 * x2
y <- rnorm(n, signal, 0.5)
ts.plot(cbind(signal, y), col = 1:2)