# How to determine the cut off value of an hyperellipsoid in order to retrieve a single quantile of a multivariate normal distribution?

Introduction

My goal is to retrieve the $\alpha$ quantile of a N(0, H) (multivariate normal) random variable $X$ where H is a known d-dimensional positive definite matrix (with $d >3$). In other words, my goal is to retrieve a single Value-at-Risk (VaR) of X at the $\alpha \in (0,1)$ confidence level, i.e.

$$\text{VaR}_{\alpha}(X) = \text{inf} \ \{x \in \mathbb{R}: F_X(x) > \alpha\}$$

I aim to retrieve the $\text{VaR}_\alpha$ for $\alpha \in \{0.1, 0.05, 0.01\}$ in the context of a portfolio comprising $d>3$ financial log returns where we assume that the portfolio log return follows a multivariate normal distribution, i.e. N(0, w^T H w). We define w as an equally weighted vector of length $d$ meaning that each asset is equally represented in our portfolio. In this context VaR is the maximum loss with 1-$\alpha$ probability. Given the definition of VaR, VaR is expected to increase with a decrease in the $\alpha$. The idea is to use VaR estimates in a backtesting analysis in order to evaluate different models used for forecasting the covariance matrix H.

Approach

In case of unimodal symmetric distributions such as the multivariate normal distribution considered in the question, quantiles correspond to Highest Density Regions. As pointed out here, the highest density region of an N(0,H) random variable is an ellipsoid centered at its mean, 0, and oriented per the covariance matrix H:

$$x: x^TH^{-1}x \le y$$

In order to find a single quantile that satisfies the definition of Value-at-Risk we proceed as follows: the set of solutions that satisfy the condition of the ellipsoid equation, the level set, is retrieved by randomly generating many points on this ellipsoid. Then, I evaluate the portfolio return at each solution satisfying the condition: for example, let's suppose $s = (x_1, x_2, \dots, x_d)$ is a solution/ point on the ellipsoid, then the portfolio return at this point will be $w\cdot s$, which corresponds to the mean of the solution vector because we stated that each asset is equally represented in the portfolio and $x_i \in s$ corresponds to the return of asset $i$. Consequently, the lowest return, the minimum, on the ellipsoid is the maximum loss.

Question

Now, in order to ensure that the maximum loss on the ellipsoid is the maximum loss with 1-$\alpha$ probability, hence is the Value-at-Risk of X at the $\alpha$ confidence level, I need to define the cutoff value for the ellipsoid, i.e. $y$. At this point, I am a bit confused concerning the definition of the cutoff value for the ellipsoid.

Here it is stated that the cutoff value for the ellipsoid can be determined from the Chi-square with d degrees of freedom. For instance, the highest density region capturing 0.95 probability of the N(0,H) is found with y= value such that Chi-square with d degrees of freedom 0.95. When defining the cutoff value from the chi-square with d degrees of freedom and following the procedure explained in the approach section, I would interpret the lowest return on the ellipsoid as the maximum loss with 0.95 probability, i.e. $\text{VaR}_{0.05}$. Is it correct to state that these maximum loss values are then the maximum loss with 1-α probability, hence is the Value-at-Risk? VaR does increase with a decrease in $\alpha$ but at the same time VaR seems a bit a large for the context which makes me doubt my interpretation and definition of the cutoff value for the ellipsoid determined from the Chi-square with d degrees of freedom and its link to the $\alpha$ quantile of the distribution (Value-at-Risk).

Any comments that could point me in the right direction are more than welcome.

Python code for obtaining random points on the ellipsoid

import numpy as np
from scipy.linalg import sqrtm
from scipy.stats import ortho_group
from scipy.stats._continuous_distns import chi2

dim = 10
# Create a positive definite matrix H with shape (dim, dim).
# evals is the vector of eigenvalues of H.  This can be replaced
# with any vector of length dim containing positive values.
evals = (np.arange(1, dim + 1)/2)**2
# Use a random orthogonal matrix to generate H.
R = ortho_group.rvs(dim)
H = R.T.dot(np.diag(evals).dot(R))

# y determines the level set to be computed.
y = chi2.ppf(q=0.95, df=dim)
# Generate random points on the ellipsoid
nsample = 100000
r = np.random.randn(H.shape, nsample)
# u contains random points on the hypersphere with radius 1.
u = r / np.linalg.norm(r, axis=0)
xrandom = sqrtm(H).dot(np.sqrt(y)*u)
# Compute maximum loss on the ellipsoid
xrandom_min = np.min(np.array([np.mean(x) for x in xrandom.T]))
print("Maximum loss i.e. Value-at-Risk:")
print(xrandom_min)

