#######################################################
##### Program: extreme_returns_tail_modelling_SP500_index.txt #####
#######################################################

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

# Version: 06/11/2023

# Objective of this file: To model the extreme behavior of the S&P 500 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

# S&P 500 index

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

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

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

SP500 <- SP500[-1,]


p_sp500_price <- ggplot(data = SP500, aes(x = Date, y = SPX.Adjusted, group = 1)) 
p_sp500_price + geom_line()

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

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

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

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

thres_SP500_up <- quantile(SP500$return,0.975)
thres_SP500_down <- quantile(SP500$return,0.025)

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


## Generalized Pareto distribution 

fit_gpd_SP500_up <- fevd(as.vector(SP500$return),method="MLE",
                         type="GP",threshold=thres_SP500_up)
plot(fit_gpd_SP500_up)
summary(fit_gpd_SP500_up,silent=FALSE)

fit_gpd_SP500_down <- fevd(as.vector(-SP500$return),method="MLE",
                         type="GP",threshold=-thres_SP500_down)
plot(fit_gpd_SP500_down)
summary(fit_gpd_SP500_down,silent=FALSE)