####WinBUGS model specification for a bivariate multilevel autoregressive model with a Wishart prior for the
covariance matrix of the random parameters (scale matrix is calculated in the accompanying Rscript).####

model{
   
###model for timepoint 2 to nt######

	for (j in 1:np)
	{ 
      		for (t in ((j-1)*tt+2):(j*tt))
		{
           		f[t,1:2]~dmnorm(muf[t,1:2],Qpre[1:2,1:2]) ###the data is multivariate normal distributed
           		 muf[t,1]<-b[j,1] + b[j,3]*z[t-1,2] + b[j,5]*z[t-1,1] ###for Y1
				muf[t,2]<-b[j,2] + b[j,6]*z[t-1,2] + b[j,4]*z[t-1,1]  ###for Y2
           		 z[t,1]<- f[t,1] - b[j,1]
			     z[t,2]<- f[t,2] - b[j,2]
           		 
		}

 

	}


####Model specification for time point 1###
	for (j in 1:np)
	{ 
      	
  	
	f[((j-1)*tt+1),1:2]~dmnorm(muf1[j,1:2],Qpre[1:2,1:2]) ###timepoint 1
	z[((j-1)*tt+1),1]<- f[((j-1)*tt+1),1] - b[j,1]  ###timepoint 1 z
	z[((j-1)*tt+1),2]<- f[((j-1)*tt+1),2] - b[j,2]  ###timepoint 1 z
        	 
	muf1[j,1]<-b[j,1] + b[j,3]*z0[j,2] + b[j,5]*z0[j,1] 
	muf1[j,2]<-b[j,2] + b[j,6]*z0[j,2] + b[j,4]*z0[j,1]

	z0[j,1] ~dnorm(0,.5)
	z0[j,2] ~dnorm(0,.5)
	
          
	}   	 


####priors###########
###for the 2x2 covariance matrix of the innovations, we can avoid the Wishart using uniform priors on the variances and correlation###     
        
Qpre[1:2,1:2]<- inverse(Qvar[1:2,1:2]) ###Qpre is the precision matrix for the innovations
Qvar[1,1] ~dunif(0,10)
Qvar[2,2] ~dunif(0,10)
Qvar[1,2] <- Qcor*sqrt(Qvar[1,1])*sqrt(Qvar[2,2])
Qcor ~ dunif(-1,1)
Qvar[2,1] <- Qvar[1,2]


###prior for random parameters###
	for (j in 1:np)
	{
     b[j,1:6]~dmnorm(bmu[1:6],bpre[1:6,1:6]) ## the means and regression coefficients are multivariate normal distributed
    }

####priors for the fixed effects###
bmu[1]~dnorm(0,.000000001) 
bmu[2]~dnorm(0,.000000001)
bmu[3]~dnorm(0,.000000001)
bmu[4]~dnorm(0,.000000001)
bmu[5]~dnorm(0,.000000001)
bmu[6]~dnorm(0,.000000001)

###priors for the precision matrix, calculate the covariance matrix and correlations
bpre[1:6,1:6]~dwish(W[,],df) ###the scale matrix W is calculated in R and then fed to winbugs.
df<- 6
bcov[1:6,1:6] <- inverse(bpre[1:6,1:6])


	for (d in 1:6)
	{
		for (g in 1:6){
		bcor[d,g] <- bcov[d,g] / ( sqrt(bcov[d,d]) * sqrt(bcov[g,g]) )
		}
	}
}