This report is a summary of lesson by Kris Boudt, Data Camp

library(tidyverse)
library(xts)
library(PerformanceAnalytics)
library(rugarch)

options(scipen = 999)

load("sp500prices.Rdata")
load("sp500ret.Rdata")
load("ret.Rdata")
load("msftret.Rdata")
load("EURUSDret.Rdata")

1 The standard GARCH model as the Workhorse model

1.1 Analyzing volatility

1.1.1 Calculating returns

  • CalculateReturns(): compute stock returns
rtn <- CalculateReturns(sp500prices)
plot(rtn, main = "Daily S&P 500 returns")

1.1.2 How to estimate return volatility

  • sd(): compute the (daily) standard deviation:
sd(sp500ret)
## [1] 0.01099357
  • Corresponding formula for \(\hat{\sigma}\) computed using \(T\) daily returns:

\[ \hat{\sigma} = \sqrt{\frac{1}{T-1}\sum_{t=1}^T(R_t - \hat{\mu})^2} \]

  • where \(\hat{\mu}\) is the mean return.

1.1.3 Annualized volatility

  • sd(sp500ret) is daily volatility
  • Annualized volatility = \(\sqrt{252}~\times\) daily volatility
sqrt(252) * sd(sp500ret)
## [1] 0.1745175

1.1.4 Rolling volatility estimation

  • Rolling estimation windows:

  • Window width? Multiple of 22(trading days)
  • chart.RollingPerformance() function
chart.RollingPerformance(R = sp500ret,
                         width = 22,
                         FUN = "sd.annualized",
                         # sd.annualized 함수 인수 scale: 1년 일수
                         scale = 252,
                         main = "Rolling 1 month volatility")

1.2 The GARCH equation for volatility prediction

GARCH(Generalized Autoregressive Conditional Heteroskedasticity) 모델은 조건부 분산(Conditional Variance)을 예측하는 모델

1.2.1 Notation

  • Input: Time series of returns

Notation Ⅰ

  • Each return is observed at a regular frequency, like daily or weekly

Notation Ⅱ

  • At time \(t-1\), you make the prediction about the future return \(R_t\), using the information set available at time \(t-1:I_{t-1}\)

Notation Ⅲ

  • Predicting the mean return: what is the best possible prediction of the actual return?
  • Expected return: \(\mu_t = E~[~R_t~|~I_{t-1}]\)
  • Prediction error: \(e_t = R_t-\mu_t\) ~ \(N(0, \sigma^2_t)\)
    • called shocks or unexpected returns
    • There is a large positive autocorrelation in the absolute value - presence of volatility clusters
    • 즉 절대값 기준으로 prediction error는 서로 뭉쳐있다!!!
# Compute the mean daily return
m <- mean(sp500ret)

# Define the series of prediction errors
e <- sp500ret - m

# Plot the absolute value of the prediction errors
par(mfrow = c(2,1),mar = c(3, 2, 2, 2))
plot(abs(e))

# Plot the acf of the absolute prediction errors
acf(abs(e))

prediction error들은 cluster를 형성하며 대부분의 lag에서 0.2 수준의 autoacorr을 보이고 있다.


Notation Ⅳ

  • We then predict the variance: how far off the return can be from its mean?

\[ \sigma_t^2 = var~(~R_t~|~I_{t-1}) \\ = E[(R_t-\mu_t)^2~|~I_{t-1}] \\ = E[e_t^2~|~I_{t-1}] \\ = \sqrt{\sigma_t^2} \]

From theory to practice: Models for the mean

  1. Rolling mean model
  • \(\mu_t = \frac{1}{M}\sum_{i=1}^MR_{t-i}\) 그냥 rolling기간동안의 산술평균임
  1. ARMA model

From theory to practice: Models for the variance

  1. Rolling variance model
  • \(\sigma_t^2 = \frac{1}{M}\sum_{i=1}^Me^2_{t-i}\)
  • 모든 \(t\)의 가중치 동일하다고 가정
  1. ARCH(p) model
  • \(\sigma_t^2 = \omega + \sum_{i=1}^p\alpha_ie^2_{t-i}\)
  • 최근 관측치에 가중치를 주는 방식(미래의 분산은 최근 사건의 영향을 더 많이 받기 때문에)

1.2.2 GARCH(1,1) model: Generalized ARCH

실무에서는 주로 GARCH(1,1)을 주로 사용함

\[\sigma_t^2 = \omega + \alpha_1 e^2_{t-1} + \beta_1 \sigma^2_{t-1}\]

  • \(\sigma_t^2\): 시점 \(t\)에서의 조건부 변동성(conditional variance)
  • \(e^2_{t-1}\): 시점 \(t-1\)에서의 실제 수익률에서 예측 수익률을 뺀 값의 제곱
  • 직전 기간 prediction error의 제곱과 직전 기간 variance prediction을 각각 계수 \(\alpha\)와 \(\beta\)에 곱해서 계산된다.
    • \(\alpha\): determines the reactivity to \(e^2_t\)
      • 과거 수익률 충격이 현재 변동성에 미치는 영향으로 클수록 최근 충격이 변동성을 크게 증가 시킨다.
    • \(\beta\): the weight on the previous variance prediction
      • 이전 변동성이 현재 변동성에 미치는 영향으로 클수록 변동성이 오래 지속됨(volatility clustering)

Parameter restrictions

To make the GARCH process realistic, we need that:

  1. \(\omega,~\alpha~and~\beta\) are > 0: this ensures that \(\sigma^2_t>0\) at all times.
  2. \(\alpha~+~\beta\) < 1: this ensures that the predicted variance \(\sigma_t^2\) always returns to the long run variance:
  • The variance is therefore “mean-reverting”
  • The long run variance equals \(\frac{\omega}{1-\alpha-\beta}\)

