EconBase
← Back to paper

Virtual Historical Simulation for estimating the conditional VaR of large portfolios

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

75,189 characters · 15 sections · 0 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Virtual Historical Simulation for estimating the conditional VaR of large portfolios

quote\begin{center} Abstract \end{center} { In order to estimate the conditional risk of a portfolio's return, two strategies can be advocated. A multivariate strategy requires estimating a dynamic model for the vector of risk factors, which is often challenging, when at all possible, for large portfolios. A univariate approach based on a dynamic model for the portfolio's return seems more attractive. However, when the combination of the individual returns is time varying, the portfolio's return series is typically non stationary which may invalidate statistical inference. An alternative approach consists in reconstituting a "virtual portfolio", whose returns are built using the current composition of the portfolio and for which a stationary dynamic model can be estimated. This paper establishes the asymptotic properties of this method, that we call Virtual Historical Simulation. Numerical illustrations on simulated and real data are provided. }

{\it JEL Classification:} C13, C22 and C58.

{\it Keywords:} Accuracy of VaR estimation, Dynamic Portfolio, Estimation risk, Filtered Historical Simulation, Virtual returns.

Introduction

The quantitative standards laid down under Basel Accord II and III allow banks to develop internal models for setting aside capital. Methods that incorporate time dependence to quantify market risks are able to use knowledge of the conditional distribution. In particular, the conditional Value-at-Risk (VaR) of financial returns, with a given risk level $\alpha$ (typically, $\alpha=1\%$ or $5\%$) is nothing else, from a statistical point of view, than the negated $\alpha$-quantile of the conditional distribution of the portfolio returns. Estimating conditional quantiles, or more generally conditional risk measures, of a time series of financial returns is thus crucial for risk management.

It is also essential, for risk management purposes, to be able to evaluate the accuracy of such estimators of conditional risks. Uncertainty implied by statistical procedures in the implementation of risk measures may lead to false security in financial markets (see e.g. Farkas, Fringuellotti and Tunaru (2016) and the references therein). Estimation risk thus needs to be accounted for, in addition to market risk. However, evaluating the estimation risk for the conditional Value-at-Risk (VaR) is generally challenging for two main reasons. Firstly, because the stochastic nature of the conditional VaR does not allow in general to reduce the problem to the estimation of a parameter. Making inference on a stochastic process is obviously more intricate than on a parameter. Secondly, quantiles being obtained as the solutions of optimization problems based on non-smooth functions, establishing asymptotic properties of conditional VaR estimators may become a difficult task.

Increasing attention has been directed in the recent econometric literature to the inference of risk measures in dynamic risk models. Francq and Zako�an (2015) derived asymptotic confidence intervals (CI) for the conditional VaR of a series of financial returns driven by a parametric dynamic model. Robust backtesting procedures were developed by Escanciano and Olmo (2010, 2011), and Gouri�roux and Zako�an (2013) studied the effect of estimation on the coverage probabilities. Several articles proposed resampling methods: among others, Christoffersen and Gon�alves (2005) and Spierdijk (2016) considered using bootstrap procedures for constructing CIs for VaR; Hurlin, Laurent, Quaedvlieg and Smeekes (2017) proposed bootstrap-based comparison tests of two conditional risk measures. Beutner, Heinemann and Smeekes (2019) established the validity of a fixed-design residual bootstrap method for the two-step conditional VaR estimator of Francq and Zako�an (2015). See Nieto and Ruiz (2016) for an extensive survey of the methods for constructing and evaluating VaR forecasts that have been proposed in the literature.

Most existing studies on risk measure inference focus on the risk of a single financial asset. The aim of the present article is to estimate conditional VaR's for {\it portfolios} of financial assets. From a statistical point of view, the extension is far from trivial. First, because evaluating the quantile of a linear combination of variables may require knowledge of the complete joint distribution of such variables. When the object of interest is a {\it conditional} quantile, this approach requires specifying a dynamic model for the vector of returns of the assets involved in the portfolio. Second, portfolios compositions are generally time-varying, in particular if the agents adopt a mean-variance approach which, in a dynamic framework, requires specifying the first two conditional moments. This typically entails non-stationarity of the portfolio's return time series, as we shall see in more detail.

A natural approach for obtaining the conditional VaR of a portfolio relies on specifying a multivariate GARCH model for the vector of underlying asset returns. Rombouts and Verbeek (2009) proposed a semi-parametric approach relying on estimating the conditional density of the innovations vector and evaluating numerically the conditional VaR of a portfolio. The asymptotic properties of similar multivariate methods--with or without the assumption of sphericity of the innovations vector--were investigated by Francq and Zako�an (2018). As noted by Rombouts and Verbeek the advantage of multivariate approaches is to "take into account the dynamic interrelationships between the portfolio components, while the model underlying the VaR calculations is independent of the portfolio composition". On the other hand, for large portfolios multivariate approaches often become untractable due to the well-known dimensionality curse.

In this paper, we consider {\it univariate} procedures aiming at handling portfolios constructed with a large number of assets. We first consider a "naive" approach in which a standard volatility model is directly fitted to the portfolio returns. Despite its empirical relevance, we will see that the naive approach is not amenable to asymptotic statistical inference (due to the inherent non stationarity of the observed portfolio's time series). We study the asymptotic properties of an alternative procedure relying on a "virtual portolio" constructed with the {\it current composition} of the portfolio, on which a univariate model is fitted. This procedure--which we call Virtual Historical Simulation (VHS)--is amenable to asymptotic statistical inference. From a numerical point of view, it allows to avoid difficulties caused by the dimensionality curse in estimation of multivariate volatility models for vectors of asset returns.

The VHS method is related to other approaches introduced in Finance. The Basel Committee and European Union directives (UCITS) recommend that banks backtest their VaR measures against both "clean" and "dirty" P&L's of their trading portfolios (see Holton, 2014). Dirty P&L's are the actual P&L's reported at the end of the time horizon. They can be impacted by changes in the composition of the portfolio that occur during the VaR horizon. Since such position changes may have exogenous causes that can not be anticipated, it is relevant to backtest the VaR with the so-called clean P&L, which is the hypothetical P&L that would occur if the composition of the portfolio remained unchanged and if market moves were the only source of P&L change (see P�rignon, Deng and Wang, 2008). Clean P&L thus excludes P&L arising from intra-day trading, new trades, changes in reserves, fees and commissions. Clean P&L's are often used in the backtesting process, but they can also be used for VaR estimation. The VHS method exploits the idea of cleaning the P&L's for computing the VaR. At each past date $t$, one can compute a virtual return--the opposite of a clean (or hypothetical) P&L--that would occur if day $t$ positions were exactly those of the current date. Even if each bank uses its own internal VaR model, most financial institutions compute VaR through filtered or simple historical simulations on plain or hypothetical (virtual) returns (see Laurent and Omidi Firouzi, 2017). This is the aim of the present paper to study the asymptotic properties of such VaR evaluation methods.

The paper is organized as follows. Section (ref) defines the conditional VaR of a portfolio whose composition at the current date may depend on the historical prices, and presents the naive and VHS estimation methods. In Section (ref) we derive the asymptotic properties of the VHS procedure based on the Gaussian Quasi-Maximum Likelihood (QML) criterion, under general assumptions on the volatility model. Section (ref) presents some numerical illustrations based on Monte Carlo experiments and real financial data. Proofs are collected in the Appendix.

Estimating the conditional VaR

Conditional VaR of a dynamic portfolio

Let $\boldsymbol{p}_t=(p_{1t}, \ldots, p_{mt})'$ denote the vector of prices of $m$ assets at time $t$. Let $\boldsymbol{y}_t=(y_{1t}, \ldots, y_{mt})'$ denote the corresponding vector of log-returns, with $y_{it}=\log(p_{it}/p_{i,t-1})$ for $i=1, \ldots, m$.

Let $V_t$ denote the value at time $t$ of a portfolio composed of $\mu_{i,t-1}$ units of asset $i$, for $i=1, \ldots, m$:

equation[equation omitted — 123 chars of source]

where the $\mu_{i,t-1}$ are measurable functions of the prices up to time $t-1$, and the $\mu_{i}$ are constants. The return of the portfolio over the period $[t-1,t]$ is, for $t\geq 1$, assuming that $V_{t-1}\ne 0,$ $$\frac{V_{t}}{V_{t-1}}-1=\sum_{i=1}^m a_{i,t-1}e^{y_{it}}-1\approx \sum_{i=1}^m a_{i,t-1}y_{it}+a_{0,t-1}$$ where $$a_{i,t-1}=\frac{\mu_{i,t-1}p_{i,t-1}}{\sum_{j=1}^m \mu_{j,t-2}p_{j,t-1}}, \quad i=1, \ldots, m \quad \mbox{and}\quad a_{0,t-1}=-1+\sum_{i=1}^ma_{i,t-1}.$$ We assume that, at date $t$, the investor may rebalance his portfolio under a "self-financing" constraint.

itemize• The portfolio is rebalanced in such a way that $\sum_{i=1}^m \mu_{i,t-1}p_{it}=\sum_{i=1}^m \mu_{i,t}p_{it}.$

In other words, the value at time $t$ of the portfolio bought at time $t-1$ equals the value at time $t$ of the portfolio bought at time $t$. An obvious consequence of the self-financing assumption {\bf SF}, is that the change of value of the portfolio between $t-1$ and $t$ is only due to the change of value of the underlying assets: $$V_{t}-V_{t-1}=\sum_{i=1}^m \mu_{i,t-1}(p_{i,t}-p_{i,t-1}).$$ Another consequence is that the weights $a_{i,t-1}$ sum up to 1, that is $a_{0,t-1}=0$. Thus, under {\bf SF} we have $\frac{V_{t}}{V_{t-1}}-1\approx r_{t}$, where

equation[equation omitted — 178 chars of source]

for $i=1, \ldots, m,$ and $\mathbf{a}_{t-1}=(a_{1,t-1}, \ldots, a_{m,t-1})'$. A portfolio is usually called {\it crystallized} when the number of units of each asset is time independent, that is $\mu_{i,t-1}=\mu_i$ for each $i=1, \ldots, m$ and for all $t$. We will call {\it static} a portfolio with fixed proportion in value of each return, that is when $a_{i,t-1}=a_i$ for each $i=1, \ldots, m$ and for all $t$.

The {\it conditional} VaR of the portfolio's return process $(r_t)$ at risk level $\alpha\in (0,1)$, denoted $\mbox{VaR}_{t-1}^{(\alpha)}(r_t)$, is characterized by

equation[equation omitted — 110 chars of source]

where $P_{t-1}$ denotes the historical distribution conditional on the information $I_{t-1}$ available at time $t-1$. The specification of $I_{t-1}$ will depend on the approach used. Multivariate approaches use full information, that is all past prices of all assets. In the next approach we describe a univariate approach which only uses the past returns of the portfolio.

The naive approach

A natural approach for evaluating the conditional VaR in ((ref)) when $I_{t-1}=\sigma(r_{s}, s<t)$ is to estimate a univariate GARCH model, or any time series model, on the series of portfolio returns. We will see that this approach, which can be called "naive", may be misleading due to the fact that the return's portfolio is a time-varying combination of the individual returns.

For simplicity, we consider a crystallized portfolio, with weight $\mu_i$ and initial price $p_{i0}$ for the asset $i\in\{1, \ldots, m\}$. The composition $\mathbf{a}_{t-1}$ of such a portfolio is non stationary in general. Indeed, we have $$\log \left(\frac{a_{i,t}}{a_{j,t}}\right)=\log\left(\frac{\mu_ip_{i,0}}{\mu_jp_{j,0}}\right)+\sum_{k=1}^t\Delta_{i,j,k}, \qquad \Delta_{i,j,k}=y_{i,k}-y_{j,k},$$ and $(\sum_{k=1}^t\Delta_{i,j,k})_{t\geq 1}$ is a non stationary integrated process of order 1 under general assumptions.\footnote{By the Chung-Fuks theorem, this is the case when $\boldsymbol{y}_t$ is iid with zero mean and a non-singular covariance matrix $\boldsymbol{\Sigma}$. The non stationarity of the process also holds, for instance, if the sequence $(\Delta_{i,j,k})_k$ is mixing and nondegenerated.} More precisely, a consequence of the following lemma is that, with probability tending to one, the composition $\mathbf{a}_{t-1}$ of the portfolio converges to the set of the vectors $\boldsymbol{e}_i$ of the canonical basis (corresponding to single-asset portfolios): $P(\mathbf{a}_{t-1}\in\{\boldsymbol{e}_1,\dots,\boldsymbol{e}_m\})\to 1$ as $t\to\infty$.

lemConsider a process $(\Delta_k)_{k\geq 1}$. Assume that there exist real sequences $a_n>0$ and $b_n$, both tending to zero, such that \begin{equation} Z_n:=a_n\sum_{k=1}^n\Delta_k+b_n\stackrel{{\cal L}}{\to} Z as n\to\infty, \end{equation} for some random variable $Z$, whose cdf is continuous at 0 and such that $p=P(Z>0)\in (0,1)$. Then, for any $c>0$, we have $P(\sum_{k=1}^n\Delta_k>c)\to p$ and $P(\sum_{k=1}^n\Delta_k<-c)\to 1-p$ as $n\to\infty$.

Note that a generalized central limit theorem of the form (ref) holds for any iid sequence $(\Delta_k)$ whenever the distribution of $\Delta_k$ belongs to the domain of attraction of $Z$, which then follows a stable distribution. If the assumptions of Lemma (ref) hold with $\Delta_k=\Delta_{i,j,k}$ for any pair $(i,j)$, with $i\neq j$, then all the ratios $a_{i,t}/a_{j,t}$ are arbitrarily close to either 1 or 0 with probability tending to 1 as $t\to\infty$. In that case, the composition $\mathbf{a}_{t-1}$ tends to be totally undiversified, but is not always close to the same single-asset composition $\boldsymbol{e}_i$. If the dynamics of the individual returns $y_{it}$ are not identical, the dynamics of the return $r_t$ will be time-varying, and the naive method based on a fixed stationary GARCH model is likely to produce poor results.

Simulation experiments reported in Section (ref) confirm that for crystallized portfolios, the naive approach may behave badly due to the non stationarity of the univariate returns $r_t$. Of course, for static portfolios the non stationarity issue vanishes, but such portfolios may be considered as artificial\footnote{as it would require rebalancing every day in order to maintain a fixed composition in percentages. In general, portfolios are rebalanced to minimize risk, leading to a non-static composition.}. The next section studies a remedy to the non stationarity issue, while keeping the univariate framework.

The VHS approach

An alternative to the naive approach consists in reconstituting a "virtual portfolio", whose returns are built using the current composition of the portfolio. Given the portfolio composition $\mathbf{a}_{t_0-1}=\boldsymbol{x}$, say, at time $t_0$, we construct a process of virtual returns $$r_{t}^*(\boldsymbol{x})=\boldsymbol{x}'\boldsymbol{y}_t, \qquad t\in \mathbb{Z}$$ and we consider the information set $I_{t_0-1}=\sigma(r_{s}^*(\boldsymbol{x}), s<t_0)$. Note that, in general, $r_{t}^*(\boldsymbol{x})\ne r_t$ because the composition of the (non virtual) portfolio is time varying ($\mathbf{a}_{t-1}\ne \boldsymbol{x}$, in general, for $t\ne t_0$). Given the stationarity of $(\boldsymbol{y}_t)$, it is clear that the series of virtual returns $\{r_t^*(\boldsymbol{x})\}$ is also stationary, with conditional moments $$E_{t-1}[r_t^*(\boldsymbol{x})] =:\mu_t(\boldsymbol{x}), \qquad \mbox{var}_{t-1}[r_t^*(\boldsymbol{x})] =:\sigma_t^2(\boldsymbol{x}),$$ where $E_{t-1}(X)=E(X\mid r_{s}^*(\boldsymbol{x}), s<t)$ for any variable $X$, and the variance is defined accordingly. Thus, $r_{t}^*(\boldsymbol{x})$ follows a model of the form

equation[equation omitted — 170 chars of source]

Noting that $r_{t_0}=r^*_{t_0}(\mathbf{a}_{t_0-1})$, the conditional VaR at time $t_0$ thus satisfies

equation[equation omitted — 180 chars of source]

where $\mbox{VaR}_{t-1}^{*(\alpha)}(X)$ is the VaR of $X$ at level $\alpha$ conditional on $I_{t-1}.$

Note that the martingale difference $(u_t)$ may not be iid, as the following example illustrates.

exemple{\rm Consider the bivariate ARCH(1) process, defined as the stationary non anticipative solution of the model $$\boldsymbol{y}_t=\left(\begin{array}{c}y_{1t}\\y_{2t}\end{array}\right)= \boldsymbol{\Sigma}_t\boldsymbol{\eta}_t, \quad \boldsymbol{\Sigma}_t=\left(\begin{array}{cc}\sigma_{1t}^2:= \omega_1+\alpha_{11}y_{1,t-1}^2+\alpha_{12}y_{2,t-1}^2&0\\0&\sigma_{2t}^2:= \omega_2+\alpha_{21}y_{1,t-1}^2\end{array}\right),$$ where $\boldsymbol{\eta}_t\mbox{is iid }(\boldsymbol{0}, \boldsymbol{I})$, and assuming that the components $\eta_{1t}$ and $\eta_{2t}$ are independent. Let a portfolio which is fully invested in the first asset, hence $r_t=(1,0) \boldsymbol{y}_t=y_{1t}$ for all $t$. Denote by ${\cal F}_{1t}$ the $\sigma$-field generated by $\{y_{1u}, u\leq t\}.$ We have $E(r_t| {\cal F}_{1,t-1})=0$ and \begin{eqnarray*}E(r_t^2| {\cal F}_{1,t-1})&=&E(\sigma_{1t}^2| {\cal F}_{1,t-1})=\omega_1+\alpha_{11}y_{1,t-1}^2+\alpha_{12}E(y_{2,t-1}^2| {\cal F}_{1,t-1})\\ &=&\omega_1+\alpha_{11}y_{1,t-1}^2+\alpha_{12}\sigma_{2,t-1}^2 E(\eta_{2,t-1}^2| {\cal F}_{1,t-1})\\ &=&\omega_1+\alpha_{11}y_{1,t-1}^2+\alpha_{12}\sigma_{2,t-1}^2 :=\sigma_t^2.\end{eqnarray*} It follows that $(r_t)$ satisfies the model $r_t=\sigma_tu_t$, where $$u_t=\frac{\sigma_{1t}}{\sigma_t}\eta_{1t}=\left(1+ \frac{\alpha_{12}\sigma_{2,t-1}^2(\eta_{2,t-1}^2-1)}{\omega_1+\alpha_{11}y_{1,t-1}^2+\alpha_{12}\sigma_{2,t-1}^2}\right)^{1/2}\eta_{1t}.$$ It is then clear that $(u_t, {\cal F}_{1,t})$ is a martingale difference but $(u_t)$ is generally not iid (except when $\alpha_{12}=0$ or $\eta_{2,t}^2$ is degenerated). }

Even in the simple previous example, the conditional quantile $\mbox{VaR}_{{t_0}-1}^{*(\alpha)}(u_{t_0})$ popping up in ((ref)) cannot be explicitly computed. Whether or not this quantity could be estimated nonparametrically is beyond the scope of this paper. Instead, we consider a "hybrid" VaR defined by

equation[equation omitted — 167 chars of source]

where $\mbox{VaR}^{(\alpha)}(u)$ is the marginal VaR of $u_{t}$ at level $\alpha$. An estimator of $\mbox{VaR}_{H,{t_0}-1}^{(\alpha)}(r_{t_0})$ is obtained as follows: given $\mathbf{a}_{{t_0}-1}=\boldsymbol{x}$,

{\sc Step 1:} Compute the virtual historical returns $r_{t}^*(\boldsymbol{x})$ for $t=1, \ldots, n$.

{\sc Step 2:} Estimate $\mu_t(\boldsymbol{x})$ and $\sigma_t(\boldsymbol{x})$. Denote by $\hat{\mu}_t(\boldsymbol{x})$ and $\hat{\sigma}_t(\boldsymbol{x})$ the resulting estimators, and by $\hat{u}_t=\{r_{t}^*(\boldsymbol{x})-\hat{\mu}_t(\boldsymbol{x})\}/\hat{\sigma}_t(\boldsymbol{x})$ the residuals.

{\sc Step 3:} Compute the $\alpha$-quantile $\xi_{n,\alpha}$ of $\{\hat{u}_s, 1\leq s \leq n\}$ and define an estimator of $\mbox{VaR}_{H,{t_0}-1}^{(\alpha)}(r_{t_0})$ as

equation[equation omitted — 168 chars of source]

Step 2 can be implemented by estimating a parametric model. This approach will be developed in Section (ref).

This procedure is particularly appropriate for large portfolios, when the large dimension of the vector of underlying assets precludes--or at least formidably complicates-- estimation of multivariate volatility models. Moreover, the following example shows that for large portfolios a univariate GARCH model is a reasonable assumption for the virtual returns.

exemple{\rm Suppose that $m$ is large and that the vector of log-returns is driven by a vector $\boldsymbol{f}_t$ of $K$ factors (with $K\ll m$) as $$\boldsymbol{y}_t=\boldsymbol{\beta}\boldsymbol{f}_t+\boldsymbol{u}_t$$ where $\boldsymbol{\beta}$ is a $m\times K$ matrix, $\mbox{var}_{t-1}(\boldsymbol{f}_t)=\boldsymbol{F}_t$ is a full-rank matrix and $\mbox{var}_{t-1}(\boldsymbol{u}_t)=\boldsymbol{\Sigma}_t$. With a composition fixed to $\boldsymbol{x}$, the virtual portfolio's returns thus satisfy $$r_t^*=\boldsymbol{x}'\boldsymbol{\beta}\boldsymbol{f}_t+\boldsymbol{x}'\boldsymbol{u}_t.$$ Suppose that the portfolio is well-diversified, so that the components $x_i$ of $\boldsymbol{x}$ satisfy $x_{i}=O(1/m)$ for $i=1, \ldots, m$ as $m\to \infty$. Under appropriate assumptions $\mbox{var}_{t-1}(\boldsymbol{x}'\boldsymbol{u}_t)$ converges to 0 as $m\to \infty$. On the other hand, $\mbox{var}_{t-1}(\boldsymbol{x}'\boldsymbol{\beta}\boldsymbol{f}_t)=\boldsymbol{x}'\boldsymbol{\beta} \boldsymbol{F}_t\boldsymbol{\beta}'\boldsymbol{x}=O_P(1)$ and does not vanish as $m$ increases under appropriate assumptions. \footnote{For instance if the matrix $\boldsymbol{\beta}$ does not contain too many zeroes or, more precisely, if at least one column $\boldsymbol{\beta}_j$ of $\boldsymbol{\beta}$ is such that $\lim \inf |\boldsymbol{x}'\boldsymbol{\beta}_j|>0$ as $m\to \infty$.} It follows that $r_t^*\approx \boldsymbol{x}'\boldsymbol{\beta}\boldsymbol{f}_t$. If now $K=1$ and the (real-valued) factor $f_t$ is the solution of a GARCH model, the process $\boldsymbol{x}'\boldsymbol{\beta} \bf_t$ will follow the same model up to a change of scale. It is therefore natural to fit a GARCH model for the virtual returns under these assumptions when $m$ is large. }

Asymptotic properties of the VHS approach

To obtain asymptotic properties of the VHS procedure, we make the following parametric assumptions on Model ((ref)). For simplicity, we consider the model without conditional mean, that is $\mu_t(\boldsymbol{x})=0.$ For some (known) function $\sigma \;:\; \mathbb{R}^{\infty}\times \Theta \rightarrow (0,\infty)$, let

equation[equation omitted — 153 chars of source]

where $\boldsymbol{\theta}_0=\boldsymbol{\theta}_0(\boldsymbol{x})$ is the true value of the finite dimensional parameter $\boldsymbol{\theta}$, belonging to some compact subset $\Theta$ of $\mathbb{R}^d$. To alleviate notations, we will denote the virtual returns by $\epsilon_t:=r_{t}^*(\boldsymbol{x})$ and replace $\sigma_t(\boldsymbol{x}; \boldsymbol{\theta})$ by $\sigma_t(\boldsymbol{\theta})$. Model (ref) thus writes

equation[equation omitted — 249 chars of source]

Given (virtual) observations $\epsilon_1,\dots, \epsilon_n$, and arbitrary initial values $\tilde{\epsilon}_i$ for $i\leq 0$, we define $$ \tilde{\sigma}_t(\boldsymbol{\theta})=\sigma(\epsilon_{t-1}, \epsilon_{t-2}, \ldots, \epsilon_1, \tilde{\epsilon}_0, \tilde{\epsilon}_{-1}, \ldots ; \boldsymbol{\theta}). $$ The Gaussian QML criterion is defined by

equation[equation omitted — 263 chars of source]

Let the QML estimator of $\boldsymbol{\theta}_0$,

equation[equation omitted — 131 chars of source]

which has the particularity of being based on virtual returns rather than on observations.

To study the asymptotic properties of the VHS estimator,

equation[equation omitted — 138 chars of source]

we introduce the following additional assumptions. Recall that iidness of the sequence $(u_t)$ is not a natural assumption in our framework (see Example (ref)). Let $\boldsymbol{D}_t(\boldsymbol{\theta})=\sigma_t^{-1}(\boldsymbol{\theta})\partial \sigma_t(\boldsymbol{\theta})/\partial\boldsymbol{\theta}$, $\boldsymbol{D}_t=\boldsymbol{D}_t(\boldsymbol{\theta}_0) $. Let also ${\cal F}_{t}$ the sigma-field generated by $\{u_k,k\le t\}$, and ${\cal F}_{t:t-s}$ the sigma-field generated by $\{u_k,t-s\le k\le t\}$ with $s>0$.

itemize• The sequence $(u_t)$ is stationary, with $E|u_t|^{4+\nu}<\infty$ for some $\nu>0$, and mixing coefficients $\left\{\alpha(h)\right\}_{h\geq 0}$ satisfying, for $\epsilon\in (0,\nu)$, $$\sum_{h=1}^{\infty}h^{r^*}\{\alpha(h)\}^{\frac{\nu-\epsilon}{4+\nu-\epsilon}}<\infty\mbox{ for some }r^*>\frac{2\kappa\{2+\nu-\epsilon\}}{\nu-\epsilon-2\kappa\{2+\nu-\epsilon\}} \mbox{ and }\kappa\in \left(0,\frac{\nu-\epsilon}{4(2+\nu-\epsilon)}\right).$$ Suppose that $E(u_t\mid {\cal F}_{t-1})=0$ and $E(u^2_t\mid {\cal F}_{t-1})=1$. Let $\xi_{\alpha}$ the $\alpha$-quantile of $u_t$. Assume that the conditional distribution of $u_t$ given ${\cal F}_{t-1}$ has a density $f_{t-1}$ such that $f_{t-1}(\xi_{\alpha})>0$ a.s. and $E\sup_{\xi\in V(\xi_{\alpha})}f^4_{t-1}(\xi)<\infty$ for some neighborhood $V(\xi_{\alpha})$ of $\xi_{\alpha}$. Assume also that this density is continuous at $\xi_{\alpha}$ uniformly in ${\cal F}_{t-1}$, in the sense that for sufficiently small $\varepsilon> 0$, there exists a stationary and ergodic sequence $(K_t)$ such that $K_{t-1}\in {\cal F}_{t-1}$ and $$\sup_{x\in[\xi_{\alpha}-\varepsilon,\xi_{\alpha}+\varepsilon]}\left|f_{t-1}(x)-f_{t-1}(\xi_{\alpha})\right|\leq K_{t-1}\varepsilon$$ with $EK^4_{t}<\infty$ a.s. • $(\epsilon_t)$ is a strictly stationary and ergodic solution of ((ref)), and there exists $s_0>0$ such that $E|\epsilon_1|^{s_0}<\infty$. • There exists a sequence $\boldsymbol{D}_{t,T_n}$ such that $\boldsymbol{D}_t=\boldsymbol{D}_{t,T_n}+o_P(1)\mbox{ as } n\to\infty$, where $T_n\to\infty$ and $T_n=O(n^{\kappa})$ (with $\kappa$ defined in {\bf A1}), $\boldsymbol{D}_{t,T_n}\in {\cal F}_{t-1:t-T_n}$ and for any $r\geq 0$ $$E\|\boldsymbol{D}_{t}\|^r<\infty,\qquad \sup_{n\geq 1}E\|\boldsymbol{D}_{t,T_n}\|^r<\infty. $$ • For some $\underline{\omega}>0$, almost surely, $\sigma_t(\boldsymbol{\theta})$ and $\tilde{\sigma}_t(\boldsymbol{\theta})$ belong to $(\underline{\omega}, \infty]$ for any $\boldsymbol{\theta}\in \Theta$ and any $t\ge 1$. For $\boldsymbol{\theta}_1, \boldsymbol{\theta}_2\in \Theta$, we have $\sigma_t(\boldsymbol{\theta}_1)=\sigma_t(\boldsymbol{\theta}_2)\; a.s.$ if and only if $\boldsymbol{\theta}_1=\boldsymbol{\theta}_2$. Moreover, for any $\boldsymbol{x}\in \mathbb{R}^m$, $\boldsymbol{x}'\boldsymbol{D}_t(\boldsymbol{\theta}_0)=0$ a.s. entails $\boldsymbol{x}=0$. • There exist a random variable $C$ which is measurable with respect to $\{\epsilon_t, t<0\}$ and a constant $\rho\in(0,1)$ such that $\sup_{\boldsymbol{\theta}\in\Theta}|\sigma_t(\boldsymbol{\theta})-\widetilde{\sigma}_t(\boldsymbol{\theta})|\leq C\rho^t$. • The function $\boldsymbol{\theta} \mapsto \sigma (x_1, x_2, \ldots ; \boldsymbol{\theta})$ has continuous second-order derivatives, and $$\sup_{\boldsymbol{\theta}\in\Theta}\left\|\frac{\partial \sigma_t(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} -\frac{\partial\widetilde{\sigma}_t(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} \right\| +\left\|\frac{\partial^2 \sigma_t(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}\partial \boldsymbol{\theta}'} -\frac{\partial^2\widetilde{\sigma}_t(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}\partial \boldsymbol{\theta}'} \right\| \leq C\rho^t,$$ where $C$ and $\rho$ are as in {\bf A5}. • There exists a neighborhood $V(\boldsymbol{\theta}_0)$ of $\boldsymbol{\theta}_0$ and $\tau>0$ such that $$\sup_{\boldsymbol{\theta} \in V(\boldsymbol{\theta}_0)} \left\|\frac1{\sigma_t(\boldsymbol{\theta})}\frac{\partial \sigma_t(\boldsymbol{\theta})}{\partial\boldsymbol{\theta}}\right\|^4, \quad \sup_{\boldsymbol{\theta} \in V(\boldsymbol{\theta}_0)} \left\|\frac1{\sigma_t(\boldsymbol{\theta})}\frac{\partial^2 \sigma_t(\boldsymbol{\theta})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}'}\right\|^2, \quad \sup_{\boldsymbol{\theta} \in V(\boldsymbol{\theta}_0)} \left|\frac{\sigma_t(\boldsymbol{\theta}_0)}{\sigma_t(\boldsymbol{\theta})}\right|^{\frac{2(4+\nu)(1+\tau)}{2+\nu}},$$ have finite expectations, where $\nu$ is as in {\bf A1}.

