Penalized Averaging of Quantile Forecasts from GARCH Models with Many Exogenous Predictors

This study explores the multi-step ahead forecasting performance of a so-called hybrid conditional quantile method, which combines relevant conditional quantile forecasts from parametric and semiparametric methods. The focus is on lower (left) and upper (right) tail quantiles of the conditional distribution of the response variable. First, we evaluate and compare out-of-sample conditional quantile forecasts obtained from a hybrid method and from five non-hybrid methods, employing a large data set of exogenous predictors generated by various GARCH model specifications. Second, we compare the accuracy of these methods by calculating conditional quantile forecasts for the risk premium of the monthly S&P 500 index, using a data set of macroeconomic predictors. Monte Carlo and empirical application results indicate that the hybrid forecasting provides more accurate quantile forecasts than non-hybrid methods. The success of the hybrid method is most prominent when compared with results obtained by a simple equal-weighted combination of quantile forecasts.


Introduction
Combining forecasts obtained by a hybrid approach has long been known to improve forecast accuracy. Because each model whose forecasts are to be combined may consider different predictors and make different assumptions about the underlying data generating process (DGP) and distributions, averaging the individual forecasts broadens the information embedded and may offset individual model biases as well; see, e.g., Zhang (2003) for additional reasons. In a one-dimensional (univariate) setting the hybridization of forecasts obtained from time series models/methods has been well researched; see, e.g., Aijaz and Agarwal (2020), Aydin and Isci Güneri (2015), Pai and Lin (2005), and Valenzuela et al. (2008) among others.
Recently, however, the enhanced availability of large databases with many time series variables as potential predictors has stimulated interest in high-dimensional forecasting. Excellent reviews of the state-of-the-art in economics and finance are given by Fan et al. (2011) and Lee (2011). Furthermore, the introductory section of the paper by Uematsu and Tanaka (2019) provides an extensive survey of the current literature on high-dimensional forecasting and variable selection. Along this direction, conditional quantile averaging procedures in conjunction with dimension reduction methods have been considered. Examples include quantile forecasting the S&P 500 equity risk premium (see, e.g., Konzen and Ziegelmann, 2016;Lima and Meng,2017;Meligkotsidou et al.,2014;De Gooijer and Zerom, 2020) inflation forecasting (Garcia et al., 2017;De Gooijer and Zerom, 2019), quantile forecasting of macroeconomic time series (Manzan, 2015;Jiang et al., 2018), and realized volatility forecasting (Meligkotsidou et al., 2019).
Most of these studies are limited by the assumption that the high-dimensional data set comes from a linear DGP, and forecasts are obtained from finite-dimensional parametric models. One exception is the study by De Gooijer and Zerom (2020). Using a large data set of predictors, involving both macroeconomic predictors and technical indicators, these authors showed that combining quantile forecasts from parametric and semiparametric methods (called hybrid quantile averaging) can be useful in practice. Semiparametric models are infinite-dimensional. As a result, the quantile forecasts are less prone to model misspecification as may happen with forecasts obtained from parametric models. In terms of quantile forecast performance the hybrid method works well in identifying relevant predictors, and more importantly results in improved combined one-step ahead forecasts over alternative (no hybridization) conditional quantile methods.
While the empirical findings by De Gooijer and Zerom (2020) are interesting, they are sample-specific which makes it difficult to generalize the results to novel situations. In addition, the focus on one-step ahead out-of-sample prediction is somewhat restrictive. Multi-step ahead out-of-sample quantile forecasting results can provide more insight in the relative performance of the hybrid quantile averaging method over a longer time period. Unfortunately, a general theoretical comparison of quantile forecasts obtained from hybrid and non-hybrid methods is not feasible due to complicated interactions of nonlinear parameter estimation methods, sparse modelling, and correlated forecasts. Indeed, further insights regarding the performance of the hybrid quantile averaging method can only be obtained via a Monte Carlo simulation study. In the first half of this paper, we provide such a study.
In the second half of the paper, we evaluate the out-of-sample multi-step ahead forecasting performance of the hybrid conditional quantile method and five alternative, non-hybrid, forecasting methods via an empirical application. More specifically, we report out-of-sample conditional quantile forecasts for the risk premium of the monthly S&P 500 index using a large data set of macroeconomic predictors. A simple equal-weighted combination of parametric and semiparametric conditional quantile forecasts is adopted as a benchmark. Our main finding is that that hybridization can be an effective way to improve quantile forecasts as compared to non-hybrid methods.
The rest of this paper unfolds as follows. First, for ease of reference, Sect. 2 summarizes the main features of the semiparametric and hybrid conditional quantile averaging methods. Section 3 describes six quantile forecasting methods used in the Monte Carlo experiment. Section 4 introduces the large data set of exogenous predictors via semiparametric and parametric GARCH model specifications. This section also discusses the evaluation of quantile forecasts. Section 5 contains the simulation results, while Sect. 6 presents the empirical results. Finally, Sect. 7 contains some concluding remarks.

