Introduction

The volatility of financial markets is a widely studied area. Understanding a market’s volaility gives investors insight into how risky their investments may be. The Chicago Board Option Exchange (CBOE) publishes the Volatility Index (VIX), a financial benchmark that provides an estimate of the stock market’s volatilty[1]. It is calculated based on the prices of S&P 500 index options. The CBOE also offers futures and options for VIX itself. For example, an investor can buy VIX calls if they belive the market will soon become more volatile than it currently is. These VIX instruments provide investors with the opportunity to manage the risk of their other investments.

In this project, I will apply existing volatility models for stock prices on the value of the VIX benchmark. We will see whether the volatility of stock prices act in a similar manner to the volatility of the volatilty of the market. I will experiement with two classes of models: partially observed Markov processes (POMP) and GARCH.

Data

The data is of the VIX index from the beginning of 2018 to the end of 2019 downloaded from the CBOE website[1]. While there was data available up to the current day (April 2020), the dramatic increase in market volatilty due to the COVID-19 pandemic may be difficult to model due to its anomalous nature. This leaves us with a little over 500 days of data.

We can see that the index has an average value between 15 and 20 with very sharp peaks at the beginning and ending of 2018. We will now look at demeaned log returns of the index as they are more likely to form a stationary process without trend that will be easier to analyze.

Here we can see that there is no easily indentifiable trend. Outside of the sharp peak at the beginning of 2018, there is no clear trend in the volatilty of the returns either. This will be the time series the rest of our analysis will be performed on.

POMP Model

Leverage Model

Leverage refers to the phenomenon that large losses dealt to a stock market index are followed by increases in market volatility [2]. Leverage, \(R_n\), on day \(n\) can be formally quantified as the correlation between the index return on day \(n-1\) and increase in the log volatilty from day \(n-1\) to \(n\). We will be using a model introduced by Bretó which models \(R_n\) as: \[R_n = \frac{\exp\{2G_n\}-1}{\exp\{2G_n\}+1} \] where \(G_n\) is a Gaussian random walk [3].

Stochastic Volatility Model

Using Bretó’s model, we can model volatilty as: \[Y_n = \exp(\frac{H_n}{2})\epsilon_n \\ H_n = \mu_h(1-\phi)+\phi H_{n-1} + \beta_{n-1}R_n\exp(\frac{-H_{n-1}}{2})+W_n \\ G_n = G_{n-1} + V_n\]

where: \[H_n = \log(\sigma^2_n) = 2\log(\sigma_n) \\ \beta_n = Y_n\sigma_\eta\sqrt{1-\phi^2} \\ \epsilon_n \sim \text{i.i.d } \mathcal{N}(0,1) \\ W_n \sim \text{i.i.d } \mathcal{N}(0,\sigma^2_\omega) \\ V_n \sim \text{i.i.d } \mathcal{N}(0,\sigma^2_v) \\ \sigma_\omega = \sigma_\eta\sqrt{(1-\phi^2)(1-R^2_n)}\] In the context of our POMP model, \(G_n\) is unobservable; \(Y_n\), the return, is observable; and \(H_n\), the log volatilty, is unobservable.

1Computing the POMP model

We can now estimate the parameters for the POMP model.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   526.4   534.0   538.4   537.5   540.8   545.6

We can see that our iterated filtering processes end us with a median log likelihood of 571.4. We can see the diagnostics of the model.

From the pairwise correlation graph, we can see that the likelihood is maximized over a range of most parameters. From the convergence diagnostics, about half of the particles achieved a log likelihood of about 570 while the other half only achieved about 540. While some of the parameters seem to have converged, others have not. Instead of randomly initializing parameters we can instead define a range from which they should be selected.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   520.0   526.9   530.5   530.5   534.9   537.3

These likelihoods are slightly greater than before. Again, we can look at the diagnostics.

We can see that all particles converged to about the same likelihood unlike last time. The parameters also converge better, but \(\mu_h\) and \(\sigma_\eta\) still do not converge well. While this model is better than the previous model, we cannot be very certain that it is a good model.

