The Relationship Between ESG Scores and Value-at-Risk: A Vine Copula-GARCH Based Approach

1. Introduction

In recent years, dependence patterns among various stocks have become increasingly complex, particularly following the global financial crisis. Numerous methodologies have been developed to capture this complexity with increased flexibility. A significant body of literature has focused on multivariate copulas as models for capturing dependence due to their ability to offer flexibility in selecting marginal distributions. According to [28], any multivariate distribution can be decomposed into its marginal distributions and a copula that models the dependence structure. This foundational result allows copulas to model complex dependence structures separately from the marginal behavior of individual assets, making them powerful tools in finance and risk management. This is especially valuable in the financial context, where assets often exhibit asymmetric and tail-dependent behavior that standard correlation measures fail to capture.

Despite their advantages, traditional multivariate copulas, including both elliptical copulas (such as Gaussian and Student’s t copulas) and Archimedean copulas, have limitations. They often imply uniform dependence structures across all pairs of variables, which does not always align with the intricacies of financial markets. This uniformity can result in an inadequate representation of symmetry and tail dependence in multivariate data. Additionally, while there is a vast selection of bivariate copulas, the set of available higher-dimensional copulas is relatively limited, restricting the ability to model more complex dependencies in multivariate settings.

To address these limitations, Vine copulas have been proposed. Regular Vines, initially constructed by [16] using mixtures of conditional distribution functions, were later given a graphical structure by [2] to better understand the sequence of conditioning variables. [7] discuss how Vine copulas offer a flexible approach to modeling multivariate dependencies by constructing higher-dimensional copulas from bivariate building blocks. In fact, Vine copulas allow for the combination of copulas from different families, such as the t-copula (elliptical) and the Clayton copula (Archimedean), each of which is particularly suited for modeling specific aspects of dependency between financial variables. This allows for a more accurate representation of the dependencies between variables, particularly in cases where pairwise dependencies vary significantly.

Elliptical copulas, like the t-copula, are effective in capturing symmetric dependencies and handling heavy tails, i.e., situations with extreme events that have low probabilities but significant impacts. This is crucial in the context of Value-at-Risk (VaR), where accurately considering the possibility of extreme losses is essential. On the other hand, Archimedean copulas, such as the Clayton copula, provide excellent representation of asymmetric dependencies, especially useful in scenarios with unidirectional heavy tails. This allows for modeling situations where the dependence between variables is stronger in one tail compared to the other. The integration of these heterogeneous copulas within vine structures enables the construction of highly customized models that better reflect the reality of complex dependence structures.

There are two main types of Vine copulas: Regular Vine (R-Vine) and Canonical Vine (C-Vine). R-Vine copulas offer the most flexibility, allowing for the specification of complex dependence structures through a sequence of linked trees, as discussed by [17] and further developed by [19]. Each tree level corresponds to conditional dependencies, enabling the model to capture intricate interactions between variables. C-Vine copulas, on the other hand, are more structured, with one central variable at each level of the Vine acting as a hub that connects to other variables. This structure is useful when one variable is expected to have a dominant influence over the others, but may be less flexible in other scenarios. [6] provides a step-by-step introduction to Vines, serving as a practical guide.

Given the increasing complexity of stock return dependencies, integrating copulas with GARCH models has become a popular approach. GARCH models are widely used to model the conditional heteroskedasticity of financial time series, capturing the dynamic nature of volatility over time. By integrating GARCH models with Vine copulas, one can account for both the time-varying volatility of individual assets and the complex, potentially non-linear dependencies between them. This integration provides a more comprehensive framework for risk management, particularly for the estimation of risk measures such as VaR.

The work of [26] first introduced copula-based GARCH modeling, and [20] highlighted the benefits of using these models over dynamic conditional correlation models. The univariate GARCH model can filter the univariate serial dependence in each marginal time series, while the Vine copula model captures the cross-sectional dependence across residuals. This combined approach was used by [3] in risk management and further explored by [21] for non-parametric estimation of the marginal residual distribution. In the context of sustainable finance, [8] studied the dependence between Environmental, Social, and Governance (ESG) indices and the market, marking a significant application of Vine copulas in the ESG context.

