Dynamic linear model tutorial

Statistical analysis of time series data is usually faced with the fact that we have only one realization of a process whose properties we might not fully understand. In analysis of correlation structures, for example, we need to assume that some distributional properties of the process generating the variability stay unchanged in time. In linear trend analysis, we assume that there is an underlying change in the background that stays approximately constant over time. Dynamic regression with state space approach tries to avoid some of the problems. By explicitly allowing for variability in the regression coefficients we let the system properties change in time. Furthermore, the use of unobservable state variables allows direct modelling of the processes that are driving the observed variability, such as seasonality or external forcing, and we can explicitly allow for some modelling error.

Dynamic regression can, in very general terms, be formulated using state space representation of the of the observations and the state of the system. With a sequential definition of the processes, having conditional dependence only on the previous time step, the classical Kalman filter formulas can be used to estimate the states given the observations. When the operators involved in the definition of the system are linear we have so called dynamic linear model, DLM.

A basic model for many climatic time series consists of four elements: slowly varying background level, seasonal component, external forcing of known processes modelled by proxy variables, and stochastic noise. The noise components might contain an autoregressive correlation structure to account for correlated model residuals. However, autocorrelated errors can be problematic as they might be caused by a long memory of the process as well as by some systematic features not included in the model.

This is a short tutorial on using dynamic linear models in time series analysis. It is based on my attempts to analyse some climatic time series. For that purposes I wrote the Matlab code described here. The examples deal with univariate time series, but the code can handle multivariate data, to some extent. Here we are mostly interested in extracting the components, and using these to infer about some interesting features of the system, but not to produce predictions about the behaviour of the system in the future, although understanding the system is a first step to be able to make predictions.

The use of DLMs in time series analysis is well documented in statistical literature, see the references at the end.

General dynamic linear model can be written with a help of observation equation and model equation as

Above \(y_{t}\) are the \(p\) observations at time \(t\), with \(t=1,\dots,n\). Vector \(x_{t}\) of length \(m\) contains the unobserved states of the system that are assumed to evolve in time according to a linear system operator \(G_{t}\) (a \(m\times m\) matrix). In time series settings \(x_t\) will have elements corresponding to various components of the time series process, like trend, seasonality, etc. We observe a linear combination of the states with noise and matrix \(F_{t}\) (\(p\times m\)) is the observation operator that transforms the model states into observations. Both observation end system equations can have additive Gaussian errors with covariance matrices \(V_{t}\) and \(W_{t}\).

This formulation is quite general and flexible and it allows handling of many time series analysis problems in a single framework. Moreover, an unified computational tool can be used, i.e. a single DLM computer code can be used for various purposes. Below we give examples of different analyses. As we are dealing with linear models, we assume that the operators \(G_{t}\) and \(F_{t}\) are linear, i.e. matrices. However, they can change with the time index \(t\). The state space framework can be extended to non linear model and non Gaussian errors, and to spatial-temporal analyses as well, see, e.g., Cressie.

The DLM formulation can be seen as a special case of a general hierarchical statistical model with three levels: data, process and parameters (see e.g. Cressie). In terms of statistical distributions we have, first, the observations uncertainty \(p(y_{t}|x_{t},\theta)\) described by the observation equation. Second, the process uncertainty of the unknown states \(x_{t}\) and their evolution given by the process equations as \(p(x_{t}|\theta)\), and third, the uncertainty related to model parameters \(p(\theta)\). These conditional formulations allows both efficient description of the system and computational tools to estimate its components. It also combines different statistical approaches, as we can have full prior probabilities for the unknowns (Bayesian approach), estimate them by maximum likelihood and plug them back (frequentistic approach), or even fix the model parameters by expert knowledge (a trivial non statistical approach).

Recall that \(y_{t}\) are the observations and \(x_{t}\) the hidden system states for \(t=1,\dots,n\). In addition, we have (static) model parameter \(\theta\) that contains auxiliary parameters needed in defining the model and observation errors \(W_{t}\) and \(V_{t}\) and the system matrices \(G_{t}\) and \(F_{t}\). For dynamic linear models we have efficient and well founded computational tools for all the relevant statistical distributions.

• \(p(x_{t+1}|x_{t},y_{1:t},\theta)\) by Kalman filter
• \(p(x_{t}|y_{1:t},\theta)\) by Kalman filter
• \(p(x_t|y_{1:n},\theta)\) by Kalman smoother
• \(p(x_{1:n}|y_{1:n},\theta)\) by simulation smoother
• \(p(y_{1:t}|\theta)\) by Kalman filter likelihood
• \(p(x_{1:n},\theta|y_{1:n})\) by MCMC
• \(p(x_{1:n}|y_{1:n})\) by MCMC

