scripts/Q4_Financial_Time_Series.R

In 3schwartz/SpecialeScrAndFun: Script and functions used for master thesis.

library(rugarch)
library(qrmtools)
#
library(zoo)
library(ADGofTest) # for ad.test()
library(moments) # for skewness(), kurtosis()
library(qrmdata)
library(qrmtools)
#### GARCH_simulation.R ####
### 1 Specify the GARCH(1,1) model ##
## GARCH(1,1) model specification
(uspec <- ugarchspec(variance.model = list(model = "sGARCH", # standard GARCH
                                           garchOrder = c(1, 1)), # GARCH(1,1)
                     mean.model = list(armaOrder = c(0, 0), # no ARMA part
                                       include.mean = FALSE), # no mean included
                     distribution.model = "norm", # normal innovations
                     fixed.pars = list(omega = 0.02, # cond. var. constant (= alpha_0)
                                       alpha1 = 0.15, # alpha_1
                                       beta1 = 0.8))) # beta_1

### 2 Generate paths from a specified model ###
## Generate two realizations of length 2000
m <- 2 # number of paths
n <- 2000 # sample size
set.seed(271) # set seed (for reproducibility)
(paths <- ugarchpath(uspec, n.sim = n, n.start = 50, # burn-in sample size
                     m.sim = m))
str(paths) # S4 object of class 'uGARCHpath'
str(paths@path) # simulated processes X_t, volatilities sigma_t and unstandardized residuals epsilon_t

## We can use getClass() to see more information about such objects
getClass("uGARCHpath")

## We can also use getSlots() to see the composition
getSlots("uGARCHpath")

### 3 Plotting uGARCHpath objects ##
### Extract the simulated series
X <- fitted(paths) # simulated process X_t = mu_t + epsilon_t for epsilon_t = sigma_t * Z_t
sig <- sigma(paths) # volatilities sigma_t (conditional standard deviations)
eps <- paths@path$residSim # unstandardized residuals epsilon_t = sigma_t * Z_t
## Note: There are no extraction methods for the unstandardized residuals epsilon_t
## Sanity checks (=> fitted() and sigma() grab out the right quantities)
stopifnot(all.equal(X,   paths@path$seriesSim, check.attributes = FALSE),
          all.equal(sig, paths@path$sigmaSim,  check.attributes = FALSE))

## Plot the simulated paths (X_t)
plot(X[,1], type = "l", ylab = expression(X[1]))
plot(X[,2], type = "l", ylab = expression(X[2]))

## Plot the corresponding volatilities (sigma_t)
plot(sig[,1], type = "h", ylab = expression(sigma[1]))
plot(sig[,2], type = "h", ylab = expression(sigma[2]))

## Plot the corresponding unstandardized residuals (epsilon_t)
plot(eps[,1], type = "l", ylab = expression(epsilon[1]))
plot(eps[,2], type = "l", ylab = expression(epsilon[2]))

## Compute the standardized residuals (Z_t)
Z <- eps/sig
plot(Z[,1], type = "p", ylab = expression(Z[1]))
plot(Z[,2], type = "p", ylab = expression(Z[2]))
## Check their N(0,1) distribution
qq_plot(Z[,1])
qq_plot(Z[,2])
## Check their joint distribution
plot(Z, xlab = expression(Z[1]), ylab = expression(Z[2]))
## Note: By doing all of this backwards, one can put different joint distributions
##       on the standard residuals. This will become clear after Chapters 6 and 7.

## There are also special plotting functions for 'uGARCHpath' objects:
plot(paths, which = 1) # plots sig[,1] (no matter how many 'paths')
plot(paths, which = 2) # plots X[,1]
plot(paths, which = 3) # plots kernel density estimate of sig[,1]
plot(paths, which = 4) # plots kernel density estimate of X[,1]
## How to see the documentation of the plot function
showMethods("plot")
getMethod("plot", signature = c(x = "uGARCHpath", y = "missing"))


### 4 Does the simulated series conform to our univariate stylized facts? ##

## We (only) consider the first path here
acf(X[,1], main = expression(X[1])) # => no serial correlation
acf(abs(X[,1]), main = expression(group("|",X[1],"|"))) # => profound serial correlation
qq_plot(X[,1], method = "empirical", main = "Normal Q-Q plot") # => heavier tailed than normal
shapiro.test(X[,1]) # Normal rejected

#### GARCH_estimation ####

### 1 S&P 500 data ##
## Load S&P 500 data
data("SP500")
plot.zoo(SP500, xlab = "Time")