This integration of advanced modeling techniques, such as copula-based GARCH models, with ESG considerations highlights the increasing need for holistic risk management frameworks that account not only for financial dependencies but also for sustainability factors. ESG criteria enable investors to assess both the financial and non-financial risks associated with companies by evaluating their long-term impact on ESG practices [14]. Companies with strong ESG profiles tend to demonstrate greater resilience during periods of financial stress [30]. Such developments gained significant momentum as regulatory bodies, particularly within the European Union, have introduced frameworks like the Sustainable Finance Disclosure Regulation (SFDR) [12] to enhance transparency and consistency in ESG disclosures, promoting a more comprehensive risk management approach.

Building on this literature, this paper aims to investigate the use of the Vine copula approach for estimating the one-year VaR of a stock portfolio, focusing on the relationship between VaR and ESG scores. By applying the model to the daily log-returns of multiple stocks in the London financial market, we seek to detect interdependencies among them. Specifically, the VaR is computed for an equal-weighted portfolio and compared with results obtained by increasing the proportion of the best ESG stocks while reducing the proportion of the worst ESG stocks.

Preliminary results indicate that increasing the proportion of top-performing ESG stocks in a portfolio can lead to a reduction in overall market risk, particularly during periods of market stress. This suggests that sustainable investments not only align with ethical considerations but also contribute to enhanced portfolio stability and risk mitigation. In particular, the VaR is calibrated over a one-year period with a confidence level of 99.5%, in accordance with the Solvency II Directive. This research can serve as a basis for future works, considering the GARCH-Vine copula as a possible internal model to calculate the Solvency Capital Requirement (SCR) for equity risk.

The structure of this article is organized as follows: Section 2 explains the methodology, Section 3 presents the numerical application and finally Section 4 provides a brief conclusion.

2. Research Methods

This study focuses on VaR as a measure of the total risk to which a stock portfolio is exposed. For a given portfolio, VaR is defined as the maximum loss that is likely to occur over a given period with a specified confidence level of $(1-\alpha )$. In other words, there is a probability equal to $\alpha$ of experiencing losses greater than the VaR prediction. Formally, the one-step-ahead VaR at the $\alpha$ level is expressed as:

$$\begin{array}{c}\hfill Pr(r_t\le Va{R}_{\alpha ,t}|{\mathit{F}}_{t-1})=\alpha ,\end{array}$$

(1)

where $r_t$ represents the one-period return from time $t-1$ to time t, $\alpha$ denotes the quantile level, and $F_t$ is the information set at time t. Thus, VaR can be interpreted as the $\alpha -th$ quantile of the conditional distribution of returns.

The VaR forecasting can be enhanced using a copula approach, which captures the complex dependencies among assets and provides more accurate VaR estimates than traditional methods. While it is established that copulas generally outperform conventional methods, our approach employs Vine copulas to further refine this accuracy. Vine copulas offer superior flexibility and precision by enabling the combination of copulas from different families, which improves both the prediction and calculation of VaR beyond the capabilities of standard multivariate copulas.

This approach not only improves the accuracy of VaR estimates but also enhances the overall robustness of risk predictions, overcoming the limitations of multivariate copulas and providing more reliable financial risk management [22]. Specifically, Vine copula models are widely employed in finance to capture cross-sectional dependence among residuals obtained from applying univariate time series models to financial returns, as will be detailed in the following sections.

2.1. GARCH Model Estimation and Standardization of Residuals

A significant body of literature has established the effectiveness of GARCH models in capturing volatility dynamics in financial markets. [11] introduced the ARCH model, which [2] extended to the GARCH framework, allowing for more flexible modeling of time-varying volatility. These models have become fundamental in the analysis of return series and the forecasting of conditional volatility, playing a crucial role in risk management and the estimation of VaR [4,18].

Further applications of GARCH models have been explored in various contexts, including portfolio optimization [29], asset pricing [5], and the assessment of market risk [10]. These models have proven particularly effective in capturing the volatility clustering and heavy-tailed behavior often observed in financial time series.

To implement the GARCH-Vine approach, the first step is to fit a GARCH$(p,q)$ model to the return series of each asset, estimating the parameters using Maximum Likelihood Estimation. The GARCH$(p,q)$ model for asset i is defined as:

$$\begin{array}{ccc}\hfill r_{i,t}& =& {\mu }_i+{ϵ}_{i,t}\hfill \\ \hfill {ϵ}_{i,t}& =& {\sigma }_{i,t}·z_{i,t}\hfill \\ \hfill {\sigma }_{i,t}^2& =& {\alpha }_0+\sum _{j=1}^q{\alpha }_j{ϵ}_{i,t-j}^2+\sum _{k=1}^p{\beta }_k{\sigma }_{i,t-k}^2,\hfill \end{array}$$

(2)

where ${\sigma }_{i,t}^2$ denotes the conditional variance of asset i’s return at time t, with $r_{i,t}$ representing the return , ${ϵ}_{i,t}$ denoting the residual and $z_{i,t}$ is a realization of an independent and identically distributed (i.i.d.) random variable with mean zero and unit variance, often assumed to follow a standard normal or t-Student distribution. The process is weakly stationary if and only if ${\sum }_{j=1}^q{\alpha }_j+{\sum }_{k=1}^p{\beta }_k<1$. Specifically, in this study, we will utilize the GARCH(1,1) model in the numerical application, as it is particularly effective for capturing the volatility clustering commonly observed in financial time series data, while also ensuring that we do not complicate the model beyond what is necessary for our analysis.

Once the GARCH model is estimated, the next step involves standardizing the residuals of the univariate return series by dividing the residuals by the estimated conditional volatilities:

$$\begin{array}{cc}\hfill z_{i,t}& ={\displaystyle \frac{{ϵ}_{i,t}}{{\sigma }_{i,t}}},\hfill \end{array}$$

(3)

where ${ϵ}_{i,t}$ represents the residual of asset i at time t, and ${\sigma }_{i,t}$ is the estimated conditional volatility at time t.

Subsequently, these standardized residuals are transformed into a series of univariate pseudo copula values using the Probability Integral Transform (PIT) [13]. This step is crucial, as (Vine-)copula models require variables to be uniformly distributed to effectively model and analyze dependencies among multiple assets, in accordance with [28], which will be explained further below. The PIT is a fundamental concept in probability theory and statistics that ensures the uniformity of transformed variables. Specifically, if X is a continuous random variable with cumulative distribution function (CDF) $F_X$, the transformation $U=F_X\left(X\right)$ results in U being uniformly distributed on the interval [0,1]. This result is due to the fact that the CDF maps values of X to probabilities, which are uniformly distributed over the interval [0,1].

In our context, the standardized residuals $z_{i,t}$ (3) are transformed into uniform variables using their empirical CDF $F_i$, such that:

$$\begin{array}{c}\hfill u_{i,t}=F_i\left(z_{i,t}\right).\end{array}$$

(4)

Here, $u_{i,t}$ represents the pseudo copula values. However, if there is a justified reason to believe that the residuals follow a specific theoretical distribution, it is possible to use a non-empirical CDF, such as that of the standard normal distribution or another distribution chosen based on the characteristics of the data. As we will discuss further below, by converting the standardized residuals into uniformly distributed variables, we ensure that the dependency structure between the assets can be captured and analyzed effectively using copula models.

2.2. Vine Copula and Dependency Modeling

The strength of copulas comes from Sklar’s theorem, which demonstrates that the structure of a multivariate distribution is independent of the univariate marginals, with this structure being defined by the copula function. Specifically, the theorem states that:

Theorem 1 

(Sklar’s Theorem). Let F be the CDF of the vector $\mathbf{X}=(X_1,\dots ,X_d)$, and let $F_i$ be the marginal CDFs for $i\in \{1,\dots ,d\}$. Then, there exists a copula C such that for all $x_i\in \mathbb{R}$:

$$F(x_1,\dots ,x_d)=C\left(F_1\left(x_1\right),\dots ,F_d\left(x_d\right)\right),$$

(5)

with the associated density:

$$f(x_1,\dots ,x_d)=c(F_1\left(x_1\right),\dots ,F_d\left(x_d\right))f_1\left(x_1\right)\cdots f_d\left(x_d\right),$$

(6)

for some d-dimensional copula C with copula density c. If the marginals $F_i$ are continuous, then C is the unique copula associated with F. Conversely, if C is a copula and $F_i$ are univariate distribution functions, then F is the joint distribution function for $\mathbf{X}$.