Assumptions {\bf A2}, {\bf A4} and {\bf A5}, together with $E(u^2_t\mid {\cal F}_{t-1})=1$, are sufficient to prove the strong consistency of the QML estimator $\hat{\boldsymbol{\theta}}_n$. Assumption {\bf A6} allows to show that the initial values $\tilde{\epsilon}_i$ do not matter for the asymptotic normality of $\hat{\boldsymbol{\theta}}_n$. Assumptions {\bf A1}, {\bf A3} and {\bf A7} are used to apply a CLT for a mixing triangular array based on the approximation of $\boldsymbol{D}_t$ by $\boldsymbol{D}_{t,T_n}$.

For particular volatility models, some of the assumptions can be simplified as the following lemma shows.

lemFor the standard GARCH(1,1) model \begin{equation} \epsilon_t=\sigma_t u_t,\quad \sigma_t^2=\omega_0+\alpha_0\epsilon_{t-1}^2+ \beta_0\sigma_{t-1}^2,\qquad \omega_0>0, \alpha_0> 0, \beta_0>0, \end{equation} where $(u_t)$ satisfies {\bf A1}, Assumptions {\bf A2-A7} reduce to: i) $E\log (\alpha_0u_t^2+\beta_0)<0$; ii) $u_t^2$ has a non-degenerate distribution; iii) $\Theta=\{(\omega, \alpha, \beta)\}$ is a compact subset of $(0,\infty)^3$ such that, for all $\theta\in \Theta$, $\omega>\underline{\omega}$ for some $\underline{\omega}>0$ and $\beta<1$. iv) $E|\epsilon_t|^{s_0}<\infty$ for $s_0>0$, where $(\epsilon_t)$ is the strictly stationary solution implied by i).

