Introduction to QuantileModels package

Christian Jorge Carreiro

This document introduces the capabilities of this package, the available models, and how to estimate them and obtain the results of those estimates. In short, it serves as a brief tutorial and user guide.

This package enables users to estimate the models proposed by Engle & Manganelli (2004), as well as their subsequent multivariate version presented by White et al. (2015). Unlike other currently available functions or code1, this package enables estimation of these models for autoregressive quantile lags of order p, and q for values of the variable itself.

1 Conditional Autoregressive value at risk (CAViaR)

This package includes four different specifications, which are listed below:

\[f_t(\theta) = \beta_0 + \sum_{i=1}^p \beta_i f_{t-i}(\theta) + \sum_{j=1}^q \gamma_j |y_{t-j}|\]

\[f_t(\theta) = \beta_0 + \sum_{i=1}^p \beta_i f_{t-i}(\theta) + \sum_{j=1}^q \left( \gamma_{1,j} (y_{t-j})^+ + \gamma_{2,j} (y_{t-j})^- \right)\]

where \({(y_{t-j})^+ = \max(y_{t-j}, 0)}\) and \((y_{t-j})^- = \min(y_{t-j}, 0)\).

\[f_t(\theta) = \left( \beta_0 + \sum_{i=1}^p \beta_i f_{t-i}^2(\theta) + \sum_{j=1}^q \gamma_j y_{t-j}^2 \right)^{1/2}\]

\[f_t(\theta) = \beta_0 + \beta_1 f_{t-1}(\theta) + (1 - \beta_1)\left(\frac{\nu}{1-\gamma_1} I(y_{t-1} > 0)+ \frac{\nu}{\gamma_1} I(y_{t-1} \leq 0)\right) |y_{t-1} - u|\]

where \({\nu = \sqrt{\gamma_1^2 + (1 - \gamma_1)^2}}\) ,\(0<\gamma_1<1\) and \(u\) is the sample mean.

The following section illustrates how to use the CAViaR() function to estimate this model. It does so using data from Engle & Manganelli (2004) for the GM series at the 5% quantile.

data=dataCAViaR

SAV <- CAViaR(Y=data$GM[1:2892],
model.type = "SAV",p=1,q=1,band.hs = TRUE,quant.type = 7,
tau=0.05,refine.estim = FALSE)

AS <- CAViaR(Y=data$GM[1:2892],
model.type = "AS",p=1,q=1,band.hs = TRUE,quant.type = 7,
tau=0.05,refine.estim = FALSE)

INDGARCH <- CAViaR(Y=data$GM[1:2892],
model.type = "INDGARCH",p=1,q=1,band.hs = TRUE,quant.type = 7,
tau=0.05,refine.estim = FALSE)

I_CAV <- CAViaR(Y=data$GM[1:2892],
model.type = "I-CAV",p=1,q=1,band.hs = TRUE,quant.type = 7,
tau=0.05,refine.estim = FALSE)

To see the results, use the summary() or print() functions. This will show you important information about the estimation process, the estimated coefficients, standard errors, p-values, and confidence intervals for the coefficients. It will also show you coverage tests for the value at risk.

summary(SAV)
#> CAViaR estimation 
#> -------------------
#> Model specification: SAV 
#> Quantile (tau): 0.05 
#> Loss function value at estimates: 551.2903 
#> In-sample coverage: 0.05014 
#> Hall-Sheather bandwidth 
#> 
#> Estimation results: 
#> ======================================================== 
#>            Coef.     S.E   P>|t|  2.5% CI 97.5% CI
#> Beta 0  -0.15817 0.09835 0.10788 -0.35102  0.03467
#> Beta 1   0.88566 0.04309 0.00000  0.80118  0.97014
#> Gamma 1 -0.11445 0.01711 0.00000 -0.14800 -0.08091
#> ======================================================== 
#> 
#> Coverage test 
#> -------------------
#> Kupiec conditional coverage test (LRcc), p-value: 0.91283 
#> Christoffersen independence test (LRind), p-value: 0.67031 
#> Christoffersen unconditional coverage test (LRuc), p-value: 0.97279