Sklar’s theorem states that any multivariate distribution can be expressed in terms of its marginals and a copula that captures the dependency structure between the variables. The proof of Sklar’s theorem can be found in [24] (pages 45-60).

To understand how copulas capture these dependencies, it is essential to understand Kendall’s $\tau$ and its interpretation. Given two pairs $(X_{11},X_{21})$ and $(X_{12},X_{22})$ drawn from a random vector $(X_1,X_2)$, Kendall’s $\tau$ is defined as:

$$\tau (X_1,X_2)=Pr[(X_{11}-X_{12})(X_{21}-X_{22})>0]-Pr[(X_{11}-X_{12})(X_{21}-X_{22})<0].$$

(7)

Here, $Pr[(X_{11}-X_{12})(X_{21}-X_{22})>0]$ denotes the probability of concordance, while $Pr[(X_{11}-X_{12})(X_{21}-X_{22})<0]$ denotes the probability of discordance. Therefore, Kendall’s $\tau$ measures the relative difference between the probability of concordance and discordance, providing insights into the strength and direction of the association between the variables.

However, while there is a wide range of parametric copula families, they tend to lack flexibility as the dimensionality of the vector $\mathbf{X}$ increases. Moreover, they do not allow for different dependence structures between subgroups of random variables. In order to overcome the limitations of traditional multivariate copulas, Vine Copulas extend the concept of capturing dependencies in multivariate distributions by offering flexibility through a sequence of linked bivariate copulas arranged in a tree structure.

Fundamentally, the concept underlying Vine copulas is the Pair Copula Construction introduced by [16] which allows for the factorization of the multivariate density using bivariate copulas. Specifically, this approach begins with the decomposition of the multivariate density into conditional densities:

$$\begin{array}{ccc}\hfill f(x_1,\dots ,x_d)& =& f_d\left(x_d\right)·f_{d-1\mid d}(x_{d-1}\mid x_d)·\hfill \\ & & f_{d-2\mid d-1,d}(x_{d-2}\mid x_{d-1},x_d)\cdots f_{1\mid 2,\dots ,d}(x_1\mid x_2,\dots ,x_d).\hfill \end{array}$$

(8)

For the bivariate case, the previous equation simplifies to:

$$f(x_1,x_2)=f_2\left(x_2\right)·f_{1\mid 2}(x_1\mid x_2).$$

(9)

The second factor of the equation above can be further expressed by applying the definition of conditional density in conjunction with equation (6):

$$\begin{array}{cc}{f}_{1\mid 2}(x_1\mid x_2)\hfill & ={\displaystyle \frac{f(x_1,x_2)}{f_2\left(x_2\right)}}\hfill \\ \hfill & ={\displaystyle \frac{c_{12}(F_1\left(x_1\right),F_2\left(x_2\right))f_1\left(x_1\right)f_2\left(x_2\right)}{f_2\left(x_2\right)}}\hfill \\ \hfill & =c_{12}(F_1\left(x_1\right),F_2\left(x_2\right))·f_1\left(x_1\right),\hfill \end{array}$$

(10)

where $c_{12}$ is the pair-copula density. This leads to the formulation of the joint density in terms of pair-copulas and marginal densities:

$$f(x_1,x_2)=c_{12}(F_1\left(x_1\right),F_2\left(x_2\right))·f_1\left(x_1\right)·f_2\left(x_2\right).$$

(11)

Next, we consider the case where $d=3$, to better understand how to decompose more complex joint density within the multivariate framework. The trivariate joint density is given by:

$$f(x_1,x_2,x_3)=f_{3\mid 1,2}(x_3\mid x_1,x_2)·f_{2\mid 1}(x_2\mid x_1)·f_1\left(x_1\right).$$

(12)

We factorize this formula by considering each part separately. In particular, the conditional density $f_{2\mid 1}(x_2\mid x_1)$ follows the structure of equation (10) presented earlier:

$$f_{2\mid 1}(x_2\mid x_1)=c_{12}(F_1\left(x_1\right),F_2\left(x_2\right))·f_2\left(x_2\right).$$

(13)

Additionally, the factor $f_{3\mid 1,2}(x_3\mid x_1,x_2)$ is the conditional density of $X_3$ given that $X_1=x_1$ and $X_2=x_2$ and it can be expressed as:

