#1)
paretoDens=function(a,b,x){
  if(x<a){
    return(0)
  }else{
    return(b*(a^b)/(x^(b+1)))
  }
}

paretoInv=function(a,b,y){
  if(y>=1 || y<0){
    return(NaN)
  }else{
    return(a/((1-y)^(1/b)))
  }
}

paretoSample = function(a,b,n){
  return(paretoInv(a,b,runif(n)))
}

ML_est = function(sample){
  a=min(sample)
  n=length(sample)
  b=-n/(n*log(a)-sum(log(sample)))
  return(c(a,b))
}

empVaR=c()
ML_VaR=c()

for(i in 1:10000){
  S=paretoSample(1,3,50)
  empVaR=c(empVaR, quantile(S, 0.995))
  ab=ML_est(S)
  ML_VaR=c(ML_VaR, paretoInv(ab[1], ab[2], 0.995))
}

mybreaks=seq(0,max(empVaR)+2,2)
hist(empVaR, breaks=mybreaks)
mybreaks=seq(0,max(ML_VaR)+2,2)
hist(ML_VaR, breaks=mybreaks)
hist(empVaR-ML_VaR)

#2)
n=c(10,100,1000)
sigma2=c(100,10,1)
mu=5
lower=c()
upper=c()
for(i in 1:3){
  for(j in 1:3){
    x=rnorm(n[i], mu, sqrt(sigma2[j]))
    lower=c(lower, mean(x)-qnorm(0.975)*sqrt(sigma2[j]/n[i]))
    upper=c(upper, mean(x)+qnorm(0.975)*sqrt(sigma2[j]/n[i]))  
  }
}
lab=c()
for(i in 1:3){
  for(j in 1:3){
    text=paste("(", n[i], sep="")
    text=paste(text, ",", sep="")
    text=paste(text, sigma2[j])
    text=paste(text, ")")
    lab=c(lab, text)
    par(mar=c(8, 4, 4, 4))
    if(i==1 & j==1){
      plot(3*(i-1)+j, upper[3*(i-1)+j], col="red", pch=19, xlim=c(0,10), ylim=c(min(lower), max(upper)), lwd=2, cex=1.2, xaxt="n", xlab="", ylab="")
      lines(3*(i-1)+j, lower[3*(i-1)+j], col="red", pch=19, type="p", lwd=2, cex=1.2)
      lines(c(3*(i-1)+j, 3*(i-1)+j),c(lower[3*(i-1)+j], upper[3*(i-1)+j]), col="red")
    }else{
      lines(3*(i-1)+j, upper[3*(i-1)+j], col="red", pch=19, type="p", lwd=2, cex=1.2)
      lines(3*(i-1)+j, lower[3*(i-1)+j], col="red", pch=19, type="p", lwd=2, cex=1.2)
      lines(c(3*(i-1)+j, 3*(i-1)+j),c(lower[3*(i-1)+j], upper[3*(i-1)+j]), col="red")
    }
  }
}
axis(1, at=c(1:9), labels=lab, las=2)
abline(5,0,col="blue")