library(mvtnorm)
library(rgl)
source("Filled.contour2.R")  
# obtain source file from http://wiki.cbr.washington.edu/qerm/sites/qerm/images/4/44/Filled.contour2.R  

# get the data set and find the estimator
?trees
attach(trees)

fit <- lm(log(Volume) ~ log(Girth)+ log(Height), x=TRUE)
summary(fit)


# Sigma Estimate

Sigma <- (sd(fit$residuals)*sqrt(30/28))^2*solve(t(fit$x) %*% fit$x )
sqrt(Sigma[1,1])
sqrt(Sigma[2,2])
sqrt(Sigma[3,3])


###################################################################
# estimating the quantile

zet2 <- function(x){
	v1 = -c(1, log(5), log(50))
	v2 = -c(1, log(25), log(50))
	v3 = -c(1, log(5), log(100))
	v4 = -c(1, log(25), log(100))
    max = max((t(x)%*%v1)[1,1], (t(x)%*%v2)[1,1], (t(x)%*%v3)[1,1], (t(x)%*%v4)[1,1])    	  
    min = min((t(x)%*%v1)[1,1], (t(x)%*%v2)[1,1], (t(x)%*%v3)[1,1], (t(x)%*%v4)[1,1])    	  
	return(list(max=max, min=min))
	}

nn <- 100000

set.seed(12345)
sim <- rmvnorm(n=nn, mean=rep(0,3), sigma=31*Sigma, method="chol")
max.sim <- rep(0,nn)
for(i in 1:nn) max.sim[i] <- zet2(sim[i,])$max

hist(max.sim, col="orange")
q	<-	sort(max.sim)[0.95*nn]


###################################################################
# picture of the result

x	<-	seq(5, 25, length.out=100)
y	<-	seq(50, 100, length.out=100)

fhat	<-	outer(x, y)
for(i in 1:length(x)){
	for(j in 1:length(y)){
		fhat[i,j]	<-	6.631617-1.982650*log(x[i])-1.117123*log(y[j])
	}	
}

c	<-	-log(30)


par(pty="m", mar=c(2.5,2.5,0.5,0.5))
filled.contour2(x, y, fhat, levels=c(min(fhat), c+q/sqrt(31), max(fhat)), col=c(gray(0.75), "white"), lwd=2, lty=2)
contour(x, y, fhat, levels=c, col="black", lwd=2, drawlabels=FALSE, add=TRUE)
points(Girth[Volume<30], Height[Volume<30])
points(Girth[Volume>=30], Height[Volume>=30], pch=16)