#######################################################
##### Program extreme_returns_tail_modelling.txt #####
#######################################################

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

# Version 23112023

# Objective of this file To model the extreme behavior of the FTSE 100 index based on the Extreme Value Theory 

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


# STEP 1 Preparation

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) # R package to downlaod data from Yahoo! Finance
library(ggplot2)
library(dplyr)
library(extRemes) # R package to estimate extreme value statistics
library(moments)

# STEP 2 Data extraction and preliminary analysis

# FTSE 100 index for the UK market
DATA - getSymbols(^FTSE, from = 2015-04-01, to   = 2023-04-01, env = NULL) 

# Nikkei 225 index for the Japanese market
#DATA - getSymbols(^N225, from = 2015-04-01, to   = 2023-04-01, env = NULL) 

# S&P 500 index for the US market
#DATA - getSymbols(^GSPC, from = 2015-04-01, to   = 2023-04-01, env = NULL) 

# EURO STOXX 50 index for the Japanese market
#DATA - getSymbols(^STOXX50E, from = 2015-04-01, to   = 2023-04-01, env = NULL) 

# CSI 300 index for the Japanese market
#DATA - getSymbols(000300.SS, from = 2015-04-01, to   = 2023-04-01, env = NULL) 


#remove NA values (if there be)
DATA - DATA[!is.na(DATA[,6]),] 

# Ashuted index 
DATA$AjustedIndex = DATA[,6]

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

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

DATA - DATA[-1,]


p_DATA_price - ggplot(data = DATA, aes(x = Date, y = AjustedIndex, group = 1)) 
p_DATA_price + geom_line()

p_DATA_rtn - ggplot(data = DATA, aes(x = Date, y = return, group = 1))
p_DATA_rtn + geom_line()

# Descriptive statistics of returns
mean(DATA$return)
median(DATA$return)
sd(DATA$return)
skewness(DATA$return)
kurtosis(DATA$return)

# Top 10 negative returns
DATA %% slice_min(return, n = 10)

# Top 10 positive returns
DATA %% slice_max(return, n = 10)


# STEP 3 Modelling of the tails with the peak-over-threshold (POT) method

# Definition of the threshold to select negative tail returns
threshold_down - quantile(DATA$return,0.025)

# Definition of the threshold to select positive tail returns
threshold_up   - quantile(DATA$return,0.975)

# Estimation of the Generalized Pareto distribution(GPD) for the left tail 
fit_gpd_left_tail - fevd(as.vector(-DATA$return),method=MLE,
                          type=GP,threshold=-threshold_down)

plot(fit_gpd_left_tail, type = c(hist), xlim = c(0, 10), 
    main=Histogram of empirical exceedings and modeled GPD density n for the left tail)
title(xlab=Negative return exceedances, line=4)

summary(fit_gpd_left_tail,silent=FALSE)

# Estimation of the Generalized Pareto distribution(GPD) for the right tail 
fit_gpd_right_tail - fevd(as.vector(DATA$return),method=MLE,
                           type=GP,threshold=threshold_up)

plot(fit_gpd_right_tail, type = c(hist), xlim = c(0, 10),
     main=Histogram of empirical exceedings and modeled GPD density n for the right tail)
title(xlab=Positive return exceedances, line=4)

summary(fit_gpd_right_tail,silent=FALSE)