Implementation

# Set parameter values
alpha <- 0.1
beta <- 0.8
omega <- var(sp500ret) * (1- alpha - beta)
# Then: var(sp500ret) = omega / (1 - alpha - beta)

# Set series of prediction error
e <- sp500ret - mean(sp500ret) # Constant mean
e2 <- e ^ 2

nobs <- length(sp500ret)
predvar <- rep(NA, nobs)

predvar[1] <- var(sp500ret)

# Loop starting at 2 because of the lagged predictor
for (t in 2:nobs) {
  predvar[t] <- omega + alpha * e2[t-1] + beta * predvar[t-1]
}

# Volatility is sqrt of predicted variance
predvol <- sqrt(predvar)
predvol <- xts(predvol, order.by = time(sp500ret))

# We compare with the unconditional volatility
uncvol <- sqrt(omega / (1 - alpha - beta))
uncvol <- xts(rep(uncvol, nobs), order.by = time(sp500ret))

# Plot
p <- plot(predvol)
p <- addSeries(uncvol, col = "red", lwd = 2, on = 1)
plot(p)

1.3 The rugarch package

핵심 parameter \(\alpha\)와 \(\beta\)를 추정하는 패키지

The normal GARCH model

\[ R_t = \mu + e_t \\ e_t \sim N(0, \sigma_t^2) \\ \sigma_t^2 = \omega + \alpha_1 e^2_{t-1} + \beta_1 \sigma^2_{t-1} \]

  • Four parameters: \(\mu,\omega,\alpha,\beta\)
  • Estimation by maximum likelihood: find the parameter values for which the GARCH model is most likely to have generated the observed return series.
    • 즉, 현재 관측된 데이터가 가장 잘 설명되도록 하는 파라미터 조합을 찾는 것!!!

1.3.1 Workflow

Three steps:

  1. ugarchspec(): Specify which GARCH model you want to use(mean \(\mu_t\), variance \(\sigma^2_t\), distribution of \(e_t\))
  2. ugarchfit(): Estimate the GARCH model on your time series with returns \(R_1,...,R_T\)
  3. ugarchforecast(): Use the estimated GARCH model to make volatility predictions for \(R_{T+1},...\)

Exercise

  1. ugarchspec() specifies which GARCH model you want to use
garchspec <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
                        variance.model = list(model = "sGARCH"),
                        distribution.model = "norm")
  • mean.model: ARMA model 사용하지 않으므로 c(0,0) 설정
  1. ugarchfit() estimates the GARCH model
  • an object that contains all the results related to the estimation of the GARCH model
  • Methods coef, uncvariance, fitted and sigma
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)

coef(garchfit)
##           mu        omega       alpha1        beta1 
## 0.0005728888 0.0000012202 0.0779146377 0.9111522954
uncvariance(garchfit)
## [1] 0.0001116064
fitted(garchfit)
##            m.c.seq.row..seq.n...seq.col..drop...FALSE.
## 1989-01-04                                0.0005728888
## 1989-01-05                                0.0005728888
## 1989-01-06                                0.0005728888
## 1989-01-09                                0.0005728888
## 1989-01-10                                0.0005728888
## 1989-01-11                                0.0005728888
## 1989-01-12                                0.0005728888
## 1989-01-13                                0.0005728888
## 1989-01-16                                0.0005728888
## 1989-01-17                                0.0005728888
##        ...                                            
## 2017-12-15                                0.0005728888
## 2017-12-18                                0.0005728888
## 2017-12-19                                0.0005728888
## 2017-12-20                                0.0005728888
## 2017-12-21                                0.0005728888
## 2017-12-22                                0.0005728888
## 2017-12-26                                0.0005728888
## 2017-12-27                                0.0005728888
## 2017-12-28                                0.0005728888
## 2017-12-29                                0.0005728888
sigma(garchfit)
##            m.c.seq.row..seq.n...seq.col..drop...FALSE.
## 1989-01-04                                 0.010994656
## 1989-01-05                                 0.011291637
## 1989-01-06                                 0.010842927
## 1989-01-09                                 0.010420726
## 1989-01-10                                 0.010009278
## 1989-01-11                                 0.009647594
## 1989-01-12                                 0.009389706
## 1989-01-13                                 0.009084604
## 1989-01-16                                 0.008757774
## 1989-01-17                                 0.008432999
##        ...                                            
## 2017-12-15                                 0.005064492
## 2017-12-18                                 0.005485433
## 2017-12-19                                 0.005515829
## 2017-12-20                                 0.005483466
## 2017-12-21                                 0.005363765
## 2017-12-22                                 0.005252574
## 2017-12-26                                 0.005142110
## 2017-12-27                                 0.005051692
## 2017-12-28                                 0.004947337
## 2017-12-29                                 0.004862675

Formula:

\[ R_t = 5.7 \times 10^{-4} + e_t \\ e_t \sim N(0, \hat{\sigma}_t^2) \\ \hat{\sigma}_t^2 = 1.2 \times 10^{-6} + 0.08 e^2_{t-1} + 0.91 \hat{\sigma}^2_{t-1} \]

  • \(\alpha~+~\beta = 0.989\)로 1미만이므로, the estimated volatility is mean reverting with long run value(square root of its unconditional variance)
sqrt(uncvariance(garchfit))
## [1] 0.01056439
  • S&P500 일일수익률의 장기 표준 편차는 1% 수준
garchvol <- sigma(garchfit)
plot(garchvol)

tail(garchvol, 1)
##                   [,1]
## 2017-12-29 0.004862675
  • 시계열의 마지막 변동성은 장기 변동성 1%의 절반 수준인 48bp 수준임, 장기 평균을 향해 점차 증가할 것이라고 예상 가능
  • What about the volatility for the days following the end of the time series?
  1. ugarchforecast() forecasts the volatility of the future returns
  • sigma
  • fitted