The following table shows the results of the estimation for the various specifications, which can be compared with those in the original study by Engle & Manganelli (2004). As is evident, the results are virtually identical. The only noticeable difference lies in the standard errors, which differ for probably two reasons:

  1. Engle & Manganelli (2004) use an “exact formula” to calculate the gradient (Jacobian) of the quantile process function, which is obviously preferable. However, in this specific function, these are calculated using numerical approximations with central finite differences via numDeriv (Gilbert & Varadhan, 2019).

  2. In the original work, the bandwidth is calculated using k-nearest neighbors. However, in this implementation, it is calculated based on the multivariate posterior proposal (White et al., 2015), whereas they use the one proposed by Machado (2013).

Estimation results for GM VaR at 5% under different CAViaR specifications.
SAV AS INDGARCH I-CAV
\(\beta_0\) -0.1582 -0.0760 0.3336 -0.0384
(0.0984) (0.0325) (0.1174) (0.0226)
\(\beta_1\) 0.8857 0.9326 0.9042 0.9419
(0.0431) (0.0142) (0.0320) (0.0143)
\(\gamma_1\) -0.1145 - 0.1220 0.3690
(0.0171) - (0.0812) (0.0747)
\(\gamma_1^+\) - -0.03977 - -
- (0.0210) - -
\(\gamma_1^-\) - 0.1218 - -
- (0.0145) - -
\(RQ\) 551.29 548.30 552.12 550.10
\(\hat{\tau}\) 0.0501 0.0498 0.0501 0.0498

Here are the results reported in Engle & Manganelli (2004).

Parameter estimation and results from Engle & Manganelli (2004)

In addition, we can visualize the conditional quantile process for GM and compare it with Figure 1 from the original paper.

graphic_opts <- par(no.readonly = TRUE)
par(mfrow=c(2,2))

plot(-SAV$VaR,.by="quarter",ylim=c(1.5,10),main="SAV",main.timespan=FALSE)
plot(-AS$VaR,.by="quarter",ylim=c(1.5,10),main="AS",main.timespan=FALSE)
plot(-INDGARCH$VaR,.by="quarter",ylim=c(1.5,10),main="INDGARCH",main.timespan=FALSE)
plot(-I_CAV$VaR,.by="quarter",ylim=c(1.5,10),main="I-CAV",main.timespan=FALSE)
GM VaR at 5% under different specifications
GM VaR at 5% under different specifications
par(graphic_opts)

Note that in this R package, the adaptive specification is not implemented; therefore, you should compare only panels a) through c) and ignore the I-CAV, as it is shown only for visualization and comparison with the other specifications2. Furthermore, for all specifications, an order greater than 1 can be used for both the autoregressive quantile and the lagged value of the series, except for the I-CAV, which currently only allows for estimating I-CAV(1,1).

“Estimated CAViaR Plots for GM: (a) Symmetric Absolute Value; (b) Asymmetric Slope; (c) GARCH; (d) Adaptive”. Source: Engle & Manganelli (2004)
“Estimated CAViaR Plots for GM: (a) Symmetric Absolute Value; (b) Asymmetric Slope; (c) GARCH; (d) Adaptive”. Source: Engle & Manganelli (2004)

Furthermore, you can use the S3 plot() method for the object returned by the CAViaR() function to visualize the result graphically.

plot(SAV,titl="GM Var at 5%")

2 VAR for VaR or MVMQ-CAViaR

As well as the univariate model estimation, it is also possible to estimate its multivariate extension proposed by White et al. (2015), which has the following specification:

\[\boldsymbol{f_t}(\boldsymbol{\theta}) = \boldsymbol{c} + \sum_{i=1}^p \boldsymbol{A_i} \boldsymbol{f_{t-i}(\boldsymbol{\theta})} + \sum_{j=1}^q \boldsymbol{B_j} |\boldsymbol{Y_{t-j}|}\]

This model can be estimated using the MVMQ_CAViaR() function. Once again, using the data from White et al. (2015), the goal is to reproduce Figure 1 and part of Table 3 (the results for Barclays and Goldman Sachs; the rest is left for the user to verify) as follows:

Barclays <- MVMQ_CAViaR(MVMQ[,c(6,1)],tau =c(0.01,0.01),band.hs = TRUE)

Goldman <- MVMQ_CAViaR(MVMQ[,c(7,4)],tau =c(0.01,0.01),band.hs = TRUE)

HSBC <- MVMQ_CAViaR(MVMQ[,c(5,3)],tau =c(0.01,0.01),band.hs = TRUE)

