# Bayesian estimation of Dynamic Linear Models with RStan

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 = $$ $$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  m0;
cov_matrix  C0;

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

m0  <- 0;
m0  <- 0;
m0  <- 0;
m0  <- 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> W_diag;

}
transformed parameters {
vector  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  <- 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);

}