Introduction to Dynamic Linear Models for Time Series Analysis
Marko Laine


Marko Laine Finnish Meteorological Institute, Helsinki, Finland, marko.laine@fmi.fi

[
    [
    December 30, 2023
=====================




Dynamic linear models (DLM) offer a very generic framework to analyse time series data. Many classical time series models can be formulated as DLMs, including ARMA models and standard multiple linear regression models. The models can be seen as general regression models where the coefficients can vary in time. In addition, they allow for a state space representation and a formulation as hierarchical statistical models, which in turn is the key for efficient estimation by Kalman formulas and by Markov chain Monte Carlo (MCMC) methods. A dynamic linear model can handle non-stationary processes, missing values and non-uniform sampling as well as observations with varying accuracies. This chapter gives an introduction to DLM and shows how to build various useful models for analysing trends and other sources of variability in geodetic time series.








*[1]

*[1]




§ INTRODUCTION TO DYNAMIC LINEAR MODELS


Statistical analysis of time series data is usually faced with the fact that we have only one realization of a process whose properties might not be fully understood. We need to assume that some distributional properties of the process that generate the observations do not change with time. In linear trend analysis, for example, we assume that there is an underlying change in the background mean that stays approximately constant over time. Dynamic regression avoids this by explicitly allowing temporal variability in the regression coefficients and by letting some of the system properties to change in time. Furthermore, the use of unobservable state variables allows direct modelling of the processes that are driving the observed variability, such as seasonal variation or external forcing, and we can explicitly allow some modelling error.

Dynamic regression can  be formulated in very general terms by using a state space representation of the observations and the hidden state of the system. With sequential definition of the processes, having conditional dependence only on the previous time step, the classical recursive Kalman filter algorithms can be used to estimate the model 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 time series in geodetic or more general environmental applications consists of four elements: a slowly varying background level, a seasonal component, external forcing from known processes modelled by proxy variables, and stochastic noise. The noise component might contain an autoregressive structure to account for temporally correlated model residuals. As we see, the basic components have some physical justification and we might be interested in their contribution to the overall variability and their temporal changes. These components are hidden in the sense that we do not observe them directly and each individual component is masked by various other sources of variability in the observations.

Below, we briefly describe the use of dynamic linear models in time series analysis. The examples deal with univariate time series, i.e.  the observation at a singe time instance is a scalar, but the framework and the computer code can handle multivariate data, too. All the model equations are written in way that support multivariate observations. In the presented applications we are mostly interested in extracting the components related to the trends and using these to infer about their magnitude and the uncertainties involved. However, these models might not be so good for 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, but they might go by different terminology and notation. In <cit.> they are called structural time series, <cit.> uses the state space approach, and the acronym DLM is used in <cit.>.



§ STATE SPACE DESCRIPTION


The state space description offers a unified formulation for the analysis of dynamic regression models. The same formulation is used extensively in signal processing and geophysical data assimilation studies, for example. A general dynamic linear model with an observation equation and a model equation is

    y_t   =   H_tx_t+ϵ_t,   ϵ_t∼
      N(0,R_t), 
    
      x_t   =   M_tx_t-1+E_t,    E_t∼ N(0,Q_t).


Above y_t is a vector of length k of observations at time t, with t=1,…,n. Vector x_t of length m contains the unobserved states of the system that evolve in time according to a linear system operator M_t (a m× 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 ϵ_t, and matrix H_t (m× k) is the observation operator that transforms the model states into observations. Both observations and the system states can have additive Gaussian errors with covariance matrices R_t (k× k) and Q_t (m× m), respectively. In univariate time series analysis we will have k=1. With multivariate data, the system matrices M_t, H_t, R_t and Q_t can be used to define correlations between the observed components.

This formulation is quite general and flexible as it allows handling of many time series analysis problems in a single framework. Moreover, a 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 M_t and H_t are linear. However, they can change with the time index t and we will drop the time index in the cases where the matrices are assumed static in time. 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., <cit.>. However, as can be seen in the following example, already the dynamic linear Gaussian formulation provides a large class of models for time series trend analyses.



 §.§ Example: spline smoothing


A simple local level and local trend model can be used as a basis for many trend related studies. Consider a mean level process μ_t which is changing smoothly in time and which we observe with additive Gaussian noise. We assume that the change in the mean, μ_t+1-μ_t, is controlled by a trend process α_t and the temporal change in these processes is assumed to be Gaussian with given variances σ^2_level and σ^2_trend. This can be written as

    y_t    = μ_t + ϵ_obs,    ϵ_obs∼ N(0,σ^2_obs),      observations,
    μ_t   = μ_t-1 + α_t-1 + ϵ_level,
         ϵ_level∼ N(0,σ^2_level),     local level,
    α_t   =α_t-1 +ϵ_trend,
        ϵ_trend∼ N(0,σ^2_trend),
         local trend,

which in state space representation transfers simply into

    x_t = [ μ_t α_t ]^T,  H= [ 1 0;  ],  M=[ 1 1; 0 1;   ],


    Q=[ σ^2_level         0;         0 σ^2_trend;           ],  and  R= diag([ σ^2_obs       … σ^2_obs ]),

with three parameters for the error variances

    θ =
      [   σ^2_obs σ^2_level σ^2_trend ]^T.

We have dropped the time index t from those elements that do not depend on time.

It is interesting to note, that if we set σ^2_level = 0, we have a second difference process for μ_t as

    Δ^2μ_t = μ_t-1 - 2μ_t + μ_t+1 = Δα_t = ϵ_trend,

and it can be shown <cit.> that this is equivalent to cubic spline smoothing  with smoothing parameter λ = σ^2_trend/σ^2_obs > 0.

Figure <ref> shows simulated observations with a true piecewise trend and the fitted mean process μ_t, t=1,…,n together with its 95% uncertainty limits. In this example, the observation uncertainty standard deviation (σ_obs=0.3) as well as the level and trend variability standard deviations (σ_level=0.0, σ_trend=0.01) are assumed to be known. In the later examples these values are estimated from the data. 





§ DLM AS HIERARCHICAL STATISTICAL MODEL


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. <cit.>), with corresponding conditional statistical distributions. First, the observation uncertainty p(y_t| x_t,θ) described by the observation equation and forming the statistical likelihood function, second, the process uncertainty of the unknown states x_t and their evolution given by the process equations as p(x_t|θ) or p(x_t| x_t-1,θ), and third, the unconditional prior uncertainty for the model parameters p(θ). This formulation allows both an efficient description of the system and computational tools to estimate the components. It also combines different statistical approaches, as we can have full prior probabilities for the unknowns (the Bayesian approach), estimate them by maximum likelihood and plug them back (frequentistic approach), or even fix the model parameters by expert knowledge (a non-statistical approach). By the Bayes formula, we can write the state and parameter posterior distributions as a product of the conditional distributions

    p(x_t,θ| y_t) ∝ p(y_t| x_t,θ) p(x_t|θ) p(θ),

which is the basis for full Bayesian estimation procedures. Next we will describe the steps needed for Bayesian DLM estimation of model states, parameters and their uncertainties.



§ STATE AND PARAMETER ESTIMATION


To recall the notation, y_t are the observations and x_t are the hidden system states for time indexes t=1,…,n. In addition, we have a static vector θ that contains auxiliary parameters needed in defining the system matrices M_t and H_t and the model and observation error covariance matrices Q_t and R_t. For dynamic linear models we have efficient and well founded computational tools for all relevant statistical distributions of interest. For the state estimation assuming a known parameter vector θ the assumptions on linearity and Gaussian errors allows us to estimate the model states by classical recursive Kalman formulas. The variance and other structural parameters appear in non-linear way and their estimation can be done either by numerical optimization or by Markov chain Monte Carlo (MCMC) methods. MCMC allows for a full Bayesian statistical analysis for the joint uncertainty in the dynamic model states and the static structural parameters <cit.>. Table <ref> relates the different statistical distributions to the algorithms, which are outlined later. The notation y_1:t, x_1:t, etc.  means the collection of observations or states from time 1 to time t.





§ RECURSIVE KALMAN FORMULAS


Below we give the relevant parts of the recursive formulas for Kalman filter and smoother to estimate the conditional distributions of DLM states given the observations and static parameters. For more details, see <cit.>. A notable feature of the linear Gaussian case is that the formulas below are exact and easily implemented in computer as long as the model state dimension or the number of observations at one time is not too large. 

To start the calculations, we assume that the initial distribution of x_0 at t=0 is available. The first step in estimating the states is to use Kalman filter forward recursion to calculate the distribution of the state vector _t  given the observations up to time t,
p(_t|y⃗_1:t,θ⃗) = N(_t,C_t), which is Gaussian by the linearity assumptions.
At each time t this step consists of first calculating, as prior, the mean and covariance matrix of one-step-ahead predicted states p(_t|_t-1,y⃗_1:t-1,θ⃗) = N(_t,C_t) and the covariance matrix of the predicted observations C_y,t as 

    _t   = M_t_t-1   prior mean for x_t,
    C_t   = M_t C_t-1M_t^T + Q_t    prior covariance for x_t,
    C_y,t   = H_t C_tH_t^T + R_t   covariance for predicting y_t.

Then the posterior state mean and its covariance are calculated using the Kalman gain matrix G_t as

    G_t    =  C_tH_t^TC_y,t^-1   Kalman gain,
    r⃗_t    = y⃗_t - H_t_t   prediction residual,
    _t   = _t + G_t r⃗_t   posterior mean for x_t,
    C_t   = C_t-G_tH_tC_t   posterior covariance for x_t.

These equations are iterated for t=1,…, and the values of x_t and C_t are stored for further calculations. As initial values, we can use _0=0⃗ and C_0=κI, i.e.  a vector of zeros and a diagonal matrix with some large value κ in the diagonal. Note that the only matrix inversion required in the above formulas is the one related to the observation prediction covariance matrix C_y,t, which is of size 1×1 when we analyse univariate time series.

The Kalman filter provides distributions of the states at each time t given the observations up to the current time. As we want to do retrospective time series analysis that accounts for all of the observations, we need to have the distributions of the states for each time, given all the observations  y⃗_1:. By the linearity of the model, these distributions are again Gaussian,
p(_t|y⃗_1:,θ⃗) = N(_t,C_t).
Using the matrices generated by the Kalman forward recursion,
the Kalman smoother backward recursion gives us the smoothed states for t=,-1,…,1. There are several equivalent versions of the backward recursion algorithm. Below we show the Rauch-Tung-Striebel recursion <cit.> for illustration. For alternatives, see <cit.>:


    C_t^+   = M_t C_tM_t^T + Q_t   propagated covariance,
    G_t   = C_tM_t^T(C_t^+)^-1   smoother gain,
    x_t-1   = x_t + G_t(x_t - x_t)    smoothed state mean,
    C_t-1   = C_t - G_t( C_t - C_t^+)G_t^T   smoothed state covariance.


In smoother recursion we are dealing with several matrix vector computations and one matrix inversion of size m× m and these formulas can be implemented quite efficiently in any general numerical analysis software. As a note, we see that the algorithms work with missing observations, too. If some observations at a time t are missing, the corresponding columns of the gain matrix Eq. (<ref>) will be zero. If all are missing, the filter posterior will be equal to the prior. Note that the above smoother recursion does not refer to the observations. All the Kalman formulas given above are for observations with uniform sampling in time, for non-uniform temporal sampling, the propagation of uncertainty to the next observation time has to be handled differently, see <cit.>.


§ SIMULATION SMOOTHER


The Kalman smoother algorithm provides the distributions p(x_t|y_1:n,θ) for each t, which 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,θ). Note that we are still conditioning on the unknown parameter vector θ and will account for it later. This 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 and use the sampling distribution for statistical analysis. This has several benefits. One important is that by comparing simulated realizations to the observations, we see how realistic the model predictions are, which can reveal if the modelling assumptions are not valid. Also, we can study the distributions of model outputs directly from the samples and do not need to resolve to approximate statistics.

A simple simulation algorithm by <cit.> is the following. The state space system equations provide a direct way to recursively sample realizations of both the states x_1:n and the observations y_1:n, but the generated states will be independent of the original observations. However, it can be shown <cit.>  that the distribution of the residual process of generated against smoothed state does not depend on y_1:n. This means that if we add simulated residuals over the original smoothed state x_1:n, we get a new realization that is conditional on the original observations y_1:n. A procedure to sample a realization x^*_1:n∼ p(x_1:n|y_1:n,θ) is thus:



  * Generate a sample using the state space system equations, Eqs. (<ref>) and (<ref>) to get x̌_1:n and y̌_1:n.

  * Smooth y̌_1:n to get x̆_1:n according to formulas in Section <ref>.

  * Add the residuals from step 2 to the original smoothed state:


    x^*_1:n = x̌_1:n - x̆_1:n + x_1:n.


This simulation smoother can be used in trend studies and as a part of more general MCMC simulation algorithm that will sample from the joint posterior distribution p(x_1:n,θ|y_1:n) and by marginalization argument also from p(x_1:n|y_1:n) where the uncertainty in θ has been integrated out <cit.>.



§ ESTIMATING THE STATIC STRUCTURAL PARAMETERS


In the first examples, the variance parameters defining the model error covariance matrix Q_t were assumed to be known. In practice we need some estimation methodology for them. Basically there are three alternatives. The first one uses subject level knowledge with trial and error to fix the parameters without any algorithmic tuning. The second one use the marginal likelihood function with a numerical optimization routine to find the maximum likelihood estimate of the parameter θ and plug the estimate back to the equations and re-fit the DLM model. The third one use MCMC algorithm to sample from the posterior distribution of the parameters to estimate the parameters and to integrate out their uncertainty.

To estimate the free parameters θ in the model formulation by optimization or by MCMC we need the marginal likelihood function p(y_1:n|θ). By the assumed Markov properties of the system, this can be obtained sequentially as a byproduct of the Kalman filter recursion <cit.>,

    -2log(p(y_1:n|θ)) = constant+ ∑_t=1^n[(y_t-H_tx_t)^TC_y,t^-1(y_t-H_tx_t)
            + log(|C_y,t|)].

On the right hand side, the parameter θ will appear in the model predictions x_t as they depend on the matrix M_t as well as on the model error Q_t. For the same reason we need the determinant of the model prediction covariance matrix |C_y,t|. A fortunate property is that this likelihood can be calculated along the DLM filter recursion without much extra effort. 

The scaled one-step prediction residuals 

    r^*_t = C_y,t^-1/2(y_t-H_tx_t)

can be used to check the goodness of fit of the model. In order of the DLM model to be consistent with the observations these residuals should be approximately independent, N(0,I) Gaussian and without serial autocorrelation. Later in the GNSS time series example, we will do model diagnostics by residual quantile-quantile and autocorrelation function plots.



§ ANALYSING TRENDS


In general terms, trend is a 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 the known sources of variability, such as seasonality, has 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. However, many standard smoothing methods do not provide statistical estimates of the smoothness parameters or asses the uncertainty related to the level of smoothing.

In typical DLM trend analyses, a slowly varying (relative to the time scale we are interested in) background level of the system is modelled as a first or higher order random walk process with variance parameters that determine the time wise smoothness of the level. These variance parameters must be estimated and their uncertainty accounted for proper uncertainty quantification. In an optimal case, the data provides information on the smoothness of the trend component, but typically we need to use subject level prior information to decide the time scale of the changes we want to extract. 
In the GNSS application example in Section <ref> we assume a global linear trend and the local non-stationary fluctuations are modelled using a local random walk model with autocorrelated residuals.
A Bayesian DLM model offers means to provide qualitative prior information in the form of the model equations and quantitative information by prior distributions on the variance parameters, see e.g.  <cit.>.

For statistical analysis we need to estimate the full state as either p(x_1:n|y_1:n,θ̂), where we plug in some estimates of the auxiliary parameters θ̂, (the maximum likelihood approach) or by p(x_1:n|y_1:n) = ∫ p(x_1:n,θ|y_1:n) dθ where the uncertainty of auxiliary parameter θ is integrated out. The latter is  the Bayesian approach and calculations can be done, e.g., by Markov chain Monte Carlo (MCMC) simulation <cit.>. 
A procedure to account the uncertainty in a DLM model and its structural parameters and to study DLM output will contain the following steps:



  * Formulate the DLM model and its marginal likelihood p(y_1:n|θ) by Kalman filter.

  * Use MCMC to sample from the posterior distribution p(θ| y_1:n) with a suitable prior distribution p(θ) for the structural parameters and with the likelihood of step 1.

  * Generate a sample from the marginal posterior p(x_1:n| y_1:n) using the simulation smoother (Section <ref>) and for each sample use a different θ from the MCMC chain from the previous step.

  * For each state realization, x^*_1:n∼ p(x_1:n| y_1:n), from step 3., calculate a trend related or any other statistics of interest and use this sample for the estimates and their uncertainties.




§ EXAMPLES OF DIFFERENT DLM MODELS


In the following, we give several useful DLM formulations for model components that are typically used in geodetic or in more general environmental analyses. They have been used in existing applications for stratospheric ozone <cit.>, ionosonde analysis <cit.> and for station temperature records <cit.>. In Section <ref>, we will show analysis for synthetic GNSS station positioning time series.



 §.§ The effect of level and trend variance parameters


In the first example in Section <ref> the variance σ^2_trend was assumed to be known and fixed. Altering the variance affects the smoothness of the fit. In Figure <ref> the effect of different variance parameters are shown for the same data. Note that by setting both σ^2_level and σ^2_trend to zero results in classical linear regression without dynamical evolution of the regression components. In this case, the 95% probability limits for the level obtained from the smoother covariance matrix C_t coincide with the classical confidence intervals for the mean. In classical non-dynamic linear regression the modelling error is included in the residual term, whereas in DLM we can include it in the model definition by allowing temporal change in model parameters.

If we estimate the parameters by the likelihood approach and MCMC outlined in Section <ref>, we get the values in the last panel of Figure <ref> corresponding to the posterior mean. Figure <ref> shows MCMC chain histograms together with estimated marginal posterior densities. It also has the point values obtained by likelihood optimization. Note by optimization we get an estimate for σ_level which is very close to zero relatively to the MCMC solution, which tries to find all values of the parameter that are consistent with the data.









 §.§ Seasonal component


Seasonal variability can be modelled by adding extra state components for the effect of each season. A common description of seasonality uses trigonometric functions and is achieved by using two model states for each harmonic component. Monthly data with annual and semiannual cycles would use four state components and the following model and observation matrices


    M_seas=[  cos(π/6)  sin(π/6)         0         0; -sin(π/6)  cos(π/6)         0         0;         0         0  cos(π/3)  sin(π/3);         0         0 -sin(π/3)  cos(π/3);           ]

and

    H_seas= [ 1 0 1 0;  ].


In addition, a corresponding part or the model error covariance matrix Q_seas has to be set up to define the allowed variability in the seasonal amplitudes. A simple approach is to use a diagonal matrix with equal values for each component as diag(Q_seas) = [σ^2_seas,σ^2_seas,σ^2_seas,σ^2_seas]^T. If we set these variances to zero, the DLM algorithm will fit a temporally fixed seasonal amplitude.


For illustration we use a simulated monthly data with yearly variation that has some randomness in the amplitude. The observations have a piecewise linear trend similar to example in Section <ref> and some values as missing to see the effect on the uncertainties. We fit a seasonal component with one harmonic function, but we allow some variability in the amplitude and trend, with σ_trend = 0.005 and σ_seas = 0.4. Figure <ref> shows the generated data together with both the fitted mean process and the fitted seasonal component. A similar example was also used in <cit.>.






 §.§ Autoregressive process


Autoregressive processes have serial dependence between the observations. A general AR(p) process is defined by p coefficients [ρ_1,…,ρ_p] and an independent innovation term ϵ as

    y_t = ρ_1y_t-1 + ρ_2y_t-2 + … + ρ_py_t-p + ϵ,     ϵ∼ N(0,σ^2_AR)


For including an autoregressive component into the state space formulation we need to use state variables that "remember" their previous values. This can be achieved by suitable evolution operator M_AR. For example, AR(3) process with coefficients [ρ_1,ρ_2,ρ_3], will have three extra states with

    M_AR =
      [ ρ_1   1   0; ρ_2   0   1; ρ_3   0   0;     ],  H_AR =
      [ 1 0 0;   ],  Q_AR =
      [ σ^2_AR      0      0;      0      0      0;      0      0      0 ].


A pure AR(3) process would then be obtained by setting the observation error σ^2_obs in Eq. (<ref>) to zero and the model error component equal to the innovation variance σ^2_AR. 
If we, in addition, have σ^2_obs>0, it will result to an ARMA. In fact all ARMA and ARIMA models can be represented as DLM models <cit.>  and many ARIMA estimation software implementations use the Kalman filter likelihood Eq. (<ref>) to formulate the cost function for estimation.



 §.§ Regression covariates and proxy variables


In many applications the variability in the observations is affected by some known external factors, such as temperature, air pressure or solar activity. Sometimes these variables can be measured directly, as for the temperature, and sometimes their effect is modelled via a proxy, such as a radio fluxes for the solar effect. As an example, assume an observations model 

    y_t = μ_t + γ_t + β_tZ_t + ϵ_obs,

where μ_t and γ_t are the mean level and the seasonal components, Z_t is a row matrix of the values of the regression variables at time t, and β_t is a vector of time-varying regression coefficients.
The effect of the covariates can be formulated by having the coefficients as extra states, x_proxy,t = β_t, using an identity model operator, and by adding the covariate values to the observation operator H_t as

    H_proxy(t)   = Z_t = [ Z_t,1,…,Z_t,p;  ],
    M_proxy   = I_p= diag(1,…,1),
    Q_proxy   = diag([ σ^2_proxy,1,…,σ^2_proxy,2;  ]).

The DLM model for equation Eq. (<ref>) is then build up as diagonal block matrix combination of the components:

    x_t   = [ x_trend,t  x_seas,t x_proxy,t ]^T,
    M_t   =[ M_trend       0       0;       0  M_seas       0;       0       0 M_proxy ],
    H_t    = [    H_trend     H_seas H_proxy(t);  ],
    Q_t   =[ Q_trend       0       0;       0  Q_seas       0;       0       0 Q_proxy ].


The covariate variances σ^2_proxy control the allowed temporal variability in the coefficients β_t and their values can be estimated or set to some prior value. By setting the variances to zero, turns this model into classical multiple linear regression. 



§ SYNTHETIC GNSS EXAMPLE


Next we estimate trends in synthetic GNSS time series provided by Machiel Bos and Jean-Philippe Montillet. In this application, the trend estimated in the GNSS time series represents the tectonic rate on the East and North components and the vertical land motion on the Up coordinate. The characteristics of the GNSS time series are discussed in details in Chapter 1 and 2. We select data for one of the stations (labeled station n:o 3 in the figures) with the three components (East, North, Up) shown in Figure <ref>, top left panel. The time series are simulated using linear trend, yearly seasonal variation and a combination of coloured and i.i.d Gaussian noise. We assume that we do not know the noise structure a priori. We are interested in the (non-local) linear trend and we need a model component for the local fluctuations seen in the data.  This chosen data sets does not contain any sudden jumps in the measured position. Modelling offset changes would require a different strategy, with some iterative estimate of the jump locations, which we will not consider here.
We use a DLM approach, where we assume that the non-stationary part can be modelled by local polynomials and the stochastic stationary part can be described as an AR or ARMA process in addition to the i.i.d.  Gaussian observation uncertainty.
See <cit.> for a somewhat similar approach, which uses state space representation and Kalman filter likelihood to model flicker and random walk type noise in several stations at the same time.

So, in contrast to the spline smoothing example in Section <ref>, which had σ^2_level = 0 and σ^2_trend > 0, we will extract a non-local linear trend, σ^2_trend = 0, and model the local non-stationary fluctuations as a local level model with σ^2_level > 0. In addition, we use a yearly seasonal component for the daily observations and an autoregressive AR(1) noise component to account for the possible residual correlation. The observation error is assumed Gaussian and to have known standard deviation, σ_obs = 1 mm for components "East" and "North" and σ_obs = 4 mm for the "Up" component. The AR(1) innovation variance σ_AR as well as the AR coefficient ρ_AR will be estimated from the data. We use Kalman filter likelihood to estimate the 2 variance parameters and the AR(1) coefficient by MCMC. We analyse the three components (East, North, Up) separately. 

The true trend coefficients used in the simulation for the three data sets were give as 12.59, 17.64, and 2.778 mm/yr. The estimates obtained for them were 12.62 ± 0.61, 17.76 ± 0.69 and 2.22 ± 1.00 mm/yr, with one-sigma posterior standard deviations after ±.
Table <ref> shows the parameter estimates obtained by combination of Kalman simulation smoother for the linear slope and seasonal amplitude, and MCMC for θ = [σ_level, σ_AR, ρ_AR]^T. Figures <ref> and <ref> visualise the results graphically. There is a hint of negative autocorrelation in the ACF plot for the East components in Figure <ref>, but otherwise the residuals, obtained from the scaled prediction residuals, equation Eq. (<ref>), look very Gaussian. In overall, the selected DLM model seems to provide statistically consistent fit and reproduce the true trends within the estimated uncertainty.













§ COMPUTER IMPLEMENTATION



The examples and code to fit DLM models described here are available from a Github repository at <https://github.com/mjlaine/dlm>. The code is written in Matlab and it contains a reference implementation of Kalman filter, smoother and simulator algorithms as well as optimization and MCMC for the structural parameters. Other software implementations for DLM include state space models toolbox for Matlab described in <cit.>, a R package  described in <cit.> and python implementations in the  package <cit.>.



§ CONCLUSIONS



DLM provides a general framework for modelling many kinds of environmental time series, including geodetic ones. Some features of GNSS time series, such as the often assumed flicker noise and handling of offsets and data jumps might still require more special treatments. However, the DLM approach provides a very useful  generalization to the ordinary linear regression model. Its strengths include the ability to model non-stationary processes by allowing temporal change in the model coefficients and the direct modelling of the processes that generate the observed variability. By guiding the analysis in terms of the generating processes and their uncertainties it provides a good basis for Bayesian statistical inference. If there is prior knowledge about the changes, such as known change points in the data, they can be included in the model. By using simulation based Bayesian DLM analysis, your prior and posterior model simulations can be checked to be consistent with physical constrains and the observations.



§ REFERENCES




spbasic

