Modeling Multivariate Financial Series and Computing Risk Measures via Gram–Charlier-Like Expansions

: This paper develops an approach based on Gram–Charlier-like expansions for modeling ﬁnancial series to take in due account features such as leptokurtosis. A Gram–Charlier-like expansion adjusts the moments of interest of a given distribution via its own orthogonal polynomials. This approach, formerly adopted for univariate series, is here extended to a multivariate context by means of spherical densities. Previous works proposed the Gram–Charlier of the multivariate Gaussian, obtained by using Hermite polynomials. This work shows how polynomial expansions of an entire class of spherical laws can be worked out with the aim of obtaining a wide set of leptokurtic multivariate distributions. A Gram–Charlier-like expansion is a distribution characterized by an additional parameter with respect to the parent spherical law. This parameter, which measures the increase in kurtosis due to the polynomial expansion, can be estimated so as to make the resulting distribution capable of describing the empirical kurtosis found in the data. An application of the Gram–Charlier-like expansions to a set of ﬁnancial assets proves their effectiveness in modeling multivariate ﬁnancial series and assessing risk measures, such as the value at risk and the expected shortfall.


Introduction
There are plenty of examples in the literature on the non-normality of asset returns and its implications for pricing and measuring financial risk. Several papers (see, e.g., Boudt et al. 2007;Jondeau and Rockinger 2005, 2006a, 2006bJurczenko et al. 2004;Mandelbrot 1997) have analyzed the consequences arising from the misspecification of models related to the normality assumption and the importance of correctly accounting for stylized facts that are typically exhibited by financial series, such as leptokurtosis, volatility clustering, and asymmetry. They have also highlighted the convenience of modeling the joint portfolio distribution by assuming multivariate leptokurtic specifications. The attraction of referring to heavy tailed distributions for modeling and assessing risk measures is also well documented in the literature (Eun et al. 2004;Harmantzis et al. 2006;Marinelli et al. 2012;Rachev 2003).
Leptokurtic distributions can be easily obtained by reshaping a given law through its own orthogonal polynomials. The resulting distribution is called Gram-Charlier (GCE) or Gram-Charlier-Like Expansion financial series and in computing risk measures, such as the value at risk and the expected shortfall, is investigated. Section 4 draws some conclusions. Appendix A provides the proofs of the methodological results presented in Section 2. Zoia (2010) and Bagnato et al. (2015) proved that the kurtosis of a given univariate even density f (x) can be adjusted via a binomial of the form:

Spherical Extensions of Gram-Charlier-Like Expansions
where β is a positive coefficient and p 4 (x) is the fourth-degree orthogonal polynomial: with coefficients b j , j = (0, 2, 4) specified in terms of the (even) moments m j of f (x) as follows: Under the condition: the binomial (1) is positive, and the Gram-Charlier-Like Expansion (GCLE): is a density whose fourth moment, m 4 , is increased by a quantity equal to β, with respect to the same moment m 4 of its parent density f (x), that is: while its lower order moments remain unchanged (see Bagnato et al. 2015). This approach, which has been successfully applied to univariate distributions, is here extended to spherical densities.
Spherical distributions can be viewed as multivariate generalization of univariate symmetrical densities and prove useful in several financial studies (see, e.g., Landsman and Valdez 2003;Szegö 2004).
As is well known, any spherically distributed random vector x admits a stochastic representation of the form: where R = ||x|| = (x x) 1/2 is a non-negative random variable, called the generating variate, independent of U, which is uniformly distributed on the (n − 1)-dimensional unit hyper sphere. In Equation (7), U produces elliptically contoured surfaces, while R determines the shape and, in particular, the kurtosis of the distributions, as shown in the next theorem, which introduces some results concerning these variables that play a crucial role in the forthcoming analysis.
Theorem 1. The density of an n-dimensional spherical variable x is: Here: k n = Γ(n/2) 2π n/2 (9) where Γ(·) denotes the gamma function and g(·) is a non-negative Lebesgue measurable function, named the density generator, with finite Mellin transform (see Zayed 1996, p. 201): The density generator is related to the density function of R as follows: where r = (x x) 1/2 . Mardia's kurtosis index (Mardia 1970) of the spherical vector x, K SP , is given by: where m 2j is the j-th moment of R 2 : that tallies with the 2j-th norm-moment of the spherical density g n (x): with dx standing for ∏ n i=1 dx i .

