######################################################
#
#  Unit 2 - Example ii.
#
#  Monte Carlo integration 
#  Monte Carlo via Importance function integration
#
######################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
######################################################
set.seed(123456)

x = seq(-10,10,length=1000)
par(mfrow=c(3,2))
plot(x,dt(x,1),xlab="",ylab="",main="X ~ CAUCHY(0,1)",type="l")
plot(x,dt(x,1),xlab="",ylab="",main="P(X>2)",type="l")
for (i in seq(2,10,by=0.1))
  segments(i,0,i,dt(i,1))
plot(x,dt(x,1),xlab="",ylab="",main="P(X>2)+P(X<-2)",type="l")
for (i in seq(2,10,by=0.1)){
  segments(i,0,i,dt(i,1))
  segments(-i,0,-i,dt(i,1))
}
u = seq(0,0.5,length=1000)
plot(u,1/(2*pi*(1+u^(-2)))/4,xlab="",ylab="",main="Transformation",type="l")
for (i in seq(0,0.5,by=0.005))
  segments(i,0,i,1/(2*pi*(1+i^(-2)))/4)

Ms = c(seq(100,1000,by=100),seq(2000,10000,by=1000),100000,1000000)
stats = NULL
for (M in Ms){
  theta = rt(M,1)
  h1    = rep(0,M)
  h1[theta>2]=1
  p1    = mean(h1)
  var1  = mean((h1-p1)^2)/M
  sd1   = sqrt(var1)
  h2 = rep(0,M)
  h2[theta>2]=1/2
  h2[theta<(-2)]=1/2
  p2    = mean(h2)
  var2  = mean((h2-p2)^2)/M
  sd2   = sqrt(var2)
  u    = 0.5*runif(M)
  h3   = u^(-2)/(2*pi*(1+u^(-2)))
  p3   = mean(h3)
  var3 = mean((h3-p3)^2)/M
  sd3   = sqrt(var3)
  stats = rbind(stats,c(c(p1,p2,p3),c(var1,var2,var3),c(sd1,sd2,sd3)))
}
stats = cbind(Ms,stats)

plot(log10(stats[,1]),stats[,2],xlab="log10(n)",ylab="",main="Monte Carlo estimators",type="l")
lines(log10(stats[,1]),stats[,3],col=3,lty=2)
lines(log10(stats[,1]),stats[,4],col=4,lty=3)
abline(h=1-pt(2,1),col=2)
plot(log10(stats[,1]),stats[,8],xlab="log10(n)",ylab="",main="Standard deviations of \n Monte Carlo estimators",ylim=range(stats[,8:10]),type="l")
lines(log10(stats[,1]),stats[,9],col=3,lty=2)
lines(log10(stats[,1]),stats[,10],col=4,lty=3)

stats[,c(1,8,9,10)]