## com13a of october 2023, 
## for kiosk and for stk4011, exercise 3.6 of PartOne  

#  investigations with bigger sim than 100 and smaller n than 500: 
#  partly: \hat\gamma_4 underestimates \gamma_\true = 6.000
#  partly: too low coverage 

inte <- function(x)
{
(x-1)^4*dexp(x)
}

integrate(inte,0,Inf) # is 9 

kurt <- 6 

sim <- 10^2 # low for figure, big for estimating coverage 
nn <- 500 # i think i choose n = 500 for the figure in kiosk 
rootn <- sqrt(nn) 

sigmahatsim <- 0*(1:sim) 
gamma4hatsim <- 0*(1:sim)
lowsim <- 0*(1:sim)
upsim <- 0*(1:sim)

for (ss in 1:sim)
{
xx <- rexp(nn) 
xbar <- mean(xx) 
sigmahatsim[ss] <- sd(xx)
# 
ww <- (xx-xbar)/sigmahatsim[ss]
gamma4hatsim[ss] <- mean(ww^4)-3
# 
lowsim[ss] <- sigmahatsim[ss]-1.645*sqrt(0.50 + 0.25*gamma4hatsim[ss])/rootn
upsim[ss] <- sigmahatsim[ss]+1.645*sqrt(0.50 + 0.25*gamma4hatsim[ss])/rootn
}

plot(density(gamma4hatsim))

print(c(nn,sim,mean(gamma4hatsim)))

# quite a bit too low: \E\,\kurt < \kurt for moderate n 

# largesample_figure13B.pdf for kiosk Ch 2 : 

par(bty="l")
par(cex.lab=1.4142) 

matplot(1:sim,cbind(lowsim,upsim),type="l",
   xlab="simulations",
   ylab="lower, upper",lwd=3)
matlines(1:sim,1+0*(1:sim),col=4,lty=2,lwd=3)

print( c(nn,sim,mean(1*(1 <= upsim)*(lowsim <= 1))) )

# n    sim   coverage mean(\hat\gamma) 
# 100, 10^6, 0.821    4.020 
# 500, 10^6, 0.875    5.472 

# so \hat\gamma_4 underestimates the true \gamma_4 = 6 

# 10^4: 0.8731 
# 10^5: 0.8761