Unified Bayesian Conditional Autoregressive Risk Measures using the Skew Exponential Power Distribution

TL;DR

This paper introduces a unified Bayesian framework for risk measurement using the Skew Exponential Power distribution, incorporating non-linearity with P-splines and demonstrating effectiveness on stock market data.

Contribution

It proposes a novel unified Bayesian approach for risk measures with a flexible non-linear extension using P-splines, adapting MCMC for inference.

Findings

Effective risk measurement on stock data

Flexible modeling of non-linearity

Improved asymmetry handling

Abstract

Conditional Autoregressive Value-at-Risk and Conditional Autoregressive Expectile have become two popular approaches for direct measurement of market risk. Since their introduction several improvements both in the Bayesian and in the classical framework have been proposed to better account for asymmetry and local non-linearity. Here we propose a unified Bayesian Conditional Autoregressive Risk Measures approach by using the Skew Exponential Power distribution. Further, we extend the proposed models using a semiparametric P-spline approximation answering for a flexible way to consider the presence of non-linearity. To make the statistical inference we adapt the MCMC algorithm proposed in Bernardi et al. (2018) to our case. The effectiveness of the whole approach is demonstrated using real data on daily return of five stock market indices.

Tables

Table 1. Table 1: Summary statistics.

	Mean	Median	Std.	Min	Max	Q1	Q3	Skewness	Normality test	Kurtosis
Nasdaq	0.0488	0.1135	1.6542	-11.1149	17.2030	-0.6883	0.8280	0.0400	0.0010	9.3400
Hang Seng	0.0320	0.0586	1.5878	-24.5202	17.2471	-0.6616	0.7974	-0.5469	0.0010	19.5871
Corea SE	0.0182	0.0340	1.6457	-12.8047	11.2844	-0.7179	0.7682	-0.1386	0.0010	8.5365
AEX	0.0241	0.0668	1.2922	-9.5903	10.0283	-0.5532	0.6523	-0.1804	0.0010	9.7683
STI	0.0173	0.0225	1.1981	-10.5446	12.8738	-0.5255	0.5856	-0.0888	0.0010	12.0541

Table 2. Table 2: Summary statistics for CAViaR models. The columns, denoted by LR uc , LR cc and DQ, report the p–values of unconditional and conditional coverage tests of Kupiec Kupiec, ( 1995 ) and Christoffersen Christoffersen, ( 1998 ) and those of the Dynamic Quantile (DQ) test of Engle and Manganelli Engle and Manganelli, ( 2004 ) .

	$\tau =0.01$							$\tau =0.05$
Name	A/E	AD Viol.		Viol.	${\mathrm{LR}}_{uc}$	${\mathrm{LR}}_{cc}$	DQ	A/E	AD Viol.		Viol.	${\mathrm{LR}}_{uc}$	${\mathrm{LR}}_{cc}$	DQ
		Mean	Max						Mean	Max
NASDAQ
BNL-CAViaR	1.1689	7.4829	15.5215	38.0000	0.3450	0.1402	0.0000	1.0028	4.9176	14.6782	163.0000	0.9679	0.9554	0.0004
SAV	1.0151	7.2093	18.0225	33.0000	0.9299	0.1355	0.0000	1.0643	4.7129	14.6590	173.0000	0.4028	0.6791	0.0011
AS	0.9843	6.5240	16.6643	32.0000	0.9296	0.6133	0.4089	1.0643	4.7200	15.2859	173.0000	0.4028	0.7028	0.0090
T-CAViaR	1.4149	6.0632	16.4689	46.0000	0.0251	0.0319	0.0001	1.1566	4.6330	15.3398	188.0000	0.0450	0.1288	0.0001
IG	1.0151	6.9912	16.2478	33.0000	0.9299	0.1355	0.0000	1.0520	4.6630	15.1494	171.0000	0.4974	0.7479	0.0001
HANG SENG
BNL-CAViaR	0.7221	8.3603	23.8684	23.0000	0.0976	0.2145	0.3949	1.1177	5.0657	23.2981	178.0000	0.1332	0.1884	0.1971
SAV	0.6279	7.0596	17.5615	20.0000	0.0236	0.0680	0.1193	0.9294	5.3574	21.8372	148.0000	0.3570	0.4801	0.1216
AS	0.5965	6.5154	11.3844	19.0000	0.0134	0.0420	0.3085	1.1366	5.1932	24.4616	181.0000	0.0825	0.0349	0.0402
T-CAViaR	0.7849	8.5494	18.0606	25.0000	0.2056	0.1916	0.5401	1.1680	5.2449	25.9787	186.0000	0.0336	0.0045	0.0052
IG	0.7849	7.2136	17.7435	25.0000	0.2056	0.3682	0.2197	1.1115	5.1525	21.5238	177.0000	0.1547	0.2192	0.0692
KOREA SE
BNL-CAViaR	0.8459	9.0809	17.3114	27.0000	0.3696	0.3230	0.0000	1.0276	4.8097	16.3244	164.0000	0.7189	0.0436	0.0540
SAV	0.7519	8.5358	21.1337	24.0000	0.1411	0.2822	0.0922	1.0965	4.6451	17.1827	175.0000	0.2163	0.0013	0.0000
AS	0.8459	7.5243	17.5978	27.0000	0.3696	0.5310	0.0136	1.0965	4.6490	17.9097	175.0000	0.2163	0.0609	0.0831
T-CAViaR	0.8772	7.7273	21.7983	28.0000	0.4774	0.6064	0.0088	1.2907	4.1103	17.8722	206.0000	0.0003	0.0002	0.0000
IG	0.6266	9.6109	21.4799	20.0000	0.0229	0.0662	0.0035	1.0902	4.6719	17.6154	174.0000	0.2471	0.0072	0.0000
AEX
BNL-CAViaR	1.1202	6.6605	17.4447	37.0000	0.4947	0.5207	0.7106	1.0718	4.5112	15.4441	177.0000	0.3474	0.4552	0.6968
SAV	1.1807	6.6082	17.6856	39.0000	0.3092	0.3741	0.1622	1.1565	4.5291	15.6997	191.0000	0.0434	0.0157	0.0003
AS	1.0294	6.5699	18.0154	34.0000	0.8646	0.6919	0.8774	1.1262	4.5025	16.5794	186.0000	0.1016	0.0982	0.2513
T-CAViaR	1.1202	6.3613	17.7032	37.0000	0.4947	0.5207	0.8304	1.1444	4.5050	17.1307	189.0000	0.0619	0.0769	0.1426
IG	1.3321	6.4102	18.6286	44.0000	0.0676	0.1039	0.1099	1.1262	4.4897	16.4462	186.0000	0.1016	0.0982	0.0443
STI
BNL-CAViaR	1.0255	6.0364	16.9641	33.0000	0.8836	0.6402	0.5338	1.0876	3.9728	15.6347	175.0000	0.2588	0.1270	0.0291
SAV	1.0566	5.7375	12.4956	34.0000	0.7480	0.6605	0.7565	1.0193	3.8977	15.5282	164.0000	0.7995	0.9429	0.7876
AS	1.0255	5.7592	13.6595	33.0000	0.8836	0.6402	0.7410	1.0690	3.9577	16.9818	172.0000	0.3722	0.4894	0.1332
T-CAViaR	1.0876	5.7188	12.9416	35.0000	0.6210	0.6189	0.6659	1.0752	3.9313	17.0519	173.0000	0.3313	0.4417	0.1417
IG	1.0876	5.8194	13.6148	35.0000	0.6210	0.6022	0.1269	1.0379	3.9325	15.0095	167.0000	0.6209	0.6401	0.8215

