* 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) } } }