## Extract the data we work with and build log-returns
SPdata <- SP500['2006-01-01/2009-12-31'] # 4 years of data
X <- returns(SPdata)
plot.zoo(X, xlab = "Time")

### 2 Fit an AR(1)--GARCH(1,1) model with normal innovations ##
## Model specification (without fixed.pars, so without keeping any parameter
## fixed during the optimization)
uspec.N <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
                      mean.model = list(armaOrder = c(1,0), # AR(1) part
                                        include.mean = TRUE), # with mean
                      distribution.model = "norm") # normal innovations
(fit.N <- ugarchfit(spec = uspec.N, data = X))
## Note: Fit to 'X', not 'SPdata'!

## The fit contains a lot of information
## - The parameter estimates are the "Optimal Parameters" (at the top)
## - plot(fit.N) will take us into a menu system (below are the most important ones)

## Series with conditional quantiles
plot(fit.N, which = 2)
layout(matrix(1:4, ncol = 2, byrow = TRUE)) # specify (2,2)-matrix of plots
plot(fit.N, which = 6) # ACF of absolute data |X_t| (shows serial correlation)
plot(fit.N, which = 9) # Q-Q plot of standardized residuals Z_t (shows leptokurtosis of Z_t; normal assumption not supported)
plot(fit.N, which = 10) # ACF of standardized residuals Z_t (shows AR dynamics do a reasonable job of explaining conditional mean)
plot(fit.N, which = 11) # ACF of squared standardized residuals Z_t^2 (shows GARCH dynamics do a reasonable job of explaining conditional sd)
layout(1) # restore layout

### 3 Fit an AR(1)--GARCH(1,1) model with Student t innovations ##
## Now consider t innovations
uspec.t <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
                      mean.model = list(armaOrder = c(1,0), include.mean = TRUE),
                      distribution.model = "std") # Student t innovations
fit.t <- ugarchfit(spec = uspec.t, data = X)

## The pictures are similar, but the Q-Q plot looks "better"
layout(matrix(1:4, ncol = 2, byrow = TRUE)) # specify (2,2)-matrix of plots
plot(fit.t, which = 6) # ACF of |X_t|
plot(fit.t, which = 9) # Q-Q plot of Z_t against a normal
plot(fit.t, which = 10) # ACF of Z_t
plot(fit.t, which = 11) # ACF of Z_t^2
layout(1)

### 4 Likelihood-ratio test for comparing the two models ###
## The normal model is "nested within" the t model (appears from the t as
## special (limiting) case). The likelihood-ratio test tests:
## H0: normal innovations are good enough
## HA: t innovations necessary
## Decision: We reject the null (in favour of the alternative) if the
##           likelihood-ratio test statistic exceeds the 0.95 quantile of a
##           chi-squared distribution with 1 degree of freedom (1 here as that's
##           the difference in the number of parameters for the two models)

## The likelihood ratio statistic is twice the difference between the log-likelihoods
LL.N <- fit.N@fit$LLH # log-likelihood of the model based on normal innovations
LL.t <- fit.t@fit$LLH # log-likelihood of the model based on t innovations
LRT <- 2*(fit.t@fit$LLH-fit.N@fit$LLH) # likelihood-ratio test statistic
LRT > qchisq(0.95, 1) # => H0 is rejected
1-pchisq(LRT, df = 1) # p-value (probability of such an extreme result if normal hypothesis were true)

### 5 Using Akaike's information criterion to compare the two models ##
## Akaike's information criterion (AIC) can be used for model comparison.
## Favour the model with smallest AIC.
## Note: AIC = 2*<number of estimated parameters> - 2 * log-likelihood
(AIC.N <- 2*length(fit.N@fit$coef) - 2*LL.N)
(AIC.t <- 2*length(fit.t@fit$coef) - 2*LL.t)
## => t model is preferred


### 6 Exploring the structure of the fitted model object ##
## Basic components
class(fit.t)
getClass("uGARCHfit")
getSlots("uGARCHfit")
fit.t.fit <- fit.t@fit
names(fit.t.fit)
fit.t.model <- fit.t@model
names(fit.t.model)
fit.t.model$modeldesc
fit.t.model$pars

## To find out about extraction methods, see ?ugarchfit
## Then follow a link to documentation for "uGARCHfit" object

## Some extraction methods
(param <- coef(fit.t)) # estimated coefficients
sig <- sigma(fit.t) # estimated volatility
VaR.99 <- quantile(fit.t, probs = 0.99) # estimated VaR at level 99%
Z <- residuals(fit.t, standardize = TRUE) # estimated standardized residuals Z_t

