#R code for computer practical 4

library(MASS)
#set the mean parameter
mu <- c(0,0)
#set the mean parameter
sigma <- matrix(c(4,1,1,4), nrow=2)

#generate a bivariate normal random variable by the above method
sigma.chol <- chol(sigma)
z <- rnorm(2)
y <- mu + t(sigma.chol) %*% z

#generate a sample of size 1000 using the R function
n <- 1000
#we sample directly from the bivariate Normal distribution; 
#in the practical, sampling is done approximately by the Gibbs sampling algorithm.
x <- mvrnorm(n, mu, sigma)

#plot the last 100 values, with subsequent values linked by a line
plot(x[901:1000,], type='b', xlab="X_1", ylab="X_2")

#look at sample plots of both variables
par(mfrow=c(2,1))
plot(x[,1], type="l", xlab="t", ylab="X_1")
plot(x[,2], type="l", xlab="t", ylab="X_2")

#compute the estimate of P(X_1 >=0, X_2 >=0)
prop <- rep(0, times=n)
prop[1] <- as.numeric((x[1,1] >= 0) && (x[1,2] >= 0))
for(i in 2:n){
  prop[i] <- prop[i-1] + as.numeric((x[i,1] >= 0) && (x[i,2] >= 0))
}

prop <- prop/1:n
plot(prop, ylim=c(0.1,0.4), type='l', xlab="t", ylab="Estimate of P(X_1>=0, X_2>=0)")