garchforcast <- ugarchforecast(fitORspec = garchfit, n.ahead = 5)

garchforcast
## 
## *------------------------------------*
## *       GARCH Model Forecast         *
## *------------------------------------*
## Model: sGARCH
## Horizon: 5
## Roll Steps: 0
## Out of Sample: 0
## 
## 0-roll forecast [T0=2017-12-29]:
##        Series    Sigma
## T+1 0.0005729 0.005035
## T+2 0.0005729 0.005127
## T+3 0.0005729 0.005217
## T+4 0.0005729 0.005305
## T+5 0.0005729 0.005390
  • 역시나 향후 sigma가 증가하는 모습을 보이고 있음

1.3.2 Application to tactical asset allocation

A portfolio that invests a percentage \(w\) in a risky asset(with volatility \(\sigma_t\)) and keeps \(1-w\) on a risk-free bank deposit account has volatility equal to

\(\sigma_p = w\sigma_t\)

How to set \(w\)? One approach is volatility targeting: \(w\) is such that the predicted annualized portfolio volatility equals a target level, say 5%. Then:

\(w^* = 0.05/\sigma_t\)

Since GARCH volatilities change, the optimal weight changes as well

Volatility targeting in tactical asset allocation

# Compute the annualized volatility
annualvol <- sqrt(252) * sigma(garchfit)

# Compute the 5% vol target weights  
vt_weights <- 0.05 / annualvol

# Compare the annualized volatility to the portfolio weights in a plot
plot(merge(annualvol, vt_weights), multi.panel = TRUE)

2 Improvements of the Normal GARCH model

2.1 Non-normality of standardized returns

일반적으로 주가 수익률이 normal하다는 것은 비현실적이다.

To account for this, ugarchspec 의 distribution.model 인수 설정을 변경해준다.

  • distribution.model = "sstd"
    • sstd: skew-student t

2.1.1 Testing the normally assumption

# Obtain standardized returns
stdret <- rugarch::residuals(garchfit, standardize = TRUE)

chart.Histogram(sp500ret, methods = c("add.normal", "add.density"),
                colorset = c("grey", "red", "blue"))

chart.Histogram(sp500ret, methods = c("add.normal", "add.density"),
                colorset = c("grey", "red", "blue"),
                xlim = c(-0.05, 0),
                ylim = c(0, 20))

실제 정규분포와 비교해보면 fat tails과 skewness가 관측된다.


2.1.2 Parameters of the skewed student t distribution

  • Compared to the normal distribution, the skewed student t distribution has two extra parameters:
    • Degrees of freedom parameter \(\nu\) (in rugarch: shape): the lower is \(\nu\) the fatter the tails
    • Skewness parameter \(\xi\) (in rugarch: skew):
      • When \(\xi=1\): symmetry.
      • When \(\xi<1\): negative skewness.
      • When \(\xi>1\): positive skewness.
    • Special cases:
      • When \(\nu = \infty\) and \(\xi = 1\): normal distribution
      • When \(\xi = 1\): student t distribution
# Specify the garch model to be used
garchspec <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
                       variance.model = list(model = "sGARCH"),
                        distribution.model = "sstd")

# Estimate the model
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)

# Inspect the coefficients
coef(garchfit)
##              mu           omega          alpha1           beta1            skew 
## 0.0005679523577 0.0000006335867 0.0750402373140 0.9219854622407 0.9437430299527 
##           shape 
## 6.3134746435264

2.2 Leverage effect

The standard GARCH model uses the squared prediction error to forecast the return variance using one single equation. It does not distinguish between positive and negative prediction errors. In reality, the sign of the prediction error matters. Variance shoots up more after a large negative unexpected return than after a large positive unexpected return.

즉, 현실세계에서는 상승할 때 보다 떨어질 때 변동성이 커진다!!!

2.2.1 Two equations

그래서 우리는 2개의 식이 필요하다.

Use seperate equations for the variance following negative and positive unexpected return \(e_t = R_t-\mu_t\):

Case when \(e_{t-1}>0\)

\(\sigma_t^2 = \omega + \alpha e^2_{t-1} + \beta \sigma^2_{t-1}\)

Case when \(e_{t-1}\le0\) * The predicted variance should be higher than after a positive surprise. * This means a higher coefficient multiplying the squared prediction error, namely \(\alpha + \gamma\) instead of \(\alpha\) with \(\gamma \le0\)

\(\sigma_t^2 = \omega + (\alpha+\gamma) e^2_{t-1} + \beta \sigma^2_{t-1}\)

This is the GJR model proposed Glosten, Jagannathan and Runkle.

How?

Change the argument variance.model of ugarchspec() from model = "sGARCH" to model = "gjrGARCH"

garchspec <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
                       variance.model = list(model = "gjrGARCH"),
                        distribution.model = "sstd")

garchfit <- ugarchfit(data = msftret, spec = garchspec)

coef(garchfit)
##             mu          omega         alpha1          beta1         gamma1 
## 0.000604527494 0.000002007641 0.034232579523 0.936331587157 0.055317064482 
##           skew          shape 
## 1.061973467137 4.414858007628
out <- newsimpact(garchfit)
plot(out$zx, out$zy, xlab = "prediction error", ylab = "predicted variance")

  • newsimpact()

prediction error가 0보다 낮을 때 variance가 더 커진다.


Estimation of GJR garch model

# Specify the GJR GARCH model
garchspec <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
                 variance.model = list(model = "gjrGARCH"),
                 distribution.model = "sstd")

# Estimate the model and compute volatility
gjrgarchfit <- ugarchfit(data = msftret, spec = garchspec)
gjrgarchvol <- sigma(gjrgarchfit)