## Plots
plot.zoo(sig,    xlab = "Time", ylab = expression(hat(sigma)[t]))
plot.zoo(VaR.99, xlab = "Time", ylab = expression(widehat(VaR)[0.99]))
plot.zoo(Z,      xlab = "Time", ylab = expression(hat(Z)[t]))

## More residual checks
(mu.Z <- mean(Z)) # ok (should be ~= 0)
(sd.Z <- as.vector(sd(Z))) # ok (should be ~= 1)
skewness(Z) # should be 0 (is < 0 => left skewed)
hist(Z) # => left skewed
kurtosis(Z) # should be 6/(nu-4)
nu <- param[["shape"]] # estimated degrees-of-freedom
6/(nu-4) # => sample kurtosis larger than it should be
pt.hat <- function(q) pt((q-mu.Z)/sd.Z, df = nu) # estimated t distribution for Z
## Z ~ t_{hat(nu)}(hat(mu), hat(sig)^2) -> F_Z(z) = t_{hat(nu)}((z-hat(mu))/hat(sig))
ad.test(as.numeric(Z), distr.fun = pt.hat) # AD test
## => Anderson--Darling test rejects the estimated t distribution

## Possible things to try:
## 1) A GARCH model with asymmetric dynamics - GJR-GARCH
## 2) A GARCH model with asymmetric innovation distribution
## 3) Different ARMA specifications for the conditional mean

### 7 Forecasting/predicting from the fitted model ##
(forc <- ugarchforecast(fit.t, n.ahead = 10))
plot(forc, which = 1) # prediction of X_t
plot(forc, which = 3) # prediction of sigma_t

#### GARCH_EWMA ####
## Original idea by Alexios Ghalanos


### Setup ##

library(rugarch)
library(zoo)
library(qrmdata)
library(qrmtools)


### 1 S&P 500 data ##

## Load S&P 500 data
data("SP500")
plot.zoo(SP500, xlab = "Time")

## Extract the data we work with and build log-returns
SPdata <- SP500['2006-01-01/2009-12-31'] # 4 years of data
X <- returns(SPdata)


### 2 Model specifications and fitting ##

## Specify models
GARCH.spec <- ugarchspec(mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                         variance.model = list(model = "sGARCH",  garchOrder = c(1,1)),
                         distribution.model = "std")
EWMAfixed.spec <- ugarchspec(mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                             variance.model = list(model = "iGARCH"),
                             fixed.pars = list(alpha1 = 1-0.94, omega = 0))
EWMAest.spec <- ugarchspec(mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                           variance.model = list(model = "iGARCH"),
                           fixed.pars = list(omega = 0))

## Model fitting
mod1 <- ugarchfit(GARCH.spec, X)
mod2 <- ugarchfilter(EWMAfixed.spec, X) # nothing to estimate, apply filter instead
mod3 <- ugarchfit(EWMAest.spec, X)


### 3 Plotting ##

## Plot volatility estimates
plot(sigma(mod1), main = "", auto.grid = FALSE, major.ticks = "auto", minor.ticks = FALSE)
lines(sigma(mod2), col = 2, lty = 2)
lines(sigma(mod3), col = 3, lty = 3)
legend("topleft", bty = "n", col = 1:3, lty = 1:3, cex = 0.9,
       legend = as.expression(c("GARCH", substitute("EWMA[fix."*lambda == x*"]", list(x = 0.94)),
                                substitute("EWMA[est."*lambda == x*"]", list(x = round(coef(mod3)[3],2))))))

#### Simulate_fit_predict_VaR_under_ARMA_GARCH ####
## Simulate an ARMA(1,1)-GARCH(1,1) process, fit such a process to
## the simulated data, estimate and backtest VaR, predict, simulate B paths,
## and compute corresponding VaR prediction with confidence intervals.
### 1 Simulation ##

## Model specification (for simulation)
armaOrder <- c(1,1) # ARMA order
garchOrder <- c(1,1) # GARCH order
nu <- 3 # d.o.f. of the standardized distribution of Z_t
fixed.p <- list(mu = 0, # our mu (intercept)
                ar1 = 0.5, # our phi_1 (AR(1) parameter of mu_t)
                ma1 = 0.3, # our theta_1 (MA(1) parameter of mu_t)
                omega = 4, # our alpha_0 (intercept)
                alpha1 = 0.4, # our alpha_1 (ARCH(1) parameter of sigma_t^2)
                beta1 = 0.2, # our beta_1 (GARCH(1) parameter of sigma_t^2)
                shape = nu) # d.o.f. nu for standardized t_nu innovations
varModel <- list(model = "sGARCH", garchOrder = garchOrder)
spec <- ugarchspec(varModel, mean.model = list(armaOrder = armaOrder),
                   fixed.pars = fixed.p, distribution.model = "std") # t standardized residuals

