# Set up the Gibbs sampler #

x<-0
y<-0
rho<-0.9
for (i in 1:10000){
  x<-c(x,rnorm(1,rho*y[i],sqrt(1-rho^2)))
  y<-c(y,rnorm(1,rho*x[i+1],sqrt(1-rho^2)))
}
vals<-cbind(x,y)

# Does the distribution look correct? #
#  (Naively using the entire chain.)

mean(vals)
cor(vals)
par(mfrow=c(2,2))
qqnorm(x)
qqline(x)
qqnorm(y)
qqline(y)
plot(x,y)


# Check the autocorrelation the hard way....# 

par(mfrow=c(3,3))
for (i in 1:9){
plot(x[1:(10001-i)],x[(1+i):10001])
}

cors<-rep(0,30)
for (i in 1:30){
cors[i]<-cor(x[1:(10001-i)],x[(1+i):10001])
}
cors


# Or use the built in time-series functions...#

par(mfrow=c(1,1))
acf(x)

# If we took every 20 or so and had a burn in of  #
#   1,000 that would give us....  #

x2<-x[1001+(0:450)*20]
y2<-y[1001+(0:450)*20]