Introduction
My goal is to retrieve a single Value-at-Risk (VaR) of a N(0, H) random variable $X$ at the $\alpha \in (0,1)$ confidence level where H is a known d-dimensional positive definite matrix (with $d >3$) In other words, I aim to retrieve the $\alpha$ quantile of a N(0, H) random variable $X$ , i.e.
$$\text{VaR}_{\alpha}(X) = \text{inf} \ \{x \in \mathbb{R}: F_X(x) > \alpha\}$$
I consider 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. I do not wish to consider a reduced model in which the portfolio's returns are considered as a single univariate time series, i.e. N(0, $\sigma^2$), as this ignores complicated correlation among the individual returns. In this context VaR is the maximum loss with 1-$\alpha$ probability. The idea is to use VaR estimates in a backtesting analysis (Kupiec test, Christoffersen Test) 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 α 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 α 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[0], 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)
