Volatility Modeling and Dependence Structure of ESG and Conventional Investments

: The question of whether environmental, social, and governance investments outperform or underperform other conventional ﬁnancial investments has been debated in the literature. In this study, we compare the volatility of rates of return of selected ESG indices and conventional ones and investigate dependence between them. Analysis of tail dependence is important to evaluate the diversiﬁcation beneﬁts between conventional investments and ESG investments, which is necessary in constructing optimal portfolios. It allows investors to diversify the risk of the portfolio and positively impact the environment by investing in environmentally friendly companies. Examples of institutions that are paying attention to ESG issues are banks, which are increasingly including products that support sustainability goals in their offers. This analysis could be also important for policymakers. The European Banking Authority (EBA) has admitted that ESG factors can contribute to risk. Therefore, it is important to model and quantify it. The conditional volatility models from the GARCH family and tail-dependence coefﬁcients from the copula-based approach are applied. The analysis period covered 2007 until 2019. The period of the COVID-19 pandemic has not been analyzed due to the relatively short time series regarding data requirements from models’ perspective. Results of the research conﬁrm the higher dependence of extreme values in the crisis period (e.g., tail-dependence values in 2009–2014 range from 0.4820/0.4933 to 0.7039/0.6083, and from 0.5002/0.5369 to 0.7296/0.6623), and low dependence of extreme values in stabilization periods (e.g., tail-dependence values in 2017–2019 range from 0.1650 until 0.6283/0.4832, and from 0.1357 until 0.6586/0.5002). Diversiﬁcation beneﬁts vary in time, and there is a need to separately analyze crisis and stabilization periods.