We are now in a position to state our main result.

theoAssume $\xi_{\alpha}<0$. Let {\bf A1-A7} hold. Then $\hat{\boldsymbol{\theta}}_n\to\boldsymbol{\theta}_0$ a.s. as $n\to \infty$, and \begin{eqnarray*} \left(\begin{array}{c}\sqrt{n}\left(\hat{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_0\right) \\ \sqrt{n}(\xi_{\alpha}-\xi_{n, \alpha}) \end{array}\right) &\stackrel{{\cal L}}{\to}&{\cal N}(0,\boldsymbol{\Sigma}_{\alpha}), \qquad \boldsymbol{\Sigma}_{\alpha}=\left(\begin{array}{cc} \boldsymbol{J}^{-1}\boldsymbol{S}^{11}\boldsymbol{J}^{-1}& \boldsymbol{\Lambda}_{\alpha}\\ \boldsymbol{\Lambda}_{\alpha}' &\zeta_{\alpha} \end{array}\right), \end{eqnarray*} where $\boldsymbol{J}=E(\boldsymbol{D}_t\boldsymbol{D}_t')$ with $\boldsymbol{D}_t=\boldsymbol{D}_t(\boldsymbol{\theta}_0)$, and \begin{eqnarray*} \boldsymbol{\Lambda}_{\alpha}&=& \frac{\xi_{\alpha}}{4E[f_{t-1}(\xi_{\alpha})]}\boldsymbol{J}^{-1}\boldsymbol{S}^{11}\boldsymbol{J}^{-1}\boldsymbol{\Psi}_{\alpha}+\frac1{2E[f_{t-1}(\xi_{\alpha})]}\boldsymbol{J}^{-1}\boldsymbol{S}_{\alpha}^{12}, \qquad \boldsymbol{\Psi}_{\alpha}=E[f_{t-1}(\xi_{\alpha})\boldsymbol{D}_t],\\ \zeta_{\alpha}&=&\frac1{\{E[f_{t-1}(\xi_{\alpha})]\}^2}\left(\frac14\xi_{\alpha}^2\boldsymbol{\Psi}_{\alpha}'\boldsymbol{J}^{-1}\boldsymbol{S}^{11}\boldsymbol{J}^{-1}\boldsymbol{\Psi}_{\alpha}+ \xi_{\alpha}\boldsymbol{\Psi}_{\alpha}'\boldsymbol{J}^{-1}\boldsymbol{S}_{\alpha}^{12}+\boldsymbol{S}_{\alpha}^{22}\right),\\ \boldsymbol{S}^{11}&=&E\left[\left(u_t^2-1\right)^2\boldsymbol{D}_t\boldsymbol{D}_t'\right], \qquad \boldsymbol{S}_{\alpha}^{22}=\sum_{h=-\infty}^{\infty}cov\left(\mathbf{1}_{\{u_t<\xi_{\alpha}\}},\mathbf{1}_{\{u_{t-h}<\xi_{\alpha}\}}\right),\\ \boldsymbol{S}_{\alpha}^{12}&=&\boldsymbol{S}_{\alpha}^{21 '}=\sum_{h=0}^{\infty}cov[(u_t^2-1)\boldsymbol{D}_t,\mathbf{1}_{\{u_{t+h}<\xi_{\alpha}\}}]. \end{eqnarray*}
remWhen the errors $u_t$ are iid, the expression of the asymptotic covariance matrix simplifies considerably. More precisely, we have \begin{eqnarray*} \boldsymbol{S}^{11}&=&\left(\kappa_4-1\right)\boldsymbol{J}, \qquad \boldsymbol{S}_{\alpha}^{22}=\alpha(1-\alpha),\quad \boldsymbol{S}_{\alpha}^{12}=\boldsymbol{S}_{\alpha}^{21 '}=[E\left(u_t^2\mathbf{1}_{\{u_{t}<\xi_{\alpha}\}}\right)-\alpha]E(\boldsymbol{D}_t):= p_{\alpha}E(\boldsymbol{D}_t),\\ \boldsymbol{\Lambda}_{\alpha}&=& \left(\frac{\xi_{\alpha}(\kappa_4-1)}{4}+\frac{p_{\alpha}}{2f(\xi_{\alpha})}\right) \boldsymbol{J}^{-1}\boldsymbol{\Omega}, \qquad \boldsymbol{\Omega}=E(\boldsymbol{D}_t),\\ \zeta_{\alpha}&=&\frac{\kappa_4-1}4\xi_{\alpha}^2\boldsymbol{\Omega}'\boldsymbol{J}^{-1}\boldsymbol{\Omega}+\frac{\xi_{\alpha}p_{\alpha}}{f(\xi_{\alpha})}\boldsymbol{\Omega}'\boldsymbol{J}^{-1}\boldsymbol{\Omega}+\frac{\alpha(1-\alpha)}{f^2(\xi_{\alpha})}= \frac{\kappa_4-1}4\xi_{\alpha}^2+\frac{\xi_{\alpha}p_{\alpha}}{f(\xi_{\alpha})}+\frac{\alpha(1-\alpha)}{f^2(\xi_{\alpha})}. \end{eqnarray*} where $f$ denotes the density of $u_t$, and $\kappa_4=Eu_t^4$. For the last equality, we used the relation $\boldsymbol{\Omega}'\boldsymbol{J}^{-1}\boldsymbol{\Omega}=1$ (see (8) in Francq and Zako�an (2013)). Hence, when $(u_t)$ is iid we have $$\boldsymbol{\Sigma}_{\alpha}=\left(\begin{array}{cc} \left(\kappa_4-1\right)\boldsymbol{J}^{-1}& \left(\frac{\xi_{\alpha}(\kappa_4-1)}{4}+\frac{p_{\alpha}}{2f(\xi_{\alpha})}\right) \boldsymbol{J}^{-1}\boldsymbol{\Omega}\\ \left(\frac{\xi_{\alpha}(\kappa_4-1)}{4}+\frac{p_{\alpha}}{2f(\xi_{\alpha})}\right)\boldsymbol{\Omega}'\boldsymbol{J}^{-1} &\frac{\kappa_4-1}4\xi_{\alpha}^2+\frac{\xi_{\alpha}p_{\alpha}}{f(\xi_{\alpha})}+\frac{\alpha(1-\alpha)}{f^2(\xi_{\alpha})} \end{array}\right).$$ In particular, the asymptotic variance of $\sqrt{n}(\xi_{\alpha}-\xi_{n,\alpha})$ only depends on the errors distribution, not on the volatility parameter $\boldsymbol{\theta}_0$. Nevertheless, the estimation of $\boldsymbol{\theta}_0$ affects the asymptotic accuracy, since the asymptotic variance would be equal to $\frac{\alpha(1-\alpha)}{f^2(\xi_{\alpha})}$ if the $u_t$'s were observed.
remConsistent estimation of ${\boldsymbol{\Sigma}}_{\alpha}$ is crucial for evaluating the estimation risk. Some matrices involved in the asymptotic covariance matrix have the form of expectations, and can therefore be straightforwardly estimated by their empirical counterpart. This is the case of $\boldsymbol{J}$ and $\boldsymbol{S}^{11}$. To estimate $\boldsymbol{S}^{22}_{\alpha}$ and $\boldsymbol{S}^{12}_{\alpha}$, a classical HAC (heteroskedasticity and autocorrelation consistent) estimator can be used (see for instance Andrews (1991), Newey and West (1987)). Estimation of the matrix $\boldsymbol{\Psi}_{\alpha}$ is more tricky. We propose two approaches: i) noting that \begin{eqnarray*} \boldsymbol{\Psi}_{\alpha}&=&\lim_{h\to 0} \frac 1h E\left[\{P_{t-1}(u_t<\xi_{\alpha}+h)-P_{t-1}(u_t<\xi_{\alpha})\}\boldsymbol{D}_t\right]\\ &=&\lim_{h\to 0} \frac 1h E\left[E_{t-1}\{\mathbf{1}_{\{u_t<\xi_{\alpha}+h\}}-\mathbf{1}_{\{u_t<\xi_{\alpha}\}}\}\boldsymbol{D}_t\right]= \lim_{h\to 0} \frac 1h E\left[(\mathbf{1}_{\{u_t<\xi_{\alpha}+h\}}-\mathbf{1}_{\{u_t<\xi_{\alpha}\}})\boldsymbol{D}_t\right], \end{eqnarray*} a natural estimator for $\boldsymbol{\Psi}_{\alpha}$ is \begin{eqnarray*} \widehat{\boldsymbol{\Psi}}_{\alpha}&=& \frac1{nh_n} \sum_{t=1}^n(\mathbf{1}_{\{\hat{u}_t<\xi_{n,\alpha}+h_n\}}-\mathbf{1}_{\{\hat{u}_t<\xi_{n,\alpha}\}})\hat{\boldsymbol{D}}_t, \end{eqnarray*} where $\hat{\boldsymbol{D}}_t=\tilde{\sigma}_t^{-1}(\hat{\boldsymbol{\theta}}_n)\partial \tilde{\sigma}_t(\hat{\boldsymbol{\theta}}_n)/\partial\boldsymbol{\theta}$ and $h_n$ is a bandwidth parameter; ii) alternatively, one may approximate $f_{t-1}$ by the density of $u_t$ conditional on $u_{t-1}$ (instead of the infinite past). A standard kernel estimator for this density at point $\xi_{\alpha}$ is $$\hat{f}_{t-1}(\xi_{\alpha})= \frac{\sum_{s=1}^nK_{h_2}(\hat{u}_{t-1}-\hat{u}_{s-1})K_{h_1}(\xi_{\alpha}-\hat{u}_{s})}{\sum_{s=1}^nK_{h_2}(\hat{u}_{t-1}-\hat{u}_{s-1})},$$ where $h_1, h_2$ are bandwidths, $K_h(x)=h^{-1}K(x/h)$ for any bandwidth $h$, and $K$ is a bounded symmetric kernel function (see e.g. Hansen (2004)). Matrix $\boldsymbol{\Psi}_{\alpha}$ can thus be estimated by \begin{eqnarray*} \widehat{{\boldsymbol{\Psi}}}^{\circ}_{\alpha}&=&\frac 1n\sum_{t=1}^n \hat{f}_{t-1}(\xi_{\alpha})\hat{\boldsymbol{D}}_t. \end{eqnarray*} The asymptotic properties under Assumptions {\bf A1-A7} of both estimators, $\widehat{\boldsymbol{\Psi}}_{\alpha}$ and $\widehat{{\boldsymbol{\Psi}}}^{\circ}_{\alpha}$, are left for further research.

We have the Taylor expansion $$F_{t_0}(\hat{\boldsymbol{\theta}}_n, \xi_{n,\alpha}):=\tilde{\sigma}_{t_0}(\hat{\boldsymbol{\theta}}_n)\xi_{n,\alpha}=F_{t_0}({\boldsymbol{\theta}}_0, \xi_{\alpha})+ \frac{\partial}{\partial (\boldsymbol{\theta}', \xi)}{F}_{t_0}({\boldsymbol{\theta}}_n^*, \xi_{n}^*) \sqrt{n}\left(

array[array omitted — 96 chars of source]

\right),$$ where $({\boldsymbol{\theta}}_n^*, \xi_n^*)$ is between $({\boldsymbol{\theta}}_n, \xi_{n,\alpha})$ and $(\boldsymbol{\theta}_0, \xi_{\alpha})$. Given a fixed information set $I_{t_0-1}$, it can be shown using {\bf A5-A6} that the observations available at time $t_0-1$ have no effect on the asymptotic distribution of Theorem \ref{hard}. Thus, an approximate {\it conditional} $(1-\alpha_0)%$ CI for the hybrid VaR (ref) has bounds given by

equation[equation omitted — 251 chars of source]

where $\widehat{\mbox{VaR}}_{VHS, t_0-1}^{(\alpha)}(r_{t_0})$ is defined in (ref), $\widehat{\boldsymbol{\Sigma}}_{\alpha}$ is a consistent estimator of ${\boldsymbol{\Sigma}}_{\alpha}$ and

equation*[equation* omitted — 228 chars of source]

In other words, $$\lim_{n\to \infty} P_{t_0-1}\left[\left|{\mbox{VaR}}_{H, t_0-1}^{(\alpha)}(r_{t_0})- \widehat{\mbox{VaR}}_{VHS, t_0-1}^{(\alpha)}(r_{t_0})\right| \le \frac1{\sqrt{n}} \Phi^{-1}(1-\alpha_0/2) \left\{ \boldsymbol{\delta}_{{t_0}-1}'\widehat{\boldsymbol{\Sigma}}_{\alpha} \boldsymbol{\delta}_{{t_0}-1} \right\}^{1/2}\right]= 1-\alpha_0. $$

remIt is worth noting that the previous CI has to be understood conditionally on $I_{t_0-1}$ for fixed $t_0$. In particular, it cannot be directly applied when $t_0-1=n$. Indeed, the bounds in ((ref)) are no longer random in this case. \footnote{except if one assumes, as in Gong, Li and Peng (2010), that the number of conditioning values is finite as in the ARCH$(q)$ case.} See Beutner, Heinemann and Smeekes (2019) for an asymptotic justification of this type of CI for conditional objects.

Numerical illustrations

An alternative to the univariate approaches is a multivariate strategy, requiring a dynamic model for the vector of risk factors $\boldsymbol{y}_t$. We describe in Appendix (ref) two multivariate approaches -- the spherical method and the Filtered Historical Simulation (FHS) method. Such methods might perform better than univariate ones because they incorporate information stemming from individual returns rather than an aggregate information stemming from portfolio returns. However, at least for large portfolios, multivariate approaches can be very challenging, when at all possible.

The first section is devoted to numerical comparisons of the different methods. We first study the reliability of the naive approach in a static framework where the returns are iid and the portfolio is crystallized. In a dynamic framework, we then compare the performance of the univariate and multivariate approaches when the dynamic multivariate model is well specified and the dimension is small. Next, we consider misspecified GARCH models and, possibly, higher dimensions. The second section concerns real data examples based on large sets of US stocks.\footnote{ The code and data used in the paper are available on the web site\\ \url{http://perso.univ-lille3.fr/ cfrancq/Christian-Francq/VaRPortfolio.html} }

Monte-Carlo experiments

The theoretical conditional VaR in ((ref)) depends on the information set $I_{t-1}$. In our first two sets of experiments, the vector of individual returns will be simulated using multivariate GARCH and the different estimators will be compared to the same target VaR obtained by taking the full information set (i.e. including the returns of the individual assets, not only the returns of the portfolio). In the third set of experiments, we will estimate a misspecified multivariate model for which the true conditional VaR is no longer available.

Static model and crystallized portfolio

For a simple illustration of Lemma (ref), we consider a crystallized equally weighted portfolio of 3 assets (of initial price $p_{i0}=1000$): $V_t=\sum_{i=1}^3 p_{it}$ and $\mu_{i,t}=1$ for $i=1,2,3$. Thus, the return portfolio composition is time varying, with coefficients $\mathbf{a}_{t-1}=(a_{1,t-1}, a_{2,t-1}, a_{3,t-1})'$ and $a_{i,t-1}=p_{i,t-1}/\sum_{j=1}^3 p_{j,t-1}$. Assume that the vector of the log-returns is iid, Gaussian, centered, with variance $\mbox{Var}(\boldsymbol{y}_t)=\boldsymbol{\Sigma}^2=\boldsymbol{D}\boldsymbol{R}\boldsymbol{D},$ with $$\boldsymbol{D}=\left(

array[array omitted — 47 chars of source]

\right),\quad \boldsymbol{R}=\left(

array[array omitted — 65 chars of source]

\right).$$ The composition $\mathbf{a}_{t-1}$ of the portfolio is plotted in Figure \ref{Composition1}. As we have seen in Section \ref{secInval}, this vector is non stationary. More precisely, as discussed in Section \ref{secInval}, with increasing probability the composition $\mathbf{a}_{t-1}$ of the portfolio is arbitrarily close to one of the three single-asset portfolios $(1,0, 0), (0,1,0)$ and $(0,0,1)$. Experiments conducted with different parameters (in particular different correlation matrices) led to the same conclusion.

It is thus non surprising to see that the univariate return series $r_t$ plotted in Figure (ref) exhibits some nonstationarity features, in particular unconditional heteroscedasticity. The increased variance in the second part of the sample reflects the fact that the portfolio tends to be less and less diversified (see Figure (ref)).

When $\boldsymbol{\Sigma}_t$ and $\boldsymbol{m}_t$ are constant, the FHS estimator in (ref) reduces to the opposite of the quantile of the portfolio's virtual returns: $\widehat{\mbox{VaR}}^{(\alpha)}_{FHS,t-1}(r_t) :=-q_{\alpha}\left(\{\mathbf{a}_{t-1}'\boldsymbol{y}_1,\dots,\mathbf{a}_{t-1}'\boldsymbol{y}_{t-1}\}\right)$ and this estimator coincides with the VHS estimator. These empirical quantiles were computed starting from $t=100$. The naive estimator is the opposite of the quantile of the portfolio's returns: $\widehat{\mbox{VaR}}^{(\alpha)}_{N,t-1}(r_t) =-q_{\alpha}\left(\{\mathbf{a}_{1}'\boldsymbol{y}_1,\dots,\mathbf{a}_{t-1}'\boldsymbol{y}_{t-1}\}\right)$. The spherical method, based on the estimation of $\boldsymbol{\Sigma}$, was computed on the same range of observations. In view of ((ref)), we have $\widehat{\mbox{VaR}}^{(\alpha)}_{S,t-1}(r_t) =\sqrt{\mathbf{a}_{t-1}'\widehat{\boldsymbol{\Sigma}}_{t-1}\mathbf{a}_{t-1}}q_{1-2\alpha}\left(\{\widehat{\boldsymbol{\Sigma}}_{t-1}^{-1}\boldsymbol{y}_1,\dots,\widehat{\boldsymbol{\Sigma}}_{t-1}^{-1}\boldsymbol{y}_{t-1}\}\right)$, where $\widehat{\boldsymbol{\Sigma}}_{t-1}$ is the empirical covariance matrix of the observations $\boldsymbol{y}_1,\ldots,\boldsymbol{y}_{t-1}$. Figure (ref) displays the sample paths of the true conditional VaR as well as the estimated VaRs. It can be seen that the spherical method converges faster to the true value than the FHS=VHS method. Unsurprisingly, the univariate naive method fails to converge to the theoretical conditional VaR based on the full information set. Now we turn to a less artificial setting.

figure[figure omitted — 235 chars of source]
figure[figure omitted — 236 chars of source]
figure[figure omitted — 251 chars of source]

Well-specified multivariate GARCH models

In this section, we simulate the process of log-returns $\boldsymbol{y}_t=\boldsymbol{\Sigma}_t\boldsymbol{\eta}_t$ from the corrected Dynamic Conditional Correlation (cDCC) GARCH model of Aielli (2013). For the multivariate approaches, we estimate the same cDCC-GARCH(1,1) model. For the univariate approaches we estimate GARCH(1,1) models, which are generally misspecified (see Example (ref)).

We consider the minimum variance portfolio variance given by

equation[equation omitted — 273 chars of source]

We simulated $N$ independent trajectories of length $n$ for the cDCC-GARCH(1,1) model. On each simulation, the first $n_1$ observations are used (i) to obtain an estimator $\widehat{\boldsymbol{\vartheta}}_{n_1}$ of the parameters involved in $\boldsymbol{\Sigma}_t$ by the three-step estimator defined in Appendix C of Francq and Zakoian (2018), and (ii) to estimate the quantiles required for the VaR estimator. On the last $n-n_1$ simulations, {\it i.e.} for $t=n_1+1,\dots,n$, we compared the theoretical VaR of the portfolio with the three estimates obtained from the spherical, FHS and VHS methods. We considered portfolios of $m=2$ assets. The different designs, displayed in Appendix (ref), correspond to spherical (designs A-H) or non spherical (designs A$^*$-H$^*$) distributions.

We took $N=100$ independent replications, and $n-n_1=1000$ out-of-sample predictions for each simulation. In each design, we then compared the corresponding $100, 000$ theoretical values of the conditional VaR with their estimates obtained by the spherical, FHS and VHS methods. Denote by $\mathrm{MSE}_S$, $\mathrm{MSE}_{FHS}$ and $\mathrm{MSE}_{VHS}$ the mean square errors (MSE) of prediction of the three methods. Table (ref) displays the relative efficiency (RE) of the spherical method with respect to the FHS and the VHS methods, as measured by the ratios $\mathrm{MSE}_{FHS}/\mathrm{MSE}_S$ and $\mathrm{MSE}_{VHS}/\mathrm{MSE}_S$. It should be underlined that all MSE's are computed with respect to the full information VaR's, which a priori favours the multivariate methods. Let us first briefly compare the two multivariate approaches: in designs A-H (with spherical distributions) the spherical method is as expected more efficient than the FHS method (for Designs C and D, the spherical method can be twice more efficient than the other multivariate method). On the contrary, the bottom panel of Table (ref) reveals that, when the density is strongly asymmetric, the FHS method can be much more efficient than the spherical method.

One question of interest is whether the VHS estimator (targeted to estimate the hybrid VaR) can be used to approximate this full-information VaR. Let us first compare the VHS method to the multivariate approaches in the spherical case. The univariate VHS method is apparently dominated by the multivariate methods in designs A-H but one has to recall that the reference VaR to which the methods are compared is designed for the multivariate framework. Better results are obtained for the VHS method in designs E-H (identically distributed returns) as opposed to A-D, and for designs $\{$A,B,E,F$\}$ as opposed to $\{$C,D,G,H$\}$ (independent returns). This can be intuitively explained as follows. Univariate methods are expected to behave better when the trajectories of the underlying returns are close, which is the case when the two assets have similar dynamic models and are strongly dependent. It is thus non surprising to note that the worst results for the VHS (by comparison with the multivariate methods) occur for designs C-D (independent returns with very different dynamics), and the best results occur with designs E-F (dependent identically distributed returns).

In the case of non-spherical distributions (designs A*-H*) the same conclusions hold: the more dependent the assets and the closer their dynamics, the better results for the VHS. This method clearly outperforms the spherical approach in designs E*-H*. Unexpectedly, the results reveal that the univariate VHS method can even outperform the FHS approach (design E* with $n_1=1000$).

From these experiments, it appears that the accuracy of the approximation provided by the VHS approach is very dependent from the model parameters. However, these first two examples clearly favour the multivariate methods by assuming that the dynamic model is well specified. In the next example, we consider a data generating process that does not belong to the GARCH class, and we will compare the different methods using backtests.

center[center omitted — 3,009 chars of source]

Misspecified GARCH models

We simulated $m$-multivariate factor models, with two GARCH factors of the form $$f_{1t}=\sigma_{1t}\eta_{1t}, \quad f_{2t}=\sigma_{2t}\eta_{2t},$$ where $(\eta_{1t})_t$ and $(\eta_{2t})_t$ are two independent sequences of iid ${\cal N}(0,1)$-distributed random variables. The volatilities follow standard GARCH$(1,1)$ equations of the form $$\sigma_{it}^2=\omega_i+\alpha_if_{i,t-1}^2+\beta_i\sigma_{i-1,t}^2.$$ We took $(\omega_1,\alpha_1,\beta_1)=(1,0.09,0.87)$ and $(\omega_2,\alpha_2,\beta_2)=(0.1,0.7,0.01)$, so that the dynamics of the two factors be quite distinct. The even and odd components of our simulated factor model are respectively of the form $$y_{2k,t}=f_{2t}+e_{2k,t},\qquad y_{2k+1,t}=f_{1t}+e_{2k+1,t},$$ where $(e_{k,t})_t$, for $k=1,\dots,m$, are idiosyncratic independent iid noises with law ${\cal N}(0,0.1^2)$. To obtain a graphical comparison of the VaR estimates, we first simulated a trajectory of size $1, 100$ of the factor model with $m=4$. A crystallized portfolio of composition $(1/m,\dots,1/m)$ at time $t=1$ has been considered. For the multivariate approaches we estimated cDCC-GARCH(1,1) models, while for the univariate approaches we estimated GARCH(1,1) models. The four competing estimators of the 5% $\mbox{VaR}_{t-1}$ at time $t=1001$ were estimated on the basis on the first $1, 000$ simulated values $\boldsymbol{y}_1,\dots,\boldsymbol{y}_{t-1}$. Then, VaR at time $t=1, 002$ was estimated based on the past $1, 000$ simulations $\boldsymbol{y}_2,\dots,\boldsymbol{y}_{t-1}$. We continued this way until we obtained the last VaR estimations at time $t=1,100$. Figure (ref) shows that the estimates obtained by the Spherical, FHS and VHS methods are very close (actually, they are not distinguishable on the figure), whereas the estimates obtained by the naive method behave differently. This can be explained by the fact that the portfolio is crystallized but not static. In other words, even if the portfolio is constituted of an equal quantity of the $m$ simulated assets, the return $r_t$ is not a fixed average of the individuals returns $y_{kt}$ (see Figure (ref)).

To compare the methods by using formal backtests, we considered the same framework of GARCH estimations on rolling windows of length $1, 000$, but the methods have been backtested on a longer period of length $2, 000$. Moreover, in order to obtain a clearcut comparison between the naive method and the VHS method, the composition of the portfolio has to be highly time-varying. We thus simulated portfolios whose composition alternates as follows: we take an equal proportion of the returns of the even assets $\epsilon_{2k,t}$ during a period of length 100, and then we switch to an equal proportion of the odd assets $\epsilon_{2k+1,t}$ during another period of length 100. Table (ref) summarizes the results of the 4 VaR estimation methods for $m=2,4, 8, 100$. This simulation exercise is intensive since 2000 DCC-GARCH models were estimated for each of the two multivariate methods, and 2000 univariate GARCH(1,1) models were estimated for each of the univariate methods. The spherical and FHS methods become rapidly too time consuming when the number $m$ of returns increases, because multivariate $m$-GARCH models have to be estimated. Interestingly, the numerical complexity of the univariate methods does not increase much with $m$, so that Table (ref) reports results on portfolios of $m=8$ and $m=100$ assets for the univariate methods only.

In this table, the column Viol gives the relative frequency of violations (in $\%$), while the columns LRuc, LRind and LRcc give respectively the $p$-values of the unconditional coverage test that the probability of violation is equal to the nominal $5\%$ level, the independence test that the violations are independent and the conditional coverage test of Christoffersen (2003). Conclusions drawn from those backtests, which solely focus on the violations, are that all methods are validated on these experiments. In particular, it is interesting to notice that the naive method does not behave so poorly in terms of backtests. It is necessary to introduce alternative statistics to differentiate the different approaches. From a portfolio's manager point of view, it is interesting to minimize the average VaR (denoted $\overline{\mbox{VaR}}$ in the table) in order to minimize the reserves. For both regulator and manager, it is also important to minimize the amount of violation. The column AV displays the average amount of violation, and the column ES gives the expected shorfall, that is the average loss when the VaR is violated: for each estimator $\widehat{VaR}_t$ of the conditional VaR, let $$\mbox{AV}=\frac{\sum_{t=1}^n-(r_t+\widehat{VaR}_t)\mathbf{1}_{\{r_t<-\widehat{VaR}_t\}}}{\sum_{t=1}^n\mathbf{1}_{\{r_t<-\widehat{VaR}_t\}}}, \qquad \mbox{ES}=\frac{\sum_{t=1}^n-r_t\mathbf{1}_{\{r_t<-\widehat{VaR}_t\}}}{\sum_{t=1}^n\mathbf{1}_{\{r_t<-\widehat{VaR}_t\}}}. $$ These statistics clearly show that the naive approach is inefficient compared to its competitors. With this method, the amount of violation tends to be higher whatever the size $m$ of the portfolio. For these statistics AV and ES, the VHS approach appears comparable to the multivariate methods when comparison is possible, that is when $m$ is not too large. Alternative comparisons are provided by introducing the loss function $$\mbox{Loss}=\frac1n \sum_{t=1}^n-(r_t+\widehat{VaR}_t)(\alpha-\mathbf{1}_{\{r_t<-\widehat{VaR}_t\}}),$$ also considered by Giacomini and Komunjer (2004), and Gneiting (2011) in the context of forecast evaluation. The last column of Table (ref) reports, for each of the three non-naive methods, $p$-values of the Diebold-Mariano (1995) test for the null that the naive method produces the same loss against the alternative that it induces higher loss. The null is rejected in each situation, leading to the same conclusion as before: the naive method is outperformed by its three competitors when $m$ is small, and by the VHS method when $m$ is large.

center[center omitted — 273 chars of source]
center[center omitted — 298 chars of source]
table[table omitted — 2,319 chars of source]

Real data

We start by plotting the returns of an actual crystallized portfolio, showing the same kind of behaviour as the simulated portfolio of Figure (ref). The portfolio is obtained by equally weighting 3 stocks, ADM (Advanced Micro Devices), BFB (Brown-Forman Corporation) and AMZ (Amazon). on the period 1997-05-15 to 2018-09-05 (5363 observations). The composition $\mathbf{a}_{t-1}$ of the portfolio is plotted in Figure (ref). It is seen that the third asset becomes preponderant at the end of the period. The plot of the portfolio's returns in Figure (ref) shows that the changes in the portfolio composition induce apparent non-stationarities. Contrary to the simulated portfolio of Figure (ref), the volatility decreases when the portfolio becomes more concentrated, which is explained by the fact the third asset is also the less volatile.

figure[figure omitted — 260 chars of source]
figure[figure omitted — 263 chars of source]

Next, we will compare the univariate methods on portfolios whose composition is strongly time varying. Then, we will consider portfolios that are regularly rebalanced so that their compositions do not vary too much.

Estimating the conditional VaR of portfolios of US stocks

We now consider portfolios built from a set of $m=49$ US stocks covering 2,489 trading days, from January 4, 1999 to December 31, 2008. The data have been kindly provided to us by S�bastien Laurent, and are described in Laurent, Lecourt and Palm (2016). The top panel of Figure (ref) displays the returns of a crystallized portfolio which was fully diversified at the beginning of the period, {\em i.e} with composition $\mathbf{a}_{0}=(1/m,\dots,1/m)$ at time $t=1$, and for which the number of units of each asset $\mu_{i,t}=1$ is time-invariant. The bottom panel of this figure displays $\max_{i\in\{1,\dots,49\}}a_{i,t}$ as function of $t$. This figure shows that the composition of the portfolio is time-varying, and that the portfolio tends to become more and more concentrated. At the beginning of the period, the return of the crystallized portfolio is an equi-weighted average of the individual returns, but at the end of the period, one of the individual returns tends to have a prominent weight (this individual return is the UNILEVER stock from February 20, 2003 onwards).

figure[figure omitted — 212 chars of source]

Given the large number of assets, we did not implement multivariate approaches to estimate the VaR of this portfolio. Estimating a multivariate GARCH(1,1) model by QML in this setting would require inverting very large correlation/covariance matrices at each step of the optimization algorithm. There exists multivariate approaches--either based on constrained models or using alternatives to the full QML (e.g. the composite likelihood as in Engle, Ledoit and Wolf (2017), or the Equation-by-Equation method of Francq and Zakoian (2016)), or using intraday data (e.g. Boudt, Laurent, Quaedvlieg, and Sauri (2017))--which do not fit into our semiparametric GARCH framework.

Figure (ref) displays the estimates of the 5%-VaR obtained from the Naive and VHS methods. Starting from $t=1, 001$, the estimates are computed from all the previous observations $r_1,\dots,r_{t-1}$. As can be seen from the figure, the naive and VHS methods provide very similar results and the backtests used in the previous section are not able to distinguish them. At first sight, this result is quite surprising, but it can be explained by the fact that the composition of the portfolio varies relatively slowly and, even if the composition is changing, the dynamics does not change drastically because most of the individual returns follow similar GARCH models. The interesting conclusion is that, even if the naive method is not supported by rigorous theoretical results, it may work surprisingly well in practice.

In a second experiment, we considered a portfolio whose composition changes every week, between an equi-weighted average of the stocks MS, F, GM and an equi-weighted average of the stocks CVX, XOM, ECX. At the beginning of the week, the portfolio is thus composed of the same amount of the three assets, then the portfolio is crystallized until the end of the week. The two sets of stocks have been chosen because of their different dynamic behaviors. Results are displayed in Table (ref). Once again, the naive and VHS methods are not much different, but there are some discrepancies in favor of the VHS method (as indicated by the DM test).

center[center omitted — 257 chars of source]
table[table omitted — 988 chars of source]

Comparing naive and VHS methods on rebalanced portfolios

In practice, portfolios are often periodically re-balanced. Intuitively, the naive and VHS approaches should behave similarly in this situation. To verify this intuition, we now build portfolios with the $m=19$ stocks that illustrate Chapter 17 of Boyd and Vandenberghe (2018). The data set covers the period from 2004-01-02 to 2013-12-31 (2517 values). We consider a portfolio which is equally diversified at the beginning of the period ({\it i.e.} we take $\mu_{i,t}=\mu_{i}=V_0/(mp_{i,0})$ at $t=0$, corresponding to the beginning date 2004-01-02), and we re-balance the portfolio every $T$ periods (such that $\mu_{i,t}=V_t/(mp_{i,t})$ at $t=kT$ for all $k\in \mathbb{N}$). If the portfolio is re-balanced at any time, {\it i.e.} $T=1$, then $a_{i,t}=1/m$ for all $t$, and thus the naive and VHS methods coincide. To limit transaction costs, the investor can maintain the same asset allocation for an extended period of time. The so-called lazy portfolios or permanent portfolios are re-balanced every year or at any time its asset allocation strays too far from its initial state. Figure (ref) shows that, when the portfolio is re-balanced every year, the composition of the portfolio can deviate much from the fully diversified portfolio (for which $a_{i,t}=1/m$ for all $i\in\{1,\dots,m\}$). This does not necessarily entail a huge difference between the VaR's estimated by the naive and VHS methods. Actually, for the nominal risk level $\alpha=1\%$, the maximum difference between the estimated VaR has been observed on 2008-12-19, with a naive VaR of 0.09163292 and a VHS VaR of 0.09775444. As illustrated by Figure (ref) the two estimated VaR are hardly distinguishable. This is interesting because it entails that the asymptotic theory built for the VHS method should also apply for the naive method. In other words, the naive method is not so naive if the portfolio is re-balanced from time to time, as recommended by finance professionals.

figure[figure omitted — 311 chars of source]
figure[figure omitted — 239 chars of source]

Conclusion

This paper developed a method for estimating the conditional VaR of a portfolio of asset returns, without relying on a joint dynamic modelling of the vector of returns. For large portfolios, using optimization routines in multivariate approaches often entails formidable numerical difficulties. By circumventing the dimensionality curse, univariate methods provide an operational alternative.

The naive method discussed in this paper has no theoretical grounds because it implicitly, and erroneously when the composition is time varying, relies on the stationarity of the returns process. In many cases, however, it behaves satisfactorily as our numerical experiments revealed. For the VHS method, we developed an asymptotic theory for a general class of dynamic models, which are not directly estimated on observations but rather on reconstituted returns. The obtained asymptotic results allow to quantify the estimation risk that should be taken into account in risk management. From our numerical experiments, the multivariate methods can be recommended when the size of the portfolio is small and the estimated multivariate GARCH model is likely to be well-specified. When the number of underlying assets is large, or when finding an appropriate multivariate specification is difficult, the univariate methods offer a valuable alternative. The VHS method requires computing virtual returns, which is negligible computational burden. Thus, the two univariate methods display similar numerical complexity. However, when the naive and VHS estimators differ substantially, the former should not be considered as reliable. When they provide similar results, the asymptotic results obtained for the VHS method could also be used to evaluate the accuracy of the naive method.