Homepage

Econometrics - USD/EUR exchange rate, an EVT approach (R walkthrough)

29/12/2020

TAGS: econometrics, evt, garch, r

Purpose

"Modeling Exchange Market Volatility Risk in Rwanda Using GARCH-EVT Approach" from Rafiki Murenzi, Kigabo Thomas and Joseph K. Mung’atu propose a case study based on McNeil & Frey EVT approach for time series. Here we’ll study the USD/EUR exchange rate and develop the programming process using R that allows to obtain comparative results on different time series. The result analysis will not be expanded in this post.

Sources for an up-to-date dataset are available HERE.

Dataset

Data must be loaded and pre-processed in following fashion (where data.csv is the downloaded time series dataset):

d <- read.csv("data.csv", header = TRUE, skip = 4)
# rename columns
colnames(d) <- c("date", "rate")
# switch data to time series
d$date <- as.Date(d$date, format = "%Y-%m-%d")
d$rate <- as.double(d$rate)

econometrics_us_eur_exchange_rate

We can transform the dataset to log-returns using the following chunk of code:

# transform data into log-returns
d$rate_logr = c(diff(log(d$rate)), NA)
# trick: copy dataframe so we can keep date
dlogr <- d
# trick: remove useless ’rate’ column
dlogr$rate <- NULL
# there should be only one NA, remove it
which(is.na(dlogr))
dlogr <- na.omit(dlogr)

This transformation looks that way:

econometrics_us_eur_log_returns

Stationarity

The original time series is not stationary from the Augmented Dickey Fuller test. Stationarity can be achieved on log-returns:

# check for stationarity: null hypothesis should be
# rejected now (log-returns series is stationary)
adf.test(dlogr$rate_logr)

GARCH models

Using the rugarch R package, several GARCH variants (APARCH, EGARCH, GJR-GARCH and standard GARCH) can be modelled with the following chunk of code:

# ARMA(1,1)-GARCH(1,1)
spec_garch <- ugarchspec(variance.model = list(model = "sGARCH",
    garchOrder = c(1, 1), submodel = NULL, external.regressors = NULL,
    variance.targeting = FALSE), mean.model = list(armaOrder = c(1,1),
	external.regressors = NULL), distribution.model = "std")

fit_garch <- ugarchfit(spec = spec_garch, data = dlogr$rate_logr,
    solver.control = list(trace = 0), solver = "hybrid")

# ARMA (1,1)-GJRGARCH(1,1)
spec_gjrgarch <- ugarchspec(variance.model = list(model = "gjrGARCH",
    garchOrder = c(1, 1), submodel = NULL, external.regressors = NULL,
    variance.targeting = FALSE), mean.model = list(armaOrder = c(1,1),
	external.regressors = NULL), distribution.model = "std")

fit_gjrgarch <- ugarchfit(spec = spec_gjrgarch, data = dlogr$rate_logr,
    solver.control = list(trace = 0), solver = "hybrid")

# ARMA (1,1)-APARCH(1,1)
spec_aparch <- ugarchspec(variance.model = list(model = "apARCH",
    garchOrder = c(1, 1), submodel = NULL, external.regressors = NULL,
    variance.targeting = FALSE), mean.model = list(armaOrder = c(1,1),
	external.regressors = NULL), distribution.model = "std")

fit_aparch <- ugarchfit(spec = spec_aparch, data = dlogr$rate_logr,
    solver.control = list(trace = 0), solver = "hybrid")

# ARMA (1,1)-EGARCH(1,1)
spec_egarch <- ugarchspec(variance.model = list(model = "eGARCH",
    garchOrder = c(1, 1), submodel = NULL, external.regressors = NULL,
    variance.targeting = FALSE), mean.model = list(armaOrder = c(1,1),
	external.regressors = NULL), distribution.model = "std")

fit_egarch <- ugarchfit(spec = spec_egarch, data = dlogr$rate_logr,
    solver.control = list(trace = 0), solver = "hybrid")

Analyzing Akaike, Bayes, Shibata, Hannan-Quinn leads to choose APARCH as the best model for Auto Regressive Conditional Hesteroskesticity.