# Specify the GJR GARCH model
sgarchspec <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
                 variance.model = list(model = "sGARCH"),
                 distribution.model = "sstd")

# Estimate the model and compute volatility
sgarchfit <- ugarchfit(data = msftret, spec = sgarchspec)
sgarchvol <- sigma(sgarchfit)

# Compare volatility
plotvol <- plot(abs(msftret), col = "grey", main = "Standard vs. GJR garch model")
plotvol <- addSeries(gjrgarchvol, col = "red", on=1)
plotvol <- addSeries(sgarchvol, col = "blue", on=1)
plotvol

2.3 Mean model

2.3.1 GARCH-in-mean model

  • Quantify the risk-reward trade-off
  • Risk: \(\sigma^2_t\). Reward: \(\mu_t\)
  • GARCH-in-mean_model:

\(\mu_t = \mu+\lambda \sigma_t^2\)

  • \(\lambda>0\) is the risk/reward parameter indicating the increase in expected return per unit of variance risk

How?

Change the argument mean.model in ugarchspec() from list(armaOrder = c(0,0)) to list(armaOrder = c(0,0), archm = TRUE, archpow = 2):

garchspec <- ugarchspec(mean.model = list(armaOrder = c(0,0), archm = TRUE, archpow = 2),
                 variance.model = list(model = "gjrGARCH"),
                 distribution.model = "sstd")
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)
round(coef(garchfit)[1:2], 4)
##     mu  archm 
## 0.0002 2.0079

Predicted mean returns

\(\hat{\mu}_t = 0.0002 + 2.0079 \hat{\sigma}_t^2\)

2.3.2 Today’s return predicts tomorrow’s return

1) AR(1)

  • The GARCH-in-mean uses the financial theory of a risk-reward trade-off to build a conditional mean model.
  • Let’s now use statistical theory to make a mean model that exploits the correlation between today’s return and tomorrow’s return.
  • The most popular model is the AR(1) model:
    • AR(1) stands for autoregressive model of order 1
    • It predicts the next return using the deviation of the return from its long term mean value \(\mu\):

\(\mu_t = \mu + \rho(R_{t-1}-\mu)\)

Autoregressive coefficient: \(\rho\)

  • \(\rho>0\):
    • A higher(resp.lower) than average return is followed by a higher(resp.lower) than average return.
    • Possible explanation: markets underreact to news and hence there is momentum in returns
  • \(|\rho|<1\):
    • Mean reversion: The deviations of \(R_t\) from \(\mu\) are transitory(일시적인)
  • \(\rho<0\)
    • A higher(resp.lower) than average return is followed by a lower(resp.higher) than average return.
    • Possible explanation: markets overreact to news and hence there is reversal in returns

Application to daily S&P 500 returns

Specification and estimation of AR(1)-GJR GARCH with sst distribution

garchspec <- ugarchspec(mean.model = list(armaOrder = c(1,0)),
                 variance.model = list(model = "gjrGARCH"),
                 distribution.model = "sstd")
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)
round(coef(garchfit)[1:2], 4)
##      mu     ar1 
##  0.0003 -0.0292

ar1 = -0.0292로 음수이기에 hints towards overreaction and thus a reversal on the next day


2) MA(1)

The Moving Average model of order 1 uses the deviation of the return from its conditional mean:

\(\mu_t = \mu + \theta(R_{t-1}-\mu_{t-1})\)

3) ARMA(1,1) combines AR(1) and MA(1):

\(\mu_t = \mu + \rho(R_{t-1}-\mu) + \theta(R_{t-1}-\mu_{t-1})\)

Effect of mean model on volatility predictions

Modeling the mean dynamics(평균 역학) typically has a major effect on the obtained predicted returns, but only a minor effect on the volatility predictions. The effect on volatility predictions is so small that, if the interest is only in the volatility dynamics, usually one can ignore the mean dynamics and just assume the most simple specification, namely the constant mean model.

즉, 변동성을 예측하는게 목표라면 mean model은 어떤 걸 사용해도 무방하다!! 일일 수익률의 평균은 대체로 0에 수렴하기때문이다.

2.4 Avoid unnecessary complexity

2.4.1 Use your intuition to avoid unneeded complexity

Restrict the parameter estimates

  • If you know that the parameters
    • are equal to a certain value
    • or, are inside an interval
  • Then you should impose this in the specification using the methods
    • setfixed()
    • setbounds()
garchspec <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
                 variance.model = list(model = "sGARCH"),
                 distribution.model = "std")

setfixed(garchspec) <- list(alpha1 = 0.05, shape = 6)
setbounds(garchspec) <- list(alpha1 = c(0.05, 0.2), beta1 = c(0.8, 0.95))

garchfit <- ugarchfit(data = EURUSDret, spec = garchspec)
coef(garchfit)
##               mu            omega           alpha1            beta1 
## -0.0000410595829  0.0000002069766  0.0500000000000  0.9489742490984 
##            shape 
##  6.0000000000000

2.4.2 Variance targeting

  • Mathematically, this means that the unconditional variance implied by the GARCH models equals the sample variance \(\hat{\sigma}^2\)
    • \(\omega\)를 장기 변동성 수준과 동일하게 맞춰놓고 나머지 parameter를 추정하여 모델을 더 안정적으로 만들고, 수렴 문제를 줄이는 효과가 있다.
    • GARCH model이 적절하게 설정되면, 장기적인 변동성 수준이 실제 데이터에서 측정한 변동성과 같아짐
  • variance.targeting = TRUE in variance.model
garchspec <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
                 variance.model = list(model = "sGARCH",
                                       variance.targeting = TRUE),
                 distribution.model = "std")
garchfit <- ugarchfit(data = EURUSDret, spec = garchspec)
# 허용오차 0.0001 적용
all.equal(uncvariance(garchfit), sd(EURUSDret)^2, tol = 1e-4)
## [1] TRUE