Here are the relevant parts of the recursive formulas for Kalman filter and smoother to estimate the marginal distributions of DLM states given the observations. We assume that the initial distributions at \(t=1\) are available. First, we perform Kalman filter forward recursion for the predicted states \(p(x_{t+1}|x_{t},y_{1:t},\theta) = N(\hat{x}_{t+1},\hat C_{t+1})\), \(t=1,2,\dots,n-1\)

Then, we apply Kalman smoother backward recursion to obtain the smoothed states \(p(x_t|y_{1:n},\theta) = N(\tilde{x}_t,\tilde C_t)\), for \(t=n,n-1,\dots,2,1\)

Both these recursion are implemented in dlmsmo.m Matlab code.

Kalman smoother algorithm provides the marginal distributions \(p(x_{t}|y_{1:n},\theta)\) for each \(t\). These are all Gaussian. However, for studying trends and other dynamic features in the system, we are interested in the joint distribution spanning the whole time range \(p(x_{1:n}|y_{1:n},\theta)\). Note that we are still conditioning on the unknown parameter vector \(\theta\) and will account for it later. This very high dimensional joint distribution is not easily accessible directly. As in many cases, instead of analytic expressions, it is more important to be able to draw realizations from the distribution.

The system equations provide a direct way to recursively produce realizations of both the states \(x_{1:n}\) and the observations \(y_{1:n}\). However, the generated states will be independent of the original observations. It can be shown (Durbin & Koopman, Section 4.9) that the residual process of generated vs. smoothed state is independent of the \(\tilde{x}_{1:n}\) and \(\tilde{y}_{1:n}\). This means that if we add these residuals over the smoothed state \(x_{1:n}\), we get a new realization that is conditional on \(y_{1:n}\), the original observations. So to produce \({x}^*_{1:n} \sim p(x_{1:n}|y_{1:n},\theta)\), we:

1. Sample from the system equations to get \(\tilde{x}_{1:n}\) and \(\tilde{y}_{1:n}\).
2. Smooth \(\tilde{y}_{1:n}\) to get \(\tilde{\tilde{x}}_{1:n}\).
3. Add the residuals to the original smoothed state, \({x}^*_{1:n} = \tilde{x}_{1:n} - \tilde{\tilde{x}}_{1:n} + {x}_{1:n}\).

This simulation smoother will be used in trend studies and as a part of more general simulation algorithm that will sample from the joint posterior distribution \(p(x_{1:n},\theta|y_{1:n})\), and by marginalization argument also from \(p(x_{1:n}|y_{1:n})\) where the uncertainty in \(\theta\) has been integrated out. Simulation smoother is implemented in dlmsmo and dlmsmosam functions. As a note, the simulation smoother can be seen as a special case of the randomize-then-optimize (RTO) method in Bardsley et al. (2104).

In general terms, trend is some change in the distributional properties, such as in the mean, of the process generating the observations. We are typically interested in slowly varying changes in the background level, i.e. in the mean process after known sources of variability, such as seasonality, have been accounted for. A common way to explore trends is to fit some kind of a smoother, such as a moving average, over the time series. Many smoothing methods do not provide statistical ways to estimate the smoothness parameters or asses the uncertainty related to the level of smoothing.

In DLM trend analysis, a slowly varying (relative to the time scale we are interested in) background level of the system is modelled as a random walk process with variance parameters that determine the time wise smoothness of the level. These variance parameters must be estimated. In an optimal case, the data provides information also on the smoothness of the trend component, but sometimes we need to use subject level prior information to decide the time scale of the changes we want to extract. DLM models offer means of providing qualitative prior information in the form of the model equations and quantitative information by prior distributions on variance parameters, see (Gamermann).

Let \(x_{\mathrm{level},t}\) be the model state that defines the background level of the process. For statistical analysis, we need to estimate the whole state, as either \(p(x_{1:n}|y_{1:n},\hat{\theta})\), where we plug in some estimates of the auxiliary parameters \(\theta\), (maximum likelihood approach) or by \(p(x_{1:n}|y_{1:n}) = \int p(x_{1:n},\theta|y_{1:n})\,d\theta\) where the uncertainty of auxiliary parameters \(\theta\) are integrated out (Bayesian approach). The latter is typically done by Markov chain Monte Carlo (MCMC) simulation.

In examples above, the variance parameters defining the model error covariance matrix \(W_t\) were assumed to be known. In practice we need some estimation methodology. Basically there are three alternatives. One. Use subject level knowledge with trial and error to fix the parameters without any algorithmic tuning. Two. Use the likelihood function with a numerical optimization routine to find maximum likelihood estimates of the parameters and plug the estimates back to the equations and re-fit the DLM model. Three. Use MCMC to sample from the posterior distribution of the parameters to estimate the parameters or to integrate out their uncertainty.

To estimate the free parameters \(\theta\) in the model formulation we need the marginal likelihood function \(p(y_{1:n}|\theta)\). By the assumed Markov properties of the system, this is obtained sequentially as a by product of the Kalman filter recursion,