$$\begin{array}{cc}{f}_{3\mid 1,2}(x_3\mid x_1,x_2)\hfill & =c_{13\mid 2}(F_{1\mid 2}(x_1\mid x_2),F_{3\mid 2}(x_3\mid x_2))·f_{3\mid 2}(x_3\mid x_2)\hfill \\ \hfill & =c_{13\mid 2}(F_{1\mid 2}(x_1\mid x_2),F_{3\mid 2}(x_3\mid x_2))·c_{23}(F_2\left(x_2\right),F_3\left(x_3\right))·f_3\left(x_3\right),\hfill \end{array}$$

(14)

which results from a recursive process (see [6]).

Thus, the joint density function (12) can be expressed in terms of bivariate copulas as follows:

$$\begin{array}{cc}\hfill f(x_1,x_2,x_3)& =f_3\left(x_3\right)·f_2\left(x_2\right)·f_1\left(x_1\right)·\hfill \\ \hfill & c_{12}(F_1\left(x_1\right),F_2\left(x_2\right))·c_{23}(F_2\left(x_2\right),F_3\left(x_3\right))·\hfill \\ \hfill & c_{13\mid 2}(F_{1\mid 2}(x_1\mid x_2),F_{3\mid 2}(x_3\mid x_2)).\hfill \end{array}$$

(15)

However, this decomposition is not unique. There are two other possible copula constructions in three dimensions by changing the order of the variables. For a large number of variables, the decomposition can become intricate and challenging to interpret. For this reason, [2] extended Joe’s work [16] by representing these decompositions as interconnected tree structures within the context of graph theory, providing a more comprehensible understanding.

Following [2], we define a parametric R-Vine distribution as follows:

Definition 1.(Parametric R-Vine distribution). A parametric R-Vine distribution $R(\mathcal{V},\mathcal{C},\theta )$ for a random vector ( $\mathbf{X}$ ) consists of three main components:

• Tree Structure: A set of linked trees $V_j$, where j ranges from 1 to $d-1$, that is defined as $\mathcal{V}=\{V_1,V_2,\dots ,V_{d-1}\}$.
• Parametric Bivariate Copulas: A set of bivariate copulas $\mathcal{C}=\mathcal{C}\left(\mathcal{V}\right)$ assigned to each edge in the tree structure. These copulas, known as pair copulas, model the pairwise dependencies between the variables connected by the respective edges.
• Corresponding Parameters: The parameters $\theta =\theta \left(\mathcal{C}\right(\mathcal{V}\left)\right)$ associated with the corresponding copulas in $\mathcal{C}$. These parameters define the specific dependency structure captured by each copula.

The density of an R-Vine for a random vector ( $\mathbf{X}$ ) is given by:

$$f\left(\mathbf{x}\right)=\left(\prod _{i=1}^d{f}_i\left(x_i\right)\right)·\left(\prod _{k=1}^{d-1}\prod _{(i,j)\in E_k}{c}_{ij;D}\left(F_{i|D}(x_i\mid x_D),F_{j|D}(x_j\mid x_D)\right)\right),$$

(16)

where $f_i\left(x_i\right)$ is the marginal density; $c_{ij;D}$ is the bivariate copulas modeling the dependence between the variables $x_i$ and $x_j$ conditionally on a set D; $F_{i|D}$ and $F_{j|D}$ represent the conditional distribution functions of the variables $x_i$ and $x_j$, respectively, conditioned on the set D; $E_k$ is the set of pairs $(i,j)$ corresponding to the edges in the k-th tree of the R-Vine structure.

Building on this foundation, the selection of the R-Vine model involves three key choices: determining the structure of all its trees, selecting the copula family for each of the ${\displaystyle \frac{d(d-1)}{2}}$ pair copulas, and specifying the parameters for each pair copula. The selection strategy employed here is heuristic and bottom-up. For further details, please refer to the paper by [9].

However, the generality of R-Vines can make them complex to specify and interpret, especially as the number of variables increases. To address this complexity, more structured Vine copulas, such as the C-Vine, have been introduced to simplify the specification of the dependency structure. The C-Vine imposes a specific structure where, at each level, there is a central variable. This structure reduces complexity and facilitates the interpretation of dependencies, particularly when some variables are thought to play a more significant role in determining the overall dependence. Graphically, as illustrated in the Figure 1 for the case $d=4$, each tree in the C-Vine has a unique node that is connected to $d-j$ edges, where j denotes the level of the tree.