| | Akaike | Bayes | Shibata | Hannan-Quinn | |--------------------------|---------|---------|---------|--------------| | ARMA(1,1)-GARCH(1,1) | -7.5171 | -7.5087 | -7.5171 | -7.5142 | | ARMA (1,1)-GJRGARCH(1,1) | -7.5170 | -7.5073 | -7.5170 | -7.5136 | | ARMA(1,1)-APARCH(1,1) | -7.5014 | -7.4905 | -7.5014 | -7.4976 | | ARMA(1,1)-EGARCH(1,1) | -7.5206 | -7.5110 | -7.5206 | -7.5173 |

This APARCH model’s residuals can now be analyzed as innovations for an EVT approach. The can be displayed as such:

dlogr$resid_aparch <- residuals(fit_aparch, standardize=TRUE)

# plot residuals
ggplot(dlogr, aes(x = date, y = resid_aparch)) +
  geom_line() +
  labs(x = "Temps", y = "Résidus", title = "Résidus du modèle APARCH(1,1)") +
  theme_minimal()

econometrics_aparch_residuals

ACF and PACF plots can be displayed in following fashion:

p1 <- ggAcf(dlogr$resid_aparch) +
  labs(x = "Lag", y = "ACF", title = "ACF des résidus") +
  theme_minimal()

p2 <- ggPacf(dlogr$resid_aparch) +
  labs(x = "Lag", y = "PACF", title = "PACF des résidus") +
  theme_minimal()

grid.arrange(p1, p2, ncol=2)

econometrics_aparch_acf_pacf

A Student distribution can be applied to APARCH residuals. We display next the QQ plot that shows a good fit for conditional distribution of residuals.

econometrics_aparch_student_qq_plot

A Normal distribution would be inadequate, as the next QQ plot testifies:

econometrics_aparch_normal_qq_plot

We can now move to EVT modelling (GPD, GEV).

Generalized Pareto Distribution

Mean Residual Plot can be displayed the following way:

resid_mrl <- mrl(dlogr$resid_aparch)
resid_grf <- gpdRangeFit(dlogr$resid_aparch)
p1 <- ggplot(resid_grf)
p2 <- ggplot(resid_mrl)

grid.arrange(p1[[1]] + ggtitle("Stability plot, scale parameter"),
             p1[[2]] + ggtitle("Stability plot, shape parameter"),
             p2 + ggtitle("MRL plot"),
             ncol=2)

econometrics_mrl_plot

If we define an arbitrary threshold of 2.0, we can display the selected points above it.

threshold <- 2.0
ggplot(data = dlogr, aes(x = date, y = resid_aparch)) +
  geom_point(color = ifelse(dlogr$resid_aparch > threshold, "red", "black"), size=0.5) +
  geom_hline(yintercept=threshold, color = "red") +
  scale_color_identity() +
  labs(title = paste0("Peak Over Thresold at ", threshold, sep = " ")) +
  theme_minimal()

econometrics_peak_over_threshold

A threshold range can be supplied for GPD adjustment.

gpd.fitrange(dlogr$resid_aparch, 1.5, 3.5, nint = 50, show = FALSE)

And the fitted model can be adjusted for a given threshold with the following command.

fpot(dlogr$resid_aparch, threshold=2.0)

The computed parameters below show that the attraction domain covers a Frechet and a Weibull. We may need to do a deviance test in order to evaluate the correct fit.

| sigma | xi | |-------------------------|---------------------------| | 0.58015 | 0.06596 | | [0.452946 ; 0.707354] | [-0.0813144 ; 0.2132344] |

The GDP densities for Gumbel and Weibull can be displayed as such:

econometrics_gev_attraction_domains

The deviance test implies a call to the anova function, as such:

fre_gpd <- fpot(dlogr$resid_aparch, threshold=2.0)
gum_gpd <- fpot(dlogr$resid_aparch, threshold=2.0, shape=0)
anova(fre_gpd, gum_gpd)

Generalized Extreme Value distribution

Similarly to GPD, a GEV can be fitted using the following call:

gev_fit <- fgev(dlogr$resid_aparch)
gev_fit

We should obtain the same shape result as the GDP. If it not the case we might consider evaluating the quality of fit visually. Which can be achieved with this simple call:

par(mfrow=c(2,2))
plot(gev_fit)

econometrics_aparch_gev_adjustment

Risk metrics

If we choose Expected Shortfall over VaR as a risk measure, the quantile can be computed with the following command:

esgumbel(c(0.9, 0.95, 0.975, 0.99, 0.995, 0.999), mu=2.0, sigma=0.58015)
esdweibull(c(0.9, 0.95, 0.975, 0.99, 0.995, 0.999), mu=-0.423, c=-0.1707, sigma=1.0836)