On the right hand side, the parameter \(\theta\) will appear in the model predictions \(\hat{x}_{t}\) as they depend on the model formulation, on the matrix \(G_t\), as well as on the model error \(W_t\). For the same reason we need the determinant of the model prediction covariance matrix \(C_{t}^y\). A good thing is that the above likelihood can be calculated along the DLM filter and smoother recursions without much extra effort.

The dlmfit Matlab function has code for both maximum likelihood and MCMC calculations of parameters in the diagonal of the model error matrix \(W\), autoregressive coefficients in \(G\), as well as for estimating one extra multiplicative factor for the observation error matrix \(V\). Other parameter estimation tasks can be programmed, with the help of the likelihood returned by the dlmsmo function, also. See the AR example below.

out = dlmsmo(y,F,V,x0,G,W,C0,X, sample);

The main workhorse for the DLM calculations. Given the definition of DLM model as input, return Kalman filter and Kalman smoother estimates of the states, an optional sample from the state and observations, and statistics related to the fit.

Input:

Above \(m\) is the dimension of the state vector, \(n\) is the number of observation times, \(p\) is the dimension of the observation at each time and \(q\) is the number of regression covariates in the model. If covariates are given in matrix X, the observation operator and the model operator are appended with suitable entries for X, hence their input dimensions are of size \(m-q\). However, the initial values x0 and C0 and model error covariance matrix W should already contain the elements corresponding to columns of X. An empty X as input means no covariates.

Output is a structure that contains the element below, and many more (to be added).

The function dlmsmo is called by dlmfit that contains code for constructing the system matrices and estimating for some of the model parameters. The ploting functions described below can be used to plot the model output.

out = dlmfit(y,s,w,x0,C0,X,options);

The function dlmfit uses dlmsmo to fit univariate DLM time series model to observations. It generated the suitable model structure according to options given. Optionally it performs optimization or MCMC sampling analysis for parameters in the diagonal of model error matrix \(W\).

The inputs are as the following.

y

The observations to be fitted.

s

The observation uncertainties as standard deviations, scalar or matrix of same size as y.

w

Square roots in the diagonal of the model error matrix \(W\), up to the last non zero element.

x0

The initial values for the model states, if not given, the code tries to use some reasonable first guess (typically 0).

C0

Covariance matrix for the initial model state uncertainty. Defaults to some large values in the diagonal.

X

Covariate proxy time series. Defaults to empty matrix [].

options

System matrix generation and optimization options as a structure.

This function uses dlmgensys to generate the system matrices. So, all the options defined for it can be used. In addition, the options structure can have the following options related to the estimation of variance parameters

As an output the function returns the result of dlmsmo after optional estimation of model parameters. Even if no iterative optimization is performed, the function dlmsmo is run at least twice to iterate the initial values x0 and C0. Also, some additional elements are added to the output structure. (These are to be documented later, see the example and source.)

Currently, the optimization and MCMC simulation can both have target from the diagonal of the model error matrix \(W\) and one scale parameter for the observation error \(V\). For some models, there are predefined optimization targets. In general you should use winds to define which indexes in diag(\(W\)) are to be estimated by maximum likelihood or sampled by MCMC.

To use winds, you need to know the order of the model components in \(W\): trend, seasonal, AR parameters, covariates. For example if the model has trend of order 1, and two harmonic seasonal components (trig=2) then having one parameter for the local trend variability and one common parameter for the seasonal variability would be achieved by setting options.winds = [0 1 2 2 2 2].

If you need to estimate something else than the diagonal of \(W\), then you can build the target function by using the lik output from dlmsmo. See the AR example.

[G,F] = dlmgensys(options);

This function generates system evolution matrix and corresponding observation operator matrix for several common dynamic linear models used in time series analyses. The input variable options is a structure having elements that define the components. The following elements are possible at the moment:

Example

[G,F] = dlmgensys(struct('order',1,'fullseas',1,'ns',4))

G =

     1     1     0     0     0
     0     1     0     0     0
     0     0    -1    -1    -1
     0     0     1     0     0
     0     0     0     1     0


F =

     1     0     1     0     0

sample = dlmsmosam(dlm,nsam);

Sample of states from DLM fit. Input: dlm fit object, nsam number of samples to draw. Output: matrix of size nstate*ntime*nsam. Uses the simulation smoother, i.e. samples from joint distribution of states at all times given the full set of observations.

dlmplotfit(dlm, time, yscale, yind);

Plots the observations with the estimated level component of the model. Input: dlm is output from dlmfit, time is optional x axis for the plots, yscale optional scaling factor for y, and yind the column of y to be plotted.

dlmplotdiag(dlm);

dlmplotdiag(dlm);

Plot residual diagnostics of DLM fit performed by dlmsmo or dlmfit.

dlmplotmcmc(dlm);

Plot the MCMC chain used in estimation of model parameter in function dlmfit.