###################
# Simulation code #
###################

# Assumed true parameters

beta0<-0
beta1<-1
beta2<-0.5
sigma<-1
N<-100		#1700
nburn<-10		#A number of the first simulated obs to be burnt=200

truevalue<-c(alpha0, alpha1, alpha2, sigma2)
library(boot)


# Bayesian simulation using WinBugs

B<-1   #200
BUGSest<-matrix(0,nrow=B,ncol=4)  #ncol=4 parameters of the distribution
out<-matrix(0,nrow=B,ncol=3)
time<-matrix(1,nrow=N,ncol=1)
X<-matrix(0,nrow=(N-nburn),ncol=2)
Y<-matrix(0,nrow=(N-nburn),ncol=1)
Y1<-matrix(0,nrow=(N-nburn),ncol=1)

# Data simulation
for (j in 1:B){
	y<-matrix(0,nrow=N,ncol=1)
	ylag<-matrix(0,nrow=N,ncol=1)
	mu<-matrix(0,nrow=N,ncol=1)
	t<-1
	mu[t]<-beta0+beta1*t+beta2*0		#y[t-1] is unknown
	y[t]<-rnorm(1,mu[t],sigma)

		for (t in 2:N){
			mu[t]<-beta0+beta1*t+beta2*y[t-1]
			y[t]<-rnorm(1,mu[t],sigma)
			ylag[t]<-y[t-1]
			time[t]<-t
		}
    
	Y<-y[(nburn+1):N]
	Y1<-ylag[(nburn+1):N]

# BUGS regression estimation
T<-N-nburn
data<-list("T","Y", "Y1")		
inits<-function()
list(tau.y=rgamma(1,0.001,0.001), 
beta0.c=rnorm(1,0,1000),
beta0.tau=rgamma(0.001, 0.001), 
beta1.c=rnorm(1,0,1000), 
beta1.tau=rgamma(0.001, 0.001),
beta2.c=rnorm(1,0,1000),
beta2.tau=rgamma(0.001, 0.001)) 

parameters<-c("beta0", "beta1", "beta2", "sigma.y", "sigma.beta0", 
"sigma.beta1", "sigma.beta2")

inits1<-list(beta0=0, beta1=0, beta2=0, sigma.y=1, sigma.beta0=1, sigma.beta1=1, sigma.beta2=1) 
inits2<-list(beta0=0, beta1=0, beta2=0, sigma.y=10, sigma.beta0=10, sigma.beta1=10, sigma.beta2=10) 
inits3<-list(beta0=0, beta1=0, beta2=0, sigma.y=100, sigma.beta0=100, sigma.beta1=100, sigma.beta2=100) 
inits<-list(inits1, inits2, inits3)

# Bugs input in R#
library("arm")
AR1bugs.sim <- bugs(data, inits, parameters, "D:AR1GammaModel.txt", n.chains=3, n.iter=1000,    bugs.directory="c:/Program Files/WinBUGS14/",
    working.directory=NULL, clearWD=TRUE, debug=TRUE)

attach(AR1bugs.sim)
print (AR1bugs.sim)
plot (AR1bugs.sim)
}