References

  1. Modeling Exchange Market Volatility Risk in Rwanda Using GARCH-EVT Approach - Rafiki Murenzi, Kigabo Thomas, Joseph K. Mung’atu

  2. Inférence statistique des modèles conditionnellement hétéroscédastiques avec innovations stables, contraste non gaussien et volatilité mal spécifiée - Guillaume Lepage

  3. Les tests de non-stationnarité (racine unité) - Arthur Charpentier

  4. Finance Empirique - Ecole des HEC, Université de Lausanne - Eric Jondeau

  5. Introduction to the rugarch package - Alexios Ghalanos

  6. A comparison of volatility models: Does anything beat a GARCH(1,1) ? - P. Reinhard Hansen and A. Lunde

  7. Why log returns? - Cathy O’Neil

  8. La volatilité boursière : des constats empiriques aux difficultés d’interprétation - Marie-Hélène Grouard, Sébastien Lévy, Catherine Lubochinsky

  9. Asymmetric GARCH type models for asymmetric volatility characteristics analysis and wind power forecasting

  10. Forecasting Value-at-Risk using GARCH and Extreme-Value-Theory Approaches for Daily Returns

  11. Estimation of tail-related risk measures for heteroscedastic financial time series: an extreme value approach - Alexander J. McNeil, Rüdiger Frey

  12. Evaluating GARCH models - Stefan Lundbergh & Timo Teräsvirta

Programming references

  1. ACF Plot with ggplot2: Setting width of geom_bar

  2. How to fit ARMA+GARCH Model In R?

  3. Plotting Time Series Statistics with ggplot2 and ggfortify

  4. Basic Time-Series Analysis: Modeling Volatility (GARCH)

  5. How to Model Volatility with ARCH and GARCH for Time Series Forecasting in Python

  6. Value at Risk estimation using GARCH model

  7. VaR with GARCH(1,1)

  8. Time Series Analysis of Apple Stock Prices Using GARCH models

  9. ARCH-GARCH Tutorial with rugarch package

  10. Toulouse School of Economics M2 Statistics and Econometrics, Univariate extreme value theory: Lab session

  11. Univariate Extreme Value Modelling using R

  12. Toulouse School of Economics, M2 Statistics and Econometrics - Univariate extreme value theory: Lab session

Appendix: full sources

require(knitr)
require(ggplot2)
require(tseries)
require(forecast)
require(rugarch)
require(WeightedPortTest)
require(fGarch)
require(evd)
require(texmex)
require(gridExtra)
require(extRemes)
require(xts)
require(ismev)
require(fExtremes)
require(VaRES)
library(reshape2)

# https://sdw.ecb.europa.eu/browseSelection.do?org.apache.struts.taglib.html.TOKEN=dfb25b6a53a7f0373e21297d35c17527&df=true&ec=&dc=&oc=&pb=&rc=&DATASET=0&removeItem=&removedItemList=&mergeFilter=&activeTab=EXR&showHide=&MAX_DOWNLOAD_SERIES=500&SERIES_MAX_NUM=50&node=qview&SERIES_KEY=120.EXR.D.USD.EUR.SP00.A
# load data
d <- read.csv("data.csv", header = TRUE, skip = 4)
# rename columns
colnames(d) <- c("date", "rate")
# switch data to time series
d$date <- as.Date(d$date, format = "%Y-%m-%d")
d$rate <- as.double(d$rate)
# check for NAs
which(is.na(d))
# remove NAs
d <- na.omit(d)
d <- d[order(as.Date(d$date, format="%Y-%m-%d")),]

# plot time series
p <- ggplot(d, aes(x = date, y = rate)) +
  geom_line() +
  labs(x = "Date", y = "Rate", title = "EUR/USD exchange rate 1999-2020") +
  theme_minimal()

# check for stationarity
adf.test(d$rate)

# ACP/PACF

ggAcf(d$rate)
ggAcf(dlogr$rate_logr)
ggPacf(d$rate)
ggPacf(dlogr$rate_logr)

# transform data into log-returns
d$rate_logr = c(diff(log(d$rate)), NA)
# trick: copy dataframe so we can keep date
dlogr <- d
# trick: remove useless "rate" column
dlogr$rate <- NULL
# there should be only one NA, remove it
which(is.na(dlogr))
dlogr <- na.omit(dlogr)

