if( !require('bootstrap') ){
    install.packages('bootstrap', dependencies=TRUE, repos='http://cran.rstudio.com/')
}
library(bootstrap)

Kurtosis = function(x){
    mean( ( (x-mean(x)) / sd(x) )^4 )
}

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
}

set.seed(3751)
niter = 500
nboot = 400
n = 50
nu = 10
trueKurtosis = 3 + 6/ (nu-4)
fn_to_bootstrap = Kurtosis
#fn_to_bootstrap = quantKurt # used for the second part of the question 
correct = matrix(nrow=niter, ncol=5)
width = matrix(nrow=niter, ncol=5)
error = matrix(nrow=niter, ncol=1)
t1 = proc.time()
for( i in 1:niter ){
    y = rt(n, nu)

    # Method 1:
    # 
    int1 = boott(y, Kurtosis, nboott=nboot, nbootsd=50)$confpoints[c(3, 9)]
    width[i, 1] = diff(int1)
    correct[i, 1] = as.numeric((int1[1] < trueKurtosis) & (trueKurtosis < int1[2]))

    # Method 2:
    # 
    int2 = bcanon(y, nboot, Kurtosis)$confpoints[c(1, 8), 2]
    width[i, 2] = diff(int2)
    correct[i, 2] = as.numeric((int2[1] < trueKurtosis) & (trueKurtosis < int2[2]))

    # Method 3:
    # 
    boot = bootstrap(y, nboot, Kurtosis)$thetastar
    int3 = Kurtosis(y) + 1.96 * c(-1, +1) * sd(boot)
    width[i, 3] = diff(int3)
    correct[i, 3] = as.numeric((int3[1] < trueKurtosis) & (trueKurtosis < int3[2]))

    # Method 4:
    # 
    int4 = quantile(boot, c(0.025, 0.975))
    width[i, 4] = diff(int4)
    correct[i, 4] = as.numeric((int4[1] < trueKurtosis) & (trueKurtosis < int4[2]))

    # Method 5:
    # 
    int5 = 2*Kurtosis(y) - quantile(boot, c(0.975, 0.025))
    width[i, 5] = diff(int5)
    correct[i, 5] = as.numeric((int5[1] < trueKurtosis) & (trueKurtosis < int5[2]))

    error[i] = mean(boot) - Kurtosis(y) 
}
t2 = proc.time()

options(digits=3)
print((t2 - t1)/60)
print(colMeans(width))
print(colMeans(correct))
print(mean(error))
print(mean(error^2))

# P 10:
#
MSE = mean(error)^2 + mean(error^2) # bias^2 + s^2
print(MSE)