2.4.3 Modeling choices

A GARCH model is a collection of assumptions about the process that generates the returns for which we wish to model the volatility dynamics. According to the French economist Malinvaud (1966) the art is in “trying to find the right set of assumptions which are sufficiently specific, yet realistic to enable us to make the best possible advantage of the available data.”

GARCH모델링의 핵심은 현실적인 가정을 세우는 것이고, 이 가정이 충분히 구체적이면서도 데이터에 적합해야 한다.

3 Performance Evaluation

3.1 Statistical significance

  • Use of statistical tests to decide whether the magnitude of the estimated parameter is significantly large enough to conclude that the true parameter in not zero.
  • How? By comparing the estimated parameter to its standard error

Standard error is the standard deviation of the parameter estimate

3.1.1 t-statistic

t-statistic = estimated parameter / standard error

  • Its absolute value is a distance measure: It tells you how many standard errors the estimated parameter is separated from 0.
  • The larger the distance, the more evidence the parameter is not 0
    • t 통계량이 클수록 통계적으로 유의미하지 않을 가능성이 높으며 0이 아니기에 중요한 변수일 수 있다.

Rule of thumb for deciding:

If |t-statistic| > 2
Then the estimated parameter is statistically significant
Therefore we conclude that true parameter \(\neq\) 0.

3.1.2 p-value

  • The p-value measures how likely it is that the parameter is zero:

Rule of thumb for deciding:

If p-value < 5%
Then the estimated parameter is statistically significant
Therefore we conclude that true parameter \(\neq\) 0.

Example for MSFT returns

flexgarchspec <- ugarchspec(mean.model = list(armaOrder = c(1,0)),
                            variance.model = list(model = "gjrGARCH"),
                            distribution.model = "sstd")
flexgarchfit <- ugarchfit(data = msftret, spec = flexgarchspec)
round(flexgarchfit@fit$matcoef, 6)
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.000610    0.000189   3.220758 0.001279
## ar1    -0.037798    0.013718  -2.755425 0.005862
## omega   0.000002    0.000001   2.617974 0.008845
## alpha1  0.034573    0.003396  10.180935 0.000000
## beta1   0.935928    0.007162 130.687742 0.000000
## gamma1  0.055485    0.009772   5.678128 0.000000
## skew    1.059964    0.020677  51.264038 0.000000
## shape   4.392674    0.256755  17.108435 0.000000

t 통계량과 p-value가 모든 parameter에서 통계적으로 유의미하다.

  • POINT
    • skew parameter는 통계적으로 유의미한 수치로 나옴
    • 근데 skew는 사실상 0이 아니라 1을 기준으로 검정을 실시해야함
    • t 통계량을 다시 계산해보면
      • t-value = \((1.059964-1) / 0.020677 = 2.900034\)
    • 즉, 기준값 1과 표준오차의 2.9배 떨어져 있으므로 여전히 유의미한 통계치임
    • 반대로 2보다 낮아질 경우 분포가 사실상 대칭일 수도 있음


3.2 Goodness of fit

  • Depends on what you want to evaluate
    • the predicted mean
    • the predicted variance
    • the predicted *distribution of the returns**

3.2.1 Goodness of fit for the mean prediction

Based on the estimated GARCH model, we have:

Predicted mean: \(\hat{\mu}_t = E\hat{[R_t|I_{t-1}]}\)

Predicted error: \(e_t = R_t - \hat{\mu}_t\)

Mean squared prediction error: \(\frac{1}{T}\sum_{t=1}^T e_t^2\) (the lower the better)

  • residuals(): 잔차 \(e_t\) 반환하며 이는 표본 수익률과 조건부 평균을 뺀 값
e <- rugarch::residuals(flexgarchfit)
mean(e ^ 2)
## [1] 0.0003857876

3.2.2 Goodness of fit for the variance prediction

The GARCH model leads to:

Predicted variance: \(\hat{\sigma}_t^2 = var~\hat{(~R_t~|~I_{t-1})}\)
\(= E\hat{[(R_t-\mu_t)^2~|~I_{t-1}]}\)
\(= E\hat{[e_t^2~|~I_{t-1}]}\)

Prediction error: \(d_t = e^2_t - \hat{\sigma}_t^2\)

Mean squared prediction error: \(\frac{1}{T}\sum_{t=1}^T d_t^2\) (the lower the better)

  • sigma(): GARCH모델의 조건부 표준편차(\(\sigma_t\))를 반환
e <- rugarch::residuals(flexgarchfit)
d <- e ^ 2 - sigma(flexgarchfit) ^ 2
mean(d ^ 2)
## [1] 0.000001577261

3.2.3 Goodness of fit for the distribution

  • The GARCH model provides a predicted density for all the returns in the sample
    • The higer the density, the more likely the return is under the estimated GARCH model
    • The likelihood of the sample is based on the product of all these densities. It measures how likely it is to that the observed returns come from the estimated GARCH model
    • The higher the likelihood, the better the model fits with your data
    • likelihood() function
      • 값 자체보다는 다른 모델의 값과 비교해서 더 높으면 good
likelihood(flexgarchfit)
## [1] 13126.92

3.2.4 Risk of overfitting

We use an in-sample evaluation approach where the estimation sample and evaluation sample coincides.

Danger of overfitting:

Overfitting consists of choosing overly complex models that provides a good fit for the returns in the estimation sample used but not for the future returns that are outside of the sample

Solution: Balance goodness of fit with a penalty for the complexity

  • A GARCH model is parsimonious(인색한, 검소한) when it has:
    • a high likelihood
    • and a relatively low number of parameters

Information criteria

  • Parsimony is measured information criteria

information criteria = - likelihood + penalty(number of parameters)

  • The lower, the better