Proof. See Appendix A.
The next theorem shows how to obtain GCLEs of spherical laws and provides their Mardia kurtosis indexes. The latter highlight the rise of kurtosis, which is introduced in spherical laws by the polynomial expansion.
Theorem 2. A GCLE, based on the fourth-degree orthogonal polynomial, of an n-dimensional spherical variable x with density function g n (x) is defined as:g Here, q(x, β) is the binomial defined in (1) with p 4 (x) specified as: with coefficients, b 2 and b 4 , determined as in (3) from the (even) norm-moments m 2j of the spherical distribution defined as in (14): Mardia's kurtosis index of the GCLE as defined in (15) is: Proof. See Appendix A.
Spherical distributions can be built from univariate laws. In the following, we focus on spherical distributions ensuing from the so-called Power-raised Hyperbolic Secant densities (PHS). This is a family of bell-shaped symmetric univariate distributions arising from the hyperbolic-secant law raised to a positive power λ, which in standard form can be specified as follows: where: In (20), ψ (1) (·) is the trigamma function and B(·, ·) is the beta function. It can be proven that this family embraces both the Gaussian law, for λ → ∞, and the Laplace law, for λ → 0, as limit distributions (see Faliva and Zoia 2017). The logistic and hyperbolic-secant distributions arise when λ = 2 and λ = 1, respectively. Thus, this family includes distributions whose kurtosis ranges from three (corresponding to the case of the Gaussian law) to six (corresponding to the case of the Laplace distribution).
The next theorem explains how to build Spherical laws from PHS densities (SPHS).
Theorem 3. An n-dimensional SPHS engendered by the power-raised hyperbolic-secant family has a density of the form: where k n is as in (9) and M(n, λ) is the (finite) Mellin transform: Mardia's kurtosis index of an SPHS variable is given by: In what follows, we consider some SPHS laws that are characterized by different kurtosis levels, namely the Gaussian, the logistic, and the hyperbolic-secant law.
Corollary 1. The n-variate spherical hyperbolic-secant has a density specified as follows: Here, (n) k is the Pochhammer symbol and: with ξ(·, ·) denoting the Hurwitz zeta function. The n-variate spherical logistic has a density specified as: where: with ζ(·) denoting the Riemann zeta function. Finally, the n-variate spherical Gaussian distribution has a density specified as follows: Proof. See Appendix A.
Spherical laws built from PHS distributions are, besides the multinormal, leptokurtic distributions. Even so, the Mardia kurtosis index that characterizes these distributions does not always resemble the value of the empirical kurtosis typically exhibited by multivariate financial series. As happens in the univariate case, the target of increasing the kurtosis index of SPHS laws can be accomplished via their polynomial expansions. The upcoming theorem is concerned with the problem of building GCLEs of SPHS laws and the kurtosis increase attainable with GCLE distributions. The corollary that follows provides closed-form expressions for the GCLE of the Gaussian, the logistic, and the hyperbolic-secant laws considered in Corollary 1.
Theorem 4. The density of a GCLE of the SPHS law as defined in (21) is: where q(x, β) is as in (16) with coefficients determined as in (3) from (even) moments, m 2j defined as: where M(.) is as in Equation (22), with a = 1 for the Gaussian case. Mardia's kurtosis index of the GCLE, as specified in (29), is given by: Proof. See Appendix A.