Deutsche <- MVMQ_CAViaR(MVMQ[,c(6,2)],tau =c(0.01,0.01),band.hs = TRUE)

Using the summary() function

summary(Barclays)
#> MVMQ CAViaR estimation 
#> Loss function at estimates: 324.0218 
#> Hall-Sheather bandwidth 
#> ======================================================== 
#> Equation: yEU_index 
#> Quantile (tau): 0.01 
#> In sample coverage 0.01049 
#> Estimation results: 
#> -------------------------------------------------------- 
#>                               Coef.     S.E   P>|t|  2.5% CI 97.5% CI
#> Const.                     -0.13209 0.04568 0.00387 -0.22166 -0.04251
#> q.yEU_index,1               0.80808 0.04854 0.00000  0.71290  0.90326
#> q.BARCLAYS - PRICE INDEX,1 -0.01093 0.00695 0.11565 -0.02455  0.00269
#> |yEU_index|,1              -0.51111 0.12095 0.00002 -0.74827 -0.27395
#> |BARCLAYS - PRICE INDEX|,1 -0.05018 0.01083 0.00000 -0.07142 -0.02894
#> 
#> Coverage test 
#> -------------------
#> Kupiec conditional coverage test (LRcc), p-value: 0.70412 
#> Christoffersen independence test (LRind), p-value: 0.42513 
#> Christoffersen unconditional coverage test (LRuc), p-value: 0.79796 
#> -------------------------------------------------------- 
#> 
#> Equation: BARCLAYS - PRICE INDEX 
#> Quantile (tau): 0.01 
#> In sample coverage 0.00976 
#> Estimation results: 
#> -------------------------------------------------------- 
#>                               Coef.     S.E   P>|t|  2.5% CI 97.5% CI
#> Const.                     -0.09075 0.04920 0.06519 -0.18722  0.00571
#> q.yEU_index,1              -0.12595 0.06035 0.03697 -0.24428 -0.00762
#> q.BARCLAYS - PRICE INDEX,1  0.95959 0.01066 0.00000  0.93868  0.98049
#> |yEU_index|,1              -0.33446 0.12710 0.00855 -0.58369 -0.08524
#> |BARCLAYS - PRICE INDEX|,1 -0.14545 0.07610 0.05608 -0.29467  0.00378
#> 
#> Coverage test 
#> -------------------
#> Kupiec conditional coverage test (LRcc), p-value: 0.08638 
#> Christoffersen independence test (LRind), p-value: 0.02713 
#> Christoffersen unconditional coverage test (LRuc), p-value: 0.90074 
#> --------------------------------------------------------

summary(Goldman)
#> MVMQ CAViaR estimation 
#> Loss function at estimates: 338.3715 
#> Hall-Sheather bandwidth 
#> ======================================================== 
#> Equation: yNA_index 
#> Quantile (tau): 0.01 
#> In sample coverage 0.01013 
#> Estimation results: 
#> -------------------------------------------------------- 
#>                      Coef.     S.E   P>|t|  2.5% CI 97.5% CI
#> Const.            -0.03359 0.02090 0.10820 -0.07458  0.00740
#> q.yNA_index,1      0.92515 0.03924 0.00000  0.84821  1.00208
#> q.GOLDMAN SACHS,1 -0.02328 0.02567 0.36470 -0.07362  0.02707
#> |yNA_index|,1     -0.20184 0.10772 0.06109 -0.41306  0.00939
#> |GOLDMAN SACHS|,1 -0.08306 0.06843 0.22497 -0.21724  0.05113
#> 
#> Coverage test 
#> -------------------
#> Kupiec conditional coverage test (LRcc), p-value: 0.74159 
#> Christoffersen independence test (LRind), p-value: 0.44108 
#> Christoffersen unconditional coverage test (LRuc), p-value: 0.94677 
#> -------------------------------------------------------- 
#> 
#> Equation: GOLDMAN SACHS 
#> Quantile (tau): 0.01 
#> In sample coverage 0.00976 
#> Estimation results: 
#> -------------------------------------------------------- 
#>                      Coef.     S.E   P>|t|  2.5% CI 97.5% CI
#> Const.            -0.01989 0.03179 0.53153 -0.08223  0.04244
#> q.yNA_index,1     -0.00081 0.04793 0.98655 -0.09479  0.09317
#> q.GOLDMAN SACHS,1  0.94425 0.03314 0.00000  0.87928  1.00923
#> |yNA_index|,1     -0.00140 0.15193 0.99267 -0.29931  0.29652
#> |GOLDMAN SACHS|,1 -0.16411 0.08959 0.06710 -0.33979  0.01157
#> 
#> Coverage test 
#> -------------------
#> Kupiec conditional coverage test (LRcc), p-value: 0.08638 
#> Christoffersen independence test (LRind), p-value: 0.02713 
#> Christoffersen unconditional coverage test (LRuc), p-value: 0.90074 
#> --------------------------------------------------------