Table 3. Table 3: Summary statistics for CARE models. The columns, denoted by LR uc , LR cc and DQ, report the p–values of unconditional and conditional coverage tests of Kupiec Kupiec, ( 1995 ) and Christoffersen Christoffersen, ( 1998 ) and those of the Dynamic Quantile (DQ) test of Engle and Manganelli Engle and Manganelli, ( 2004 ) .

	$\tau =0.01$							$\tau =0.05$
Name	A/E	AD Viol.		Viol.	${\mathrm{LR}}_{uc}$	${\mathrm{LR}}_{cc}$	DQ	A/E	AD Viol.		Viol.	${\mathrm{LR}}_{uc}$	${\mathrm{LR}}_{cc}$	DQ
		Mean	Max						Mean	Max
NASDAQ
BNL-CARE	0.7075	8.5479	16.6762	23.0000	0.0772	0.0764	0.0275	0.8920	5.3056	15.0000	145.0000	0.1517	0.3496	0.0002
SAV	0.8920	8.0395	17.9063	29.0000	0.5297	0.0703	0.0005	1.0889	4.8234	14.9587	177.0000	0.2497	0.5032	0.0012
AS	0.7075	8.0443	16.5115	23.0000	0.0772	0.1780	0.1135	1.0151	4.8593	15.5370	165.0000	0.8409	0.6495	0.0191
T-CARE	0.6152	7.9794	16.3227	20.0000	0.0177	0.0530	0.0437	1.0458	4.8962	15.8722	170.0000	0.5490	0.7937	0.0006
IG	0.7075	8.0717	16.3479	23.0000	0.0772	0.0764	0.0580	1.0581	4.8623	15.4656	172.0000	0.4486	0.7153	0.0000
HANG SENG
BNL-CARE	0.6279	10.0314	23.1171	20.0000	0.0236	0.0680	0.0116	1.0801	5.3312	22.1913	172.0000	0.3039	0.0725	0.0730
SAV	0.6593	6.8444	17.8665	21.0000	0.0396	0.1046	0.1401	0.9922	5.1728	23.5544	158.0000	0.9222	0.2856	0.0689
AS	0.5651	6.0058	10.1083	18.0000	0.0073	0.0246	0.1670	1.1303	5.1605	25.2195	180.0000	0.0973	0.0832	0.1786
T-CARE	0.5024	8.3517	12.6012	16.0000	0.0018	0.0071	0.2169	1.0612	5.3918	15.7235	169.0000	0.4299	0.1099	0.0259
IG	0.7849	7.2146	17.4942	25.0000	0.2056	0.3682	0.2180	1.0487	5.2958	22.4331	167.0000	0.5291	0.6650	0.5228
KOREA SE
BNL-CARE	0.6892	8.6710	19.1366	22.0000	0.0618	0.1500	0.1472	0.9837	4.9451	16.0867	157.0000	0.8355	0.0464	0.0507
SAV	0.7832	8.5001	20.7419	25.0000	0.2013	0.3628	0.0525	1.0276	4.8199	17.1703	164.0000	0.7189	0.0882	0.0178
AS	0.7519	7.8578	17.6468	24.0000	0.1411	0.2822	0.0075	1.0526	4.7782	18.0365	168.0000	0.4960	0.0552	0.0405
T-CARE	0.8145	7.7491	17.6432	26.0000	0.2773	0.4477	0.0153	1.2218	4.3582	17.7446	195.0000	0.0053	0.0068	0.0193
IG	0.7832	8.7727	21.1006	25.0000	0.2013	0.3628	0.0172	1.0526	4.8380	17.3711	168.0000	0.4960	0.0268	0.0002
AEX
BNL-CARE	1.0294	7.5506	17.3401	34.0000	0.8646	0.6919	0.1884	1.1262	4.6867	17.7054	186.0000	0.1016	0.1425	0.0469
SAV	1.2716	7.8593	16.8078	42.0000	0.1318	0.1870	0.0000	1.2231	4.5865	14.9509	202.0000	0.0044	0.0047	0.0001
AS	1.1202	7.6155	16.0772	37.0000	0.4947	0.5207	0.0095	1.1505	4.6591	15.7784	190.0000	0.0519	0.1227	0.0066
T-CARE	1.2110	7.8290	17.3910	40.0000	0.2374	0.3046	0.0000	1.2171	4.5207	15.4260	201.0000	0.0055	0.0185	0.0016
IG	1.2716	7.5170	16.7142	42.0000	0.1318	0.2719	0.0003	1.1565	4.5903	15.7695	191.0000	0.0434	0.0268	0.0085
STI
BNL-CARE	0.9944	7.1228	16.9036	32.0000	0.9759	0.6198	0.0001	1.2927	3.7855	14.9127	208.0000	0.0003	0.0011	0.0001
SAV	0.8701	7.4723	16.3825	28.0000	0.4499	0.3822	0.0000	0.9323	4.2518	14.3585	150.0000	0.3749	0.6746	0.7608
AS	0.9012	7.2291	15.0321	29.0000	0.5678	0.4559	0.0009	1.0068	4.3077	15.3054	162.0000	0.9260	0.9940	0.0005
T-CARE	0.8080	7.4438	17.3445	26.0000	0.2581	0.2383	0.0028	0.9944	4.3194	15.4082	160.0000	0.9451	0.9975	0.0012
IG	0.9012	7.0910	16.3101	29.0000	0.5678	0.6525	0.0308	0.9447	4.2920	14.0801	152.0000	0.4701	0.4451	0.2950