2.3. Forecasting and VaR Calculation

With the dependency structure in place, the next stage involves forecasting returns and computing the VaR for the portfolio over the chosen time horizon. This process can be broken down into key steps:

1) The first step involves fitting the R-Vine copula to the standardized residuals obtained from the GARCH model. Before fitting the copula, the residuals are standardized and then transformed into uniform pseudo copula values using the chosen CDF, which can be either empirical or theoretical, as discussed above. This fitting process is crucial for accurately capturing the dependency structure among the assets.

2) Using the fitted R-Vine model, we generate a random sample of pseudo copula values $u_{t+1},u_{t+2},\dots$, simulating joint realizations over the chosen time horizon to estimate the cumulative returns.

3) The simulated pseudo copula values are then transformed into standardized residuals by applying the inverse CDFs of the residual series. This step ensures that the simulated values match the distributional properties (e.g., t-distribution) of the original residuals. Mathematically, for each asset, this transformation is expressed as:

$$\begin{array}{cc}\hfill {\widehat{z}}_{i,t+1}& =F_i^{-1}\left(u_{i,t+1}\right),\hfill \end{array}$$

(17)

where $F_i^{-1}$ is the inverse CDF corresponding to the distribution of the standardized residuals in the context of one-step-ahead forecasting.

4) After transforming the pseudo copula values into standardized residuals, the estimated GARCH model is employed to forecast the returns for each asset. The return forecast for asset i is given by:

$$\begin{array}{cc}\hfill r_{i,t+1}& ={\mu }_i+{\sigma }_{i,t+1}·{\widehat{z}}_{i,t+1},\hfill \end{array}$$

(18)

where ${\sigma }_{i,t+1}$ is the forecasted conditional volatility from the GARCH model. The model simulates a number of time steps based on the chosen time horizon.

5) The cumulative return of the portfolio is then derived from the simulated returns for each component. This simulation process is repeated multiple times, producing a distribution of potential portfolio returns over the specified time horizon. From this distribution, the VaR is calculated, representing the maximum expected loss over the chosen period at a given confidence level.

3. Numerical application

In this section, the previously described methodology is applied to an empirical dataset. The aim is to compute the one-year VaR for different compositions of the porfolio by using the GARCH-Vine copula model. The R-Vine model is implemented using the VineCopula package developed by [27]. The forecasted conditional volatility is derived from the GARCH(1,1) model, and the simulation process is executed over 252 steps to estimate cumulative returns for the one-year time horizon. By generating 10,000 simulations, we obtain a distribution of simulated one-year cumulative returns. From this distribution, the VaR is calculated with a confidence level of $(1-\alpha )=99.5\%$, in linee with the requirements of the Solvency II framework for determining the SCR for equity risk.

3.1. Data and Descriptive Statistics

The analysis is conducted on close-to-close daily log-returns of sixteen companies listed on the London Stock Exchange. Raw stock data was sourced from finance.yahoo.com website for the period spanning from January 1, 2018, to December 31, 2022. Relative ESG scores were obtained from Bloomberg to select the top eight companies and the bottom eight companies based on their average ESG scores over the five-year period considered.

Table 1 presents the ESG scores and sector classifications of the portfolio’s components. The data reveals a diverse range of industries, from Beverages and Real Estate Investment Trusts to Oil & Gas Producers and Mining. Notably, companies like Diageo (7.20) and British Land Company (7.10) exhibit high ESG scores, indicating strong sustainability practices. In contrast, companies such as Frasers Group (1.96) and Howden Joinery Group (2.45) have lower scores, reflecting more significant ESG challenges. In addition, it is worth noting that the Financial Services sector exhibits a wide range of ESG scores, reflecting diverse approaches to sustainability within the industry. For example, the London Stock Exchange Group, with an ESG score of 6.38, demonstrates a strong commitment to ESG practices. In contrast, Beazley Group, with a significantly lower ESG score of 2.58, highlights the sector’s variability and suggests that some companies may have more significant ESG-related challenges. Overall, Table 1 underscores the significant variability in ESG scores across different industries, offering valuable insights for investors and stakeholders interested in sector-specific ESG practices.