It can be seen that the estimates are very similar to those obtained by the authors, although there are differences in certain parameters and their respective significance. In this regard, future improvements could be made to ensure that the results are more consistent.

However, from a graphical perspective, these results generate virtually the same quantile process, as can be seen in the following graphs, which are virtually identical to those in Figure 1 of White et al. (2015).

dates <- as.Date(zoo::index(MVMQ))

#Barclays
plot(dates,as.vector(MVMQ[,1]),type = "p",ylim = c(-60,60),cex=0.6,xaxt="n",cex.axis=0.8,col=2,xlab="",ylab = "",main = "Barclays")
lines(dates,Barclays[[5]][,2],type = "l")
axis.Date(side=1,at=seq(dates[1],dates[2765],by="year"),cex.axis=0.8,las=2)
legend("topleft",legend = c("Daily Returns","1% VaR"),col = 2:1,pch = c("o","-"))


#Deutsche
plot(dates,as.vector(MVMQ[,2]),type = "p",ylim = c(-60,60),cex=0.6,xaxt="n",cex.axis=0.8,col=2,xlab="",ylab = "",main = "Deutsche")
lines(dates,Deutsche[[5]][,2],type = "l")
axis.Date(side=1,at=seq(dates[1],dates[2765],by="year"),cex.axis=0.8,las=2)
legend("topleft",legend = c("Daily Returns","1% VaR"),col = 2:1,pch = c("o","-"))


#Goldman
plot(dates,as.vector(MVMQ[,4]),type = "p",ylim = c(-60,60),cex=0.6,xaxt="n",cex.axis=0.8,col=2,xlab="",ylab = "",main = "Goldman")
lines(dates,Goldman[[5]][,2],type = "l")
axis.Date(side=1,at=seq(dates[1],dates[2765],by="year"),cex.axis=0.8,las=2)
legend("topleft",legend = c("Daily Returns","1% VaR"),col = 2:1,pch =c("o","-"))


#HSBC
plot(dates,as.vector(MVMQ[,3]),type = "p",ylim = c(-60,60),cex=0.6,xaxt="n",cex.axis=0.8,col=2,xlab="",ylab = "",main = "HSBC")
lines(dates,HSBC[[5]][,2],type = "l")
axis.Date(side=1,at=seq(dates[1],dates[2765],by="year"),cex.axis=0.8,las=2)
legend("topleft",legend = c("Daily Returns","1% VaR"),col = 2:1,pch = c("o","-"))

Of course, just as in the univariate case, the resulting object can also be plotted using plot()

plot(Barclays,rows=2,columns=1)

plot(Goldman,rows=2,columns=1)

References

Engle, R. F., & Manganelli, S. (2004). CAViaR. Journal of Business & Economic Statistics, 22(4), 367–381. https://doi.org/10.1198/073500104000000370
Huang, D., Yu, B., Fabozzi, F. J., & Fukushima, M. (2009). CAViaR-based forecast for oil price risk. Energy Economics, 31(4), 511–518. https://doi.org/10.1016/j.eneco.2008.12.006
Machado, J. A. F. (2013). Quantile regression and heteroskedasticity.
White, H., Kim, T.-H., & Manganelli, S. (2015). VAR for VaR: Measuring tail dependence using multivariate regression quantiles. Journal of Econometrics, 187(1), 169–188. https://doi.org/10.1016/j.jeconom.2015.02.004

  1. To the best of my knowledge↩︎

  2. Please also disregard the dates, as the data they have available does not include dates; therefore, dummy dates have simply been included to create the graph. Additionally, note that in the original graphs, they show the VaR for the in-sample and out-of-sample all together, whereas in the graphs presented here, the graph is plotted only within the estimation sample.↩︎