Table 4. Table 4: Bootstrap test p-values for zero mean of the standardardized residuals. Test based on 1000 post sample estimates of the conditional 1% and 5% ES.

	$\tau =0.01$					$\tau =0.05$
Name	BNL-CARE	SAV	AS	T-CARE	IG	BNL-CARE	SAV	AS	T-CARE	IG
Nasdaq	0.5260	0.9400	0.8590	0.2220	0.7300	0.0760	0.1750	0.5370	0.6780	0.3840
Hang Seng	0.0000	0.0000	0.0000	0.8340	0.0000	0.0000	0.0000	0.0000	0.0200	0.0010
Corea SE	0.1310	0.0080	0.0080	0.0440	0.0150	0.1090	0.0000	0.0000	0.0080	0.0050
AEX	0.0080	0.0000	0.0020	0.0010	0.0000	0.0060	0.0490	0.0050	0.0090	0.4720
STI	0.0000	0.0000	0.0000	0.0000	0.0000	0.0130	0.0010	0.0000	0.0000	0.0050

Equations

$$y_t$$

$$q_{t,\tau }(\theta )$$

$$\displaystyle\ell\left(\mbox{\boldmath$\beta$},y_{t-1}\right)$$

$$\displaystyle\ell\left(\mbox{\boldmath$\beta$}y_{t-1}\right)$$

$$q_{t,\tau}(\theta)=\left(\omega+\gamma q^{2}_{t-1,\tau}(\theta)+\ell\left(\mbox{\boldmath$\beta$},y_{t-1}\right)\right)^{1/2}$$

$$\ell\left(\mbox{\boldmath$\beta$},y_{t-1}\right)=\beta y^{2}_{t-1}$$

$$\underset{\theta }{\min }\frac{1}{T}\sum _{i=1}^T{\rho }_{\tau }(y_t-q_{t,\tau }(\theta ))$$

$$y_t$$

$${\mu }_{t,\nu }(\theta )$$

$$\underset{\theta }{\min }\frac{1}{T}\sum _{i=1}^T{r}_{\tau }(y_t-{\mu }_{t,\tau }(\theta ))$$

$${\mathrm{ES}}_t\left(\tau \right)=\left(1+\frac{\nu }{\left(1-2\nu \right)\tau }\right){\mu }_{t,\nu }(\theta )-\frac{\nu }{\left(1-2\nu \right)\tau }\mathbb{E}\left(y_t\right).$$

