#
# 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) 

source("./util_fns.R")

# P 1: EPage 527
#
library("fEcofin")
library("mnormt")
Berndt = berndtInvest[,5:6]
names(Berndt)

source('../Chapter7/ml_fit_multivariate_t.R')
result = ml_fit_multivariate_t(Berndt)


# P 2:
# 

# Method I: Use fitdistr to estimate the t-distribution of the portfolio returns:
#
w = as.vector(c(0.3,0.7))
r_port = as.matrix(Berndt) %*% w
params = fitdistr( r_port, densfun="t" )
params = as.vector( params$estimate )
mu_1 = params[1] 
lambda_1 = params[2] # estimate of lambda = sqrt( (nu-2)/nu ) standard_deviation
nu_1 = params[3]
var_1 = ( lambda_1^2 ) * ( nu_1 / (nu_1-2) ) 

# Method II: Use projection: 
#
mu_2 = sum( result$fit$center * w )
nu_2 = result$nu
scale_2 = result$fit$cov 
var_2 = t(w) %*% ( ( nu_2 / (nu_2-2) ) * scale_2 ) %*% w # variance 
lambda_2 = sqrt( ( (nu_2-2)/nu_2 ) * var_2 )

list( m=mu_2, s=lambda_2, df=nu_2 ) # matches that from method I

# Estimate VaR and ES:
#
S = 100000
alpha = 0.05

VaR_1 = VaR_t(alpha,mu_1,lambda_1,nu_1,S)
VaR_2 = VaR_t(alpha,mu_2,lambda_2,nu_2,S)
print(sprintf("VaR_1= %10.6f; VaR_2= %10.6f", VaR_1, VaR_2))

ES_1 = ES_t(alpha,mu_1,lambda_1,nu_1,S)
ES_2 = ES_t(alpha,mu_2,lambda_2,nu_2,S)
print(sprintf("ES_1= %10.6f; ES_2= %10.6f", ES_1, ES_2))

# Check the formulas above with example 19.3 from the book 
#
mu_hat = 0.000689
lambda_hat = 0.007164
nu_hat = 2.984

VaR_t(0.05,mu_hat,lambda_hat,nu_hat,20000) # ~ 323.42
ES_t(0.05,mu_hat,lambda_hat,nu_hat,20000) # ~ 543.81 

# P 3
#
# Draw samples with replacement, use fitdistr to estimate paramters of the t-distribution, and compute the VaR^t
#
library(boot)

boot_fn = function(orig_data,inds){
  # Fit a t-distribution to this data and estimate VaR and ES:
  #
  data = unique( orig_data[inds] ) # duplicate values seem to give us trouble 
  params = fitdistr( data, densfun="t", start=list(m=0.015, s=0.077, df=5.38) )
  params = as.vector( params$estimate )
  mu = params[1] 
  lambda = params[2] # estimate of lambda = sqrt( (nu-2)/nu ) standard_deviation
  nu = params[3]

  VaR = VaR_t(0.05,mu,lambda,nu,100000)
  ES = ES_t(0.05,mu,lambda,nu,100000)
  
  c( mu, lambda, nu, VaR, ES )
}
boot_results = boot(r_port,boot_fn,R=250)


VaR_boots = boot_results$t[,4]
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter19/chap_19_prob_3_VaR_density_plot.eps", onefile=FALSE, horizontal=FALSE) }
plot( density(VaR_boots), lwd=2, main="VaR (kde)" )
if( save_plots ){ dev.off() }

shapiro.test(VaR_boots)


df_boots = boot_results$t[,3]
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter19/chap_19_prob_3_df_density_plot.eps", onefile=FALSE, horizontal=FALSE) }
plot( density(df_boots), lwd=2, main="df (kde)" )
if( save_plots ){ dev.off() }

shapiro.test(df_boots)


# P 4
#
# Use the returns from DEC to estimate the left tail index using the Hill estimator
#
DEC_rets = Berndt[,2]

sorted_rets = sort( DEC_rets )
Ys = abs(sorted_rets)
a_hill = rep( 0, length(Ys) ) 
for( nc in 1:length(Ys) ){
  c = Ys[nc]
  a_hill[nc] = nc / sum( log(Ys[1:nc]/c) )
}

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter19/chap_19_prob_4_hill_plot.eps", onefile=FALSE, horizontal=FALSE) }
plot( 1:length(Ys), a_hill, xlim=c(2,55), ylim=c(-0.1,4), xlab="nc", ylab="a", main="a^{Hill} estimate" )
grid()
if( save_plots ){ dev.off() }