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