## Simulate (X_t)
n <- 1000 # sample size (= length of simulated paths)
x <- ugarchpath(spec, n.sim = n, m.sim = 1, rseed = 271) # n.sim length of simulated path; m.sim = number of paths
## Note the difference:
## - ugarchpath(): simulate from a specified model
## - ugarchsim():  simulate from a fitted object

## Extract the resulting series
X <- fitted(x) # simulated process X_t = mu_t + epsilon_t for epsilon_t = sigma_t * Z_t
sig <- sigma(x) # volatilities sigma_t (conditional standard deviations)
eps <- x@path$residSim # unstandardized residuals epsilon_t = sigma_t * Z_t
## Note: There are no extraction methods for the unstandardized residuals epsilon_t

## Sanity checks (=> fitted() and sigma() grab out the right quantities)
stopifnot(all.equal(X,   x@path$seriesSim, check.attributes = FALSE),
          all.equal(sig, x@path$sigmaSim,  check.attributes = FALSE))

## Plots
plot(X,   type = "l", xlab = "t", ylab = expression(X[t]))
plot(sig, type = "h", xlab = "t", ylab = expression(sigma[t]))
plot(eps, type = "l", xlab = "t", ylab = expression(epsilon[t]))


### 2 Fitting and checks ##

## Now we do as if we don't know the data-generating model. One would then
## normally fit models of various orders and decide which one to take. This
## step is omitted here for the sake of simplicity.

## Fit an ARMA(1,1)-GARCH(1,1) model
spec <- ugarchspec(varModel, mean.model = list(armaOrder = armaOrder),
                   distribution.model = "std") # without fixed parameters here
fit <- ugarchfit(spec, data = X) # fit

## Extract the resulting series
mu. <- fitted(fit) # fitted hat{mu}_t (= hat{X}_t)
sig. <- sigma(fit) # fitted hat{sigma}_t

## Sanity checks (=> fitted() and sigma() grab out the right quantities)
stopifnot(all.equal(as.numeric(mu.),  fit@fit$fitted.values),
          all.equal(as.numeric(sig.), fit@fit$sigma))

## Plot data X_t and fitted hat{mu}_t
plot(X, type = "l", xlab = "t",
     ylab = expression("Data"~X[t]~"and fitted values"~hat(mu)[t]))
lines(as.numeric(mu.), col = adjustcolor("blue", alpha.f = 0.5))
legend("bottomright", bty = "n", lty = c(1,1),
       col = c("black", adjustcolor("blue", alpha.f = 0.5)),
       legend = c(expression(X[t]), expression(hat(mu)[t])))

## Plot the unstandardized residuals epsilon_t
resi <- as.numeric(residuals(fit))
stopifnot(all.equal(fit@fit$residuals, resi))
plot(resi, type = "l", xlab = "t", ylab = expression(epsilon[t])) # check residuals epsilon_t

## Q-Q plot of the standardized residuals Z_t against their specified t
## (t_nu with variance 1); as we generated the data, we know the exact t
## innovation distribution
Z <- fit@fit$z
stopifnot(all.equal(Z, as.numeric(resi/sig.)))
qq_plot(Z, FUN = function(p) sqrt((nu-2)/nu) * qt(p, df = nu),
        main = substitute("Q-Q plot of ("*Z[t]*") against the (here: correct) standardized"~
                            t[nu.], list(nu. = nu)))


### 3 VaR estimates, checks and backtest ##

## VaR confidence level we consider here
alpha <- 0.99

## Compute VaR_alpha based on fitted mu_t and sigma_t
VaR. <- as.numeric(quantile(fit, probs = alpha))
## Build manually and compare the two
nu. <- fit@fit$coef["shape"] # extract (fitted) d.o.f. nu
VaR.. <- as.numeric(mu. + sig. * sqrt((nu.-2)/nu.) * qt(alpha, df = nu.)) # VaR_alpha computed manually
stopifnot(all.equal(VaR.., VaR.))
## => quantile(<rugarch object>, probs = alpha) provides VaR_alpha = hat{mu}_t + hat{sigma}_t * q_Z(alpha)


### 4 Backtest VaR estimates ##

## Note: VaRTest() is written for the lower tail (not sign-adjusted losses)
##       (hence the complicated call here, requiring to refit the process to -X)
btest <- VaRTest(1-alpha, actual = -X,
                 VaR = quantile(ugarchfit(spec, data = -X), probs = 1-alpha))