Data Description
The aim of this section is to evaluate the performance of the GCLEs in both fitting leptokurtic financial series and quantifying risk measures such as the Value at Risk (VaR) and the Expected Shortfall (ES). To this end, we considered a set of daily indexes: the Hang Seng Index (HSI), the 5-year U.S. Treasury Yield Index (FVX), and the CAC40index (FCHI), observed over a period that goes from 1 January 2009 to 31 December 2016, so as to cover the 2008 financial crisis and its aftermath; the data refer to three different geographical areas representing three different markets: the Chinese, the American, and the European one.
The dataset, including T = 1903 observations for each index (excluding missing cases), was split into two sub-samples of size T 1 = 952 and T 2 = 951, respectively. The former was employed to evaluate the goodness of fit of the GCLEs and their in-sample capability in assessing VaR and ES. The latter was used to evaluate the out-of-sample performance of said distributions in estimating the above risk measures.
For each series, the returns from data were computed as the log-ratio between daily closing prices at time t and t − 1. The preliminary statistics on the returns in the whole sample period are shown in Table 1. All the kurtosis indexes exceed three significantly, while the Jarque-Bera normality test statistics (Jarque and Bera 1987) are significant at the 1% level, thus providing evidence that the above returns are not normally distributed. Moreover, all the Ljung-Box Q 2 (20) statistics for the squared returns (Ljung and Box 1978) are significant at the 1% level. Finally, the D'Agostino test (D'Agostino 1970) signals skewness only at the 5% significance level for FCHI and FVX. As the returns are only mildly affected by skewness, the use of spherical distributions is suitable for their modeling.

Model Estimation and Fit
This sections explains how GCLEs can be fit to the empirical data described above. First of all, we modeled the three mentioned returns with a trivariate spherical law. More precisely, we considered a trivariate Gaussian, logistic, and hyperbolic-secant density.
According to Equations (24), (26) and (28) in Section 2, the spherical Gaussian, logistic, and hyperbolic-secant laws, denoted by g N , g L , and g HS , are given by: Taking into account Equations (32)-(34) in Corollary 2, the GCLE versions of these distributions, denoted by GCN, GCL, and GCHS respectively, are: Here, β is the excess-kurtosis of the empirical trivariate series with respect to the (Mardia) kurtosis index (see Equation (12) in Appendix A) of the spherical Gaussian, logistic, and hyperbolic-secant laws, given in Equations (35)-(37).
By denoting withK the kurtosis of the empirical trivariate series and with K the (Mardia) kurtosis index of the spherical Gaussian, logistic, and hyperbolic-secant law that have been fit, in turn, to the empirical series, the excess kurtosis β can be estimated via the Method of Moments (MM) as follows: Alternatively, Maximum Likelihood (ML) estimation can be employed to estimate µ, Σ, and β from the data, numerically maximizing (via the optim function available in the R software, (R Core Team 2020)) the log-likelihood: whereg(.) is one of the trivariate densities in (38)-(40) and z = Σ −1/2 (x − µ), thus obtaining the relatedβ ML estimate. Table 2 reports the estimates of β, the AIC information criterion values for the GCLEs, together with those of the multivariate Student t (Mt; Kotz and Nadarajah 2004) and the power exponential (MPE; Gómez et al. 1998) distributions, which have been fit to the series for comparison reasons. Furthermore, for the GCLEs, the levels of kurtosis of each parent distribution (K SPHS ) and of the associated expansions (K SPHS GCLE ) are also reported. For the MM distribution, the latter clearly corresponds to the empirical Mardia's kurtosis index of the data.
Looking at the AIC results, when these criteria are applicable, we conclude that GCLEs fit better the empirical series than the other distributions, and this is particularly true for the GCHS. In order to better investigate the goodness of fit of these polynomially adjusted densities in comparison with the popular Mt and MPE distributions, the multivariate version of the Wilcoxon rank test (Liu and Singh 1993) was also employed. To implement the test, the following procedure was devised: • An artificial population of one million data points was generated from each distribution, using the β estimates in Table 2. For the GCLE distributions, since no standard routine is available, the Hamiltonian Monte Carlo (HMC; Neal 2011) was employed; • B = 1000 samples were extracted from each population. The MW test was performed via the R package DepthProc (Kosiorowski and Zawadzki 2014) on each of the samples coming from each population. This yields a B-dimensional vector containing the binary outcome of the tests (acceptance = 0/rejection = 1); • A non-parametric bootstrap procedure was then performed on the above vector, to obtain a confidence interval on the proportion of rejections. The higher the proportion of rejections is, the worse the performance of the underlying model. Table 3 shows the average p-value of the MW test across the B replicates, along with the proportion of rejections (p R ) and its 95% confidence interval (p l , p u ), for each model. According to the test, a p-value higher than the significance threshold provides evidence that the two datasets come from the same underlying population. In our situation, one dataset represents a simulation from a given model (thus, it can be regarded as a proxy of the model itself), while the other represents the empirical data. This entails that the models with an average p-value above the significance threshold have a sufficiently adequate fit to the data. The results provide evidence that GCLEs perform well compared to Mt and MPE, with the only exception of the GCN, since they are the only distributions with an average p-value above the 10% level. Looking at the results shown in the other columns of the table, we can see that GCLEs show also a much lower rate of rejection of the null hypothesis of the MW test across the B replications, except for the GCN. Furthermore, looking at both the average MW p-value and the rejection rate, we conclude that the GCHS distributions obtained via the method of moments estimation also adequately fit the data.

