ARMA(1,1) as part of hierarchical model term

Modeling
14

sure.


data{
  int<lower=0> T;
  int<lower=0> N;
  matrix[N,T]  x;
  matrix[N,T]  y;
}

parameters{
  real                    beta;
  real<lower=-1, upper=1> phi;
  real                    theta;
  real<lower=0,upper=5>   sig2;
  real<lower=0,upper=5>   sig2v;
  vector[N]               v;
  real<lower=0>           sig_a;
  vector[T]               a_std;
}

transformed parameters {
  vector[T] eta;
  vector[T] a = a_std * sig_a; // non-centered innovations a

  eta[1] = a[1]; // Not sure for boundary
  for(t in 2:N) {
    // assume phi/beta/theta are given as data
    eta[t] = phi * eta[t-1] + theta * a[t-1] + a[t];
  }
}

model{
  for(i in 1:N){
    for(t in 1:T){
      y[i,t] ~ normal(x[i,t]*beta+v[i]+eta[t], sig2);
    }
  }
  sig2 ~ uniform(0,5);
  for(i in 1:N){    
    v[i] ~ normal(0,sig2v);
  }
  sig2v ~ uniform(0,5);
  phi ~ uniform(-1,1);
  theta ~ normal(0,2);
  a_std ~ normal(0, 1);
  sig_a ~ normal(0, 1);
}

And this is my simulation data part:

N = 10
T = 120
phi = 0.3
theta = 0.5
sig2 = 0.4
sig2v = 0.1
sig2eta = 0.2
beta = 0.3
set.seed(123)
eta = arima.sim(n = T, list(ar = phi, ma = theta),sd = sig2eta)
v = rnorm(N,mean=0,sd=sig2v)
epsilon = matrix(NA,nrow=N,ncol=T)
for(i in 1:N){
  epsilon[i,] = rnorm(T, mean=0, sd=sig2)}
x = matrix(NA,nrow=N,ncol=T)
truew = matrix(NA,nrow=N,ncol=T)
for(i in 1:N){
  for(t in 1:T){
    x[i,t] = rnorm(1,mean=1,sd=1)
    truew[i,t] = x[i,t]*beta+v[i]+eta[t]+epsilon[i,t]
  }
}

state_test = list(N=N,T=T,y=truew,x=x)
test = stan(file = "try.stan",dat = state_test,  
             pars = c("phi","theta","sig2","sig2v","beta"),
             iter=2000,warmup=500,control = list(adapt_delta=0.99,max_treedepth = 10))
test_result = extract(test,permute=TRUE)  
test