#############################################################################
#############################################################################

# Application of EVT - BM GEV & POT GPD to analyse the negative tail risk of
# the S&P 500 Index.

# Shengyu ZHENG

# ESSEC Business School 

# September 17, 2022

#############################################################################
#############################################################################


#1. Packages and libraries 

# Installs pacman ("package manager") 
if (!require("pacman")) install.packages("pacman")

# Use pacman to load add-on packages as desired
pacman::p_load(pacman, dplyr, GGally, ggplot2, ggthemes, 
               ggvis, httr, plotly, rio, shiny, stringr, 
               tidyr,Hmisc,naniar,moments, data.table, fixest, 
               magrittr,evd,quantmod,data.table,extRemes) 


#2. Data loading, processing and preliminary inspection 


#2.1. Loading SPX and calculating log-return 

## loading the SPX Index (daily from 1970/01/01 to 2022/08/31) 
getSymbols("^GSPC",from="1970-01-01",to="2022-08-31")


#2.2. Data processing 

## calculating log-returns (with respect to the daily closing price)
GSPC$SPX.rtn = diff(log(GSPC$GSPC.Adjusted))*100
GSPC$SPX.rtn[1]=0

## converting the xts object to a dataframe and processing dates 
GSPC$Date1 = rownames(GSPC)
GSPC <- data.frame(date=index(GSPC), coredata(GSPC))
GSPC$date <- as.POSIXct(GSPC$date)
GSPC$mo <- strftime(GSPC$date,"%m")
GSPC$yr <- strftime(GSPC$date,"%Y")


#2.3. Data preliminary inspection 

# Descriptive statistics 

describe(GSPC)
summary(GSPC)
sd(GSPC$SPX.rtn)
kurtosis(GSPC$SPX.rtn)
skewness(GSPC$SPX.rtn)


#3. Block Maxima - Generalized Extreme Value (BM-GEV) distribution 

#3.1. Block maxima (monthly minima in this case)

extrema <- setDT(GSPC)[, .(SPX_down = -min(SPX.rtn)), by = .(yr, mo)]


#3.2. Descriptive statistics of the monthly minima 

describe(extrema)
summary(extrema)
sd(extrema$SPX_down)
kurtosis(extrema$SPX_down)
skewness(extrema$SPX_down)


#3.3. GEV fitting 

fit_gev_SPX_down<-fevd(as.vector(extrema$SPX_down), method = "MLE", type="GEV")
plot(fit_gev_SPX_down)
summary(fit_gev_SPX_down, silent=FALSE) #to get the standard error of estimates




#4.Peak Over Threshold - Generalized Pareto distribution (POT-GPD)

#4.1. Peak over threshold 

thres <- 0.01*100
nb_exceedances <-GSPC%>%
                  filter(SPX.rtn>=thres)%>%
                  count()
1-nb_exceedances/length(GSPC$SPX.rtn)


#4.2. Generalized Pareto distribution 

fit_gpd_SPX_down <- fevd(as.vector(-GSPC$SPX.rtn),method="MLE",
                                    type="GP",threshold=thres)
plot(fit_gpd_SPX_down)
summary(fit_gpd_SPX_down,silent=FALSE)