$$f_{\mathsf{S}\mathsf{E}\mathsf{P}}\left(y_{t};g_{t},\sigma,\tau,\alpha\right)=c^{-1}\times\left\{\begin{array}[]{lc}\exp\left\{-\left(1-\tau\right)\left(\frac{g_{t}-y_{t}}{\sigma}\right)^{\alpha}\right\},&{\rm if}\quad y_{t}<g_{t}\\ \exp\left\{-\tau\left(\frac{y_{t}-g_{t}}{\sigma}\right)^{\alpha}\right\},&{\rm if}\quad y_{t}\geq g_{t},\end{array}\right.$$

$$y_t$$

$$g_{t,\nu }(\theta )$$

$$\displaystyle\ell\left(\mbox{\boldmath$\beta$},y_{t-1}\right)$$

$$\displaystyle\ell\left(\mbox{\boldmath$\beta$}y_{t-1}\right)$$

$$g_{t,\tau}(\theta)=\left(\omega+\gamma g^{2}_{t,\tau}(\theta)+\ell\left(\mbox{\boldmath$\beta$},y_{t-1}\right)\right)^{1/2}$$

$$\ell\left(\mbox{\boldmath$\beta$},y_{t-1}\right)=\beta y^{2}_{t-1}$$

$$\ell\left(\mbox{\boldmath$\beta$},y_{t-1}\right)=\sum_{\nu=1}^{k+d}\beta_{\nu}B_{\nu}\left(y_{t-1}\right).$$

$${\beta }_{\nu }=2{\beta }_{\nu -1}-{\beta }_{\nu -2}+u_{\nu },\forall \nu =1,2,\dots ,k+d.$$

$$\pi\left(\mbox{\boldmath$\beta$}\mid\phi\right)\propto\mathcal{N}_{k+d}\left(0,\phi^{2}\left(D_{2}^{\prime}D_{2}\right)^{-1}\right),$$

$${\mathcal{L}}_{\tau }(\theta ,\sigma$$

$$\prod _{t=1}^T\left[\exp \left\{-\left(1-\tau \right){\left(\frac{g_t-y_t}{\sigma }\right)}^{\alpha }\right\}\text{mathbbm}{1}_{\left(-\infty ,g_t\left(\theta \right)\right)}+\exp \left\{-\tau {\left(\frac{y_t-g_t}{\sigma }\right)}^{\alpha }\right\}\text{mathbbm}{1}_{\left(g_t\left(\theta \right),\infty \right)}\right],$$

$$\pi \left(\mathbf{\Xi }\right)=\pi \left(\mathbf{\beta }\mid {\varphi }^2\right)\pi \left({\varphi }^2\right)\pi \left(\omega \right)\pi \left(\gamma \right)\pi \left(\sigma \right),$$

$$\displaystyle\pi\left(\mbox{\boldmath$\beta$}\mid\phi^{2}\right)$$

$$\pi \left({\varphi }^2\right)$$

$$\pi \left(\gamma \right)$$

$$\pi \left(\omega \right)$$

$$\pi \left(\sigma \right)$$

$$\pi \left(\mathbf{\Xi }\mid \mathbf{y}\right)\propto {\mathcal{L}}_{\tau }\left(\theta ,\sigma \mid \alpha ,\mathbf{y}\right)\pi \left(\mathbf{\beta }\mid {\varphi }^2\right)\pi \left({\varphi }^2\right)\pi \left(\omega \right)\pi \left(\gamma \right)\pi \left(\sigma \right),$$

$$\displaystyle q\left(\mbox{\boldmath$\beta$}_{i-1},\mbox{\boldmath$\beta$}^{*}_{i}\right)$$

$$q\left({\omega }_{i-1},{\omega }_i^{\ast }\right)$$

$$q\left({\gamma }_{i-1},{\gamma }_i^{\ast }\right)$$

$$q\left({\sigma }_{i-1},{\sigma }_i^{\ast }\right)$$

$$\lambda \left({\xi }_j^{\left(i-1\right)},{\xi }_j^{\ast }\right)=\min \left\{1,\frac{\mathcal{L}\left({\xi }_j^{\ast },{\mathbf{\Xi }}_{-j}^{\left(i-1\right)}\mid \mathbf{y},\mathbf{x}\right)}{\mathcal{L}\left({\mathbf{\Xi }}^{\left(i-1\right)}\mid \mathbf{y},\mathbf{x}\right)}\frac{\pi \left({\xi }_j^{\ast }\right)}{\pi \left({\xi }_j^{\left(i-1\right)}\right)}\frac{q\left({\xi }_j^{\left(i-1\right)}\right)}{q\left({\xi }_j^{\ast }\right)}\right\},$$

$$\displaystyle\mbox{\boldmath$\mu$}_{\beta}^{(i+1)}$$

Full text

1 Introduction

After the recent financial crisis an accurate risk measurement is a primary need for financial institutions and investors. Within the instruments for market risk measurement, Value–at–Risk (VaR) (Jorion, (2007)) and Expected Shortfall (ES) (Artzner et al., (1999)) are certainly the most popular and used approaches. VaR answers the question on what is the maximum potential loss that will be exceeded with a certain probability in the next days. It can be simply understood as a specific (say $\tau$) conditional quantile of the portfolio returns given the current information, i.e., $\mathbb{P}\left(Y_{t}<-{\rm VaR}_{t}\mid\mathcal{F}_{t}\right)=\tau$, where $Y_{t}$ and $\mathcal{F}_{t}$ denote the return of a portfolio and the information set available at time $t$, respectively, while $\tau\mathcal{2}\left(0,1\right)$ denotes the quantile confidence level associated with the VaR. Even though it is widely used among financial institutions VaR has been criticized because of the absence of the sub–additivity property, namely, it does not guarantee that a diversified portfolio is less risky than a concentrated one. In addition VaR gives no information regarding possible exceedances beyond the quantile which may be quite important in evaluating risks. Artzner et al., (1999) first recognized the lack of coherency of the VaR and proposed the ES as an alternative coherent risk measure which gives more information about the returns’ distribution in the tails. In particular the ES is defined as the conditional expected loss given that the loss exceed the VaR, i.e., $\mathbb{E}\left(Y_{t}\mid y_{t}<-{\rm VaR}\right)$. Despite ES is a coherent risk measure, it is in general more difficult to calculate and to backtest it, i.e., to verify how accurately the strategy or method would have predicted actual results. Its modeling and measuring is still an ongoing problem without reaching a consensus. Moreover, it is not as simple to interpret as the VaR, and, for these reasons there is not a prevailing risk measure between VaR and ES.

In financial time series it is usually the case that the distribution of the returns typically changes over time, for this reason the recent literature focuses on modeling VaR and the ES by considering conditional autoregressive models. A recent development in the VaR literature is the Conditional Autoregressive Value at Risk (CAViaR) class of models (see Engle and Manganelli, (2004)). In particular it specifies the evolution of the quantile over time using an autoregressive approach, and estimates the parameters using a quantiles regression framework. This approach has strong appeal in that it focuses the tail of return distribution directly and does not rely on any distributional assumption.

In order to impose a dynamic evolution over ES risk measure, Taylor, (2008) elaborates the Conditional Autoregressive Expectile (CARE) class of models which delivers estimates for both VaR and ES by considering the theory of expectiles. In practice the expectile is used as the basis for estimating ES, obtained applying the equation (12) of Section 2. Moreover, the one–to–one mapping from expectiles to quantiles can be used to obtain also an estimate of the VaR by following Efron, (1991) who showed that the estimator of the $\tau$–th quantile is the corresponding $\nu$–th expectile having the proportion of observations below equal to $\tau\%$.

The inferential approach proposed in the literature to estimate the parameters of the two class of models consider both the frequentist and the Bayesian tools. In the CaViaR framework the former approach is based on quantile regression methods (Koenker, (2005), Engle and Manganelli, (2004)) consisting in minimizing the loss function introduced by Koenker and Basset, (1978). The Bayesian approach instead relies on the Asymmetric Laplace (AL) distribution tool; see for example Yu and Moyeed, (2001), Kottas and Gelfand, (2001), Kottas and Krnjajic, (2009), Sriram et al., (2013) and Bernardi et al., (2015). From a CARE model point of view the frequentist approach relies on the Asymmetric Least Squares (ALS) as in Newey and Powell, (1987) while the Bayesian paradigm consider the Asymmetric Gaussian (AG) distribution assumption (see, e.g., Gerlach and Wang, (2015), Gerlach and Chen, (2014), Wichitaksorn et al., (2014), Gerlach et al., (2016), Gerlach and Chen, (2016)).

With this paper we introduce three main innovations on the existing literature on CaViaR and CARE models. First, we develop a unified Bayesian Conditional Autoregressive Risk Measure (B–CARM) model which encompasses both the CAViaR and the CARE ones as particular cases by using the Skew Exponential Power (SEP) as working likelihood. Using the properties and the parametrization of the SEP presented in Kobayashi, (2015), we show how to estimate the CAViaR and the CARE class of Conditional Autoregressive Risk models by varying one of the parameter of the SEP distribution. A parallel idea is developed in Kobayashi, (2015) where quantiles and expectiles are calculated for a stochastic volatility model when the SEP distribution is assumed. Differently from Kobayashi, (2015) our paper uses the properties of the SEP distribution to estimate of Conditional Autoregressive models in the risk measures framework.

The second new issue of the paper consists in proposing a new nonlinear and semiparametric model specification (the BNL–CARM) of B–CARM one by using the penalized B–Splines approximation for the News Impact Curve (NIC) (see, Boor, (2001), Eilers et al., (1996), Lang and Brezger, (2004)) to estimate the relation between the quantile/expectile and the exogenous variables. The need for a model that allows for nonlinearity without assuming any particular shape restrictions is of clear interest in financial literature. Most of the existing literature on CaViaR and CARE imposes a priori the nonlinearity relation between the observed variables and the current quantile or expectile level (see i.e. Engle and Manganelli, (2004), Gerlach et al., (2011) and Gerlach et al., (2012), Chen et al., (2009) and Chen et al., (2012), Gerlach and Chen, (2014)).

Finally, a new Adaptive–Independent Metropolis–Hastings (AIMH) algorithm is implemented to estimate efficiently the model parameters both in the parametric and semiparametric versions. In particular modify and readapt the method proposed in Bernardi et al., (2018) to our case. The algorithm is in line with the Adaptive MCMC methods proposed in Faming Liang, (2010) that allows the proposal distribution to be updated at each iteration to tailor its shape to that of the target one. These methods do not require the prior specification of the proposal parameters and their theoretical properties are now well understood, see, e.g., Andrieu and Thoms, (2008) and Faming Liang, (2010), Atchad’e and Rosenthal, (2005), Atchad’e et al., (2011).

From an applied point of view we estimate the BNL–CARM and the B–CARM models for five stock market indices comparing their performances. In particular we conduct a backtesting procedure showing that, even without imposing any restrictive assumption on the relations among variables, the proposed BNL–CARM performances are in line with those of the competitors. By that we can consider the BNL–CARM as a valid and more general alternative to the existing models.

The remaining of the paper is organized as follows. Section 2 briefly reviews conditional autoregressive risk models; Section 3 introduces the SEP as working likelihood for the Bayesian Conditional Autoregressive Risk Model with the CAViaR and the CARE models being special cases, while Section 4 proposes a nonlinear and semiparametric extension. Section 5 shows and discusses the implemented Bayesian methododology; Section 6 shows results from real datasets and Section 7 provides concluding remarks.

2 CAViaR and CARE models

The most important VaR and ES dynamic models proposed in literature are the class of conditional autoregressive risk measure models known as the CAViaR and the CARE models introduced by Engle and Manganelli, (2004) and Taylor, (2008), respectively. The CAViaR class of models attempt to compute the $\tau$–th level VaR by estimating the $\tau$–th level quantile of the portfolio returns through a conditional autoregressive equations structure. More specifically, let $y_{t}$ be the return at time $t$, the CAViaR model has the following form:

[TABLE]

where $q_{t,\tau}(\theta)$ is the $\tau$ level quantile of $y_{t}$ conditional upon the information set up to time $t-1$, $F_{t-1}$. It is defined as the value that minimizes the function $\mathbb{E}\left[\left|\tau-\mathbbm{1}_{\left(-\infty,q_{t,\tau}(\theta)\right)}\left(y_{t}\right)\right|\left|y_{t}-q_{t,\tau}(\theta)\right|\right]$, where $\mathbbm{1}_{A}$ is the indicator function of the set $A$, \theta=\left(\omega,\gamma,\mbox{\boldmath\beta}\right) are the parameters of the model and $l\left(\cdot\right)$ is an unknown function of the past returns. Here, $\varepsilon_{t}$, for any $t=1,2,\dots,T$, are independent random variables which are supposed to have zero $\tau$–th quantile and constant variance. As noted by Engle and Manganelli, (2004), $\ell\left(\cdot\right)$ can be interpreted as the News Impact Curve (NIC) introduced by Engle and Ng, (1993) for ARCH–type models. The form of the function $\ell\left(\cdot\right)$ is one of the most addressed tasks in the risk modeling literature. Indeed, we can recognize different CAViaR models by considering different form of $\ell\left(\cdot\right)$ such as:

[TABLE]

where $z_{t}$ is an observed threshold variable that could be exogenous or $y_{t}$ itself, i.e., $z_{t}=y_{t}$; $r$ is the threshold value, typically set equal to zero, i.e., $r=0$ (see Gerlach et al., (2011)). Another possible configuration of the CAViaR model is the Indirect Garch(1,1) where in equation (2) we have

[TABLE]

and

[TABLE]

From an inferential point of view, the frequentist approach estimate CAViaR through quantile regression methods (Koenker, (2005), Engle and Manganelli, (2004)) by Koenker and Basset, (1978) i.e. solving the problem:

[TABLE]

where $\rho_{\tau}(u)\equiv u\left(\tau-\mathbbm{1}_{(u<0)}\right)$ is the well known quantile check function.

The Bayesian approach, instead, relies on the AL distribution assumption as inferential tool to perform the statistical analysis (see, e.g., Yu and Moyeed, (2001), Kottas and Gelfand, (2001), Kottas and Krnjajic, (2009), Sriram et al., (2013), Bernardi et al., (2015)).

The same structure considered for the CAViaR models is used to build the CARE model defined as:

[TABLE]

where $\mu_{t,\nu}(\theta)$ is the $\nu$–th expectile of $y_{t}$ defined as the value that minimizes the function $\mathbb{E}\left[\left|\nu-\mathbbm{1}_{\left(-\infty,\mu_{t,\nu}(\theta)\right)}\left(y_{t}\right)\right|\left(y_{t}-\mu_{t,\nu}(\theta)\right)^{2}\right]$. In this model $\varepsilon_{t}$, for any $t=1,2,\dots,T$, are independent random variables having zero $\nu$–th expectile and constant variance. The specifications of $\ell\left(\cdot\right)$ in equations (3)–(5) remain valid to define different CARE models with the Indirect GARCH(1,1) CARE obtained by substituting in (6) $q_{t,\nu}(\theta)$ with $\mu_{t,\nu}(\theta)$. The estimation procedure for the generic expectile is addressed in the frequentist approach by using the Asymmetric Least Square (ALS) estimator as in Newey and Powell, (1987), i.e. by minimizing

[TABLE]

where $r_{\tau}(u)\equiv u^{2}\left(\tau-\mathbbm{1}_{(u<0)}\right)$ is the expectile check function.

In the Bayesian paradigm the literature relies on the use of the AG distribution as working likelihood in order to make inference (see, e.g., Gerlach and Wang, (2015), Gerlach and Chen, (2014), Wichitaksorn et al., (2014), Gerlach et al., (2016), Gerlach and Chen, (2016)).

Finally, since the expectile does not give immediately information on quantile and ES another passage is required to obtain VaR and ES from the estimated expectile. In the first case we simply use the expectile as an estimator of the quantile by iteratively searching for the $\nu$–th expectile for which we observe $\tau\%$ observations in–sample below it. This procedure was suggested by Efron, (1991) and allows us to obtain the quantile of interest and consequently the associated VaR level. The ES is instead obtained using the one to one mapping between expectile and ES suggested in Taylor, (2008) given by:

[TABLE]

In the next Section we provide a unified risk measure framework structure which encompasse the CAViaR and CARE ones.

3 Bayesian Conditional Autoregressive Risk Measures

In this Section we develop a unified Conditional Autoregressive Risk model which encompasses both the CAViaR and the CARE as particular cases by using the Skew Exponential Power (SEP) proposed by Kobayashi, (2015). In particular we show how one of the parameter of the SEP distribution is responsable for the choise of the autoregressive risk model considered. We approach this problem from a Bayesian point of view following Bernardi et al., (2018), producing a unified framework that we call the Bayesian Conditional Autoregressive Risk Measure (B–CARM) class of models. For this reason let’s recall the SEP desity function:

[TABLE]

where $g_{t}$ is the location parameter and $\tau\in\left(0,1\right)$ is the skewness parameter. In addiction, $\sigma\in\Re^{+}$ and $\alpha\in\left(0,\infty\right)$ are the scale and shape parameters, respectively, while $c=\sigma\Gamma\left(1+\frac{1}{\alpha}\right)\left(\frac{1}{\tau^{\frac{1}{\alpha}}}+\frac{1}{\left(1-\tau\right)^{\frac{1}{\alpha}}}\right)$ and $\Gamma\left(\cdot\right)$ is the complete gamma function. It is easy to check that the SEP distribution encompass the AL and the AG distributions by setting $\alpha=1$ and $\alpha=2$, respectively. In this way by considering the generic model specification

[TABLE]

with

[TABLE]

and

[TABLE]

where

[TABLE]

and using the SEP distribution as working likelihood we retrieve the CAViaR models by fixing $\alpha=1$ and the CARE ones by fixing $\alpha=2$. Indeed if we use $\alpha=1$ in (13) than $g_{t}$ becomes the quantile, while with $\alpha=2$ the $g_{t}$ becomes the expectile. In this way we are able to build a unified class of autoregressive risk models by using the SEP distribution. As said before, since the expectile does not represent a risk measure itself, when $\alpha=2$ we estimate the VaR by iteratively searching for the expectile for which we observe a given percentage of observations below it, while the ES is obtained using the one to one mapping between expectile and ES given in equation (12).

4 Nonlinear CARM

To keep as general as possible the model specified in equations (14)–(20) we present a new Bayesian Non–linear CARM (BNL–CARM) model where in (14)–(15) we set up a nonlinear and semiparametric framework by using a spline approach to model the function $\ell\left(\cdot\right)$. In this way we avoid to impose an apriori functional form for $\ell\left(\cdot\right)$ like the ones considered in the previous sections allowing the data to decide about its shape. In particular we use the B–spline functions of order $d$ with $k$ knots as:

[TABLE]

where $B_{\nu}\left(y_{t-1}\right)$ denote the B–spline basis functions and $\beta_{\nu}$ are unknown coefficients to be estimated. As it is well known in equation (21) the value of the estimated coefficients and the shape of the fitted function depend upon the number and the position of the knots. In this paper, in absence of specific prior information we assume the equidistance of the knots. Moreover in order to capture the smoothness of the data it is important to carefully choose the number of knots to take in two account of the trade–off between too few and too many knots, which may cause underfitting or overfitting respectively. A possible solution of this problem is known as Penalized Spline (P–Spline) and proposed by (O’Sullivan, (1986) and O’Sullivan, (1988)) and generalized by Eilers et al., (1996). To ensure enough flexibility without incurring in the overfitting problem we use a relatively large number of knots jointly with a penalization term able to smooth sufficiently the fitted curve following Eilers et al., (1996). In more details we consider a penalty element on the second differences of the B–Spline coefficients which can be embedded in a Bayesian framework by using a second order random walk type prior for all the B–Spline coefficients following Lang and Brezger, (2004), Brezger and Lang, (2006), Brezger and Steiner, (2008) and Bernardi et al., (2018) i.e:

[TABLE]

Here the generic stochastic component $u_{\nu}$ has a Gaussian distribution with zero mean and variance equal to $\phi^{2}$, i.e. $u_{\nu}\sim\mathcal{N}\left(0,\phi^{2}\right)$ and $\beta_{\nu-1}$ and $\beta_{\nu-2}$ are initialized with diffuse priors (i.e. $\propto 1$). The smoothness of the fitted curve is controlled by the variance of the error term which correspond to the inverse of the penalization parameter used by Eilers et al., (1996) in the frequentist framework. We choose a conjugate Inverse Gamma prior for $\phi^{2}$, that is $\phi^{2}\sim\mathcal{IG}\left(a^{\left(\phi\right)},b^{\left(\phi\right)}\right)$ with $a^{\left(\phi\right)}=b^{\left(\phi\right)}=0.001$. Different choices of hyper parameters are allowed but they all bring to very similar results. Finally, it is possible to write the prior distribution for the B–Spline coefficients as

[TABLE]

where \mbox{\boldmath\beta}=\left(\beta_{1},\beta_{2}\dots,\beta_{k+d}\right)^{\prime}, $D_{2}$ is the difference matrix of dimension $\left(k+d-2\right)\times\left(k+d\right)$ i.e. the differential order of the random walk in equation (22) (see also Bernardi et al., (2018))

5 Bayesian inference

Bayesian inference requires the specification of the likelihood function as well as the prior distribution for all the parameters of interest. The likelihood function of the model, based on the SEP distribution showed in Section 3 is given by:

[TABLE]

where the vector $\mathbf{y}=(y_{1},\ldots,y_{T})$ is the sample of observations. The parameter $\alpha$ is fixed equal to $\alpha=1$ or $\alpha=2$, depending on the model we want to estimate, i.e. the quantile model or the expectile one, as showed in previous section. Concerning the prior specification, we assume the following hierarchical prior structure independent on the value of $\tau$:

[TABLE]

with

[TABLE]

where $\bm{\Xi}=\left(\theta,\phi^{2},\sigma\right)$, $\left(a^{\left(\phi\right)},b^{\left(\phi\right)},\sigma_{\gamma}^{2},\sigma_{\omega}^{2},a,b\right)$ are given positive hyperparameters, while $\mathcal{N}\left(\cdot\right)$ and $\mathcal{IG}\left(\cdot\right)$ denote the Normal and the Inverse Gamma distributions respectively.

5.1 The Adaptive Independent Metropolis within Gibbs sampler

The Bayesian inference is carried out using an adaptive MCMC sampling scheme, similar to the one proposed by Bernardi et al., (2018), based on the following posterior distribution

[TABLE]

where $\mathcal{L_{\tau}}\left(\theta,\sigma\mid\alpha,\mathbf{y}\right)$ indicates the likelihood function specified in equation (5).

After choosing a set of initial values for the parameter vector $\mathbf{\Xi}^{(0)}$, the block–move Independent Metropolis within Gibbs (IMG) proceeds by iteratively simulate candidate values from a given set of proposal distributions which are subsequently accepted o rejected according to the usual Metropolis–Hastings acceptance rule, to preserve the detailed balance condition. Therefore, at the $i$–th iteration, for $i=1,2,\dots$ the simulation algorithm requires a proposal distribution for the parameters $\left(\theta,\sigma\right)$. Specifically, we consider the following set of proposal distributions

[TABLE]

where the scale parameter $\tilde{\sigma}=\log\left(\sigma\right)$ is transformed on the logarithmic scale and subsequently transformed back to preserve positiveness. The Jacobian term in equation (35) is then required to account for the distribution of the deterministic logarithmic transformation of $\sigma$.

At the $i$–th iteration, the IMG algorithm draws a candidate parameter $\Upsilon^{*}=\left(\xi_{1}^{*},\xi_{2}^{*},\xi_{3}^{*},\xi_{4}^{*}\right)=\left(\bm{\beta}^{*},\omega^{*},\gamma^{*},\sigma^{*}\right)$ from (32)–(35), which is subsequently accepted or rejected, with probability

[TABLE]

for $j=1,2,3,4$, where $\lambda\left(\xi_{j}^{\left(i-1\right)},\xi^{*}_{j}\right)$denotes the probability of moving to the new state of the chain, $\pi\left(\cdot\right)$ is the prior given in equations (26)–(30) and $\mathbf{\Xi}_{-j}^{\left(i-1\right)}$ refers to the whole set of parameters at iteration $i-1$ without the $j$–th element of $\Upsilon^{*}$. The last step of the algorithm update $\phi^{2}_{j}$, for $j=1,2,\dots,J$ by simulating directly from the respective full conditional distributions which are proportional to an Inverse Gamma (Brezger and Steiner, (2008)) with shape $a^{\left(\phi\right)}_{i}=a^{\left(\phi\right)}_{i-1}+rank\left(D_{2}^{\prime}D_{2}\right)/2$ and scale $b^{\left(\phi\right)}_{i}=b^{\left(\phi\right)}_{i-1}+\beta^{\prime}D_{2}^{\prime}D_{2}\beta$. Since most of the statistical properties of the Markov chain, as well as the performance of the Monte Carlo estimators, depend crucially on the definition of the proposal distribution $q\left(\cdot\right)$ (see, e.g., Andrieu and Moulines, (2006) and Andrieu and Thoms, (2008)), we improve the basic IMG–MCMC algorithm with an additional tuning step that adapts the proposal parameters to the target using the following equations:

[TABLE]

where $\varsigma^{(i+1)}$ denotes a tuning parameter that should be carefully selected in order to ensure the convergence and the ergodicity of the resulting chain (see Andrieu and Moulines, (2006)). Roberts and Rosenthal, (2007) provide two conditions for the convergence of the chain: the diminishing adaptation condition, which is satisfied if and only if $\varsigma^{(i)}\longrightarrow 0$, as $i\rightarrow+\infty$, and the bounded convergence condition, which guarantees that all the considered transition kernels have bounded convergence time. Moreover, Andrieu and Moulines, (2006) show that both conditions are satisfied if and only if $\varsigma^{(i)}\propto i^{-d}$ where $d\in\left[0.5,1\right]$. Therefore, we choose $\varsigma^{\left(i\right)}=\frac{1}{Ci^{0.5}}$ where $C=10$. As argued by Roberts and Rosenthal, (2007), these two conditions together ensure the asymptotic convergence of the adaptive algorithm and the existence of a weak law of large numbers, respectively.

6 Empirical applications

To show the usefulness and the simplicity of the unified risk model proposed and to asses its performance on real data we treat five stock market indices: Nasdaq (US); Straits Times Index (STI, Singapore); Hang Seng Index (Hong Kong); Corea SE (Corea) and AEX (Holland). Daily closing stock prices ($P_{t}$) from January 1, 1988 to November 30, 2018 are obtained from the Thomson Reuters (Datastream) database. The percentage returns are computed as $y_{t}=(log(P_{t})-log(P_{t-1}))\times 100$ and the full sample is divided into an in-sample period (from January 1, 1988 to December 30, 2005) and a forecast (out-of-sample) period of 3250 observations (from January 1, 2006 to November 30, 2018), which covers both the 2008-2009 General Financial Crisis and the following Sovereign Debt Crisis. Summary statistics are displayed in Table 1. All the series present the typical stylized facts of financial returns such as positive kurtosis and negative skewness. Moreover, the p-value of the Jarque-Bera test (9-th column) always reject the null hypothesis of normality. In the next two subsections we evaluate the forecast ability of the BNL-CAViaR and of the BNL-CARE models i.e. the BNL-CARM when $\alpha=1$ and $\alpha=2$ respectively, for all the series considered and we compare them with the regular CAViaR and CARE specifications. Tables 2 and 3 show the results of the out-of-sample estimation exercise conduced on all the series specified above.

6.1 CAViaR forecast evaluation

As mentioned before the CAViaR models allow to obtain a dynamic VaR estimation. In this section we study the behaviour of the BNL-CAViaR model and of all the CAViaR specification considered in Section 3 comparing the accuracy of the VAR estimation. In particular we label the CAViaR specifications as SAV for the Symmetric Absolute Value, AS for the Asymmetric Slope, T-CAViaR for the Threshold CAViaR and IG for the Indirect Garch.

The VaR estimation accuracy is evaluated by using the most common backtesting methods based on the comparison between the actual returns with their 1 day ahead VaR forecast.

In Table 2 we show the results of those tests applied to all the examined series for two $\tau$ levels, i.e. $\tau=0.01$ and $\tau=0.05$. The first column reports the ratio between the actual and the expected number of violations for a given coverage level $\tau$, i.e. $A/E=\frac{1}{m\tau}\sum_{t=n+1}^{n+m}I\left(y_{t}<-VaR_{t}\right)$, where $m=2456$ is the length of the forecast window. The last three columns report the p-values for common backtesting methods: the unconditional (LRuc), the conditional (LRcc) coverage tests of Kupiec, (1995) and Christoffersen, (1998) and the CAViaR Dynamic Quantile (DQ) test of Engle and Manganelli, (2004) respectively. The first two are likelihood ratio tests based on the assumption that the hit variable $I_{t}\left(\tau\right)$ defined as

[TABLE]

has a Bernoulli distribution with probability $\tau$. Under this hypothesis, the LRuc test verifies that the violation probability is equal to the coverage rate, i.e. $P\left(I_{t}\left(\tau\right)\right)=E\left(I_{t}\left(\tau\right)\right)=\tau$, while the LRcc test, in addition to LRuc test, also examines the independence hypothesis between violations observed at two different dates. The DQ test is instead a regression type test based on the demeaned process associated to $I_{t}\left(\tau\right)$, namely $H_{t}\left(\tau\right)=I_{t}\left(\tau\right)-\tau$. This test uses a regression model to assess the hypothesis of a linear relation between $H_{t}$, its lagged values and other relevant regressors and is known to be more powerful than the LRuc and the LRcc tests.

6.2 CARE forecast evaluation

While CAViaR models provide an estimate of the $\tau$ level quantile i.e. the risk measure $\tau$-th VaR, CARE models provide an estimate of the $\nu-th$ expectile, which is not a risk measure itself, but can be used to obtain estimates of both the VaR and the ES. Specifically, as anticipated in Section 2, from the different CARE specification we derive the $\tau$ level VaR by searching the $\nu-th$ expectile for which the proportion of observations below is $\tau\%$. The value of $\nu$ that satisfies this condition is obtained from the iterative estimation of the CARE models on a grid of step $0.001$ of $\nu$. The $\nu-th$ expectile founded by this recursive searching mechanism, i.e. $\mu_{t,\nu}\left(\theta\right)$, represents an estimate of the $\tau-th$ VaR. Therefore results can be evaluated using the same techniques showed in the previous subsection, i.e. LRcc, LRuc and DQ test. The output of these tests, using the specification presented in Section 2, are reported in table 3.

Moreover, CARE models can be used to deliver estimate of the ES through equation (12) that allows to map the $\nu$-th conditional expectile, $\mu_{t,\nu}\left(\theta\right)$, to the $\tau$-th expected shortfall, $ES_{t}\left(\tau\right)$. ES results are more difficult to evaluate and optimal assessment for ES forecasts is still an issue under investigation. Here we follow the approach of McNeil et al., (2005) and Taylor, (2008) based on a direct test of the residuals, i.e. the difference between the observations and the ES level for only those observations beyond the quantile VaR prediction. The test assess whether the residuals, standardized by the conditional volatility (the conditional quantile estimate in our case), are i.i.d with zero mean. Specifically, a bootstrap test is implemented (as in Efron and Tibshirani, (1993) page 224) in order to avoid distributional assumptions. The results of this tests are showed in Table 4.

6.3 Summary of VaR and ES results

By analyzing the empirical results given in tables 2 - 4 it is clear that both the BNL-CAViaR and BNL-CARE models show a forecasting performance in line with the other competitor models considered. This can be deducted by looking at the backtesting results reveling that no model outperforms the others in terms of A/E rate, number of violations, or p-values of the LRuc, LRuc and DQ tests. By following the traffic light approach suggested by the Basel Committee, (1996), we can say that all the models considered are classified as acceptable (green zone) or at least disputable (yellow zone) and that they never appear to be seriously flawed (red zone). It is worth noting that the time lenght of the series considered covers both the 2008-2009 General Financial Crisis and the following Sovereign Debt Crisis, events that shocked the market by creating large unexpected loss making the forecast exercise more complex. Despite this situation both BNL-CAViaR and the BNL-CARE perform quite well.

In figures 1 and 4 we report the posterior estimate of the NICs from the BNL-CAViaR (BNL-CARM with $\alpha=1$) and BNL-CARE (BNL-CARM with $\alpha=2$) models respectively. For each figure, left panels exhibit the NIC (black line) along with the $95\%$ HPD regions (grey areas) at the quantile confidence levels $\tau=\left(0.05\right)$, while right panels exhibit the NIC (black line) along with the $95\%$ HPD regions (grey areas) at the quantile confidence levels $\tau=\left(0.01\right)$. By looking at the tables and at the fugures we can assert that the Spline approach has at least two important advantages despite the bigger number of parameters to be estimated.

The first advantages consists to avoid an apriori choice for the shape of the NIC function allowing the data to speak by themself. Indeed all the CAViaR and CARE specifications considered in Section 2 impose a given structure of the NIC providing only a partial estimate of its true shape. The ASV, for istance, impose a ”V” shape relation between the risk measure at a given time $t$ and the past observations $y_{t-1}$ allowing just for an estimate of the overall slope of the NIC. The AS model only marginally improve upon the ASV allowing for a different slope for positive and negative returns but still it forces a ”V” shape for the NIC function.

From figures 1 and 4 we can see that the proposed BNL-CARM provides an estimate of the entire shape of NIC giving us unconstrained information on the impact of new information on the risk measure.

Second, the P-Spline approach used to model the NIC also allows us to estimate the threshold level of the returns such that asymmetric response to the risk is observed instead of assuming it equal to zero like in some of the forementioned models (see for example model (4) and (5)). It is well known infact that a common stylized fact of financial time series is that positive and negative returns have a different impact on volatility. In our model the threshold level is naturally estimated through the NIC function by the data; this point correspond to the (local) minimum of the NIC as reported in figures 1 and 4 and as can be observed is in general slightly greater than zero. This evidence is besides consistent with the recent financial crisis which has increased the degree of risk aversion among traders.

7 Conclusion

In this paper we present a new extension of the CAViaR and CARE models unified them in a Bayesian quantile regression framework called the B–CARM models by using the SEP likelihood approach. In addiction we propose a semiparametric P–spline framework i.e. the BNL-CARM model to relax some of the shape constrains present in those models in order to have a more flexible way to account for the well known stylized facts about financial time series. The Bayesian estimation methodology is carried out using the new adaptive Independent Metropolis–Hastings within Gibbs (IMG) technique already considered in Bernardi et al., (2018). The results obtained from the model validation backtesting criteria from five stock market indexes show that the model and the estimation methodology effectively capture the nonlinear relation between the unobserved $\tau$–level quantile and its determinants.

Appendix A Tables

Appendix B Figures