# Parametric bootstrap of Law Data using fitted
# Multivariate normal distribution

# Initialisation
library(bootstrap)
library(MASS)
B<-2500

# First we calculate the sample mean vector and sample
# variance/covariance matrix for the original
# law data in `law'.

mu<-mean(law)
Sigma<-var(law)

# Then we generate B samples from the fitted bivariate normal
# and accumulate the resulting parametric bootstrap replicates
# for the correlation coefficient.

pbrep<-rep(0, B)
for(i in 1:B){                             
 pbsample<-mvrnorm(15,mu,Sigma)
 pbrep[i]<-cor(pbsample[,1],pbsample[,2])
}
hist(pbrep, xlim=c(-1.05,1.05), breaks=seq(-1.05, 1.05, by=0.05))
sd(pbrep)
mean(pbrep) - Sigma[1,2]/sqrt(Sigma[1,1]*Sigma[2,2])