Rule of thumb for deciding:
Choose the model with lowest information criterion.

  • infocriteria()
infocriteria(flexgarchfit)
##                       
## Akaike       -5.489087
## Bayes        -5.478255
## Shibata      -5.489092
## Hannan-Quinn -5.485282

3.3 Diagnosing absolute standardized returns

Check 1: Mean and standard deviation of standardized returns

Formula standardized returns:

\(Z_t = \frac{R_t-\hat{\mu}_t}{\hat{\sigma}_t}\)

  • First check of model validity:
    • Sample mean of standardized returns \(\approx 0\)
    • Sample standard deviation of standardized returns \(\approx 1\)

Check 2: Time series plot of standardized returns

  • Second check of model validity:
    • time series plot of standardized returns
    • standardized returns should have constant variability

Check 3: No predictability in the absolute standardized returns

  • Third check of model validity:
    • verify that there is no correlation between the past absolute standardized returns and the current absolute standardized returns.
    • this means:\(Corr(|Z_{t-k}|,|Z_t|) \approx 0,~for~k > 0\)
  • Why?
    • The magnitude of the absolute standardized returns should be constant -> no correlations in the absolute standardized returns.
  • acf() function

Check 4: Ljung-Box test

  • Fourth check of model validity:
    • Ljung-Box test that the first k autocorrelations in the absolute standardized returns \(|Z_t|\) are zero
    • Similar to a \(t\)-test for statistical significance of the estimated parameters, but here we want to have 0 in order to have a good model.

Rule of thumb: p-value less than 5% indicates that the model used is not valid
즉 반대로 p-value가 5%보다 더 커야 모델이 유효하다.

  • Box.test(type = "Ljung-Box") function

Application to MSFT

garchspec <- ugarchspec(mean.model = list(armaOrder = c(1,0)),
                        variance.model = list(model = "gjrGARCH"),
                        distribution.model = "sstd")
garchfit <- ugarchfit(data = msftret, spec = garchspec)
stdmsftret <- rugarch::residuals(garchfit, standardize = TRUE)

acf(abs(msftret), 22, main = "MSFT")

acf(abs(stdmsftret), 22, main = "MSFT - standardized returns")

Box.test(abs(stdmsftret), 22, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  abs(stdmsftret)
## X-squared = 25.245, df = 22, p-value = 0.2855

3.4 Backtesting using ugarchroll

Use only the data that were available at the time of prediction

Volatility prediction by applying sigma() to ugarchfit object

Look ahead bias: * Future returns are used to make the volatility estimate. * 즉 과거, 현재 그리고 미래 수익률 전부가 과거 변동성 추정치를 얻기 위해 사용됨

sigma(garchfit)
##            m.c.seq.row..seq.n...seq.col..drop...FALSE.
## 1999-01-04                                  0.01964148
## 1999-01-05                                  0.01928779
## 1999-01-06                                  0.02007005
## 1999-01-07                                  0.02042743
## 1999-01-08                                  0.01985669
## 1999-01-11                                  0.01932047
## 1999-01-12                                  0.01939895
## 1999-01-13                                  0.02189037
## 1999-01-14                                  0.02129767
## 1999-01-15                                  0.02110930
##        ...                                            
## 2017-12-15                                  0.01366397
## 2017-12-18                                  0.01405902
## 2017-12-19                                  0.01376024
## 2017-12-20                                  0.01356153
## 2017-12-21                                  0.01326542
## 2017-12-22                                  0.01291569
## 2017-12-26                                  0.01257694
## 2017-12-27                                  0.01226387
## 2017-12-28                                  0.01196220
## 2017-12-29                                  0.01166048

Solution to avoid look-ahead bias: Rolling estimation

  • Program a for loop: Iterate over the prediction times and use ugarchfit() to estimate the model and ugarchforecast() to make the volatility forecast
  • Built-in function ugarchroll() in the rugarch package
  • Options:
    • Length of the estimation sample (표본 길이)
    • Frequency of model estimation (추정 빈도)

Moving window estimation

  • Reduce the computational cost by estimating the model every \(K\) observations

3.4.1 Function ugarchroll

garchroll <- ugarchroll(garchspec, data = EURUSDret, n.start = 2500, refit.window = "moving", refit.every = 500)

Arguments to specify:

  1. GARCH specification used
  2. data: return data to use
  3. n.start: the size of the initial estimation sample
  4. refit.window: how to change that sample through time: “moving” or “expanding” (롤링 방식)
  5. refit.every: how often to re-estimate the model (재추정 빈도)
  • return값은 as.data.frame으로 변환시켜줘야함

Exercise: In-sample versus rolling sample vol

garchspec <- ugarchspec(mean.model = list(armaOrder = c(1,0)),
                        variance.model = list(model = "gjrGARCH"),
                        distribution.model = "sstd")

# Estimate the GARCH model using all the returns and compute the in-sample estimates of volatility
garchinsample <- ugarchfit(data = sp500ret, spec = garchspec)
garchvolinsample <- sigma(garchinsample)

# Use ugarchroll for rolling estimation of the GARCH model 
garchroll <- ugarchroll(garchspec, data = sp500ret, 
        n.start = 2000, refit.window = "moving", refit.every = 2500)

# Set preds to the data frame with rolling predictions
preds <- as.data.frame(garchroll)

# Compare in-sample and rolling sample volatility in one plot
garchvolroll <- xts(preds$Sigma, order.by = as.Date(rownames(preds)))
volplot <- plot(garchvolinsample, col = "darkgrey", lwd = 1.5, main = "In-sample versus rolling vol forecasts")
volplot <- addSeries(garchvolroll, col = "blue", on = 1)
plot(volplot)

변동성이 처음에는 불일치하다가 점점 일치되는 모습이다.


4 Application

4.1 Value-at-risk

  • A popular measure of downside risk: 5% VaR

  • Forward looking risk management uses the predicted quantiles from the GARCH estimation
  • How? Method quantile() applied to a ugarchroll object

4.1.1 Workflow to obtain predicted 5% quantiles from ugarchroll

# Specify which GARCH model you want to use
garchspec <- ugarchspec(mean.model = list(armaOrder = c(1,0)),
                        variance.model = list(model = "gjrGARCH"),
                        distribution.model = "sstd")
# Estimate the GARCH model on rolling estimation samples
garchroll <- ugarchroll(garchspec, data = sp500ret, 
        n.start = 2500, refit.window = "moving", refit.every = 100)
# Compute the predicted quantile
garchVaR <- quantile(garchroll, probs = 0.05)

# Plot
actual <- xts(as.data.frame(garchroll)$Realized, time(garchVaR))
VaRplot(alpha = 0.05, actual = actual, VaR = garchVaR)

4.1.2 Exceedance and VaR coverage

A VaR exceedance occurs when the actual return is less than the predicted value-at-risk:
\(R_t < VaR_t\)

The frequency of VaR exceedances is called the VaR coverage

mean(actual < garchVaR)
## [1] 0.05159143

유의수준인 5% 수준을 보이고 있다!!


  • If coverage > \(\alpha\): too many exceedances: the predicted quantile should be more negative.
    • Risk of losing money has been underestimated.
    • 예상보다 돈 더 잃는 거 같은데? 더 보수적으로 해야해
  • If coverage < \(\alpha\): too many exceedances: the predicted quantile was too negative.
    • Risk of losing money has been overestimated.
    • 너무 보수적으로 설정한거 같아. 과대평가 했어

4.2 For production and simulation

4.2.1 Use in production

production이란 개발한 모델을 실제 환경에서 운용하는 단계를 의미하며 이제부터 out of sample forecasting이기에 개발 단계에서와는 다른 데이터를 사용하게 된다.

1) New functionality

  • ugarchfilter(): for analyzing the recent dynamics in the mean and volatility
  • ugarchforecast(): applied to a ugarchspec object(instead of ugarchfit()) object for making predictions about the future mean and volatility