Some Notation
To simplify presentation, we only discuss the case of making h ¼ 1 step ahead quantile forecasts, where h denotes the forecast horizon. The extension to multihorizon ðh [ 1Þ conditional forecasting is straightforward. Let fY t g n t¼1 be a set of observations obtained from a strictly stationary time series process fY t ; t 2 Zg that depends on q y ! 1 past values of Y t , and on a q z -dimensional vector Z t ¼ ðZ 1;t ; . . .; Z q z ;t Þ T that consists of exogenous, possibly lagged, stationary time series.
. . .; Y tÀq y Þ T , and q ¼ q y þ q z . Given the observed data set fðX t ; Y t Þg n t¼1 , the one-step ahead out-of-sample sth ð0\s\1Þ conditional quantile of the unobserved random variable Y nþ1 given X t ¼ X n will be denoted by Q Y nþ1 ðsjX n Þ.
For convenience of later analysis, we define the associated process fðX t ; Y Ã t Þg 2 R q Â R where the components of the predictor vector X t ¼ ðX 1;t ; . . .; X q;t Þ T are given by ( and Y Ã t ¼ Y tþq y ; t ¼ 1; . . .; N ; ðN n À q y Þ: Furthermore, for a given quantile level s, we define as the set of relevant predictors that truly influence Q Y Ã t ðsjxÞ. We assume that sparsity of the predictors exists, i.e., jM s j ( q which means that only a small subset of predictors influences Q Y Ã t ðsjxÞ.

Semiparametric (SP) Conditional Quantile Averaging
For q large (many predictors), the semiparametric (SP) quantile prediction averaging method consists of the following three steps.
1. Using gradient boosting, obtain b h 0;j ðsjx j Þ as a nonparametric estimate of a marginal conditional quantile, defined as where c 0;j ðsÞ are weights depending on the quantile level s.
2. Using b h 0;j ðsjX j;t Þ as regressors and applying penalized quantile regression, obtain the set c f M s ¼ fj : b c 0;j ðsÞ 6 ¼ 0g. This is an estimate of the set of relevant one-step ahead predictors f M s ¼ fj : c j ðsÞ 6 ¼ 0g identified by the averaging model in (2).
Then, compute an estimate b cðsÞ ¼ ðb c 1 ðsÞ; . . .; b c q ðsÞÞ T of the q Â 1 parameter vector of optimal weights cðsÞ ¼ ðc 1 ðsÞ; . . .; c q ðsÞÞ T by considering the following weighted L 1 -penalized quantile estimator b cðsÞ ¼ arg min whereq s ðzÞ ¼ fs À Iðz\0Þgz is the quantile 'tick' loss function, IðÁÞ the indicator function, k [ 0 is a tuning (or penalization) parameter, w j ¼ jb c Ã j j Àc , c [ 0, and b c Ã j is an initial parameter estimate. Observe that for w j ¼ 1, 8j, (3) becomes the usual LASSO (L) penalty. For the adaptive LASSO (aL) penalty function used in this paper, the typical choice is given by w j ¼ ðb c ðLÞ j þ 1=N Þ À1 . 3. Finally, using all the ingredients given above, compute the one-step ahead penalized averaged (PA) and aL-based quantile forecast of Y nþ1 , i.e., We adopt a prediction-based criterion for the selection of k. For a given k, let b c k ðsÞ ¼ À b c 1;k ðsÞ; . . .; b c q;k ðsÞ Á T be the penalized estimates and j c f M k;s j be the corresponding number of non-zero estimates. We choose k by minimizing the following high-dimensional criterion where C N ¼ logðqÞ logðlogðN ÞÞ= logðN Þ.

Hybrid (H) Conditional Quantile Averaging
In addition to the forecasting approach of Section 2.2, there are also situations were a marginal conditional quantile forecast of Y tþ1 can be computed from a fully specified parametric DGP. Combining such a parametric quantile with the marginal semiparametric conditional quantile e Q Y Ã t ðsjxÞ in (2) is the first step in a so-called hybrid (H) quantile averaging method. To formalize, we assume that a parametric model is available with the corresponding estimated marginal quantile functions f b h 1 ðsjx 1;t Þ; . . .; b h q Ã ðsjx q Ã ;t Þg ðq Ã ! 1Þ. Then, we define a hybrid extension of e where fc 0;j ðsÞg q j¼1 and fc j ðsÞg q Ã j¼1 are sets of weights depending on the quantile level s, which we summarize by the set fc ðHÞ u ðsÞg qþq Ã u¼1 with q þ q Ã q y þ q z þ q Ã . The design matrix X t consists of the set of potential predictors: its uth component will be denoted by b h 3 Forecasting Procedures In Sect. 5, we report recursive prediction results for the following six forecasting methods.
(  (4), and the set of relevant predictors consists of q marginal conditional quantiles. Thus, PA-aL only considers the first component of (6). (4) Penalized Linear Quantile Regression (P-Lin-QR) with L-based penalty: This is the well-known quantile estimation method of Koenker (2005). It does not allow for nonlinearities in the predictors. In this case the forecast results are based on a set of q marginal conditional quantiles. (5) Additive QR (Ad-QR): This method is a generalization of method (3) where each additively entered predictor has a non-parametric (possibly nonlinear) effect. Unlike methods (2)-(4), we do not conduct predictor selection for Ad-QR. Instead, we use the predictor selection results of method (3) and implement a low dimensional additive model. (6) Hybrid with aL-based penalty (H-aL): The conditional quantile forecast is defined by (8). In this case the forecasting results are based on a set with ðq þ q Ã Þ marginal conditional quantiles. In Sect. 4, we extend the set of potential predictors with another set of marginal conditional quantiles representing a parametric approximation for the GARCH-type model under study.
Note, methods (1) - (3) and (6) are logically connected in the sense that they are based on the concept of averaging forecasts, assigning weights (equal or varying) to all predictors. Forecasting methods (4) and (5) are popular approaches used in the quantile regression literature. Unlike methods (4) and (5), penalized averaging does not assume that Q Y nþh ðsjX n Þ has an additive structure.

Simulation Design
In Sect. 5, we evaluate/compare the multi-step out-of-sample forecasting performance of the quantile prediction methods using a recursive forecasting strategy with an increasing window size for parameter re-estimation. More concretely, the h ¼ 1 conditional quantile forecast is based on a method with parameter estimates using fY t g 300 t¼1 observations, the second one on a parameter vector using fY t g 301 t¼1 observations, and so on. The last one is based on a parameter vector using fY t g 399 t¼1 observations. Thus, in total there are N Ã ¼ 100 one-step ahead out-of-sample conditional quantile forecasts, and where at each step of the recursive forecasting strategy the available in-sample size N is updated with one observations, starting from N ¼ 300 till N ¼ 399. No lagged values of Y t are included in the set of predictors, i.e., q y ¼ 0. This implies that at each step of the recursive strategy n N and q q z . For h [ 1, the recursive scheme starts at observation number N À h þ 1 and ends at N À h þ N Ã .
where g j ðÁÞ is an unknown real-valued function. In (9), we assume that fV tþ1 g are i.i. d. random variables from U(0, 1) and F À1 m ðÁÞ is the inverse of the cumulative distribution function (CDF) of the Student t m distribution with m ¼ 4 degrees of freedom. So, the model has heavier tails than a traditional GARCH model with Gaussian innovations.
Assuming that the root of 1 À bz ¼ 0 lies outside the unit circle, r tþ1 can be written in an infinite-dimensional ARCH(1) representation. Under some regularity conditions (Xiao and Koenker, 2009) the coefficients of this model decrease geometrically. Let p ðp ( n À hÞ denote a truncation parameter which depends on the number of in-sample observations. Then, the finite-order approximation of the ARCHð1Þ model is given by where a 0 ¼ a 0 =ð1 À bÞ, f j ðÁÞ ¼ g j ðÁÞ=ð1 À bÞ ðj ¼ 1; . . .; qÞ, and a u ¼ a 1 b uÀ1 ðu ¼ 1; . . .; pÞ. To ensure the identifiability of r tþ1 , we normalize (11) so that a 0 ¼ 1.
The truncation parameter p is small relative to the sample size n, but large enough to avoid bias of the parameter estimates. Xiao and Koenker (2009) use p ¼ 3n 1=4 in their study of conditional quantile estimators for GARCH models, which in our case gives p ¼ 12 or 13 with sample sizes n ¼ 300; . . .; 399. However, being conservative, we set p ¼ 15 and q ¼ 28 throughout the simulation study. Thus, the design matrix X t for the SP approach of Sect. 2.2 consists of d ¼ p þ q ¼ 43 predictors. Moreover, we consider the following set of parameter values: a 0 ¼ 0:1, a 1 ¼ 0:3 and b ¼ 0:5.
In each step of the recursive forecast strategy, the one-step ahead quantile forecast is given by where b h ðSPÞ j ðsjX j;N Àpþ1 Þ are estimates of the recursively obtained marginal conditional quantiles h ðSPÞ j ðsjX j;N Àpþ1 Þ. Parametric: For the parametric (P) model, we adopt the following avGARCH (1, 1) process Y tþ1 ¼ r j;tþ1 e tþ1 ; e tþ1 $ i:i:d: N ð0; 1Þ; ð15Þ r j;tþ1 ¼ n 0 þ / j Z j;t þ n 1 jY t j þ n 2 r j;t ; ðj ¼ 1; . . .; qÞ; where Z 1;t and Z 2;t are generated as in (13), and where the remaining Z j;t 's are again i. i.d. U(0, 1) random variables. Using the apARCH model specification in the R-rugarch package, we estimate the parameters in (15)-(16) for j ¼ 1; . . .; q. The conditional quantiles are given by where U À1 ðsÞ is the sth quantile of the CDF of e t , and b r j;tþ1 is an estimate of r j;tþ1 . These marginal parametric quantile estimates and the marginal semiparametric ðsjZ j;t Þ serve as "pilots" of the H-aL conditional quantile estimates. Note that in the case of the hybrid forecasting method, the design matrix X t consists of q þ d ¼ 28 þ 43 ¼ 71 potential regressors, i.e.

Expected loss
Recall that q y ¼ 0. Then, let fe be the set of h-step ahead out-of-sample quantile prediction errors (QPEs) obtained by the recursive forecasting strategy, where the superscript ðÁÞ refers to one of the six methods discussed in Sect. 3, Y Ã t ¼ Y tþhþpÀ1 , and N Ã ¼ 100 for all values of h. As a measure of performance, we calculate an average of the tick loss function q s ðÁÞ of QPE values, This average is an estimate of the expected loss L

Encompassing
The topic of combining conditional quantile forecasts can also be investigated via the principle of "encompassing" (ENC). To this end, consider the set of h-step ahead Under a sufficient weak dependence condition (and additional regularity conditions), it can be shown (Koenker, 2005) that within the context of conditional quantile autoregression and for fixed values of s and h, ffiffi ffi n p À e Q n;h ðsÞ À e Q h ðsÞ  The kernel-type matrix b R 1;h builds upon a method proposed by Powell (1991), with h n ¼ mn À1=3 ðm [ 0Þ a bandwidth parameter satisfying h n ! 0 and nh 2 n ! 1 as n ! 1. 3 For a fixed value of the forecast horizon h, it follows from (White, 2001, Thm. 4  in the left tail ðs ¼ 0:05 and s ¼ 0:1) and the right tail ðs ¼ 0:9 and s ¼ 0:95) of the h-step ahead distribution of the QPE values (h ¼ 1 and h ¼ 6) for independent (r ¼ 0) and dependent (r ¼ 0:8) exogenous predictors. The main findings are summarized as follows.

QPE Loss Values
1. When r ¼ 0 and r ¼ 0:8, the estimated QPE loss ratios are virtually the same for h ¼ 1 and h ¼ 6. 2. Relative to the EW method, the H-aL method beats all other methods in terms of the lowest QPE loss ratios, irrespective of the level of dependence r between the predictors. Compared to all other methods, the reduction in loss is less pronounced for the Ad-QR method. 3. As revealed by their lower loss ratios, there are differences between the h ¼ 1 forecasting performance of all methods in the right tail and the left tail of the distribution of the QPE values. The differences are no longer present for h ¼ 6. 4. For all methods the QPE loss ratios are much smaller at the boundaries s ¼ 0:05 and s ¼ 0:95 than at, respectively, s ¼ 0:1 and s ¼ 0:9. With the exception of H-aL, the differences in QPE loss ratios are small for all other quantile forecasting methods at s ¼ 0:5 (not shown here). In this case the set of exogenous predictors fZ j;t g consists of 28 i.i.d. U(0,1) random variables.
Equal Forecast Accuracy Table 2 shows averages (based on 50 replications) of computed p-values of the DM ðEW;MethodÞ s;h test statistic for r ¼ 0, with h ¼ 1 and h ¼ 6. Except for the PA-L method at s ¼ 0:9 and h ¼ 6, all null hypotheses are rejected at a 5% nominal significance level, implying that the use of a particular nonparametric conditional quantile method yields statistically more accurate quantile forecasts than using the EW method. Again, we note that the dependencies between the exogenous predictors have hardly any impact on the reported p-values. Overall, the PA-aL method has the lowest average p-value but the differences with the P-Lin-QR, Ad-QR, and H-aL methods are minimal. The DM ðEW;MethodÞ s;h test results are generally the same for r ¼ 0:8 (not shown here).
Encompassing Table 3 shows the average p-values of the ENC ðiÞ s;h ði ¼ 1; 2Þ test statistics for h ¼ 1 and h ¼ 6. Here we make pairwise comparisons between quantile forecasts obtained by methods 1)-5) relative to quantile forecasts obtained by the H-aL method. Two general results are evident from the table.
1. In the left tail of the quantile forecast distribution neither H 1;0 and H 2;0 are rejected at a 1% nominal significance level for h ¼ 1 and h ¼ 6. This indicates that quantile forecasts of the H-aL method encompasses those of methods 1)-5).
Overall, the largest evidence for not rejecting H i;0 ði ¼ 1; 2Þ seems to occur for quantile forecasts obtained by the PA-L method. 2. Very different conclusions emerge for the right tail of the quantile forecast distribution. In particular, at a 5% nominal significance level both ENC ðiÞ s;h tests plainly reject H i;0 ði ¼ 1; 2Þ for h ¼ 1 and h ¼ 6, indicating that quantile    Table 1 forecasts obtained by the H-aL method do not encompass quantile forecasts obtained from methods 1)-5).
6 Empirical Application

Data and Forecast Procedure
In this section we make pairwise comparisons between quantile forecasts obtained by methods 2)-6) relative to quantile forecasts obtained by the EW method. The time series of interest is the risk premium of the monthly S&P 500 index, denoted by R t . The set of potential exogenous predictors consists of q z ¼ 14 macroeconomic variables covering the time period 1951:01-2016:12 (792 observations). 4 Specifically, we consider the following predictors: dividend-price ratio (DP); dividend yield (DY); earnings-price ratio (EP); dividend-payout ratio (DE); equity risk premium volatility (RVOL); book-to-market ratio (BM); net equity expansion (NTIS); treasury bill rate (TBL); long-term yield (LTY); long-term return (LTR); term spread TMS); default yield spread (DFY); default return spread (DFR); and inflation (NFL). These series have been the subject of conditional quantile prediction studies by Lima and Meng, (2017), Meligkotsidou et al. (2014Meligkotsidou et al. ( , 2019, and De Gooijer and Zerom (2020) among others. To obtain estimates of the marginal parametric conditional quantiles of the R t series, we adopt a so-called time-varying mean model with exponential GARCH-Z volatility, i.e., R tþ1 ¼ b 0;j þ b 1;j Z j;t þ e j;tþ1 ; e j;tþ1 $ i:i:d: N ð0; r 2 j;tþ1 Þ; logðr 2 j;tþ1 Þ ¼ d 0;j þ d 1;j Z j;t þ d 2;j logðr 2 j;t Þ þ d 3;j e j;t r j;t þ d 4;j e j;t r j;t ; ðj ¼ 1; . . .; q z Þ: This is a very wide model specification from the literature on the predictability of the S&P 500 risk premium; see, e.g., Cenesizoglu and Timmermann (2012). Thus, with the 14 semiparametric marginal quantiles, the complete portfolio of predictors for the hybrid quantile averaging method consists of 28 potential variables.
We use a recursive forecasting strategy similar to the one discussed in Sect. 4. For h ¼ 1, the first forecast origin is 1965:12 and hence the first in-sample (estimation) period is from 1951:01 to 1964:12 ðN ¼ 168Þ: The second in-sample covers the period 1951:01-1965:01 ðN ¼ 169Þ. The last in-sample covers the period 1951:01-2016:11 ðN ¼ 791Þ. Thus, in total there are N Ã ¼ 624 one-step ahead out-of-sample conditional quantile forecasts. For h ¼ 6, the first forecast origin is 1964:07, and the last forecast origin is at 2016:06. for s ¼ 0:05; 0:1; 0:9, and 0.95, and h ¼ 1 and h ¼ 6. With the exception of Ad-QR, the values of the estimated QPE loss ratios are all below one, irrespective of the forecast horizon h. This indicates that there are improvements in conditional quantile forecasts over out-of-sample quantile forecasts obtained by the EW method. Overall, the best forecast performance seems to occur for P-Lin-QR and H-aL. Interestingly, for s ¼ 0:05 and s ¼ 0:95 the P-Lin-QR method has slightly lower QPE loss ratios than the H-aL method irrespective of the value of h. In this respect, recall that conditional quantiles based on the P-Lin-QR method are best for linear functionals of predictors. So the choice of s can affect the ranking of P-Lin-QR and H-aL.

Empirical Results
The bottom panel of Table 4 shows p-values of the DM ðEW;MethodÞ s;h test statistic. Note, that the null hypothesis of equal forecast accuracy is not rejected for PA-aL and Ad-QR for all values of s and h. We also see that the p-values of DM ðEW;PAÀLÞ s;h for h ¼ 1; 6, and s ¼ 0:9 are close to one, indicating that there is no significant difference between the forecast performance of EW and PA-L. Overall, the top and bottom panels confirm the good performance of the H-aL method as opposed to the EW method of combining quantile forecasts. Table 5 shows that the gain in forecasting performance of the H-aL method mainly comes from including the predictor variables DFR and RVOL in the parametric component. By contrast, the semiparametric component includes a variety of predictor variables with significant weights on the set of potential predictors. Finally, as revealed by their selection frequencies less than 5%, there is evidence that the variables DE, BM, and LTY hardly contribute to the predictability of R t .

Conclusion
In the fist half of the paper, we have explored the multi-step ahead quantile forecasting performance of a hybrid method versus quantile forecast results obtained by five alternative, non-hybrid methods, via a Monte Carlo study. The study yields two main findings. First, relative to quantile forecasts obtained from simple averaging quantile forecasts (i.e., the EW method), the hybrid H-aL method is more accurate in terms of the lowest QPE loss ratios, irrespective of the level of dependence between the predictors. This evidence is strongest in the tails of the forecast distribution of the response variable, and most likely due to the benefits of including information from multiple predictor variables. Second, quantile forecasts obtained by the H-aL method do not encompass quantile forecasts from five competitors in the right tail of the quantile forecast distribution. The second half of the paper provided an empirical application to forecast the risk premium of the monthly S&P 500 index. Very similar results are obtained as in the Monte Carlo simulation study. In particular, our findings strongly suggest that it is possible to use pre-selected macroeconomic predictors to produce better out-ofsample quantile forecasts than the EW combination approach which uses fixed, timeinvariant, weights of conditional quantile forecasts. The H-aL method on the other hand uses stochastic rather than fixed weights. As such, we showed that variable selection and random weights jointly can achieve forecast efficiency gains in the tails of the quantile return distribution. We also provided some insight as to where the gains of the hybrid method come from.
Acknowledgements The author would like to acknowledge an anonymous referee for constructive suggestions and comments.
Funding Funding not applicable. No funding available.

Declarations
Competing Interests The author declares that he has no conflict of interest. Compliance with ethical standards The author declares that he has no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.