btest$expected.exceed # number of expected exceedances = (1-alpha) * n
btest$actual.exceed # actual exceedances
## Unconditional test
btest$uc.H0 # corresponding null hypothesis
btest$uc.Decision # test decision
## Conditional test
btest$cc.H0 # corresponding null hypothesis
btest$cc.Decision # test decision

### 5 Predict (X_t) and VaR_alpha ##

## Predict from the fitted process
fspec <- getspec(fit) # specification of the fitted process
setfixed(fspec) <- as.list(coef(fit)) # set the parameters to the fitted ones
m <- ceiling(n / 10) # number of steps to forecast (roll/iterate m-1 times forward with frequency 1)
pred <- ugarchforecast(fspec, data = X, n.ahead = 1, n.roll = m-1, out.sample = m) # predict from the fitted process

## Extract the resulting series
mu.predict <- fitted(pred) # extract predicted X_t (= conditional mean mu_t; note: E[Z] = 0)
sig.predict <- sigma(pred) # extract predicted sigma_t
VaR.predict <- as.numeric(quantile(pred, probs = alpha)) # corresponding predicted VaR_alpha

## Checks
## Sanity checks
stopifnot(all.equal(mu.predict, pred@forecast$seriesFor, check.attributes = FALSE),
          all.equal(sig.predict, pred@forecast$sigmaFor, check.attributes = FALSE)) # sanity check
## Build predicted VaR_alpha manually and compare the two
VaR.predict. <- as.numeric(mu.predict + sig.predict * sqrt((nu.-2)/nu.) *
                             qt(alpha, df = nu.)) # VaR_alpha computed manually
stopifnot(all.equal(VaR.predict., VaR.predict))


### 6 Simulate B paths from the fitted model ##

## Simulate B paths
B <- 1000
set.seed(271)
X.sim.obj <- ugarchpath(fspec, n.sim = m, m.sim = B) # simulate future paths

## Compute simulated VaR_alpha and corresponding (simulated) confidence intervals
## Note: Each series is now an (m, B) matrix (each column is one path of length m)
X.sim <- fitted(X.sim.obj) # extract simulated X_t
sig.sim <- sigma(X.sim.obj) # extract sigma_t
eps.sim <- X.sim.obj@path$residSim # extract epsilon_t
VaR.sim <- (X.sim - eps.sim) + sig.sim * sqrt((nu.-2)/nu.) * qt(alpha, df = nu.) # (m, B) matrix
VaR.CI <- apply(VaR.sim, 1, function(x) quantile(x, probs = c(0.025, 0.975)))


### 7 Plot ##

## Setup
yran <- range(X, # simulated path
              mu., VaR., # fitted conditional mean and VaR_alpha
              mu.predict, VaR.predict, VaR.CI) # predicted mean, VaR and CIs
myran <- max(abs(yran))
yran <- c(-myran, myran) # y-range for the plot
xran <- c(1, length(X) + m) # x-range for the plot
## Simulated (original) data (X_t), fitted conditional mean mu_t and VaR_alpha
plot(X, type = "l", xlim = xran, ylim = yran, xlab = "Time t", ylab = "",
     main = "Simulated ARMA-GARCH, fit, VaR, VaR predictions and CIs")
lines(as.numeric(mu.), col = adjustcolor("darkblue", alpha.f = 0.5)) # hat{\mu}_t
lines(VaR., col = "darkred") # estimated VaR_alpha
mtext(paste0("Expected exceed.: ",btest$expected.exceed,"   ",
             "Actual exceed.: ",btest$actual.exceed,"   ",
             "Test: ", btest$cc.Decision),
      side = 4, adj = 0, line = 0.5, cex = 0.9) # label
## Predictions
t. <- length(X) + seq_len(m) # future time points
lines(t., mu.predict, col = "blue") # predicted process X_t (or mu_t)
lines(t., VaR.predict, col = "red") # predicted VaR_alpha
lines(t., VaR.CI[1,], col = "orange") # lower 95%-CI for VaR_alpha
lines(t., VaR.CI[2,], col = "orange") # upper 95%-CI for VaR_alpha
legend("bottomright", bty = "n", lty = rep(1, 6), lwd = 1.6,
       col = c("black", adjustcolor("darkblue", alpha.f = 0.5), "blue",
               "darkred", "red", "orange"),
       legend = c(expression(X[t]), expression(hat(mu)[t]),
                  expression("Predicted"~mu[t]~"(or"~X[t]*")"),
                  substitute(widehat(VaR)[a], list(a = alpha)),
                  substitute("Predicted"~VaR[a], list(a = alpha)),
                  substitute("Simulated 95%-CI for"~VaR[a], list(a = alpha))))