# This script tries to calculate the measures of (a) extreme tail independency in Table 2 (page  594) and (b) extreme tail dependency in Table 3 (page 597)
# of Poon, Rockinger, Tawn(2004), "Extreme Value Dependence in Financial Markets:Diagnostics, Models and Financial Implications", The Review of Financial Studies, Vol 17, No2, pp 481-610
# between DJIA and NASDAQ's left and right tails from 1/2/1990 to 12/31/2002 for both filtered and unfiltered returns.

# This code runs in R 2.10.1

# To filter the returns, TGARCH(1,1)with threshold =1 was used for the asymmetric version of GARCH(1,1) in page 606

# Expected results are:
# a. measures for the (a) tail independency must range between [0,1] and (2) tail dependency must range between (-1,1] - see pages 586 to 587 respectively.
# b. the measures for the unfiltered returns must be GREATER than for the filtered returns

# The journal is refered to in the script as PRT(2004)


library(QRMlib)
library(fBasics)
library(POT)

#GRAPHS OF PRICES AND RETURNS

prices <- read.table("yPrices.txt",header=TRUE) 	# daily closing prices
prices1 <- ConvertDFToTimeSeries(prices)

returns <- read.table("yReturns.txt",header=TRUE) 	# unfiltered daily returns
returns1 <- ConvertDFToTimeSeries(returns)

DJIA_r<-returns[,"DJIA_r"]				# DJIA's unfiltered returns. Used for the right tail.
NASDAQ_r=returns[,"NASDAQ_r"]				# NASDAQ's filtered returns. Used for the right tail.


DJIA_r_l=DJIA_r*-1 					# DJIA's negative unfiltered returns are positive. Used for the left tails.
NASDAQ_r_l=NASDAQ_r*-1					# NASDAQ's negative unfiltered returns are positive. Used for the left tails.


p=0.05							# top 5% of the total observations is considered extremes
N=length(DJIA_r_l)					# total observations
kval=N*p							# no of exceedances


descStats<-basicStats(returns1)			# descriptive stats for both variables
descStats

residuals <- read.table("yResiduals.txt",header=TRUE) # filtered daily returns
residuals1 <- ConvertDFToTimeSeries(residuals)

DJIA_rs<-residuals[,"DJIA_rs"] 			# DJIA's filtered returns from TGARCH model. Used for the righ tails.
NASDAQ_rs<-residuals[,"NASDAQ_rs"]			# NASDAQ's filtered returns from TGARCH model. Used for the right tails.


DJIA_rs_l=DJIA_rs*-1					# DJIA's negative filtered returns are positive. Used for the left tails.
NASDAQ_rs_l=NASDAQ_rs*-1				# NASDAQ's negative filtered returns are positive. Used for the left tails. 


#BIVARIATE :  A PORTFOLIO containing DJIA and NASDAQ: 
# Test for asymptotic dependence

#------------
# FILTERED LEFT TAIL (F.L)
#------------

mod1<- fit.GPD(DJIA_rs_l,nextremes=kval)		
mod2<- fit.GPD(NASDAQ_rs_l,nextremes=kval)


u1=mod1$threshold 
nu1=mod1$n.exceed 

u2=mod2$threshold
nu2=mod2$n.exceed

#Implements eq(1) on page 584 using eq(12) in page 590 which estimates GPD and uses the empirical distribution function

FnS <-ecdf(DJIA_rs_l)
FS<- ifelse(DJIA_rs_l<u1,FnS(DJIA_rs_l),1-(1-FnS(DJIA_rs_l))*(1-pgpd(DJIA_rs_l,loc=u1,shape=mod1$par.ests[1],scale=mod1$par.ests[2])))


FnE <-ecdf(NASDAQ_rs_l)
FE<- ifelse(NASDAQ_rs_l<u2,FnE(NASDAQ_rs_l),1-(1-FnE(NASDAQ_rs_l))*(1-pgpd(NASDAQ_rs_l,loc=u2,shape=mod2$par.ests[1],scale=mod2$par.ests[2])))

S <-(-1)/log(FS)
U <-(-1)/log(FE)

# Let ZA=min(S,U) 

# Implements eq(10) on page 589

ZA<-pmin(S,U)
n=length(ZA)