• @MartijnWeterings The definition of value at risk is the maximum loss with alpha percent probability. I try to obtain this from the set of (indeed not-unique) boundaries of the hyperellipsoid that captures alpha percent of the multivariate normal variable N(0,H). For every boundary, I define the loss as the mean of the coordinates specifying the boundary point. Then I take the minimum (maximum loss) of this set of losses and my question is then whether I can state if that maximum loss is the maximum loss with alpha percent probability, hence the value-at-risk definition. Hope its more clear – Paul Jul 16 '18 at 12:45
• Those hyper-ellipsoids refer to cases $x_1,x_2,...x_d$ of equal probability density, or equal Mahalanobis distance. They do not refer cases of equal loss. – Sextus Empiricus Jul 16 '18 at 14:11
• I understand that the hyper-ellipsoid refers to cases $x_1, x_2, \dots, x_d$ of equal probability, or equal Mahalanobis distance as you say. The intuition I was after was that when I find the minimum return, as specified below $w^t x$, of cases of equal Mahalanobis distance that covers 95% percent of the data. Is it then correct to state that that minimum return is the maximum loss with 95% probability (as data outside this hyper-ellipsoid will cover the remaining 5% of the data). – Paul Jul 16 '18 at 14:18
• No, you need to sample the 95% of samples that optimize the loss, not the 95% of samples that optimize the Mahalonobis distance. These two will cover two entirely different regions, or at least shapes (your loss isosurfaces, defined by $w^tx=c$ are hyperplanes), thus your Mahalobis ellipsoid should contain some points outside the region of the 95% samples with optimized loss. Therefore the calculated value in this way is lower than the VaR. – Sextus Empiricus Jul 16 '18 at 14:29
• I tried to give as much information to support the question (background and possible approach) without intent to confuse. What it boils down to is this: A mean vector of zero and covariance matrix H describe the first two moments of a set of correlated returns. I am trying to find the lower quantiles of this return distribution = higher quantiles of loss distribution explicitly taking into account the correlation between individual returns. I believe without any distributional assumptions this will be very hard therefore the normality assumption. – Paul Jul 16 '18 at 14:52

It seems that you made an error very early in your story when you say

where we assume that the portfolio log return follows a multivariate normal distribution, i.e. N(0, w^T H w).

That distribution is a univariate normal distribution and not a multivariate normal distribution.

## What you should use is

If $$R=w^T x$$ is your loss/return, then this means that you are looking for the distribution of a variable that is the linear combination of the components of the multivariate normal distributed variable

The sum of any linear combination of the components of a multivariate distributed normal distribution is a univariate normal distributed variable (it is one way how you can define a multivariate normal distribution).

$$R \sim N(0, w^THw)$$

for background about this $$w^THw$$ term see Matrix notation for the variance of a linear combination

Now you can just look for the percentiles of this univariate normal distribution.

So you look for $$\text{VaR}_\alpha(w^TX)=\text{VaR}_\alpha(R)$$ instead of $$\text{VaR}_\alpha(X)$$

## Why and how you get a difference:

• You need to look at the minimum of the 95% of samples that optimize the loss.
• Instead you use the minimum of the 95% of samples that optimize the Mahalanobis distance.

These two will cover two entirely different regions, or at least shapes (your loss isosurfaces are defined by $$w^T x=c$$, which are hyperplanes).

Thus your Mahalanobis ellipsoids should contain some points outside the region of the 95% samples with optimized loss. Therefore the calculated value in this way is lower than the VaR. Because inside the 95% ellipsoids you include some of 5% point beyond the 5% hyperplane with lowest cost.

### Graphical example (comparing two 95% boundaries)

Below is an image that demonstrates the principle (in two dimensions but the principle remains the same for more dimensions). Thousand datapoints are simulated and drawn for $$H=\begin{bmatrix} 4 & 1 \\ 1 & 2 \end{bmatrix}$$. The colours depict varying return which is calculated as $$w \cdot x$$, and with equal weights this is $$x_1+x_2$$.

You see that the the boundary around the 95% points with highest density, is something different than the boundary for the 95% with the highest return (or lowest loss). • thanks @Martijn for showing the difference between the boundary around the 95% points with highest density the boundary for the 95% with the lowest return (or highest loss). That answer the question whether it is correct to state that the maximum loss value on the hyperellipsoid coincides with the maximum loss with 1-α probability, hence is the Value-at-Risk: It is not. – Paul Jul 17 '18 at 9:56

If you're just after VaR, one way to calculate VaR is the "Delta Normal" method.

If $\Phi^{-1}$ is the inverse CDF function for a standard normal random variable, $\mathbf{w}$ is your nonrandom weight vector, and $\mathbf{H}$ is the covariance matrix of the future returns, then $$\text{VaR}_{\alpha}(X) = \Phi^{-1}(\alpha)\sqrt{\mathbf{w}^T \mathbf{H} \mathbf{w}}.$$

You can see the probability of your return being less than this lower bound is $\alpha$: \begin{align*} P\left(R \le \Phi^{-1}(\alpha)\sqrt{\mathbf{w}^T\mathbf{H} \mathbf{w}}\right) &= P\left(R(\mathbf{w}^T\mathbf{H} \mathbf{w})^{-1/2} \le \Phi^{-1}(\alpha)\right)\\ &= \Phi \left( \Phi^{-1}(\alpha)\right) = \alpha. \end{align*}

This does not require simulation, and the only place where matrices and vectors come into play is when you take the variance of the scalar return: $\operatorname{Var}(R) = \operatorname{Var}(\mathbf{w}^T \mathbf{X}) = \mathbf{w}^T \operatorname{Var}(\mathbf{X}) \mathbf{w}= \mathbf{w}^T \mathbf{H} \mathbf{w}$.

Here's some code:

import numpy as np
from scipy.stats import norm

def compute_var(alpha, cov_mat, weight_arr):
variance = np.dot(np.dot(weight_arr, cov_mat),weight_arr)
std_dev = np.sqrt(variance)
return norm.ppf(alpha)*std_dev

weight = np.repeat(1,3)/3.0 #equal weights
Sigma = np.array([[1.0,.5, .5],[.5,1.0,.5],[.5,.5,1.0]])
compute_var(.05, Sigma, weight)