#######################################################
##### Program: extreme_returns_tail_modelling_SX5E_index.txt #####
#######################################################

# Author: Shengyu ZHENG (ESSEC Business School, Grande Ecole Program, Master in Management, 2020-2024)

# Version: 08/11/2023

# Objective of this file: To model the extreme behavior of the EURO STOXX 50 index based on the Extreme Value Theory 

# Structure of this file: 
#   Step 1:  Prepration
#   Step 2:  Data extraction and preliminary analysis
#   Step 3:  Modelling of the tails


# STEP 1: Prepration

if (!require(quantmod)) install.packages('quantmod')
if (!require(ggplot2)) install.packages('ggplot2')
if (!require(dplyr)) install.packages('dplyr')
if (!require(extRemes)) install.packages('extRemes')
if (!require(moments)) install.packages('moments')

library(quantmod)
library(ggplot2)
library(dplyr)
library(extRemes)
library(moments)

# STEP 2: Data extraction and preliminary analysis

# EURO STOXX 50 index

SX5E <- getSymbols("^STOXX50E", from = "2015-04-01", to   = "2023-04-01", env = NULL) 
SX5E <- SX5E[!is.na(SX5E[,6]),] #remove NA values (if there be)

# Compute logarithmic percentage returns (with respect to the daily closing prices)
SX5E$return = diff(log(SX5E[,6]))*100

# Convert the xts object (from the QuantMod R package) to a data frame and process dates 
SX5E <- data.frame(Date=index(SX5E), coredata(SX5E))
SX5E$Date <- as.POSIXct(SX5E$Date)

SX5E <- SX5E[-1,]


p_SX5E_price <- ggplot(data = SX5E, aes(x = Date, y = STOXX50E.Adjusted, group = 1)) 
p_SX5E_price + geom_line()

p_SX5E_rtn <- ggplot(data = SX5E, aes(x = Date, y = return, group = 1))
p_SX5E_rtn + geom_line()

# Descriptive statistics
mean(SX5E$return)
median(SX5E$return)
sd(SX5E$return)
skewness(SX5E$return)
kurtosis(SX5E$return)

#top 10 negative and positive returns
SX5E %>% slice_min(return, n = 10)
SX5E %>% slice_max(return, n = 10)

#STEP 3: Modelling of the tails
## Peak over threshold 

thres_SX5E_up   <- quantile(SX5E$return,0.975)
thres_SX5E_down <- quantile(SX5E$return,0.025)

nb_exceedances  <- SX5E%>%
  filter(return>=thres_SX5E_up)%>%
  count()
thres[1]
1-nb_exceedances/length(SX5E$return)


## Generalized Pareto distribution 

fit_gpd_SX5E_up <- fevd(as.vector(SX5E$return),method="MLE",
                         type="GP",threshold=thres_SX5E_up)
plot(fit_gpd_SX5E_up)
summary(fit_gpd_SX5E_up,silent=FALSE)

fit_gpd_SX5E_down <- fevd(as.vector(-SX5E$return),method="MLE",
                           type="GP",threshold=-thres_SX5E_down)
plot(fit_gpd_SX5E_down)
summary(fit_gpd_SX5E_down,silent=FALSE)