GARCH Models

Generalized autoregressive heteroscedasticity (GARCH) models are widely used to model financial time series. A process \(X_n\) is GARCH(p,q) if: \[X_n = \mu_n + \sigma_n\epsilon_n\] where \(\epsilon_n\) is a white noise process with mean 0 and variance 1, and: \[\sigma^2_n = \alpha_0 + \sum_{i=1}^p\beta_i\sigma_{n-i}^2+\sum_{j=1}^q\alpha_j\tilde{X}_{n-j}^2\] where \(\tilde{X}_n = X_n - \mu_n\) [4]. GARCH processes are self exciting in that current volatility simulates an increase in volatility in the near future. GARCH models are good when the error in the volatilties follow an ARMA process. In practice, GARCH(1,1) is used, and it is rare to use a GARCH model of higher order. We will fit a GARCH(1,1) model to our time series of VIX returns.

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = logvix.demean, trace = F) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x7f8c4b48a570>
##  [data = logvix.demean]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##         mu       omega      alpha1       beta1  
## 5.5415e-18  1.2299e-03  2.3926e-01  6.1910e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu     5.541e-18   3.249e-03    0.000 1.000000    
## omega  1.230e-03   3.380e-04    3.639 0.000274 ***
## alpha1 2.393e-01   5.452e-02    4.388 1.14e-05 ***
## beta1  6.191e-01   6.931e-02    8.932  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  547.3731    normalized:  1.090385 
## 
## Description:
##  Sun May  3 19:01:54 2020 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  345.2637  0           
##  Shapiro-Wilk Test  R    W      0.9301414 1.492946e-14
##  Ljung-Box Test     R    Q(10)  5.852502  0.8274994   
##  Ljung-Box Test     R    Q(15)  6.370844  0.9728228   
##  Ljung-Box Test     R    Q(20)  11.32157  0.937453    
##  Ljung-Box Test     R^2  Q(10)  10.56683  0.3922456   
##  Ljung-Box Test     R^2  Q(15)  14.30063  0.502862    
##  Ljung-Box Test     R^2  Q(20)  15.93586  0.7205955   
##  LM Arch Test       R    TR^2   11.10798  0.5196878   
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -2.164833 -2.131219 -2.164959 -2.151645

We can see that the log likelihood is 547.4, which is lower than the POMP model. However, by default, the model assumes the residuals follow a Gaussian normal distribution. From the summary, we can see that the Shaprio-Wilk Test gave a p-value very close to 0, indicating that the residuals are not normally distributed. A look at the QQ plots verifies this.

We can see that the right end of the data is heavier tailed and the left end is lighter tailed than a normal distribution. Therefore, we can instead fit a model where the residuals follow a skewed Student’s t distribution. We also will explore combining the GARCH model with an ARMA model. Here, the conditional mean of the index is modeled using an ARMA process, and its conditional variance is modeled with a GARCH process [5]. We will fit ARMA(p,q)-GARCH(1,1) models for p and q from 0 to 5 and compare the AIC values of each model.

q0 q1 q2 q3 q4 q5
p0 -2.34637 -2.35345 -2.35353 -2.36235 -2.39189 -2.39405
p1 -2.35271 -2.37409 -2.38897 -2.38741 -2.39246 -2.40302
p2 -2.35081 -2.38929 -2.38824 -2.37705 -2.38993 -2.40316
p3 -2.35309 -2.38846 -2.38177 -2.43651 -2.39172 -2.40027
p4 -2.36990 -2.39361 -2.41331 -2.42940 -2.38079 -2.42001
p5 -2.36912 -2.39986 -2.39901 -2.39822 -2.41681 -2.42501

