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

# Install some needed packages if they are not already installed:
# 
if( "fEcofin" %in% rownames(installed.packages()) == FALSE ){
    install.packages("fEcofin")
}

set.seed(0)

source("./util_fns.R")

# E 1 EPage 503
# 
library(MASS)
library(fEcofin)

rets = bmwRet[,2]

S = 100000
alpha = 0.01 

# nonparametric:
#
t = quantile( rets, probs=alpha )
VaR_np = -S * t
ES_np = -S * mean( rets[ rets < t ] )
as.vector( c( VaR_np, ES_np ) )

# normal assumption:
# 
m = mean(rets)
s = sd(rets)
c( VaR_norm(alpha,m,s,S), ES_norm(alpha,m,s,S) )

# t-distribution assumption:
#
t_params = fitdistr( rets, densfun="t" )
t_params = as.vector( t_params$estimate ) 
c( VaR_t(alpha,t_params[1],t_params[2],t_params[3],S), ES_t(alpha,t_params[1],t_params[2],t_params[3],S) )


# E 2
#
252 * ( 0.05 / 0.005 )^(1/3.1) 


# E 3
#
library(tseries)
ticker = "IBM"
data = get.hist.quote(instrument=ticker,quiet=TRUE)
log_rets = diff(log(data$Close))
rets = log_rets[(length(log_rets)-999):length(log_rets)] # get a sample of 1000 returns
rets = as.vector( rets ) 

S = 10000
alpha = 0.025

# a) t-distribution assumption:
#
library(MASS)
t_params = fitdistr( rets, densfun="t" )
t_params = as.vector( t_params$estimate ) 
VaR_t_pt_a = VaR_t(alpha,t_params[1],t_params[2],t_params[3],S)

# b) nonparametric:
#
t = quantile( rets, probs=alpha )
VaR_np = -S * t
ES_np = -S * mean( rets[ rets < t ] )
VaR_np

# c) t-plot (following the examples in the Rlab in Chapter 4): 
#
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter19/chap_19_exer_3_pt_c_norm_qqplot.eps", onefile=FALSE, horizontal=FALSE) }
qqnorm(rets)
qqline(rets)
grid()
if( save_plots ){ dev.off() }

n = length(rets)
q.grid = (1:n)/(n+1)
df = t_params[3]
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter19/chap_19_exer_3_pt_c_t_qqplot.eps", onefile=FALSE, horizontal=FALSE) }
qqplot(rets, qt(q.grid,df=df), main=paste("QQ plot for", ticker, ": df=", df))
abline(lm( qt(c(0.25,0.75),df=df) ~ quantile(rets,c(0.25,0.75)) ))
grid()
if( save_plots ){ dev.off() }

# d) estimate left-tail index a (using the Hill estimator):
#
sorted_rets = sort( 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_exer_3_pt_d_hill_plot.eps", onefile=FALSE, horizontal=FALSE) }
plot( 1:length(Ys), a_hill, xlim=c(2,200), ylim=c(-0.1,4), xlab="nc", ylab="a", main="a^{Hill} estimate" )
grid()
if( save_plots ){ dev.off() }

a = mean( a_hill[50:125] )
VaR_t_pt_a * ( 0.025 / 0.0025 )^(1/a) 


# E 4
#
library("fEcofin")
library(boot)
close = msft.dat$Close
rets = diff(log(close))

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.001589855, s=0.024606088, df=3.998675607) )
  params = as.vector( params$estimate )
  mu = params[1] 
  lambda = params[2] # estimate of lambda = sqrt( (nu-2)/nu ) standard_deviation
  nu = params[3]

  alpha = 0.005
  VaR = VaR_t(alpha,mu,lambda,nu,100000)
  ES = ES_t(alpha,mu,lambda,nu,100000)
  
  c( mu, lambda, nu, VaR, ES )
}
boot_results = boot(rets,boot_fn,R=1000)


VaR_boots = boot_results$t[,4]
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter19/chap_19_exer_4_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)


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

shapiro.test(ES_boots)

# Compute the two required confidence intervals:
#
boot.ci(boot_results, conf=0.95, type="norm", index=4)

boot.ci(boot_results, conf=0.95, type="norm", index=5)


# E 6 EPage 529
#
alpha = 0.1

#dat = read.csv("../../BookCode/Data/Stock_FX_Bond.csv",header=T) # first edition 
dat = read.csv("../../BookCode/Data/Stock_Bond.csv",header=T) # second edition (but the data seems to be the same as in the previous csv file)
prices = as.matrix(dat[1:500,c(3,5,7,9,11)])

# a:
n_prices = dim(prices)[1]
rets = ( prices[2:n_prices,] - prices[1:(n_prices-1),] ) / prices[1:(n_prices-1),]

stock_means = apply( rets, 2, mean )

stock_covs = cov( rets )

# b:
P = 50*106
n_stocks = dim(rets)[2]

(P/n_stocks)/prices[500,]

# c:
#
# Compute the mean return and standard deviation of this portfolio 
portfolio_mean = sum( (1/5) * stock_means )
w = rep( 1/5, n_stocks ) 
portfolio_var = t(w) %*% stock_covs %*% w

VaR_norm(0.1,portfolio_mean,sqrt(portfolio_var),S=100000)

# d (for looking five days ahead)
#
portfolio_mean = 5 * sum( (1/5) * stock_means )
portfolio_var = 5 * sum( (1/5)^2 * diag(stock_covs) ) # assuming uncorrelated returns

VaR_norm(0.1,portfolio_mean,sqrt(portfolio_var),S=100000)