Evaluation of VaR and ES via GCLEs
So far, we have evaluated the capability of the GCLEs to fit financial series. In the following, the validity of these densities in computing risk measures such as the value at risk and the expected shortfall is also investigated (see Chen 2018;McNeil et al. 2005), both in-and out-of-sample, for a linear portfolio: To this end, the entire observation period, which ranges from 1 January 2009 to 31 December 2012 (1903 observations), was split into two sub-periods having approximately the same length: T 1 , going from 1 January 2009 to 31 December 2012 (952 observations), and T 2 , going from 1 January 2013 to 31 December 2016 (951 observations).
The observations of the period T 1 were used to obtain parameter estimates of the GCLEs, the Mt, and the MPE law. Then, the latter were employed to estimate risk measures, such as the value at risk and the expected shortfall.
Given the spherical nature of the mentioned distributions, VaR and ES were estimated by using the approach proposed by Kamdem (2005). The value at risk at a given level α, VaR α , of a linear portfolio described by an n-dimensional elliptical random variable z with vector mean µ and covariance matrix V, was evaluated as: Here, δ is a vector of weights, Σ = n m 2 V with m 2 as defined in (17), and q α,n is the unique positive solution of the following equation: In (45), g n (·) is the density of the vector z, and Γ(·) denotes the gamma function. In the case of GCLEs of spherical variables and by setting equal weights δ k = 1 n ∀k, Equation (44) simply becomes: where q α,n is the solution of (45) obtained by replacing g n (u) withg n (x, β) defined as in (15). Following Kamdem (2005) and Nesmith et al. (2017), the expected shortfall at level α, ES α , is given by: where the symbols are defined as in (44)-(46). In the case of GCLEs of spherical variables, by setting equal weights in δ and by replacing g n (u) in (47) withg n (u, β) as defined in (15), the above equation becomes: Equation (48) was used to compute ES α by using VaR α estimated as in (44). In particular, GCLEs as defined in (38)-(40) were employed to compute ES α as given in (48). The empirical counterparts of VaR and ES, VaR α,emp and ES α,emp , respectively, of the linear portfolio P z are given by: whereF P z is the empirical cumulative distribution function of the portfolio defined as in Equation (43). α,emp at the 95% confidence level are also reported, along with the absolute difference between VaR α and VaR α,emp . These intervals were obtained by extracting 1000 bootstrap samples from the whitened data of the first period. To ensure preservation of the structure of the data, the maximum entropy bootstrap of the method of (Vinod 2006) was employed, via the R package meboot (Vinod and López-de Lacalle 2009). A distribution whose VaR α falls inside the 95% confidence interval of VaR α,emp can be regarded as adequate for this risk measure. Looking at the results shown in the table, we see that GCLEs manage to fall inside the confidence intervals across all α's. In particular, the GCN is adequate for the two highest risk levels, while the GCL and the GCHS are adequate for the two lowest risk levels. Neither the Mt nor the MPE follow a particular pattern and tend to underestimate the risk for α = 0.025 and α = 0.05. The only distribution that is able to estimate VaR α,emp satisfactorily across all levels of α is the GCHS estimated via the method of moments. Finally, looking at the results reported in the last column of the table, we can conclude that the GCHS distributions provide the closest fit to VaR α,emp for α = 0.01 and α = 0.025, while the MPE the closest fit for α = 0.05 and the Mt the closest fit for α = 0.1, immediately followed by the GCHS.
A similar analysis was carried out for the expected shortfall. In particular, the last two tables highlight how the most leptokurtic distribution, the GCHS, has a great advantage in estimating losses associated with high risk levels, thanks to a tail heaviness that is more pronounced relatively to the other GCLEs The out-of sample performance of the CGLE in estimating the value at risk was evaluated by using Kupiec's Likelihood Ratio test (LR;Kupiec 1995) and the Binomial Test (BT; Lee and Su 2012). Both the LR and BT null hypotheses state that the percentage of portfolio losses that in the second part of the sample exceed VaR α is equal to α. Since VaR α is estimated via GCLEs in the first period of the sample, this analysis should provide evidence of the out-of-sample stability of the GCLEs in assessing this risk measure.
The p-values of both the LR (p-value LR) and binomial test (p-value BT) are shown in Table 6 together with two loss measures: the Average Binary Loss Function (ABLF) and the Average Quadratic Loss Function (AQLF). The former is the average of the Binary Loss (BL) function, which gives a penalty of one when each return r t exceeds the VaR, without concern for the magnitude of the excess: If a VaR model truly provides the level of coverage defined by its confidence level α, then the ABLF should be as close as possible to the latter. The AQLF (Lopez 1999) is the average of the Quadratic Loss (QL) function: and pays more attention to the magnitude of the violations, as large violations are more penalized than the small ones. Table 6 shows the results obtained by performing the LR and the BN test. They confirm the good performance offered by GCLEs, with the exception of the GCN. Looking at the loss functions, we can conclude that the GCHS outperforms all the other distributions for α = 0.025. In particular, the GCL and the Mt are equivalent for α = 0.01; the GCN (ML) and the GCHS (MM) provide the best fitting for α = 0.05 (the former for the ABLF, the latter for the AQLF); and the GCN is the only distribution not underestimating the ABLF for α = 0.1, which is better approximated by the MPE according to the AQLF.
Finally, the out-of-sample GCLE performance in estimating the ES was evaluated by performing the Z 1 and Z 2 tests of Acerbi and Szekely (2014). The null hypothesis of both tests assumes that the estimated ES α tallies with the ES calculated from the data. This entails, in the present application, that ES α should provide a good estimate of the empirical expected shortfall computed in the second part of the sample. Under the alternative, this does not happen, and ES α systematically underestimates the effective loss mean, ES emp , thus implying financial damage.
The first test is based on the following statistic: where N T = ∑ T t=1 I t is the number of losses L t , which in the second part of the sample, break the threshold VaR α , and I t is an indicator variable that is equal to one if L t > VaR α and zero otherwise.
As observed by Acerbi and Szekely (2014), this test is completely insensitive to an excessive number of exceptions as it is simply an average taken over the exceptions themselves. Under the null, the realized value Z 1 is expected to be zero, and it signals a problem when it displays values significantly greater than zero.
The second test statistic is: where T denotes the sample size of the second period, and the other symbols are defined as above.
As for the previous test, under the null, the realized value Z 2 is expected to be zero, and it signals a problem when it assumes values significantly greater than zero. To perform the test and obtain the associated p-values, the following parametric bootstrap approach was used: 1. A population of one million observations was generated from each distribution under investigation (Mt, MPE, or one of the GCLEs considered in the paper), with parameters estimated by using the data of the first sample period. For the class of the GCLEs, HMC was employed to this aim; 2. B = 5000 samples, of a size equal to the dimension of the out-of-sample dataset, were drawn from each of these populations; 3. the bootstrap distributions of the statistics Z 1 and Z 2 were obtained and used to calculate the p-values, p Boot , of the statistics z 1 and z 2 computed by using the observed sample data: The results are reported in Table 7. They allow concluding that the GCLEs adequately estimate the ES, with the exception of the GCN. According to the Z 1 statistic, the Mt distribution underestimates the expected losses for α = 0.01, while the MPE underestimates the ES according to the Z 2 statistic for α = 0.1. As in the case of the in-sample assessment, distributions that are appreciably more leptokurtic than the GCN show a better tail sensitivity and are thus particularly adequate to evaluate risk measures at high levels of α.