# Implements the eq following eq(10) on page 589
# Draphically diagnose whether the threshold of 5% exceedances for the estimated GPD is good

modZ <- fit.GPD(ZA,nextremes=kval)  	
uZF.L=modZ$threshold
nuZF.L=modZ$n.exceed

par(mfrow = c(2, 2))
xiplot(ZA, ci=0.95,auto.scale=TRUE, labels=TRUE, reverse=FALSE);
abline(v=modZ$n.exceed, lty=4); 
MEplot(ZA)
abline(v=modZ$n.exceed, lty=4); 
par(mfrow = c(1, 1))
dev.copy(png,'threshold selection plots.png')
dev.off()


Zu<-ZA[ZA>uZF.L]

# Calculate the estimates of the two TDCs in page 589

chi_barF.L<-(2*(sum(log(Zu/uZF.L))/nuZF.L)-1)
varchi_barF.L=((chi_barF.L+1)^2)/nuZF.L
se_chi_barF.L=sqrt(varchi_barF.L) 
chi_bar_calcF.L=chi_barF.L+1.96*sqrt(varchi_barF.L)

# If there is no enough evidence to reject Ho, estimate chi:

chiF.L=(uZF.L*nuZF.L)/n
varchiF.L=(uZF.L^2*nuZF.L*(n-nuZF.L))/n^3
se_varchiF.L=sqrt(varchiF.L)

# PRT(2004) results are in Table 2 on page 594

#RESULTS OF FILTERED, LEFT TAIL
uZF.L
chi_barF.L
se_chi_barF.L
chi_bar_calcF.L
chiF.L
se_varchiF.L

#------------
# UNFILTERED LEFT TAIL (U.L)
#------------

mod5 <- fit.GPD(DJIA_r_l,nextremes=kval)
mod6 <- fit.GPD(NASDAQ_r_l,nextremes=kval)

u5=mod5$threshold
nu5=mod5$n.exceed

u6=mod6$threshold
nu6=mod6$n.exceed

FnS <-ecdf(DJIA_r_l)
FS<- ifelse(DJIA_r_l<u5,FnS(DJIA_r_l),1-(1-FnS(DJIA_r_l))*(1-pgpd(DJIA_r_l,loc=u5,shape=mod5$par.ests[1],scale=mod5$par.ests[2])))

FnE <-ecdf(NASDAQ_r_l)
FE<- ifelse(NASDAQ_r_l<u6,FnE(NASDAQ_r_l),1-(1-FnE(NASDAQ_r_l))*(1-pgpd(NASDAQ_r_l,loc=u6,shape=mod6$par.ests[1],scale=mod6$par.ests[2])))


S <-(-1)/log(FS)
U <-(-1)/log(FE)

ZC<-pmin(S,U)
n=length(ZC)

modZC <- fit.GPD(ZC,nextremes=kval)
mleZC <- fitgpd(ZC, thresh = modZC$threshold, est = "mle")

uZU.L=modZC$threshold
nuZU.L=modZC$n.exceed

Zu<-ZC[ZC>uZU.L]

chi_barU.L<-(2*(sum(log(Zu/uZU.L))/nuZU.L)-1)
varchi_barU.L=((chi_barU.L+1)^2)/nuZU.L
se_chi_barU.L=sqrt(varchi_barU.L) 
chi_bar_calcU.L=chi_barU.L+1.96*sqrt(varchi_barU.L)

chiU.L=(uZU.L*
nuZU.L)/n
varchiU.L=(uZU.L^2*nuZU.L*(n-nuZU.L))/n^3
se_varchiU.L=sqrt(varchiU.L)

#RESULTS OF UNFILTERED, LEFT TAIL
uZU.L
chi_barU.L
se_chi_barU.L
chi_bar_calcU.L
chiU.L
se_varchiU.L


# Compare results for unfiltered and filtered returns for left tails
chi_bar_calcF.L; chi_bar_calcU.L;
chi_barF.L; chi_barU.L;


#-------
# FILTERRED + RIGHT TAIL  (F.R)
#-------

mod3 <- fit.GPD(DJIA_rs,nextremes=kval)
mod4 <- fit.GPD(NASDAQ_rs,nextremes=kval)


