
model{
   for (i in 1:N){
       LS[i,1:2]~dmnorm(mu[1:2], inv_sig_LS[1:2,1:2])
       Y[i,1]~dnorm(y[i,1],inv_sig_e2)

       y[i,1]<-LS[i,1]

        for (t in 2:T){
             Y[i,t]~dnorm(y[i,t], inv_sig_e2)
             d[i,t-1]<-gamma[i,t-1]*y[i,t-1]+LS[i,2]
             y[i,t]<-d[i,t-1]+y[i,t-1]
				
				## another way to deal this
				gamma[i, t-1]~dnorm(mugamma[i,t-1], inv_sig_v2)
				mugamma[i, t-1]<-beta0+beta1*X[i,t-1] + beta2*Z[i]
        }
    }

       inv_sig_e2~dgamma(.001,.001)
       inv_sig_v2~dgamma(.001,.001)
       beta0~dnorm(0,.00001)
       beta1~dnorm(0,.00001)
       beta2~dnorm(0,.00001)
      
       mu[1]~dnorm(0,.00001)
       mu[2]~dnorm(0,.00001)

      inv_sig_LS[1:2,1:2]~dwish(sigs[1:2,1:2],2)
       sigs[1,1]<-1
       sigs[2,2]<-1
       sigs[1,2]<-sigs[2,1]
       sigs[2,1]<-0
      
       sig_e2<-1/inv_sig_e2
       sig_v2<-1/inv_sig_v2
       sig_LS[1:2,1:2]<-inverse(inv_sig_LS[1:2,1:2])
       
       para[1]<-mu[1]
       para[2]<-mu[2]
       para[3]<-sig_e2
       para[4]<-sig_v2
       para[5]<-sig_LS[1,1]
       para[6]<-sig_LS[1,2]
       para[7]<-sig_LS[2,2]
       
       para[8]<-beta0
       para[9]<-beta1
       para[10]<-beta2
}