To complement this analysis, Table 2 reports the descriptive statistics for all return series associated with the portfolio.

As depicted in Figure 2, it is evident that all close-to-close daily log-returns exhibit sufficient stationarity, fluctuating around a common level.

Notably, significant return changes are followed by larger changes, while smaller changes are followed by smaller ones, indicating volatility clustering, a characteristic captured by GARCH models. In essence, the series display conditional heteroskedasticity.

Moreover, Table 2 highlights heavy tails, high kurtosis, and serial correlation across the series. The skewness values differ from zero, implying asymmetric distribution. Specifically, many series exhibit negative skewness, suggesting longer left tails than right tails. The exceptionally high kurtosis values indicate leptokurtosis, signifying heavier tails compared to a normal distribution. Consistent with these statistics, the high Jarque-Bera test statistics and low p-values reject the null hypothesis of normality for all analyzed series.

Furthermore, considering the validity of the GARCH process with correlated squared residuals, the Box-Pierce test on squared returns at lag 5 rejects the null hypothesis of no serial autocorrelation due to very low p-values, indicating non-zero autocorrelations up to the tested lags. Lastly, the Augmented Dickey-Fuller test confirms stationarity, rejecting the unit root hypothesis and confirming no need for differentiation. In the next section, the serial correlation identified in the previous test is modeled through the GARCH model, which assumes autoregressive behavior in conditional variance.

3.2. Model Estimation and Results

Firstly, the optimal time series model needs to be identified and selected. A GARCH(1,1) model is applied across all series due to the statistical significance of all coefficients, aiming to maintain model simplicity. Consistent with earlier findings, the errors are assumed to follow a t-distribution.

Following the calculation process, as described in Section 2.3, the GARCH model was calibrated in order to obtain standardized residuals, subsequently transformed into pseudo copula values. Figure 3 shows the corresponding pairwise normalized contour plots for the sixteen assets.

The bivariate correlations are measured in terms of Kendall’s tau, which is a non-linear measure of dependence that ranges between [-1,1]. The correlation among series is not very strong, except for the correlation between Shell (SHEL) and Bp plc (BP), which is 0.71. This is not surprising since they belong to the same sector of oil and gas. Most of the series do not hold independence but only an imperfect positive dependence. Additionally, bivariate distributions tend to exhibit elliptical shapes with varying tail dependence structures and central symmetry, underscoring the suitability of the t-copula in pair copula constructions.

The R-Vine model is used to describe the joint distribution among univariate residual time series. Although our fitted model includes 15 estimated trees, we limit the presentation to the first two because of space limitations. The initial trees are more representative as they describe the strongest direct dependencies between the variables, capturing the primary relationships within the system. In fact, the R-Vine model builds a hierarchy of dependencies, where the first tree captures the fundamental relationships and the subsequent trees refine these connections by adding conditional levels.

The estimation results of pairwise copula families for Tree 1 and Tree 2 are provided in Table 3 and Table 4, respectively, while Figure 4 and Figure 5 show the corresponding graphical representations, along with the fitted Kendall’s tau.

In Table 3, the "conditioned" column presents the variables directly involved in each copula pair. For our portfolio stocks, the t-student copula consistently emerges as the optimal choice, indicating symmetric dependency patterns within the British financial market. The "parameters" column provides the estimated coefficients, offering insights into the strength and nature of these dependencies. The first coefficient represents the correlation parameter, while the second corresponds to the degrees of freedom, which vary considerably, ranging from 4.06 to 14.29. This wide range suggests diverse tail behaviors among the pairs, a variability that can not be fully captured by a single t-student copula, which only allows for one parameter of degrees of freedom.

In Figure 4, the St James’s Place (STJ) stock emerges as a significant node within the dependence network. This observation is also supported by the C-Vine structure analysis, which identifies STJ as the root node. It is important to note that the selection of the root node is optimized to maximize the sum of pairwise dependencies. Therefore, STJ plays a crucial role in driving the performance of our stock portfolio, which aligns with the prominent position of the banking sector in the British economy. A comparison of Figure 1 with Figure 4 and Figure 5 reveals the differences in the tree structures arising from the application of C-Vine and R-Vine copulas.

