• Using Bugs

##Model part
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]<-beta1*X[i,t-1] + beta2*Z[i] 
				X[i,t-1]~dnorm(0,1)
			}
			Z[i]~dbern(.5)
    }

   }

## Starting values
list(N=100, T=5, beta1=-.2,beta2=0.1,inv_sig_e2=1,inv_sig_v2=100, mu=c(10,-.2), inv_sig_LS=structure(.Data=c(.5,0,0,20),.Dim=c(2,2)))


## Run OpenBUGS scripts
[zzhang4@crcfe01 data]$ /pscratch/zzhang4/Private/bugs/bin/OpenBUGS
OpenBUGS version 3.2.1 rev 781
type 'modelQuit()' to quit
modelCheck('simmodel.txt')
modelData('simdata.txt')
modelCompile()
modelGenInits()
modelUpdate(1000)
samplesSet('X')
samplesSet('Y')
samplesSet('Z')
modelUpdate(1000)
samplesCoda('X','X')
samplesCoda('Y','Y')
samplesCoda('Z','Z')
modelQuit()

• Using R

library(mvtnorm)
library(rbugs)
gendata<-function(N, T, r){
	Y<-array(NA, c(N,T))
y<-array(NA, c(N,T))
d<-array(NA, c(N,T-1))
LS<-array(NA, c(N,2))
x<-array(rnorm(N*T), c(N,T))
z<-sample(c(0,1),N, replace=T)
sig_LS<-array(c(2, 0, 0, .15), c(2,2))
mu_LS<-c(3.5, 10)
sig_e2<-1
sig_e<-sqrt(sig_e2)
sig_v2<-.1
sig_v<-sqrt(sig_v2)
beta0<- -.2
beta1<- -.2
beta2<- .1

 for (i in 1:N){
       LS[i,]<-rmvnorm(1, mu_LS, sig_LS)
       y[i,1]<-LS[i,1]
       Y[i,1]<-rnorm(1, y[i,1],sig_e)

        for (t in 2:T){
        	 beta<-beta0 + beta1*x[i,t-1] + beta2*z[i] + rnorm(1,0,sig_v)
        	 d[i,t-1]<-beta*y[i,t-1]+LS[i,2]
             y[i,t]<-d[i,t-1]+y[i,t-1]
             Y[i,t]<-rnorm(1, y[i,t], sig_e)             
        }
    }
bugdata<-list(N=N, T=T, Y=Y, X=x, Z=z)
filename<-paste('simdata-N',N,'-T',T,'-',r,'.txt',sep='')
genDataFile(bugdata, filename)
	
}


for (N in c(50,100,200)){
	for (T in c(4,5,10)){
		for (r in 1:1000){
			gendata(N,T,r)
		}
	}
}