Conclusions
In this paper, we extend the Gram-Charlier-Like Expansion (GCLE) based approach to modify the kurtosis of multivariate distributions. We obtain a class of multivariate leptokurtic distributions, called GCLEs, which prove suitable to model the excess kurtosis typically found in financial data. GCLEs are tail sensitive densities and potentially fit well the tails of the empirical distributions of multivariate financial returns. As such, they can be conveniently used to compute risks measures like the Value at Risk (VaR) and the Expected Shortfall (ES). An application of GCLEs to a portfolio of a set of financial assets provides evidence of their effectiveness in computing the mentioned risk measures both in-and out-of-sample. Furthermore, a comparison of GCLEs with other commonly employed spherical distributions such as the Mt and the MPE show that the former can yield a comparable, if not better performance with respect to the latter, both in fitting the data and in estimating risk measures. Generally, we find no systematic reason to prefer the MM estimation to ML, barring computational convenience. The paper opens up a potentially fruitful line of research providing a class of multivariate distributions that are of interest in applications involving the modeling of heavy-tailed distributions. The results thus far obtained suggest that GCLEs can be an interesting tool for risk management. This paper investigated their potential in modeling leptokurtic data that do not need any significant adjustment for skewness. In fact, the GCLEs being spherical laws, they are not suitable to account for skewness. In the context of the polynomial expansion based approach for modeling data, skewness could be dealt with by using as starting distributions other laws than the spherical ones. Thus, supplementary research could include the consideration of other moments (e.g., skewness and volatility dynamics) and the investigation of the GCLEs' performance in other financial applications such as asset pricing or risk forecasting. which, in turn, taking into account (8), implies what follows: As for (12), it must be considered that spherical variables are the standardized form of elliptical variables and that the Mardia kurtosis index, K, of elliptical variables is defined as follows (see Gómez et al. 2003): (12) in the case of spherical variables, as simple computations prove. To obtain this, let us observe that a vector y of elliptical variables has a density of the form: where µ is the mean vector and Σ is an n × n matrix related to the variance-covariance matrix of y by the the relationship: Then, according to (A6), V(x) = Σ = I, as occurs in the spherical case, if and only if n = E R 2 = m 2 . Substituting this result into (A4) yields (12).
Proof of Theorem 2. Simple computations prove that the density g n (x) can be also expressed as follows: where h(·) is the density of the generating variate R (see (11)). Now, adjusting h(·) with a binomial specified as in (1) with p 4 (.) defined as in (16) gives: h((x x) 1/2 ) = q(x, β)h((x x) 1/2 ) (A8) whereh((x x) 1/2 ) is the density of a generating variate,R, whose fourth moment,m 4 , is increased by a quantity β with respect to the same moment, m 4 , of the generating variate, R, of the parent density h((x x) 1/2 ), that is (see Nicolussi and Zoia 2020): The latter, according to (13), can be also expressed as follows: Simple computations prove that the density of the spherical variable havingh((x x) 1/2 ) as the generating variate:g n (x) = k n (x x) (1−n)/2h ((x x) 1/2 ) (A11) tallies with that of the GCLE, g n (x), given in (15). Moreover, in light of (12), Mardia's kurtosis index of g n (x) is equal to the fourth moment ofR, which proves (18).