Furthermore, Table 4 includes a "conditioning" column, representing the variables acting as conditioning factors within the Vine structure and appearing only from the second tree onward. These variables have already established a relationship at the previous level of the tree structure and are now conditioning the dependence between the "conditioned" variables. As shown in Table 4, Tree 2 includes non-elliptical pair copula families, yet the dependency strength, as indicated by Kendall’s tau values, is relatively weak and diminishes further in subsequent trees. Specifically, we observe Gaussian, Clayton, and Gumbel copulas, each of which presents a single parameter value. For the Gaussian copula, this value represents the correlation coefficient, which varies between $[-1,1]$, allowing it to capture a wide range of dependence, from perfect negative correlation to perfect positive correlation. Similarly, for the Clayton and Gumbel copulas, the single parameter reflects the dependence structure’s specific characteristics: in the case of the Clayton copula, which has a parameter $\theta >0$, the value indicates the degree of lower tail dependence. For the Gumbel copula, with a parameter constrained to $\theta \ge 1$, it characterizes the upper tail dependence. Despite the presence of these alternative copula families, even in the second tree, the t-student copula remains the most frequently selected option across the conditional pairs, with corresponding degrees of freedom reaching up to 24.83. In subsequent trees ($j>2$), even larger estimated degrees of freedom are observed, indicating a reduction in tail dependence for these conditional pairs.

The fitted R-vine model is then utilized to simulate cumulative returns in order to forecast the one-year VaR for four different portfolio compositions, following the procedure outlined in Section 2.3. The initial total value of the portfolio at the start of forecasting is assumed to be £100 million. Portfolio 1 consists of an equal-weighted selection of 16 stocks: the top 8 and bottom 8 companies ranked by ESG scores within the FTSE 100 (London Stock Exchange). Portfolio 2 increases investments in stocks with the highest ESG scores by 10%, while decreasing investments in stocks with the lowest ESG scores by 10%. Portfolio 3 follows the same strategy but with a 20% adjustment, and Portfolio 4 with a 30% adjustment. In addition, the same procedure is computed for the sub-period that goes from January 31, 2020 to the end of 2022. It allows to isolate the effects of the COVID crisis.

Table 5 shows the VaR values for the different portfolios, along with their respective percentage changes. The latter are computed using the VaR of Portfolio 1 as the baseline. As per convention, VaR is reported as a positive value. It can be observed that an increase in green assets, at the expense of those with a low ESG level, results in a reduction of the one-year VaR. It is easy to see from the last column how, during periods of crisis, investors should direct their choices towards sustainable investments, which are more resilient to market variations. In fact, when the period is characterized by stressed market conditions, an increment in the portion of the highest ESG assets implies a stronger reduction in the VaR compared to the period that also includes normal market conditions (2018-2022). For example, by looking at Portfolio 2, an increase of 10 % of the best ESG stocks (and a corresponding decrease in the share invested in the worst ESG stocks) allows a decrease in the VaR of 3.16 % for the period from 2018 to 2022. The same change in the portfolio composition results in a reduction of 5.89% in the VaR for the period from 2020 to 2022.

4. Conclusions

This study is intended to serve as a foundational exploration of the relationship between ESG scores and VaR within a Vine copula framework. While numerous articles have demonstrated the efficacy of the Vine model in predicting VaR, few have explored its potential for reducing the latter through sustainable investments.

Specifically, this paper emphasizes that increasing investments in top-performing ESG stocks can reduce overall portfolio market risk. This effect is particularly pronounced during periods of market stress, when investors should consider sustainable options that enhance portfolio resilience.

Moreover, the portfolio’s VaR is calibrated as required by the Solvency II Directive for the calculation of the SCR. It paves the way for using the GARCH-Vine as an internal model and anticipating the possible relationship between ESG scores and solvency.

In light of these results, future research could build on this study by incorporating a time-varying covariate, enabling the use of a "conditional copula" approach. This would introduce an innovative dimension to the analysis, offering deeper insights into how dynamic factors influence the relationship between ESG scores and VaR. Furthermore, expanding the analysis to include a larger dataset of stocks or financial instruments would bolster the evidence supporting the proposed thesis and provide a more comprehensive understanding of how ESG investments impact portfolio risk management.