# plot log-return time series
plog <- ggplot(dlogr, aes(x = date, y = rate_logr)) +
  geom_line() +
  labs(x = "Date", y = "Rate log-returns", title = "EUR/USD log-returns exchange rate 1999-2020") +
  theme_minimal()

# check for stationarity: null hypothesis should be rejected now (log-returns series is stationary)
adf.test(dlogr$rate_logr)

# Auto ARIMA for p and q order solving: bad idea
mdl.auto <- auto.arima(dlogr$rate_logr)
arimaorder(mdl.auto)

# ARMA(1,1)-GARCH(1,1)
spec_garch <- ugarchspec(variance.model = list(model = "sGARCH",
                                               garchOrder = c(1, 1),
                                               submodel = NULL,
                                               external.regressors = NULL,
                                               variance.targeting = FALSE),
                         mean.model = list(armaOrder = c(1, 1),
                                           external.regressors = NULL),
                         distribution.model = ’std’)

fit_garch <- ugarchfit(spec = spec_garch, data = dlogr$rate_logr, solver.control = list(trace=0), solver="hybrid")

# ARMA (1,1)-GJRGARCH(1,1)
spec_gjrgarch <- ugarchspec(variance.model = list(model = "gjrGARCH",
                                                  garchOrder = c(1, 1),
                                                  submodel = NULL,
                                                  external.regressors = NULL,
                                                  variance.targeting = FALSE),
                            mean.model = list(armaOrder = c(1, 1),
                                              external.regressors = NULL),
                            distribution.model = ’std’)

fit_gjrgarch <- ugarchfit(spec = spec_gjrgarch, data = dlogr$rate_logr, solver.control = list(trace=0), solver="hybrid")

# ARMA (1,1)-APARCH(1,1)
spec_aparch <- ugarchspec(variance.model = list(model = "apARCH",
                                                garchOrder = c(1, 1),
                                                submodel = NULL,
                                                external.regressors = NULL,
                                                variance.targeting = FALSE),
                          mean.model = list(armaOrder = c(1, 1),
                                            external.regressors = NULL),
                          distribution.model = ’std’)

fit_aparch <- ugarchfit(spec = spec_aparch, data = dlogr$rate_logr, solver.control = list(trace=0), solver="hybrid")

# ARMA (1,1)-EGARCH(1,1)
spec_egarch <- ugarchspec(variance.model = list(model = "eGARCH",
                                                garchOrder = c(1, 1),
                                                submodel = NULL,
                                                external.regressors = NULL,
                                                variance.targeting = FALSE),
                          mean.model = list(armaOrder = c(1, 1),
                                            external.regressors = NULL),
                          distribution.model = ’std’)

fit_egarch <- ugarchfit(spec = spec_egarch, data = dlogr$rate_logr, solver.control = list(trace=0), solver="hybrid")

# check parameters and fit metrics (AIC,...)
show(fit_garch)
show(fit_gjrgarch)
show(fit_aparch)
show(fit_egarch)

# Li-Mak tests
Weighted.LM.test(fit_garch@fit$residuals, fit_garch@fit$sigma^2, lag = 2,type = c("correlation", "partial"),fitdf = 1, weighted = FALSE)
Weighted.LM.test(fit_gjrgarch@fit$residuals, fit_gjrgarch@fit$sigma^2, lag = 2,type = c("correlation", "partial"),fitdf = 1, weighted = FALSE)
Weighted.LM.test(fit_aparch@fit$residuals, fit_aparch@fit$sigma^2, lag = 2,type = c("correlation", "partial"),fitdf = 1, weighted = FALSE)
Weighted.LM.test(fit_egarch@fit$residuals, fit_egarch@fit$sigma^2, lag = 2,type = c("correlation", "partial"),fitdf = 1, weighted = FALSE)

model_garch <- garchFit(formula = rate_logr ~ garch(1, 1), data = dlogr, trace=FALSE)
model_aparch <- garchFit(formula = rate_logr ~ aparch(1, 1), data = dlogr, trace=FALSE)
summary(model_garch)
summary(model_aparch)

# plot residuals
dlogr$resid_aparch <- residuals(fit_aparch, standardize=TRUE)
ggAcf(dlogr$resid_aparch)
ggPacf(dlogr$resid_aparch)

# QQ t-student plot
ggplot(dlogr, aes(sample = rate_logr)) +
  geom_qq(distribution = qstd) +
  geom_qq_line(col = "red") +
  labs(title = "Conditional t-Student QQ plot") +
  theme_minimal()

