Computing Highest Density Region given multivariate normal distribution with dimension $d$ > 3

Background:

Suppose I have a unimodal symmetric distribution of dimension $d$ > 3 such as the multivariate normal distribution ~ N(0, H), where H is a known $d$-dimensional covariance matrix.

Question:

I am interested in finding the Highest Density Region that covers 95% of the area under this distribution. Now, I understand the notion of Highest Density Regions in the univariate and bivariate case and how to proceed for obtaining these using for example the hdrcde package. However, I am interested in the following: How do I obtain the Highest Density Region given the distributional assumption of multivariate normality and a $d$-dimensional covariance matrix H where $d$=10 for example?

Idea how to proceed:

I am aware that in case of unimodal symmetric distributions such as the multivariate normal distribution considered in the question, the $\alpha$% HDR coincides with the $\alpha$% quantiles of the distribution. I could use a grid-based algorithm to approximate these quantiles for instance. The idea would be to define M grid points for each dimension such that relevant space is covered. This would would lead to M^$d$ joint grid values for $d$-dimensional distribution. Then I would evaluate the pdf of the multivariate normal density at each joint grid point. And then from this M^$d$ grid of densities I am not sure how to proceed in order to obtain $\alpha$-quantiles.

Any suggestions how to proceed are very welcome. Any other approach how to obtain High Density Regions for a $d$-dimensional normal distribution with $d>2$ are also more than welcome as I think M^$d$ function evaluations might become slow in higher dimensions. Thanking you in advance.

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$. The cutoff value for the ellipsoid can be determined from the Chi-square with $d$ degrees of freedom.

Let $y =$ value such that Chi-square with $d$ degrees of freedom $\le 0.95$. Then the highest density region capturing $0.95$ probability of the $N(0,H)$ is
$$x: x^TH^{-1}x \le y$$

For instance, when $d = 10, y = 18.307038$.

• Thank you for answer. There remains some confusion wrt determination of the cut off value which I have formulated in a detailed context here. Perhaps you could point me in the right direction? Many thanks for your time.
– Paul
Jul 15 '18 at 12:55
• The $-2log(\alpha)$ applies only to d =2, and is the chisquare result when d = 2. There is no inconsistency there. If they don't come out the same for you when d = 2, then you are making a mistake somewhere. Also note that boundaries of ellipsoid can be determined by going sqrt(y)*sqrt(eigenvalue) along corresponding eigenvector of H from the origin. Jul 15 '18 at 13:58
• -2log($\alpha$) and chisquare indeed come out the same for me when d=2. No inconsistency there. I also check all randomly solutions for x if they satisfy the condition $x^T H^{-1} x = \chi^2(\alpha)$. I am just not completely sure whether my interpretation of the minimum mean of all obtained solutions corresponds to the maximum loss (value-at-risk) with $\alpha$ percent probability. Such that for instance for y=chisquare(0.95) the obtained minimum mean solution corresponds to the maximum loss with 95% probability, i.e. VaR$_{0.05}$. At this point, I get pretty low minima for the context.
– Paul
Jul 15 '18 at 20:04
• Leaving aside whether you're making a mistake or misinterpretation, the real issue is that log returns aren't Normal. You're likely way underestimating true VaR, especially for very small $\alpha$, because actual log returns have fatter tails than Normal. VaR based on Normal might be o.k. for $\alpha = 0.2$, as I've seen used sometimes, which gets to smoothness of earnings or returns under nice conditions, but does not provide information pertaining to existential type losses, in which the tail of the distribution is paramount. VaR based on Normal with $\alpha = 0.05, 0.01, 0.001$ is bogus. Jul 15 '18 at 20:30
• I fully agree that log returns generally exhibit fat tails and a skewness that will not be adequately captured by assuming potfolio returns to follow a multivariate normal distribution. However, the assumption of normality is in the context of obtaining VaR estimates which are used for backtesting (kupiec test, christofferson markov test) different models used for estimating time-varying covariance matrix H. The peculiarity is that I am actually overestimating true VaR (overly conservative models) at this point. Hence, the reason I want to be sure whether my interpretation is correct.
– Paul
Jul 15 '18 at 20:40