2) ugarchfit vs. ugarchfilter

  • ugarchfit은 GARCH모델을 데이터에 적합(fit)하고 파라미터를 추정하는 과정
  • ugarchfilter는 이미 주어진 GARCH 모델의 파라미터를 고정한 상태에서 새로운 데이터에 대해 조건부 분산과 잔차를 계산
  • 즉 ugarchfit은 모델을 학습하는 과정, ugarchfilter는 기존 모델을 적용하여 새로운 데이터의 변동성을 계산하는 과정이다.
  • 만약 ugarchfit을 새로운 데이터가 들어올 때마다 실행한다면 과거 데이터 전체를 다시 학습하므로, 미래 정보가 반영될 위험이 있기에 이미 적합된 모델을 유지하면서 새로운 데이터에서 변동성을 업데이트 하기 위함

3) Example on MSFT returns

  • msft: 1999-2017 daily returns
  • 2010년 말 기준 데이터로 model fitting
  • 2017년 말 미래 변동성 예측하기 위해 해당 모델 사용

Step 1: Defines the final model specification

# Fit the best model
garchspec <- ugarchspec(mean.model = list(armaOrder = c(1,0)),
                        variance.model = list(model = "gjrGARCH"),
                        distribution.model = "sstd")
garchfit <- ugarchfit(data = msftret["/2010-12"], spec = garchspec)

progarchspec <- garchspec
setfixed(progarchspec) <- as.list(coef(garchfit))

Step 2: Analysis of mean and volatility dynamics

garchfiliter <- ugarchfilter(data = msftret, spec = progarchspec)
plot(sigma(garchfiliter))

Step 3: Make predictions about future returns

garchforecast <- ugarchforecast(data = msftret,
                                fitORspec = progarchspec,
                                n.ahead = 10)
cbind(fitted(garchforecast), sigma(garchforecast))
##        2017-12-29 2017-12-29
## T+1  0.0004781652 0.01124870
## T+2  0.0003610396 0.01132549
## T+3  0.0003663608 0.01140170
## T+4  0.0003661190 0.01147733
## T+5  0.0003661300 0.01155238
## T+6  0.0003661295 0.01162688
## T+7  0.0003661295 0.01170082
## T+8  0.0003661295 0.01177423
## T+9  0.0003661295 0.01184712
## T+10 0.0003661295 0.01191948

4.2.2 Use in simulation

Instead of applying the complete model to analyze observed returns, you can use it to simulate artificial log-returns:

\(r_t = \log(P_t) - \log(P_{t-1})\)

Useful to assess the randomness in future returns and the impact on prices, since the future price equals:

\(P_{t+h} = P_t~\exp(r_{t+1}~+~r_{t+2}~+~...~+~r_{t+h})\)

미래 \(t+h\) 시점 가격은 \(t\) 시점 가격에 \(t\)와 \(t+h\) 시점 사이 로그 수익률의 합을 곱한 것이다.

Step 1: Calibrate the simulation model

Use the log-returns in the estimation

sp500logret <- diff(log(sp500prices))[-1]

garchspec <- ugarchspec(mean.model = list(armaOrder = c(1,0)),
                        variance.model = list(model = "gjrGARCH"),
                        distribution.model = "sstd")
garchfit <- ugarchfit(data = sp500logret, spec = garchspec)

simgarchspec <- garchspec
setfixed(simgarchspec) <- as.list(coef(garchfit))

Step 2: Run the simulation with ugarchpath()

Simulation using the ugarchpath() function requires to choose:

  • spec: completely specified GARCH model
  • m.sim: number of time series of simulated returns you want
  • n.sim: number of observations in the simulated time series(e.g. 252)
  • rseed: any number to fix the seed used to generate the simulated series(needed for reproducibility)