From the table, we can see that the ARMA(3,3)-GARCH(1,1) model has the lowest AIC value. However, this model was found to have invertibility issues. Other models with low AIC values could not have their parameters estimated with significance. The best viable model turned out to be the ARMA(1,1)-GARCH(1,1) model. We fit this model to the data.

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(1, 1) + garch(1, 1), data = logvix.demean, 
##     cond.dist = "sstd", trace = F, grad = "numerical") 
## 
## Mean and Variance Equation:
##  data ~ arma(1, 1) + garch(1, 1)
## <environment: 0x7f8c5003bd28>
##  [data = logvix.demean]
## 
## Conditional Distribution:
##  sstd 
## 
## Coefficient(s):
##          mu          ar1          ma1        omega       alpha1  
## -5.5415e-18   6.9341e-01  -8.3863e-01   5.4019e-04   9.8310e-02  
##       beta1         skew        shape  
##  8.3008e-01   1.6377e+00   5.3387e+00  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -5.541e-18   6.357e-04    0.000  1.00000    
## ar1     6.934e-01   5.830e-02   11.893  < 2e-16 ***
## ma1    -8.386e-01   5.110e-02  -16.411  < 2e-16 ***
## omega   5.402e-04   2.185e-04    2.472  0.01343 *  
## alpha1  9.831e-02   3.761e-02    2.614  0.00895 ** 
## beta1   8.301e-01   5.177e-02   16.035  < 2e-16 ***
## skew    1.638e+00   1.382e-01   11.852  < 2e-16 ***
## shape   5.339e+00   1.265e+00    4.220 2.45e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  603.8977    normalized:  1.202983 
## 
## Description:
##  Sun May  3 19:02:20 2020 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value  
##  Jarque-Bera Test   R    Chi^2  987.5975  0        
##  Shapiro-Wilk Test  R    W      0.8989639 0        
##  Ljung-Box Test     R    Q(10)  12.81964  0.2339316
##  Ljung-Box Test     R    Q(15)  13.11432  0.5934683
##  Ljung-Box Test     R    Q(20)  18.20261  0.5740633
##  Ljung-Box Test     R^2  Q(10)  11.27493  0.3365035
##  Ljung-Box Test     R^2  Q(15)  13.77478  0.5426744
##  Ljung-Box Test     R^2  Q(20)  14.63155  0.7970877
##  LM Arch Test       R    TR^2   11.76621  0.4646339
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -2.374094 -2.306866 -2.374592 -2.347718

We can see that all the parameters are significant, and that the log likelihood is 603.9, greater than our best POMP model. We can see how well the residuals are modeled with the skewed Student’s t distribution with 5.34 degrees of freedom and a skew parameter of 1.64.

We can see that this model better captures the residuals when compared to normal distribution model.

Conclusion

In this project, we saw how the VIX can be modeled using a POMP model and a GARCH model. The log likelihood of the GARCH model was greater. There were also issues with the convergence of some parameters in the POMP model. These factors lead us to say that the GARCH model is better at modeling the volatilty of the VIX index. However, this does not mean that the POMP model is useless. Some of the parameters in the POMP model have better interpretability due to its underlying leverage model. Meanwhile, the parameters of the GARCH model are less interpretable, so while we may model the index better, we don’t learn too much about the nature of the market.

This work gives evidence that the volatility of the volatility of the market can be modeled using models meant for the volatility of markets themselves. Further work could involve running larger POMP models requiring greater computation, as well as looking other volatility indexes that may exist for other markets.

References

  1. http://www.cboe.com/vix
  2. Edward Ionides, “14. Case study: POMP modeling to investigate financial volatility”, https://ionides.github.io/531w18/14/notes14.html#arch-and-garch-models.
  3. Bretó, C. (2014). On idiosyncratic stochasticity of financial leverage efects, Statistics & Probability Letters 91: 20-26. URL: http://www.sciencedirect.com/science/article/pii/S0167715214001345
  4. Brian Thelen, “Lecture 9: Introduction to ARCH/GARCH,” University of Michigan STATS 509 Notes, Winter 2019.
  5. Brian Thelen, “Lecture 10: Advanced Topics in GARCH Time Series Analysis,” University of Michigan STATS 509 Notes, Winter 2019.

I got ideas for this project from the Winter 2018 final project titled “Investigation on Financial Volatility of NASDAQ” (no author listed). Code for this project was adapted from code from the course notes from throughout the semester.