Introduction
The pandemic has highlighted social and global inequality and spiked interests in environmental, social, and governance (ESG) investing. ESG assets reached $35.3 trillion in 2020 from around $30.7 trillion in 2018, reaching a third of current total global assets under management, according to the Global Sustainable Investment Association (http: //www.gsi-alliance.org/, accessed on 18 November 2021).
According to the 2020 Trends Report, investors are considering ESG factors across USD 17 trillion of professionally managed assets, a 42% increase since 2018. Such continued growth is expected over the long term, too. Since 1995, the value of US sustainable investment assets (USD 639 billion) has increased more than 25-fold (USD 16.6 trillion in 2020), at a compound annual growth rate of 14% (US SIF 2020).
A few terms are used interchangeably to describe environmental, social, and governance investments, e.g., socially responsible investing (SRI), responsible investing, sustainable investing, and impact investing. In this paper, we understand ESG investing in terms of the ESG factors, which enhance traditional financial analysis, making it more complete.
There are many ESG indices; in this paper, we select five ESG indices (see Table A1 for details). Conventional stock indices are represented by the Dow Jones Industry Average (DJI) and the S&P 500 (GSPC).
Because little is known about the dependence structure between ESG and conventional investments, the contribution of this paper is to fill this research gap. The pioneering character of this research-application of tail dependence to ESG investments, is highly contributing both to the knowledge of the field as well as to the practice of assets selection in portfolio construction. The paper also contributes to improving the understanding of volatility and the (tail) dependence structure between ESG and conventional investments. By quantifying the overall and lower tail risk between these two assets, we indicate that these dependencies exist, can be quantified, and are not negligible, especially in times of crisis.
The research results confirm the higher dependence of extreme values in the crisis period (declines of indices); however, this is also due to increased volatility during this period.
The paper is structured as follows: the introduction, literature review, research methodology, data description, empirical results, discussion, conclusions, and references. Managi et al. (2012) report no statistical difference in means and volatilities generated from the SRI indices and conventional indices in neither of the studied regions (the US, the UK, and Japan). Furthermore, they found strong comovements between the two indices in both regimes (bear and bull). In contrast, Ortas et al. (2014) found that in the period of the global financial crisis of 2008, social and responsible investment strategies were less risky in comparison to conventional investments.

Literature Review
There is no consensus in the literature as to whether ESG investments are characterized by very high returns and very low risks compared to conventional ones. Weber and Ang (2016) analyzed the performance of an emerging market SRI index concerning its financial performance compared to conventional indices. Their results indicate that the SRI index outperformed in terms of mean return the majority of the conventional emerging market portfolios. Similarly, Verheyden et al. (2016) found that both global and developed-markets portfolios (a 10% best-in-class ESG screening approach) show higher returns, lower (tail) risk, and no significant reduction of diversification potential.
On the other hand, Giese and Lee (2019) reported no clear consensus on whether ESG criteria have enhanced risk-adjusted returns. The empirical findings in the report on ESG indices performance (HKEX 2020), indicate that many ESG indices tended to have similar, if not better, risk-return performances then their parent indices. This implies that investment in ESG indices may provide equally good or even better returns while pursuing an ethics-focused investment strategy. Ouchen (2021) empirically verified whether the series of returns of an ESG index was less volatile than that of a conventional stock index. He concluded that the ESG index was relatively less turbulent than the stock index. Jain et al. (2019) reported there is no significant difference in the performance between sustainable indices and conventional indices. Plastun et al. (2022) investigated returns on ESG and conventional indices. They showed no significant differences between ESG and conventional indices. The types of price effects detected by them were the same for the cases of ESG and conventional indices (but their power was different in some cases). Charles et al. (2016) compared the risk-adjusted performance between ESG and conventional indices, as well as within the ESG indices, examining it based on standard and tail risk measures. They showed that the ESG screens for equities lead neither to a significant outperformance nor an underperformance compared to the benchmarks. They Risks 2022, 10, 20 4 of 25 also indicated that the weights used to construct these indices (sustainability-score weights vs. market cap-weights) seemed to impact their risk and performance. Apergis et al. (2015) employed a standard cointegration methodology and a novel time-varying quantile cointegration approach to investigate whether the US Dow Jones Sustainability Index and its conventional parent index are integrated. The results confirmed the presence of an asymmetric long-run relationship between these indices that is not detected by the standard methodology of cointegration.
Classical Markowitz portfolio theory does not consider the role currently played by the ESG investments on the market. It applies risk and return as single criteria, assuming that investors are rational and seek the highest return at the lowest level of risk, and their utility functions are convex (Markowitz 1952). Incorporating the ESG-based factor into the portfolio selection problem, Pedersen et al. (2021) proposed a hypothesis that explained how the increasingly widespread adoption of ESG affected portfolio choice and equilibrium asset prices.
In the case of dependency modeling for classical portfolio theory, linear correlation coefficients were used (assuming elliptic distributions of returns). The problem appears because the returns are not correlated strongly when they are around zero; however, the correlation increases in the tails. Then an appropriate tool for dependency modeling is the copula function (Sklar 1959).
Empirical research indicates that stock returns also display an asymmetric dependence in growing and declining markets, i.e., this dependence may be stronger in bearish markets than in bullish markets and tends to increase in the periods of violent fluctuations of prices (Ang and Bekaert 2002;Jondeau 2016;Longin and Solnik 2001). Because of the detected asymptotic dependence of random variables in tails (Patton 2006), the authors apply copula functions. This approach allows investigating the dependence in variance and in tails, which is not possible using standard dependence measures.
While more than 2000 empirical studies have been conducted analyzing the ESG factors and financial performance, still little is known about the dependence structure and the associated risks (Friede 2019). This is especially important as ESG scores are often related to investment risk (Bax et al. 2021).
Due to the imperfections of financial time series (they are not normally distributed) and correlation coefficients (which measure linear relationship and are constant in time), we applied GARCH family models and copula functions accompanied with some heavy-tailed marginal distributions.

Research Methodology
The ability to forecast volatility of assets is vital for portfolio selection and asset management, as well as for the pricing of primary and derivative assets (Engle and Ng 1993).
Early studies point to volatility clustering, leptokurtosis, and the leverage effect in stock-returns time series (Mandelbrot 1963) and (Fama and Fama 1965). The additional features of financial time series observed across different financial assets (stocks, stock indices, exchange rates) are as follows: stationarity, fat tails, asymmetry, aggregational Gaussianity, quasi-long-range dependence, and seasonality (e.g., Rydberg 2000;Taylor 1986). In the GARCH model, the variance is influenced by the square of the lagged innovation. However in the equity returns the leverage effect (higher impact of negative shocks on volatility) is observed, the simple GARCH model fails to describe it. GARCH (1,1) with a generalized residuals distribution can capture more volatility assessment than other models. On the other hand, the impact of asymmetry on stock market volatility and return analysis is beyond the descriptive power of the asymmetric GARCH models, which could capture more details.
There are several limitations to GARCH models. The most important one is the inability to capture the asymmetric performance. For that reason, EGARCH, GJR-GARCH, and APGARCH models were proposed. Furthermore, the asymmetric GARCH models can measure the effect of positive or negative shocks on stock market returns and volatility incompletely, and the GARCH (1,1) comparatively fails to accomplish this. The GJR-GARCH model performs better in the face of asymmetry, producing a predictable conditional variance during the period of high volatility. In addition, among the asymmetric GARCH models, the performance of the EGARCH model appeared to be superior.
Based on the properties of the studied time series: volatility clustering, leptokurtosis, asymmetry, leverage effects, mean-reversion, and stationarity-we apply the following models from the GARCH family: GARCH, EGARCH, GJR-GARCH, APGARCH, and AVGARCH in the study. We selected the GARCH model using the Akaike (AIC) and Bayesian (BIC) information criteria. Tables A2-A5 present only the results of estimation for the best-fitted models according to these criteria (with the minimum criteria values).
Base ARMA(p, q) model is as follows (Box and Jenkins 1983;Brockwell and Davis 1991): The autoregressive conditional heteroskedasticity (ARCH) models were introduced by Engle (1982) and their generalization, the GARCH models, by Bollerslev (1986). The standard GARCH(q, p) model (Bollerslev 1986) may be written as: where h t is the conditional variance, ω the intercept, and ε 2 t the residuals from the ARMA model. Some researchers pointed out limitations of the GARCH model. The most important one is that GARCH cannot capture asymmetric performance. Later, for improving this problem, EGARCH, GJR-GARCH, and APGARCH were proposed.
The exponential GARCH (EGARCH) is designed to model the logarithm of the variance rather than the level, and this model accounts for an asymmetric response to a shock. The exponential GARCH model of Nelson (1991) is defined as: where the coefficient α i captures the sign effect and γ i the size effect. E|z t−i | is the expected value of the absolute standardized innovation z t .
The Glosten-Jagannathan-Runkle GARCH (GJR-GARCH) as GARCH model captures features of financial time series like leptokurtic returns and volatility clustering.However, the GJR-GARCH model of Glosten et al. (1993) models positive and negative shocks on the conditional variance asymmetrically by the use of the indicator function I: where γ i now represents the 'leverage' term. The indicator function I takes on value of 1 for ε ≤ 0 and zero otherwise. The asymmetric power ARCH (APARCH) model of Ding et al. (1993) allows for both leverage and the Taylor effect, named after Taylor (1986) who observed that the sample autocorrelation of absolute returns was usually larger than that of squared returns: where δ ∈ R + , being a Box-Cox transformation of √ h t , and γ i the coefficient in the leverage term.
The absolute value GARCH (AVGARCH) model of Taylor (1986) and Schwert (1990): where η 1i and η 2i are rotations and shifts parameters respectively. This paper examines the structure of interdependence between ESG and conventional indices. In order to achieve this goal, we fit different theoretical distributions to the series of returns. Then, we assumed the best-fitted distribution describing the process of returns as a marginal distribution applied for our copula estimation. We used the two-stage maximum likelihood method to estimate the parameters of the considered two-dimensional copulas. In addition to testing the goodness of fit of alternative copulas, we also verified a set of hypotheses relating to the correlation matrix.
The concept of copula, was first introduced by Sklar (1959). The theoretical background for copulas is provided by Sklar's theorem. There exists a copula C such that: where F is two-dimensional joint distribution with the marginal distributions F 1 , F 2 of random variables (X 1 , X 2 ). If F 1 , F 2 are continuous, the copula C is unique: The proof is provided by, e.g., Nelsen (2006). One may use a wide range of parametric copula families to capture the different structures of dependence (e.g., Gaussian, Archimedean). In the paper, we applied the following copulas: Student's t, Joe-Clayton (a combination of the Joe copula and the Clayton copula), Gumbel-Clayton (a combination of the Clayton copula and the Gumbel copula), and survival. The Gaussian copula does not capture tail dependence, Student's t-copula has symmetric tail dependence in both lower and upper tails, and the Clayton and Gumbel copulas have only lower and upper tail dependence, respectively. Survival copulas correspond to rotation by 180 degrees.
For two dimensions following copulas are defined as (Patton 2006): (1) Gaussian/normal (N) copula where N is the normal joint distribution and Φ −1 is the quantile of the univariate normal distribution; (2) Student's t/t (t) copula where ρ ∈ [−1, 1], t ν,ρ is the joint Student's t distribution and t −1 ν is the univariate Student's t distribution with ν degrees of freedom; (3) Clayton-Gumbel (BB1) copula Risks 2022, 10, 20 7 of 25 (4) Joe-Clayton (BB7) copula Survival copula is the copula of (1 − u 1 ) and (1 − u 2 ) instead of u 1 and u 2 , respectively. Its function can measure the asymmetric dependence on the opposite side of the distribution as compared to the original function.
Our focus was on the extreme downside market risk, so we investigated the lower tail dependence in detail. The tail-dependence coefficients (Patton 2006) are: • Lower tail-dependence coefficient: The bivariate normal distribution is tail independent, the bivariate Student's t-distribution exhibits the same upper and lower tail dependence, the bivariate Joe-Clayton and Gumbel-Clayton distributions have both lower and upper tail dependence. The concept of tail dependence is embedded within the copula theory.
Instead of Pearson's correlation coefficient in the copula theory, we use the Kendall's τ coefficient (Nelsen 2006;Patton 2006). For each pair, the Kendall's τ is estimated and the p-values of the independence test based on the Kendall's τ were determined in the study.
In this paper, we used functions from the VineCopula library in R (Stoeber et al. 2018). We made the following assumptions during the estimation process: we took 39 copulas into consideration, we applied the Maximum Likelihood Estimation (MLE) method, and we used the AIC selection criterion to select the best-fitted copula. To measure the discrepancy between a hypothesized model and the empirical model, we used the goodness-of-fit (GoF) statistics based on the Kendall's process as proposed by Wang and Wells (2000). For computation of p-values, the parametric bootstrap described by Genest et al. (2006) was used.

Data Description
In the study, we used five selected ESG indices (presented in Table A1) and two stock indices-Dow Jones Industry and the S&P 500. Daily logarithmic rates of return for them were calculated, as a difference between logarithms of two consecutive closing prices multiplied by 100 (they can be interpreted as percentage changes). The data set was retrieved from the Thomson Reuters database. In Table 1, Table 4, Table 7 and Table  10 descriptive statistics for rates of return of selected indices were given. The period of analysis covered 3 July 2007 until 31 December 2019. We considered, based on important market events (e.g., global financial crisis, debt crisis in the EU, fall in oil prices) and data requirements from the models' perspectives, the following four subperiods: 1.

General Remarks
In order to identify the processes of the series of daily rates of return of the ESG indices and conventional indices (S&P 500, DJI) during the period between 3 July 2007 and 31 December 2019 (excluding the COVID-19 pandemic), it is essential to conduct graphical, statistical, and econometric examinations of these time series to check for their stationarity and the presence of the ARCH effect. We conduct such analysis for four subperiods.
The graphic examination of values of all indices shows that they are not stationary, but their daily rates of return are stationary (see Figures 1, 3, 5 and 7). We applied unit root tests, i.e., the Augmented Dickey-Fuller, Phillips-Perron, and KPSS to confirm stationarity.
At the same time, the distribution of the returns has tails, which are heavier than the tails of the normal distribution. To confirm this, we used the Jarque-Bera test and Q-Q plots.
Generally, the daily returns for both types of indices exhibit no significant autocorrelation, supporting the hypothesis that the returns are uncorrelated across time. To confirm this, we applied the Ljung-Box test. Finally, we checked whether the rates of return are characterized by the ARCH effect using the McLeoda and Li test. The null hypothesis is that the rates of return do not have the ARCH effect, while the alternative hypothesis is that they have the ARCH effect. In all subperiods, the ARCH effect was detected in the returns.
In the case of the standardized innovations from the ARMA-GARCH models, at the assumed 5% level of significance, the p-values for the Engle test are greater than 0.05 (see Tables A3-A5), which means that the null hypothesis was not rejected, i.e., the ARCH effect is not present in the innovations. Hence, the models are free from conditional heteroskedasticity in almost all cases (only in the first period for three indices the ARCH effect is present- Table A2). To verify the autocorrelation in the innovations, we used the Ljung-Box test. The null hypothesis that the innovations are independently distributed was not rejected in all the cases (see Tables A2-A5). It means that the models have been well-chosen and fitted.
Finally, the persistence parameters are close to one, which is high (see Tables A2-A5), meaning that the variance moves slowly through time.
We employ not only normally distributed innovations, but also the Student's t-distributions, the generalized error distribution (GED) and a skewed version of both. The reason for considering distributions other than normal is that a GARCH model with conditional normal errors has fatter tails than the normal distribution, and for many financial time series the standardized innovations still appear to be leptokurtic. Therefore, assuming a leptokurtic unconditional distribution for the innovations seems more appropriate.
Because the market returns are not normally distributed, the Gaussian copula would not capture tail dependence. Therefore, we fitted the Student's t-copula and the combination of Clayton and Gumbel copulas. At a 1% significance level for the chosen copulas, we cannot reject the null hypothesis (GoF test), i.e., our copulas are the true copulas. Due to the observed nonnormality in the returns distribution, we measured the dependence by using Kendall's τ coefficient. All results for the Kendall's τ coefficient are statistically significant. In the first two subperiods, dependencies were higher than in the next two subperiods (see Table 2, Table 3, Table 5, Table 6, Table 8, Table 9, Table 11 and Table 12). To quantify the degree of tail dependence in each pair, the Kendall's τ is estimated and the p-values of the independence test based on the Kendall's τ are determined (see Figure 2 and Tables 2 and 3). The results indicate the existence of positive, significant dependence.
We analyzed the dependence for 10 pairs of indices, namely S&P and ESG-5 pairs and DJI and ESG-5 pairs. Dependence between S&P and DJI exhibited lower tail dependence in three subperiods; only in the 2008 global financial crisis was the tail dependence symmetric. In the first and last subperiod, ESG and conventional investments showed lower and symmetric tail dependence-in the first subperiod for 5 pairs in the lower tail, for 5 pairs symmetric; in the third subperiod for 3 pairs in the lower tail and for 7 pairs symmetric. In the second period, ESG and conventional investments exhibited lower and upper tail dependence-for 6 pairs in the lower tail, for 3 pairs symmetric (see Tables 2, 3 , 5, 6, 8, 9, 11 and 12).

The Global Financial Crisis Period (3 July 2007-30 January 2009)
In the global financial crisis period, all indices are not normally distributed (high kurtosis, mostly positive asymmetry, as confirmed by the Jarque-Bera test) with negative means and similar standard deviations (lower only for SGESGSEP and TRESQ1- Table 1).
The year 2008 was characterized by high volatility. Volatility clustering was observed for all indices (Figure 1). There is no consensus in modeling the volatility of conventional indices. There were APGARCH (innovations with normal distribution-norm) for S&P 500 and EGARCH (innovations with skewed Student's t-distribution-sstd) for the fitted DJI. For ESG indices, GJR-GARCH (innovations with normal distribution) for A1SGI and for other AVGARCH (innovations with normal distribution) fitted best. From Table A2, we can see that all the parameters have very small p-values, which shows their significance. Persistence ranges from 0.9654 to 0.9818 for S&P 500 and DJI, and from 0.9649 to 0.9730 for ESG indices, indicating the slow movement of variance through time.  The highest values of Kendall's tau (see Figure 2) are between indices of the same type (S&P 500 and DJI, and SEESGSEP and SGESGSEP). High dependence is observed between S&P 500 and A1SGI and between DJI and A1SGI, low dependence between S&P 500 and SEESGSEP and between DJI and SEESGSEP. For example, for the GSPC-SEESGSEP relationship-the observed copula has the lowest dependency-depends by 34.24% on the upper and lower tail (Student's t-copula was employed). The GSPC-SXWESGP relationship has a dependency of 47.54% on the upper tail, and of 54.41% on the lower tail. It means the interaction has a greater effect in the lower tail (also for GSPC-SGESGSEP, GSPC-TRESGQ1, DJI-SXWESGP, DJI-SGESGSEP, DJI-TRESGQ1).  The dependence between the S&P 500 and ESG indices and between the DJI and ESG indices is modeled by Student's t and the Joe-Clayton copula (see Tables 2 and 3). There is one difference-the dependence between S&P 500 and SGESGSEP is described by the Joe-Clayton copula; between DJI and SGESGSEP by the t copula.  The dependence between the S&P 500 and ESG indices and between the DJI and ESG indices is modeled by Student's t and the Joe-Clayton copula (see Tables 2 and 3). There is one difference-the dependence between S&P 500 and SGESGSEP is described by the Joe-Clayton copula; between DJI and SGESGSEP by the t copula.

The Period of Debt Crisis in EU (2 December 2009-30 July 2014)
In the period of debt crisis in the EU, all indices are not normally distributed (mostly high kurtosis, negative asymmetry, as confirmed by the Jarque-Bera test) with positive means and similar standard deviations (lower only for DJI and A1SGI- Table 4). Two periods were characterized by high volatility (visible volatility clustering at all indices in 2009 and 2011- Figure 3). To model the volatility of conventional indices ARMA-EGARCH models (with innovations with skewed GED-sged) were fitted. In the case of ESG indices ARMA-EGARCH models also fitted best, but different distributions for innovation were applied (mostly skewed Student's t-distribution and skewed GED for SEESGSEP). From Table A3 we can see that all the parameters have very small p-values, which shows their statistical significance. Persistence for S&P 500 and DJI is similar (0.9410-0.9411) for ESG indices ranges from 0.9594 to 0.9910, indicating slow movement of variance through time.
The highest values of Kendall's τ (see Figure 4) are between indices of the same type (S&P 500 and DJI, and SEESGSEP and SGESGSEP). High dependence is observed between S&P 500 and A1SGI and between DJI and A1SGI; low dependence between S&P 500 and SEESGSEP and between DJI and SEESGSEP. For example, for the GSPC-SEESGSEP relationship-the observed copula has the lowest dependency (but higher comparing to previous period)-depends by 52.82% on the upper tail and 50.95% on the lower tail. The GSPC-A1SGI relationship has a dependency of 84.33% on the upper tail, and of 89.85% on the lower tail. It means the interaction has a greater effect in the upper tail for the first pair and for DJI-SEESGSEP, DJI-SGESGSEP, and GSPC-SGESGSEP. The greater effect in the lower tail is observed for the second pair and for the other indices.    The dependence between the S&P 500 and ESG indices and between the DJI and ESG indices is modeled by the Clayton-Gumbel and Joe-Clayton copulas (see Tables 5 and 6). The dependences between S&P 500 and SXWESGP, SGESGSEP, and TRESGQ1 comparing to DJI's dependence were modeled differently. For example, GSPC-SGESGSEP relationship is modeled by the Clayton-Gumbel copula, and DJI-SGESGSEP by the Joe-Clayton copula.  The dependence between the S&P 500 and ESG indices and between the DJI and ESG indices is modeled by the Clayton-Gumbel and Joe-Clayton copulas (see Tables 5 and 6). The dependences between S&P 500 and SXWESGP, SGESGSEP, and TRESGQ1 comparing to DJI's dependence were modeled differently. For example, GSPC-SGESGSEP relationship is modeled by the Clayton-Gumbel copula, and DJI-SGESGSEP by the Joe-Clayton copula.

The Russian Financial Crises, Fall in Oil Prices (3 June 2014-31 January 2017)
In the Russian financial crises and fall in oil prices period, all indices are not normally distributed (mostly high kurtosis, negative asymmetry, as confirmed by the Jarque-Bera test) with mostly positive means (for three ESG indices the mean was negative) and similar standard deviations for S&P 500 and DJI. In the case of standard deviation for the ESG indices, large differences were observed-low for A1SGI and high for SEESGSEP (Table 7). In this period, three subperiods characterized by high volatility (visible volatility clustering of all indices in 2015 and the beginning of 2016, and a large drop in the middle of 2016- Figure 5) were present. The volatility of three indices, namely S&P 500, DJI and A1SGI, was modeled by AVGARCH (innovations with skewed Student's t-distribution). There is no consensus in the modeling volatility of the ESG indices. GARCH, EGARCH, and AVGARCH (innovations with Student's t-distribution in all models) were applied. From Table A4, we can see that all the parameters have very small p-values, which shows their statistical significance. Persistence ranges from 0.9655 to 0.9749 for S&P 500 and DJI, and from 0.8972 to 0.9716 for ESG indices.
The highest values of Kendall's τ (see Figure 6) are between indices of the same type (S&P 500 and DJI, and SEESGSEP and SGESGSEP). High dependence is observed between S&P 500 and A1SGI and between DJI and A1SGI; low dependence between S&P 500 and SEESGSEP and between DJI and SEESGSEP. For example, for the DJI-SEESGSEP relationship-the observed copula has the lowest dependency (lower compared to previous periods too)-depends by 21.90% on the upper tail and 21.90% on the lower tail. The GSPC-A1SGI relationship has a dependency of 85.21% on the upper tail, and of 85.21% on the lower tail. It means that the interaction has the same effect in the upper tail and in the lower tail. The greater effects in the lower tail are observed for GSPC-TRESGQ1, DJI-TRESGQ1 and DJI-A1SGI.   The dependences between the S&P 500 and ESG indices and between the DJI and ESG indices were modeled by the Clayton-Gumbel and t copulas (see Tables 8 and 9). The dependence between S&P 500 and A1SGI compared to the DJI dependence was modeled differently. For example, the GSPC-A1SGI relationship is modeled by Student's t copula, and DJI-GESGSEP by Clayton-Gumbel copula.  The dependences between the S&P 500 and ESG indices and between the DJI and ESG indices were modeled by the Clayton-Gumbel and t copulas (see Tables 8 and 9). The dependence between S&P 500 and A1SGI compared to the DJI dependence was modeled differently. For example, the GSPC-A1SGI relationship is modeled by Student's t copula, and DJI-GESGSEP by Clayton-Gumbel copula.

The Stabilization Period (4 January 2017-31 December 2019)
In the stabilization period, all indices are not normally distributed (mostly low kurtosis comparing to previous periods, but higher than for normal distribution, negative asymmetry, as confirmed by the Jarque-Bera test) with positive means and similar standard deviations for S&P 500, DJI, and A1SGI. In the case of standard deviation for other ESG indices, large differences were observed-low for SGESGSEP and high for SEESGSEP (Table 10). In the whole period, two subperiods characterized by high volatility (visible volatility clustering in beginning of 2018 and the beginning of 2019- Figure 7) were observed for S&P 500, DJI, and A1SGI. For other ESG indices not only these clusters were observed (high volatility was during the whole of 2018 and the first half of 2019). Volatility of three indices, namely S&P 500, DJI, and A1SGI, was modeled by AVGARCH (innovations with skewed Student's t and skewed GED distributions). There is no consensus in the modeling volatility of the ESG indices. There were applied GJR-GARCH and EGARCH (innovations with Student's t and skewed normal distributions). From Table A5, we can see that all the parameters have very small p-values, which shows their statistical significance. Persistence ranges from 0.9581 to 0.9564 for S&P 500 and DJI, and from 0.9169 to 0.9628 for the ESG indices.   The highest values of Kendall's τ (see Figure 8) are between indices of the sam (S&P 500 and DJI, and SEESGSEP and SGESGSEP). High dependence is observed bet S&P 500 and A1SGI and between DJI and A1SGI, low dependence between S&P 50 The highest values of Kendall's τ (see Figure 8) are between indices of the same type (S&P 500 and DJI, and SEESGSEP and SGESGSEP). High dependence is observed between S&P 500 and A1SGI and between DJI and A1SGI, low dependence between S&P 500 and SEESGSEP, and between DJI and SEESGSEP. For example, for the GSPC-SEESGSEP relationship-the observed copula has the lowest dependency-depends by 16.50% on the upper tail and 16.50% on the lower tail (the same dependency). The pair GSPC-A1SGI has a dependency of 88.22% on the lower tail, and of 85.32% on the upper tail. It means the interaction has a greater effect in the lower tail. This interaction was also observed in the pairs DJI-A1SGI, DJI-TRESGQ1, DJI-SXWESGP, and GSPC-TRESGQ1.
Risks 2022, 10, x FOR PEER REVIEW 18 of 27 the interaction has a greater effect in the lower tail. This interaction was also observed in the pairs DJI-A1SGI, DJI-TRESGQ1, DJI-SXWESGP, and GSPC-TRESGQ1. The dependences between the S&P 500 and ESG indices and between the DJI and ESG indices were modeled by the Clayton-Gumbel, t, and Joe-Clayton copula (see Tables  11 and 12). Dependences between S&P 500 and SEESGSEP, and between S&P 500 and SGESGSEP compared to DJI were modeled by using the t copula. In the other cases, we used other copulas. For example, the GSPC-TRESGQ1 relationship is modeled by the Clayton-Gumbel copula, and DJI-SGESGSEP by the Survival Joe-Clayton copula.  The dependences between the S&P 500 and ESG indices and between the DJI and ESG indices were modeled by the Clayton-Gumbel, t, and Joe-Clayton copula (see Tables 11 and 12). Dependences between S&P 500 and SEESGSEP, and between S&P 500 and SGESGSEP compared to DJI were modeled by using the t copula. In the other cases, we used other copulas. For example, the GSPC-TRESGQ1 relationship is modeled by the Clayton-Gumbel copula, and DJI-SGESGSEP by the Survival Joe-Clayton copula.

Discussion
Our study does not confirm the outperformance of the ESG indices compared to conventional ones in terms of risk in the considered subperiods. However, generalization of the results is limited by the selection of a market index as a proxy of an investment portfolio. Other empirical studies demonstrate a strong correlation between the lower risk related to sustainability and better financial performance (Whelan et al. 2021). During the 2008 global financial crisis Fernández et al. (2019) found that German green mutual funds had risk-adjusted returns slightly better than their peers (in the noncrisis period, they were equal to conventional funds but better than the SRI funds). Similarly, ESG stock indices performed better and recovered faster after the 2008 global financial crisis (Wu et al. 2017). Other results confirm these findings, as in economic downturns, the high-rated mutual funds outperformed the low-rated funds, based on the Sharpe ratio (Das et al. 2018a(Das et al. , 2018bKhajenouri and Schmidt 2020). In line are research by Abate et al. (2021) that mutual funds investing in high ESG stocks perform better than investing in low ESG score stocks. Gil-Bazo et al. (2010) confirm that SRI funds perform better than their conventional counterparts.
In the study we evaluate whether conventional stock portfolios, including ESG companies, can effectively decrease portfolio risk, especially in times of financial distress on the market. Generally, we observe in all subperiods almost-weak to high lower-tail dependence between the ESG and conventional indices. Our findings indicate that there is low or symmetric tail dependence between ESG and conventional indices, meaning that if the ESG index decreases, the conventional one will also decrease accordingly. In the two first subperiods (economic downturn periods) lower tail dependence coefficients are higher comparing to the next two periods ('stabilization' periods). We conclude that dependencies exist, can be quantified, and are not negligible, especially in times of crisis.
In risk management of an asset portfolio, we are interested whether the decline of one (or more) assets may influence the behavior of the other assets in a portfolio. Especially, the occurrence of simultaneous extreme events on the market implies that risk diversification breaks down just when it is crucial. In case of extreme events, classical dependence analysis (e.g., linear correlation) fails, a copula approach is used. Some researchers argue that considering ESG practices when creating an equity portfolio (selecting companies with high ESG scores) can act as protection against left-tail risk; therefore, reducing ex-ante expectations of a left-tail event (Shafer and Szado 2020;De and Clayman 2015;Djoutsa Wamba et al. 2020).
The measurement of the tail dependence between ESG and conventional investment based on the copula approach, allows to monitor extreme risks between them. Understanding dependencies and risks is important for setting up adequate risk management as well as construction portfolios-ESG-diversified and resilient to crises. Bax et al. (2021) use the R-vine copula ESG risk model. They estimate all the conditional dependencies among assets as well as specify their interactions as modeled by different copulas families, but they also introduce three ESG risk measures that capture ESG risk, market risk conditionally on the ESG class, as well as an idiosyncratic risk component.

Conclusions
As interest in ESG investments grows, it is crucial to better understand the various risks and return tradeoffs between ESG and conventional stocks and the dependence structures between them. We used GARCH family models to estimate conditional volatilities and the copula approach, in particular the (tail) dependence structure between ESG and conventional investments to quantify the overall and lower tail risk between these two investments.
Hypothesis 1, that ESG investments outperform conventional ones in terms of risk was negatively verified-there are indices that underperform the conventional indices in selected subperiods. Hypothesis 2, that asymptotic dependence increases during the crisis on the market (declines of stock market indices), and it stabilizes during non-crisis periods was positively verified.
The findings indicate that there is no significant difference in daily returns of ESG indices and conventional ones.
The parameters are significant among all the estimations and comparisons, showing that GARCH family models may appropriately model the ESG and conventional index data. In most subperiods, the data for the ESG and conventional indices were negatively skewed. The EGARCH models this type of behavior, but in this study, also AVGARCH model was applied. Surprisingly, GJR-GARCH was not often used.
We found relations between the selected ESG indices. A1GSI is related in construction to S&P 500, while three other ESG indices-SXWESGP, SEESGSEP, and SGESGSEP are related to each other. These relations were visible in similar behavior of the rates of return. Volatility clustering observed at S&P 500, DJI, and A1SGI was different from that in other ESG indices (but the scale of these differences depends on the subperiod). In the first subperiod, one cluster was observed for all the indices (there is no difference visible). In the second and third subperiod, two clusters were found for all indices, but the difference between these two groups of indices was present. The most visible differences were observed in the last subperiod. Only in the second period, results of model estimation are consistent for all the indices (ARMA-EGARCH models with innovations sged and sstd).
We cannot confirm that the volatility of conventional indices should be modeled using different models than the volatility of the ESG-this depends on the time period and the market events. There are periods when volatility of all the indices may be modeled by using the same models, when there is a difference between the modeling of conventional and ESG indices, and when there is a difference between the modeling of ESG indices.
The characteristics of the time series indicated the need to apply to the dependence analysis copula approach. To capture tail dependence, Student's t-copula and the combination of the Clayton and Gumbel copulas were fitted best in this study. The choice of an appropriate copula function is crucial. Two features are important regarding the copula selection. The general structure of the chosen copula should coincide with the dependence structure of the real data. If the data show tail dependence than we must apply a copula which comprises tail dependence.
In the periods of economic downturn (declines of stock market indices), the dependencies measured by Kendell's tau coefficients were higher than in the less turbulent periods.
The lower tail dependence and symmetric dependence between ESG and conventional investments were detected. High values of low tail-dependence coefficients were observed in the economic downturn periods; low in stabilization periods. This signifies higher dependence of extreme values in the economic downturn periods and low dependence of extreme values in stabilization periods.
We conclude that when selecting the right model, the preliminary analysis of data is necessary, and the selection of the volatility model should be carried out for subperiods regarding different market-event characteristics. Results show, as in cases of systemic risk, that analysis of volatility and dependence structure should be carried out separately-in the periods of economic downturn and in less turbulent times.
To extend this research for the future, a more detailed analysis, including ESG and non-ESG companies from different markets (developed and developing) would be beneficial. In addition, application of a time-varying model (time-varying copula) would give insights into dependence structure.
Funding: This research was funded by Wroclaw University of Economics and Business.

Conflicts of Interest:
The authors declare no conflict of interest.