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

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

save_plots = F

set.seed(0) 

source("./util_fns.R")

## Problem 1:
## 
CokePepsi = read.table('../../BookCode/Data/CokePepsi.csv', header=TRUE)
price = CokePepsi[,1]
returns = diff(price)/lag(price)[-1]
ts.plot(returns)

S = 4000
alpha = 0.05
library(MASS)
res = fitdistr(returns,'t')
mu = res$estimate['m']
lambda = res$estimate['s']
nu = res$estimate['df']
print(qt(alpha, df=nu))
print(dt(qt(alpha, df=nu), df=nu))


## Problem 2:
##
v = VaR_t(alpha, mu, lambda, nu, S)
es = ES_t(alpha, mu, lambda, nu, S)
print(sprintf('VaR_t= %f; ES_t= %f', v, es))


## Problem 3:
##
library(fGarch) # for qstd() function
library(rugarch)
garch.t = ugarchspec(
    mean.model=list(armaOrder=c(0, 0)),
    variance.model=list(garchOrder=c(1, 1)),
    distribution.model='std')
KO.garch.t = ugarchfit(data=returns, spec=garch.t)
show(KO.garch.t)
plot(KO.garch.t, which=2)

pred = ugarchforecast(KO.garch.t, data=returns, n.ahead=1); pred 
fitted(pred); sigma(pred) # extraction functions to get the conditional mean and the conditional standard deviation

## Problem 5:
## 
nu = as.numeric(coef(KO.garch.t)[5])
q = qstd(alpha, mean=fitted(pred), sd=sigma(pred), nu=nu); q # the quantile of the return associated with this value of alpha 
VaR = -S * q

lambda = sigma(pred)/sqrt( nu/(nu-2) ) # the conditional scale parameter 
es = ES_t(alpha, fitted(pred), lambda, nu, S)
print(sprintf('VaR_t= %f; ES_t= %f', VaR, es))



# P 6:
#
berndtInvest = read.csv("../../BookCode/Data/berndtInvest.csv", header=TRUE)
Berndt = berndtInvest[,5:6]
names(Berndt)

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


# P 7:
# 

# 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 8
#
# 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 9
#
# 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() }