0
$\begingroup$

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);

}
$\endgroup$

1 Answer 1

3
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Thanks for you comment. As you said, getting good samples seems to be a bit problematic. I've added my Stan code, would appreciate it if you take a look at it (for obvious errors) $\endgroup$
    – Alex
    Commented Jul 1, 2016 at 18:16
  • $\begingroup$ I believe that is considered a bit off-topic for this website but if anyone else is interested, the code is being discussed on Stan users. $\endgroup$ Commented Jul 4, 2016 at 21:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.