# confirm that an assumption of conditional normality is unrealistic
ggplot(dlogr, aes(sample = rate_logr)) +
  geom_qq(distribution = qnorm) +
  geom_qq_line(col = "red") +
  labs(title = "Conditional Normality QQ plot") +
  theme_minimal()

ks.test(dlogr$resid_aparch, "qstd")

kurtosis(dlogr$resid_aparch)

# Generalized Pareto Distribution

resid_mrl <- mrl(dlogr$resid_aparch)
resid_grf <- gpdRangeFit(dlogr$resid_aparch)
p1 <- ggplot(resid_grf)
p2 <- ggplot(resid_mrl)

grid.arrange(p1[[1]] + ggtitle("Stability plot, scale parameter"),
             p1[[2]] + ggtitle("Stability plot, shape parameter"),
             p2 + ggtitle("MRL plot"),
             ncol=2)

ggplot(resid_mrl) + ggtitle("MRL plot")
mrlplot(dlogr$resid_aparch, main="Mean Residual Life Plot")

resid_gpd <- fpot(dlogr$resid_aparch, threshold=2.0, model=c("pp"))
resid_gpd <- fpot(dlogr$resid_aparch, threshold=2.0)
par(mfrow = c(2,2))
plot(resid_gpd)
fpot(dlogr$resid_aparch, 2.2)
fpot(dlogr$resid_aparch, 2.4)
fpot(dlogr$resid_aparch, 2.6)
fpot(dlogr$resid_aparch, 2.8)
fpot(dlogr$resid_aparch, 3.0)
fpot(dlogr$resid_aparch, 3.2)
fpot(dlogr$resid_aparch, 3.4)
par(mfrow=c(2,2))
plot(resid_gpd)

gpd.fitrange(dlogr$resid_aparch, 1.5, 3.5, nint = 50, show = FALSE)

evm(dlogr$resid_aparch, th=2.0, family=gpd)
evm(dlogr$resid_aparch, th=2.2, family=gpd)
evm(dlogr$resid_aparch, th=2.4, family=gpd)
evm(dlogr$resid_aparch, th=2.6, family=gpd)
evm(dlogr$resid_aparch, th=2.8, family=gpd)
evm(dlogr$resid_aparch, th=3.0, family=gpd)
evm(dlogr$resid_aparch, th=3.2, family=gpd)
evm(dlogr$resid_aparch, th=3.4, family=gpd)

threshold <- 2.0
ggplot(data = dlogr, aes(x = date, y = resid_aparch)) +
  geom_point(color = ifelse(dlogr$resid_aparch > threshold, "red", "black"), size=0.5) +
  geom_hline(yintercept=threshold, color = "red") +
  scale_color_identity() +
  labs(title = paste0("Peak Over Thresold at ", threshold, sep = " ")) +
  theme_minimal()

# we need to do a deviance test because the confidence interval
# for the modelled GPD covers a Gumbel model
fre_gpd <- fpot(dlogr$resid_aparch, threshold=2.0)
gum_gpd <- fpot(dlogr$resid_aparch, threshold=2.0, shape=0)
anova(fre_gpd, gum_gpd)


# GEV distribution
gev_fit <- fgev(dlogr$resid_aparch)
par(mfrow=c(2,2))
plot(gev_fit)

# we need to do a deviance test because the confidence interval
# for the modelled GEV covers a Gumbel model
gev_fit <- fgev(dlogr$resid_aparch)
gum_gev <- fgev(dlogr$resid_aparch, shape=0)
anova(gev_fit, gum_gev)

# density plot
Gumbel <- rgumbel(1000, 2.0, 0.58015)
Weibull <- rgev(1000, 0.423, 1.0836, -0.1707)
df <- data.frame(Gumbel, Weibull)
# Transform data from wide to long format
df <- melt(df)
# Plot density using ggplot2
ggplot(df, aes(value, color = variable)) +
  geom_density() +
  ggtitle("Weibull and Gumbel densities with estimated parameters")


# risk measures: expected shortfall
esgumbel(c(0.9, 0.95, 0.975, 0.99, 0.995, 0.999), mu=2.0, sigma=0.58015)
esdweibull(c(0.9, 0.95, 0.975, 0.99, 0.995, 0.999), mu=-0.423, c=-0.1707, sigma=1.0836)