Estimating a Bayesian TVP-VAR Model in R

This post guides through the estimation process of a Bayesian TVP-VAR model, where the algorithm proposed by Durbin & Koopman (2002) is used to draw samples of the parameter coefficients. It is intended to provide a brief, hands-on introduction to the implementation of multivariate state space models in R.

Time-varying parameter (TVP) methods become increasingly attractive in the realm of finance and macroeconomics. This seems straightforward, since macroeconomic relations are likely to evolve as a result of technological progress and changes in behavioural patterns. In empirical work changes in economic relations are frequently analysed with state space models. These models employ Kalman filtering and smoothing algorithms to obtain estimates of parameters, which are assumed to change over time. The filter and smoothing algorithm used in this post was originally proposed by Durbin & Koopman (2002) (DK), which is computationally more efficient than previously developed algorithms like Frühwirth-Schnatter (1994), Carter & Kohn (1994) or De Jong & Shephard (1995).

Each state space model consists of a measurement equation and a state equation. In the context of the TVP-VAR model considered in this post the measurement equation is y_t = Z_t b_t + \epsilon_t, where y_t is the vector of endogenous variables in period t, b_t is the k\times 1 vector of parameter coefficients and \epsilon_t \sim N(0,H). Z_t is a transformation of the vector of lagged endogenous variables and deterministic terms. Specifically, Z_t = {x_t}'\otimes I_n with x_t being a column vector. The state equation describes the evolution of the parameter values. In the present model we assume the standard form b_{t+1} = b_t + \eta_t with \eta_t \sim N(0,Q).

For the estimation of the model we use a Gibbs sampler which consists of the following blocks:

  • Draw initialisation values for the DK algorithm
  • Draw parameter coefficient values using the DK algorithm and the drawn initial values
  • Draw the covariance matrix of the state equation from an inverse Wishart distribution
  • Draw the covariance matrix of the measurement equation from an inverse Wishart distribution

Initialisation of the DK algorithm
A central challenge in estimating TVP models is the initialisation of the Kalman filter. It requires a vector of initial coefficient values a_0 from which the filter begins to run. Durbin & Koopman (2001) present methods how this can be done. The approach taken in this post follows Koop & Korobilis (2010). It starts by assuming that the values of b_t are deviations from the initial vector b_0. This leads to y_t = Z_t b_0 + Z_t b_t so that y^*_t = Z_t b_0 with y^*_t = y_t - Z_t b_t can be easily estimated. In this post I use an independent normal distribution to draw samples of the initial vector. For an introduction to multivariate Bayesian econometrics see Koop & Korobilis (2010), who I owe much of my understanding of the methods presented here.

