#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

save_plots = F

set.seed(0) 

# P 1: EPage 
#
bmwRet = read.csv('../../BookCode/Data/bmwRet.csv')
library("fGarch")  # for skewed t functions 

n = dim(bmwRet)[1]

kurt = kurtosis(bmwRet[,2],method="moment")
skew = skewness(bmwRet[,2],method="moment")
fit_skewt = sstdFit(bmwRet[,2])

q.grid = (1:n)/(n+1)
qqplot( bmwRet[,2],
       qsstd(q.grid, fit_skewt$estimate[1], fit_skewt$estimate[2], fit_skewt$estimate[3], fit_skewt$estimate[4]),
       ylab="skewed-t quantiles" )


# P 2: EPage 
#
tQuantKurt = function(nu,p1=0.025,p2=0.25){
    # for a t-distribution we have: 
    numer = qt(1-p1,nu) - qt(p1,nu)
    denom = qt(1-p2,nu) - qt(p2,nu)
    k = numer / denom     
}
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter6/rlab_prob_2.eps", onefile=FALSE, horizontal=FALSE) }
nu_values = seq( 1, 10, by=0.25 )
plot( nu_values, tQuantKurt(nu_values), type="l", xlab=expression(nu) )
if( save_plots ){ dev.off() }


# P 3: EPage 
# 
quantKurt = function(y,p1=0.025,p2=0.25){
    Q = quantile(y,c(p1,p2,1-p2,1-p1))
    k = (Q[4]-Q[1])/(Q[3]-Q[2])
    k 
}
nboot = 5000
ModelFree_kurt = rep(0,nboot)
ModelBased_kurt = rep(0,nboot)

set.seed(5640)
for( i in 1:nboot ){
    samp_ModelFree = sample( bmwRet[,2], n, replace=TRUE )
    samp_ModelBased = rsstd(n, fit_skewt$estimate[1], fit_skewt$estimate[2], fit_skewt$estimate[3], fit_skewt$estimate[4])
    ModelFree_kurt[i] = quantKurt(samp_ModelFree)
    ModelBased_kurt[i] = quantKurt(samp_ModelBased)
}

d1 = density(ModelFree_kurt)
d2 = density(ModelBased_kurt)
x_min = min(d1$x,d2$x)
x_max = max(d1$x,d2$x)
y_min = min(d1$y,d2$y)
y_max = max(d1$y,d2$y)

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter6/rlab_prob_3.eps", onefile=FALSE, horizontal=FALSE) }
plot( d1$x, d1$y, type="l", xlim=c(x_min,x_max), ylim=c(0,y_max), xlab="x", ylab="y" )
lines( d2$x, d2$y, lwd=2, lty=5 )
legend("topright",c("Model Free","Model Based"),lwd=2,lty=c(1,5))
if( save_plots ){ dev.off() }

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter6/rlab_prob_3_boxplots.eps", onefile=FALSE, horizontal=FALSE) }
df = data.frame( kurt=c( ModelFree_kurt, ModelBased_kurt ), model=c( rep( 1, length(ModelFree_kurt) ), rep( 2, length(ModelBased_kurt) ) ) )
boxplot( kurt ~ model, data=df ) 
if( save_plots ){ dev.off() }


# P 4: EPage 
#
quantile( ModelFree_kurt, c(0.05,1-0.05) )
quantile( ModelBased_kurt, c(0.05,1-0.05) )


# P 5: EPage 
#
if( !require('bootstrap') ){
    install.packages('bootstrap', dependencies=TRUE, repos='http://cran.rstudio.com/')
}
library(bootstrap)
set.seed("1234")
bc_out = bcanon( bmwRet[,2], 5000, quantKurt, alpha=c(0.05,0.95) )
print( bc_out$confpoints )