u3=mod3$threshold
nu3=mod3$n.exceed

u4=mod4$threshold
nu4=mod4$n.exceed

FnS <-ecdf(DJIA_rs)
FS<- ifelse(DJIA_rs<u3,FnS(DJIA_rs),1-(1-FnS(DJIA_rs))*(1-pgpd(DJIA_rs,loc=u3,shape=mod3$par.ests[1],scale=mod3$par.ests[2])))

FnT <-ecdf(NASDAQ_rs)
FE<- ifelse(NASDAQ_rs<u4,FnE(NASDAQ_rs),1-(1-FnE(NASDAQ_rs))*(1-pgpd(NASDAQ_rs,loc=u4,shape=mod4$par.ests[1],scale=mod4$par.ests[2])))


S <-(-1)/log(FS)
U <-(-1)/log(FE)

ZB<-pmin(S,U)
n=length(ZB)

modZB <- fit.GPD(ZB,nextremes=kval)
mleZB <- fitgpd(ZB, thresh = modZ$threshold, est = "mle")

uZF.R=modZB$threshold
nuZF.R=modZB$n.exceed

Zu<-ZB[ZB>uZF.R]

chi_barF.R<-(2*(sum(log(Zu/uZF.R))/nuZF.R)-1)
varchi_barF.R=((chi_barF.R+1)^2)/nuZF.R
se_chi_barF.R=sqrt(varchi_barF.R) 
chi_bar_calcF.R=chi_barF.R+1.96*sqrt(varchi_barF.R)

chiF.R=(uZF.R*nuZF.R)/n
varchiF.R=(uZF.R^2*nuZF.R*(n-nuZF.R))/n^3
se_varchiF.R=sqrt(varchiF.R)

#RESULTS OF FILTERED, RIGHT TAIL
uZF.R
chi_barF.R
se_chi_barF.R
chi_bar_calcF.R
chiF.R
se_varchiF.R

#-------
# UNFILTERRED + RIGHT TAIL (U.R)
#-------
mod7 <- fit.GPD(DJIA_r,nextremes=kval)
mod8 <- fit.GPD(NASDAQ_r,nextremes=kval)

u7=mod7$threshold
nu7=mod7$n.exceed

u8=mod8$threshold
nu8=mod8$n.exceed

FnS <-ecdf(DJIA_r)
FS<- ifelse(DJIA_r<u7,FnS(DJIA_r),1-(1-FnS(DJIA_r))*(1-pgpd(DJIA_r,loc=u7,shape=mod7$par.ests[1],scale=mod7$par.ests[2])))

FnT <-ecdf(NASDAQ_r)
FE<- ifelse(NASDAQ_r<u8,FnE(NASDAQ_r),1-(1-FnE(NASDAQ_r))*(1-pgpd(NASDAQ_r,loc=u8,shape=mod8$par.ests[1],scale=mod8$par.ests[2])))


S <-(-1)/log(FS)
U <-(-1)/log(FE)

ZD<-pmin(S,U)
n=length(ZD)

modZD <- fit.GPD(ZD,nextremes=kval)
mleZD <- fitgpd(ZD, thresh = modZ$threshold, est = "mle")

uZU.R=modZD$threshold
nuZU.R=modZD$n.exceed

Zu<-ZD[ZD>uZU.R]

chi_barU.R<-(2*(sum(log(Zu/uZU.R))/nuZU.R)-1)
varchi_barU.R=((chi_barU.R+1)^2)/nuZU.R
se_chi_barU.R=sqrt(varchi_barU.R) 
chi_bar_calcU.R=chi_barU.R+1.96*sqrt(varchi_barU.R)

chiU.R=(uZU.R*nuZU.R)/n
varchiU.R=(uZU.R^2*nuZU.R*(n-nuZU.R))/n^3
se_varchiU.R=sqrt(varchiU.R)

#RESULTS OF UNFILTERED, RIGHT TAIL
uZU.R
chi_barU.R
se_chi_barU.R
chi_bar_calcU.R
chiU.R
se_varchiU.R

# Compare results for unfiltered and filtered returns for right tails
chi_bar_calcF.L; chi_bar_calcU.L;
chi_bar_calcF.R;chi_bar_calcU.R