for (i in 1:t){ # t is the length of the time series
y0[,i] <- y[,i]-Z[,,i]%*%b.filter[,i+1] # b.filter is the a_t
ZHZ <- matrix(0,n.vars,n.vars) # n.vars is k
ZHy <- matrix(0,n.vars,1)
for (i in 1:t){
ZHZ <- ZHZ + t(Z[,,i])%*%H.i%*%Z[,,i] # H.i is the inverse of H
ZHy <- ZHy + t(Z[,,i])%*%H.i%*%matrix(y[,i],n)
Vpost <- solve(Vprior.i + ZHZ,tol=1e-26)
bpost <- Vpost%*%(Vprior.i%*%bprior + ZHy)
b0 <- t(rmvnorm(1,bpost,Vpost))

Drawing parameter coefficient values
Before running the DK algorithm y_t has to be transformed to \hat{y}_t = y_t - Z_t a_0 and \hat{y}_t is used as the dependent variable in the DK algorithm, where a_0 is set to zero. P_0 can be an arbitrary symmetric covariance matrix. In this example P_0 = I_{k}. The code for the DK algorithm can be obtained from my GitHub repository.

for (i in 1:t) {[,i] <- y[,i]-Z[,,i]%*%b0
b.filter <- dk.alg2(k=n,m=n.vars,,Z=Z,H=H,Q=Q,TT=t,a1=0,P1=diag(1,n.vars))
b <- matrix(b0,n.vars,t+1) + b.filter

Drawing the covariance matrix of the state equation
Since the covariance matrix of the state parameters Q is not known, it has to be estimated. In this example I use an independent Wishart distribution with t+k+2 degrees of freedom and a prior variance 0.0001 I_{k}, which allows for relatively few variation of the state parameters.

DSum <- matrix(0,n.vars,n.vars)
for (i in 1:t){
res <- matrix(b[,i+1] - b[,i],n.vars)
DSum <- DSum + res%*%t(res)
Q <- solve(rWishart(1,t+n.vars+2,solve(Qprior + DSum))[,,1])

Drawing the covariance matrix of the measurement equation
Like the covariance matrix of the state equation I use an independent Wishart distribution with t+2 degrees of freedom and a prior variance 0.0001 I_{n} to draw samples of the measurement covariance matrix.

DSum <- matrix(0,n,n)
for (i in 1:t){
res <- y[,i] - Z[,,i]%*%b[,i+1]
DSum <- DSum + res%*%t(res)
H.i <- rWishart(1,t+n+2,solve(Hprior + DSum))[,,1]

I use seasonally adjusted and log-differenced, quarterly data on U.S income and consumption expenditures from 1947 to 2015.

The posterior medians (blue) and their respective 90% highest posterior density intervals (red) of the parameter estimates of the estimated TVP-VAR model with one lag and a constant are presented in figure 1. As a reference the OLS estimate of the respective parameter is added as a green horizontal line. The estimated model suggests that the negative sensitivity of income to own shocks has increased over the past 60 years, whereas it remained relative constant with respect to changes in consumption. Similarly, consumption reacts equally to shocks in income today than in the 1950s, but less to own shocks.

Figure 1

There are many possibilities how to extend this model which are not presented here. For further reading you might be interested in Primiceri (2005) who adds stochastic volatility to the TVP-VAR specification of this post. The paper belongs to the central contributions in the field of TVP. Fabian Krüger wrote a fast and easy to use package based on Primiceri’s model and its corrigendum 2015.

The code of this example is available at my GitHub repository.


Carter, C. K., & Kohn, R.. (1994). On Gibbs Sampling for State Space Models. Biometrika, 81, 541-553.
De Jong, P., & Shephard, N.. (1995). The Simulation Smoother for Time Series Models. Biometrika, 82, 339-350.
Del Negro, M., & Primiceri, G. E. (2015). Time Varying Structural Vector Autoregressions and Monetary Policy: A Corrigendum. Review Of Economic Studies, 82(4), 1342-1345. (Working paper version)
Durbin, J., & Koopman, S. J.. (2001). Time Series Analysis by State Space Models. Oxford.
Durbin, J., & Koopman, S. J.. (2002). A Simple and Efficient Simulation Smoother for State Space Time Series Analysis. Biometrika, 89(3), 603–615.
Frühwirth-Schnatter, S.. (1994). Data Augmentation and Dynamic Linear Models. Journal of Time Series Analysis, 15, 183-202.
Koop, G., & Korobilis, D.. (2010). Bayesian Multivariate Time Series Methods for Empirical Macroeconomics, Foundations and Trends in Econometrics, 3, 267-358.
Primiceri, G. E.. (2005). Time Varying Structural Vector Autoregressions and Monetary Policy. The Review of Economic Studies, 72(3), 821–852.



  1. Thanks a lot for this example! This is really helpful.
    On a related note, are you aware of any R implementation of panel BVAR model? I tried hard to search but so far I couldn’t find anything. I’m struggling with a dataset I have which has a small T (25) and large N (~1000) and my hypotheses need a VAR model. Perhaps the only way I can have reliable estimates is by using cross-sectional variation, which means using pVAR. But I also want to incorporate time-varying parameters, which probably needs a panel BVAR. I will really appreciate any help from you.

    1. Thank you! I’m glad that my code was helpful.
      Unfortunately, I’m not aware of any R implementation of Bayesian TVP-PVAR models. The reason for this might be that it is a very recent topic. Perhaps you will find the working paper by Koop and Korobilis ( helpful, if you decide to develop the code on your own. Usually, they also publish the MATLAB code for their papers. I have not found it yet, but I’m confident that it will be put online once the paper is accepted.

      1. Thanks! I didn’t know about this paper. I will check with Gary Koop whether he can share the Matlab code.
        Keep up the good work. Your blog is very informative. I have shared it with our PhD students at ESSEC.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s