simgarch <- ugarchpath(spec = simgarchspec,
                       m.sim = 4,
                       n.sim = 10 * 252,
                       rseed = 12345)

Step 3: Analysis of simulated returns

Method fitted() provides the simulated returns:

simret <- fitted(simgarch)
plot.zoo(simret)

plot.zoo(sigma(simgarch))

Bonus: Analysis of simulated prices

Plotting 4 simulations of 10 years of stock prices, with initial price set at 1:

simprices <- exp(apply(simret, 2, "cumsum"))
matplot(simprices, type = "l", lwd = 3)

4.3 Model risk

There is always a risk of making decisions by using a wrong model

Sources of model risk and solutions

  • Sources:
    • modeling choices
    • starting values in the optimization
    • outliers in the return series
  • Solution: Protect yourself through a robust approach
    • model-averaging: averaging the predictions of multiple models
    • trying several starting values and choosing the one that leads to the highest likelihood
    • cleaning the data

4.3.1 Model averaging

library(purrr)

full_df <- expand.grid(
  variance_models = c("sGARCH", "gjrGARCH"),
  distribution_models = c("norm", "std", "sstd"),
  stringsAsFactors = FALSE
)

full_df <- full_df %>% 
  mutate(garchspec = map2(variance_models, distribution_models, ~ ugarchspec(variance.model = list(model = .x), 
                                                                             mean.model = list(armaOrder = c(0,0)),
                                                                             distribution.model = .y)),
         garchfit = map(garchspec, ~ ugarchfit(data = msftret, spec = .x)),
         sigma = map(garchfit, sigma))

msigma <- purrr::reduce(full_df$sigma, merge)
plot.xts(msigma, main = "Sigma : 6 combinations of GARCH models")

6가지 모델의 조합들은 서로 다른 조건임에도 강력한 상관관계를 보이고 있습니다. 이들의 평균을 선택합니다.


avesigma <- xts(rowMeans(msigma), order.by = time(msigma))
plot.xts(avesigma)

4.3.2 Robustness to starting values

  • GARCH estimates are the result of a complex optimization of the likelihood function.
  • Optimization is numeric and iterative: step by step improvement, which can be sensitive to starting values
  • rugarch has a default approach in getting sensible starting values
  • You can specify your own starting values by applying the setstart() method to your ugarchspec() GARCH model specification
coef(garchfit)
##              mu             ar1           omega          alpha1           beta1 
##  0.000313538306 -0.029728513352  0.000001109796  0.001577025057  0.918098634597 
##          gamma1            skew           shape 
##  0.140188716993  0.907605478105  6.867275288833
likelihood(garchfit)
## [1] 24378.18
setstart(garchspec) <- list(alpha1 = 0.05, beta1 = 0.9)
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)
coef(garchfit)
##              mu             ar1           omega          alpha1           beta1 
##  0.000340979959 -0.029162105125  0.000001071185  0.001550946695  0.918146310962 
##          gamma1            skew           shape 
##  0.143043402013  0.917627330722  6.942479840645
likelihood(garchfit)
## [1] 24377.2

setstart()로 기본 파라미터 시작점을 설정하면 새로운 추정치가 비슷한 likelihood가 출력됩니다. 시작점에 대한 민감도가 거의 없으므로 신뢰도가 상승합니다.


4.3.3 Cleaning the data

  • Avoid that outliers distort the volatility predictions
  • How? Through winsorization: reduce the magnitude of the return to an acceptable level using the function Return.clean() in the PerformanceAnalytics with method = "boudt"
    • winsorization은 극단치를 제거하는 대신 특정 백분위 값으로 치환하는 방식입니다.
library(PerformanceAnalytics)
clmsftret <- Return.clean(msftret, method = "boudt")

plotret <- plot(msftret, col = "red", main = "MSFT(Red) vs. Cleaned MSFT(Blue)")
plotret <- addSeries(clmsftret, col = "blue", on = 1)
plotret

변동성 클러스터는 유지하면서 극단적인 것들의 winsorization 된 것이 확인된다.


4.4 GARCH covariance

GARCH volatility leads to time-varying variability of the returns. When variance between assets is time-varying, their covariance changes over time.

If two asset returns \(R_{1,t}\) and \(R_{2,t}\) have correlation \(\rho\) and time varying volatility \(\sigma_{1,t}\) and \(\sigma_{2,t}\), then their covariance is:

\(\sigma_{12,t} = \rho \sigma_{1,t} \sigma_{2,t}\)

4.4.1 GARCH covariance estimation

Step 1: Use ugarchfit() to estimate the GARCH model for each return series.

msftgarchfit <- ugarchfit(data = msftret["2008-06-20/"], spec = garchspec)
retgarchfit <- ugarchfit(data = ret, spec = garchspec)

Step 2: Use residuals()` to compute the standardized returns.

stdmsftret <- residuals(msftgarchfit, standardize = TRUE)
stdret <- residuals(retgarchfit, standardize = TRUE)

Step 3: Use cor() to estimate \(\rho\) as the sample correlation of the standardized returns.

msftretcor <- as.numeric(cor(stdmsftret, stdret))
msftretcor
## [1] 0.04654836

Step 4: Compute the GARCH covariance by multiplying the estimated correlation and volatilities

msftretcov <- msftretcor * sigma(msftgarchfit) * sigma(retgarchfit)
plot.xts(msftretcov)

4.4.2 Applications of covariance in finance

  • Optimizing the variance of the portfolio
    • It depends on the:
      • portfolio weights
      • the variance of all the assets
      • the covariance between the asset returns
  • Dynamic beta
    • The estimation of a stock’s beta is the systematic risk of a stock