Bayesian estimation of Dynamic Linear Models with RStan

Question

I'm reading the Dynamic Linear Models with R book, where most of chapter 4 is devoted to bayesian estimation of parameters. They code most of it manually though, and it seems it can get quite tricky for complicated models. I wonder if this can be more easily done in Stan/Bugs. For example, on p.186, the model to be estimated consists of local linear trend and AR(2) components: \begin{align*} \mathbf{y}_t &= \mathbf{F}\boldsymbol{\Theta}_t + \boldsymbol{\nu}_t, \qquad \boldsymbol{\nu}_t \sim \mathcal{N}(0, \mathbf{V})\\ \boldsymbol{\Theta}_t &= \mathbf{G}\boldsymbol{\Theta}_{t-1} + \boldsymbol{\Omega}_t, \qquad \boldsymbol{\Omega}_t \sim \mathcal{N}(0, \mathbf{W}) \end{align*} where $$F = [1, 0, 1, 0]$$ $$G = \left[\begin{array}{ccc} 1 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & \phi_1 & 0 \\ 0 & 0 & \phi_2 & 0\end{array}\right]$$ $$ V = [0] $$ $$ W = diag(\sigma_e^2,\sigma_z^2,\sigma_u^2,0 )$$

Priors for ($\phi_1$, $\phi_2$) are $ N(0,(2/3)^2)$ and $ N(0,(1/3)^2)$ and the inverses of variances are assumed to be independent with Gamma priors $g(a^2/b, a/b)$ where a = 1 and b =1000.

My question is, can this model be estimated using Stan/Bugs, without manual tweaking of the gibbs sampler? EDIT: Here is my attempt at it using Stan (data). Has a lot of convergence problems.

data {
  int <lower = 0> N;
  matrix [1, N] y;
}
transformed data {
  matrix [4, 1] F;
  vector [4] m0;
  cov_matrix [4] C0;

  F [1, 1] <- 1;
  F [2, 1] <- 0;
  F [3, 1] <- 1;
  F [4, 1] <- 0;

  m0 [1] <- 0;
  m0 [2] <- 0;
  m0 [3] <- 0;
  m0 [4] <- 0;

  C0 <- diag_matrix(rep_vector(1.0e+7, 4));
}
parameters {
  // need to impose stationarity constraints
 real <lower=-1,upper=1> phi2; 
 real <upper=(1 - fabs(phi2))> phi1; 


  real <lower = 0> sigma ; // for V
  real<lower = 0> W1;
  real<lower = 0> W2;
  real<lower = 0> W3;

  // vector<lower = 0>[4] W_diag;

}
transformed parameters {
  vector [1] V;
  matrix [4, 4] W;
  matrix [4, 4] G;

  G [1, 1] <- 1;
  G [1, 2] <- 1;
  G [1, 3] <- 0;
  G [1, 4] <- 0;
  G [2, 1] <- 0;
  G [2, 2] <- 1;
  G [2, 3] <- 0;
  G [2, 4] <- 0;  
  G [3, 1] <- 0;
  G [3, 2] <- 0;
  G [3, 3] <- phi1;
  G [3, 4] <- 1;
  G [4, 1] <- 0;
  G [4, 2] <- 0;  
  G [4, 3] <- phi2;
  G [4, 4] <- 0;


  V [1] <- sigma  * sigma ;
  // W <- diag_matrix(W_diag);
  W [1, 1] <- W1;
  W [1, 2] <- 0;
  W [1, 3] <- 0;
  W [1, 4] <- 0;
  W [2, 1] <- 0;
  W [2, 2] <- W2;
  W [2, 3] <- 0;
  W [2, 4] <- 0;
  W [3, 1] <- 0;
  W [3, 2] <- 0;
  W [3, 3] <- W3;
  W [3, 4] <- 0;
  W [4, 1] <- 0;
  W [4, 2] <- 0;
  W [4, 3] <- 0;
  W [4, 4] <- 0;

}
model {
  sigma ~ uniform (0, 5);
  phi2 ~ normal(0,2.0/3);
  phi1 ~ normal(0,1.0/3);
  // W[1,1] ~ inv_gamma(0.001, 0.001);
  // W[2,2] ~ inv_gamma(0.001, 0.001);
  // W[3,3] ~ inv_gamma(0.001, 0.001);
  W1 ~ inv_gamma(1, 0.001);
  W2 ~ inv_gamma(1, 0.001);
  W3 ~ inv_gamma(1, 0.001);
  y ~ gaussian_dlm_obs (F, G, V, W, m0, C0);

}

Answer 1

It is perhaps easier to emphasize what models the Stan and the BUGS family cannot estimate. Stan cannot estimate models with discrete unknowns, and the BUGS family cannot estimate models that cannot be written in the form of a Directed Acyclic Graph. Thus, a DLM can, in principle, be estimated with either. Whether you can get good samples from the posterior distribution in practice is another question.