Modeling risk dependence and portfolio VaR forecast through vine copula for cryptocurrencies

Risk in finance may come from (negative) asset returns whilst payment loss is a typical risk in insurance. It is often that we encounter several risks, in practice, instead of single risk. In this paper, we construct a dependence modeling for financial risks and form a portfolio risk of cryptocurrencies. The marginal risk model is assumed to follow a heteroscedastic process of GARCH(1,1) model. The dependence structure is presented through vine copula. We carry out numerical analysis of cryptocurrencies returns and compute Value-at-Risk (VaR) forecast along with its accuracy assessed through different backtesting methods. It is found that the VaR forecast of returns, by considering vine copula-based dependence among different returns, has higher forecast accuracy than that of returns under prefect dependence assumption as benchmark. In addition, through vine copula, the aggregate VaR forecast has not only lower value but also higher accuracy than the simple sum of individual VaR forecasts. This shows that vine copula-based forecasting procedure not only performs better but also provides a well-diversified portfolio.


Introduction
In finance and insurance, one of the major and challenging issues is managing quantitative risk, specifically forecasting future risk. Risk forecast is not only important for reserving capital but also for anticipating the worse risk. Risk in finance may come from (negative) asset returns whilst, in insurance, a typical risk is a payment loss. It is often that we encounter several (dependent) risks, in practice, instead of a single risk. Dependent random risks or losses occur in many applications and have challenging statistical features, see e.g. Embrechts et al. [1,2], Gräler et al. [3], McNeil et al. [4], Naifar [5], Patton [6], Trabelsi [7], Usman et al. [8] and Zhang et al. [9].
There are several interesting topics that relate to dependent risks. The first one is construction of dependent risks either as an aggregate risk model or a multivariate risk model, see e.g. Kim and Kim [10]. The second topic lies on the fact that dependent risks have a technical problem with regard to finding an exact form of distribution function (cdf) or probability function (pdf). Both topics above eventually bring us to learn and employ a more sophisticated a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 method in dependent risks namely copula. Copula is a system that does not only accommodate either non normal or unidentical marginal distributions into a uniform distribution but also simplify number of parameters of a joint distribution. Meanwhile, a more step method than a copula is vine copula. It is also a system that constructs pair-copula with high flexibility when decomposing conditional distribution, e.g. Aas et al. [11] and Trucíos et al. [12]. Furthermore, vine copula provides us informative and complete dependence structure. It is well known that understanding dependence is an important step to make strategy on diversification.
Modeling dependence through vine copula approach is basically aimed to broaden flexibility of dependence structure among random risks, compared to classical dependence measures of Pearson's ρ and Kendall's τ or even copula. The recent use of vine copula modeling may be found in (energy) economic and finance applications. For instance, Mejdoub and Ghorbel [13] investigated conditional dependence between oil price and renewable energy stock prices and considered threshold Generalized Autoregressive Conditional Heteroscedastic (GARCH) model, see also Trabelsi [7] for tail risk dependence between oil and stocks of oil-exporting countries. Kumar et al. [14] examined conditional dependence among not only energy commodities but also agricultural and precious metals commodities. Furthermore, Ç ekin et al. [15] studied dependence structure among economic policy uncertainty (EPU) of Latin American countries. Meanwhile, Hernandez et al. [16] compared risk of portfolio: Gulf Cooperation Council (GCC) Islamic and conventional bank indices. They studied tail asymmetric dependence among Islamic banks' relationship. In addition, Usman et al. [8] explored dependence modeling between Islamic and conventional stocks through copula whilst Naifar [5] employed Archimedean copulas to model tail dependence structure between Islamic bonds and stock market.
This paper considers risk in finance defined as (negative) asset returns. In particular, we construct a dependence modeling for financial risks and form a portfolio risk. Whilst marginal risk is assumed to follow a Generalized Autoregressive Conditional Heteroscedastic (GARCH) model of order one, forecasting future risk is carried out by risk measure of Value-at-Risk (VaR). VaR, along with Conditional VaR, has be applied to both single risk and dependent risks to make use in practice, see e.g. Nieto and Ruiz [17] for latest review on VaR and its backtesting. Basically, VaR is a maximum tolerated risk or loss for either single or aggregate risk at a given level of confidence. VaR forecast is crucial for assessing the performance of financial institution. Embrechts et al. [1,2] emphasized that VaR is the industry and regulatory standard for risk capital calculation in both banking and insurance.
We pay particular attention to portfolio risk of cryptocurrencies returns; particularly, Bitcoin (BTC), Ethereum (ETH) and Litecoin (LTC) are considered. Cryptocurrency has been one of the major interests among financial practitioners, investors, academia even policy makers. Bitcoin, since introduced by Nakamoto [18], has shown a dramatic increasing value during a year period of 2017 and 2018 for about one and a half times ( [19]). It has been a prominent digital asset ever since. However, there is still debate whether cryptocurrencies are defined as currency, commodity or investment asset. Jiménez et al. [20] argued that Bitcoin is a digital asset or investment as known before (bonds and equities). Their interests are on volatility clustering and leptokurtosis which are typical in asset returns. However, we may observe an extreme event on BTC that leads to greater market instability. There are some financial implications on cryptocurrencies due to uncontrolled monitoring by monetary regulator, see. e.g. Corbet et al. [21] who recorded pricing bubbles in BTC and ETH. They also claimed that cryptocurrencies may become source of financial instability.
Boako et al. [19] and Trucíos et al. [12] are among authors recently used returns of cryptocurrencies to illustrate vine copula modeling and risk measures forecasting. It is interesting to classify several areas done by authors. The first area is cryptocurrencies testing in order to have a clear position: as a currency, commodity or asset; see e.g. White et al. [22] and Kwon [23]. The related area to the first is pricing formation of cryptocurrencies and other possible measurement such as ratio of gold to platinum (GP) prices forecast Bitcoin return, e.g. Huynh et al. [24]. The second area lies on cryptocurrencies empirical facts like other returns along with return-volatility correlation, volatility clustering or asymmetric volatility, see e.g. Bouri et al. [25], Klein et al. [26], Boako et al. [19] and Wajdi et al. [27]. Such relationship may be measured by what so called "volatility surprise" or unexpected volatility. It is the difference between squared innovation and conditional volatility, see Bouri et al. [28]. Furthermore, interdependence that accounts time and frequency and market interconnection among cryptocurrencies may be explored through wavelet-based approaches CWT, continuous wavelet transformation, and XWT, cross wavelet transform) as carried out by e.g. Qureshi et al. [29]. Interdependence may also be observed between cryptocurrency and energy or agricultural commodities. Ji et al. [30], for instance, examined information spillovers via entropy-based method among both cryptocurrencies and commodities. They found that, first, it changes over time. Secondly, unlike the spillover of cryptocurrencies, energy commodities spillover contribution to the system is dependent on their price dynamics. The third area is concerned about economic or financial implication of cryptocurrencies, see e.g. Bouri et al. [31] who studied global financial stress, whilst Trucíos et al. [12] argued that it is a shelter against economic and financial turmoil.
In this paper, we show the portfolio or aggregate VaR forecast by considering vine copulabased dependence among individual returns. Different from the work of Boako et al. [19] and Trucíos et al. [12] who found VaR forecast by using classical historical simulation (HS) method, our forecast calculation considers "estimative" VaR forecast as in Kabaila and Syuhada [32,33] and Syuhada [34]. The aggregate VaR forecast is compared to the simple sum of individual VaR forecasts that actually considers perfect dependence assumption. To evaluate their performance, we assess their accuracy by adopting several backtesting methods as recently used by Syuhada [34] and Jiménez et al. [20]. In addition, comparing the aggregate VaR and the simple sum of individual VaRs also leads us to investigate and measure benefits of the portfolio diversification.

Returns and marginal risk model
The marginal risk is (negative) returns of three cryptocurrencies defined as below where P B t , P E t and P L t denote (closing) price at time t for Bitcoin (BTC), Ethereum (ETH) and Litecoin (LTC), respectively.
We assume a marginal risk model of Generalized Autoregressive Conditional Heteroscedastic [35] of order one or GARCH(1,1) for each risk process. This is mainly due to dynamic volatility property found from each risk. Specifically, the GARCH(1,1) models for (negative) returns of Bitcoin {X t }, Ethereum {Y t } and Litecoin {Z t } are, respectively, given by where ω x , ω y , ω z > 0 and δ x , δ y , δ z , β x , β y , β z � 0, for t � 0. The restriction on persistence parameter δ x + β x < 1, δ y + β y < 1 and δ z + β z < 1 is needed to ensure the stationarity of all processes {X t }, {Y t } and {Z t }, respectively. We assume that innovation {ε x;t , t � 0} is white noise. In addition, innovation ε x;t , t � 0, and volatility σ x;t as well as ε x;t , t � 0, and information, up to time (t − 1), F x;tÀ 1 , are independent. Note that such assumptions also apply to {ε y;t , t � 0} and {ε z;t , t � 0}. The estimates for each innovation is calculated and the goodness-of-fit procedure for innovation distribution needs to be carried out. From X t = σ x;t ε x;t , the innovation ε x;t is formulated as ε x;t = X t /σ x;t . Thus, ε x;t may be estimated bŷ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Meanwhile, we assume standard Student's t distribution for such innovation. To do so, suppose that a random variable T x has Student's t distribution with degrees of freedom ν x 2 (0, where B(�, �) is beta function. Note that its mean is EðT x Þ ¼ 0, for ν x > 1, whilst its variance, VðT x Þ ¼ n x n x À 2 , is positive and finite, for ν x > 2. By defining ε x;t ¼ T x ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi VðT x Þ p ¼ T x ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi n x À 2 n x s ; we obtain Eðε x;t Þ ¼ 0 and Vðε x;t Þ ¼ 1 and the innovation ε x;t is said to have standard Student's t distribution with probability function Parameter function ffi ffi ffi ffi ffi ffi ffi n x À 2 n x q may be viewed as scale parameter, so that ε x;t � Student0s t 0; ffi ffi ffi ffi ffi ffi ffiffi Now, by applying standard Student's t distribution to GARCH(1,1) model, the conditional probability function of X t , given information F x;tÀ 1 , is f X t jF x;tÀ 1 ðxÞ ¼ 1 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi s 2 x;t ðn x À 2Þ In other words, X t jF x;tÀ 1 is Student's t distributed with parameter n x ; s x;t ffi ffi ffi ffi ffi ffi ffi The explicit form of the conditional probability function above leads us to employ the conditional maximum likelihood method, as in McNeil et al. [4], to estimate its parameters. Note that this procedure is analogous to Y t jF y;tÀ 1 and Z t jF z;tÀ 1 .

Copula-based modeling
Copula and dependence measures. Suppose that a continuous bivariate random variable (X, Y), representing a joint risk, has marginal risk distribution function F X and F Y , respectively. Suppose also that U = F X (X) and V = F Y (Y) so that they are uniformly distributed over unit interval, i.e. U; V � Uð0; 1Þ. Copula C U,V is a joint distribution function of (U, V) such that where C U,V : [0, 1] 2 ! [0, 1]. Its corresponding probability function is which is called copula density, see e.g. Nelsen [36] and McNeil et al. [4]. Furthermore, Sklar's theorem stated that joint distribution function F X,Y of (X, Y) may be determined through copula C U,V of marginal distribution functions F X and F Y , i.e. F X;Y ðx; yÞ ¼ C U;V ½F X ðxÞ; F Y ðyÞ�; ðx; yÞ 2 R 2 : Sklar's theorem has shown us that the use of copula provides many choices of joint distribution function of (X, Y). The corresponding joint probability function f X,Y of (X, Y) is given by To measure dependence between X and Y, the dependence measures of Pearson's ρ and Kendall's τ are required. According to Schweizer and Wolff [37], such dependence measures are defined as It is shown from Eq (5) that Pearson's ρ X,Y depends not only on copula but also on marginal of X and Y. Thus, Pearson's ρ is invariant only under linear transformation. Meanwhile, from Eq (6), it is shown that Kendall's τ X,Y depends only on copula. In addition, τ X,Y = τ U,V . Thus, Kendall's τ is invariant under both linear and nonlinear transformations. We may say that Kendall's τ is copula-based dependence measure. Note that for parameter θ of copula C U,V , Kendall's τ U,V is a function of θ i.e. τ U,V (θ).
There are several copulas commonly used such as Archimedean and elliptical copulas. Based on Joe [38], examples of Archimedean copulas are Clayton, Gumbel and Frank with the following function and the corresponding density: respectively, where θ 2 (−1, 1). Note that F(�) and F 2 (�, �;θ) are distribution function of standard univariate normal and standard bivariate normal random variables, respectively. Furthermore, F T,T 0 (�, �;ν, θ) is joint distribution function of identical Student's t random variables (T, T 0 ) with marginal distribution function F T (�;ν) and degrees of freedom ν 2 (0, 1). The corresponding densities of Gaussian and Student's t copulas are, respectively, given by ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 À y 2 nÞ. The dependence measures of Kendall's τ for Archimedean and elliptical copulas are provided in Table 1. It is shown that Kendall's τ may be expressed in copula parameter θ. In addition, copula parameter θ may also be expressed in Kendall's τ, except for Frank copula.
For equal value of Kendall's τ, we present in Fig 1 the contour plot for the density of Archimedean and elliptical copulas. It may be observed that Clayton, Gumbel and Student's t copulas perform tail dependence. In more detail, Clayton copula is appropriate for lower-tail dependent risks whilst upper-tail dependence may be captured by Gumbel copula. Furthermore, Student's t copula displays symmetrical lower-and upper-tail dependence. Meanwhile, Frank and Gaussian copulas have symmetrical lower-and upper-tail independence.
Copula selection. For a given risk data set, the challenging task is selecting the best copula which fits well to the data. We may do this by considering several criteria. One of them is Akaike Information Criterion (AIC) introduced by Akaike [39]. Suppose that a copula C U,V (�, �;θ) with its corresponding density c U,V (�, �;θ) has parameter θ. Based on data fðu i ; v i Þg n i¼1 of (U, V), the estimate for such parameter may be obtained through maximum likelihood (ML) method with likelihood function LðyÞ ¼ Q n i¼1 c U;V ðu i ; v i ; yÞ. By replacing θ with its estimate,ŷ, we haveĈ U;V ð�; �;ŷÞ as the parametric estimate for C U,V (�, �;θ). The AIC value from such copula is defined as where b is number of parameters in θ. Among several choices of copulas, a copula with the lowest value of AIC may be decided as the best copula.
Note: Kendall's τ is expressed as function of copula parameter θ and such θ is expressed as function of Kendall's τ.
For Frank copula, its Kendall's τ is depend on Debye function of D 1 ðyÞ ¼ 1 y R y 0 t e t À 1 dt whose inverse function has no explicit form. https://doi.org/10.1371/journal.pone.0242102.t001 In addition, we employ other criteria by considering empirical version of copula. The empirical copulaĈ n , for the data fðu i ; v i Þg n i¼1 , is defined aŝ where I A is indicator function on set A. Through graphical approach, the best fitting copula may typically be selected by comparing the surface of such empirical copula to that of To give more convenient interpretation, we consider its univariate version i.e.
Such empirical function is the nonparametric estimate of distribution function, Kðq; yÞ ¼ PðQ � qÞ, for the so-called Kendall's transform, Q = C U,V (U, V;θ), defining Kendall's τ in Eq 6. According to Genest and Rivest [40], this leads us to visualize the curve for K n ðqÞ and its corresponding function of λ, defined bŷ along with the parametric version derived fromĈ U;V ðu; v;ŷÞ. Furthermore, we may perform a goodness-of-fit (GoF) test for such copula. It is carried out by testing

Vine copula-based modeling
Pair-copula construction method. Suppose that (X, Y, Z) is a dependent risk model consisting of three continuous random variables with joint probability function where the marginal risk probability functions are f X , f Y , f Z and c X,Y, Z is a trivariate copula density. Since there is limitation on classes of trivariate copulas, we aim to find joint distribution of (X, Y, Z) through another approach. Joe [42] suggested a decomposition for f X,Y, Z as Since Now, from Eqs (9) and (10), we obtain This shows that the trivariate copula density c X,Y, Z for (X, Y, Z) model may be constructed from bivariate copula densities: c X,Y , c X,Z , c Y,Z|X . In other words, it provides flexible way to determine joint distribution of (X, Y, Z). The first two components c X,Y and c X,Z show dependence of (X, Y) and of (X, Z), respectively, with dependence measure κ X,Y and κ X,Z . Meanwhile, the component c Y,Z|X shows partial dependence of (Y, Z), given X, with partial dependence measure κ Y,Z|X . Note that the dependence measure of κ may be either Person's ρ or Kendall's τ.
Suppose that F X , F Y , F Z denote marginal distribution functions of (X, Y, Z). The copula density c X,Y, Z in Eq 11 is basically expressed as for all ðx; y; zÞ 2 R 3 , where F Y|X and F Z|X are conditional distribution function of Y and Z, given X, respectively. According to Joe [42], they may be defined as This shows that conditional copula density c Y,Z|X is determined through copula C X,Y and C Y,Z . From Eqs (12) and (13), copula density c X,Y, Z is completely given by According to Joe [38], function of C Y|X for each copula C X,Y , from Archimedean and elliptical copulas, is given in Table 2. Now, for u = F X (x), v = F Y (y) and w = F Z (z), copula density c X,Y, Z is given by Table 2. Function of C Y|X .

Copula Function of C Y|X (v|u)
Clayton Graph structure and vine (copula). Suppose that V denotes an empty and finite set. Suppose also that [V] 2 = {{v, v 0 }:v, v 0 2 V} is a set of all pairs of two unordered elements in V. According to Diestel [43], a graph is a pair (V, E) of sets for some E � [V] 2 . The element of V is called node whilst the element of E is an edge. A simple graph T ¼ ðV; EÞ that is connected and has no cycle is called a tree graph.
Based on Bedford and Cooke [44,45] and Kurowicka and Cooke [46], V m is called vine on m elements, with 3. for all j = 2, 3, . . ., m − 1, T j is a tree graph with a set of nodes V j = E j−1 and a set of edges In short, vine is a collection of nested trees where edge of the jth tree is a node of the (j + 1) In other words, two adjacent nodes in the tree graph T j are two adjacent edges in the tree graph T jÀ 1 ". There are two special cases of R-vine that are drawable vine (D-vine) and canonical vine (C-vine). R-vine V m is a D-vine if all nodes in the tree graph T 1 have maximum degrees of 2. Meanwhile, V m is called a C-vine if the tree graph T j has exactly one node with degrees of m − j, for j = 1, 2, . . ., m − 1; in tree graph T 1 , such node is called a root. For an illustration, it is shown in Fig 2 graph structure of D-vine and C-vine, for m = 4 and V 4 ¼ fT 1 ; T 2 ; T 3 g.
Now, an R-vine V 3 is employed to represent dependence model (X, Y, Z) determined through trivariate copula by pair-copula construction method. Here, node is a random variable whilst edge is a bivariate copula added by absolute dependence measure as weight. Such weight is used to find appropriate R-vine i.e. an R-vine where each tree graph has maximum of the sum of weights, according to Dißmann et al. [47].
First, define a complete graph K 3 (or a cricle graph C 3 ) with three nodes representing random variables of X, Y, Z. Let absolute value of each dependence measure κ X,Y , κ X,Z , κ Y,Z be defined as weight of edge that connects two nodes. Then, select a maximum spanning tree graph, that is a tree subgraph of the complete graph K 3 maximizing the sum of weights. We denote the resulted R-vine copula by ðF; V 3 ; C; KÞ, where F is a collection of marginal distribution functions of (X, Y, Z). Meanwhile, C and K consist of bivariate copulas and the corresponding dependence measures, respectively.

Forecasting Value-at-Risk (VaR)
Individual VaR forecast. The one-step-ahead VaR forecast is actually a forecast(ing limit) of future risk or return, X n+1 , given previous information up to time n, F x;n . At a specified confidence level of 1 − α, it may be calculated through the following formula, see e.g. McNeil et al. [4] and Syuhada [34], x;nþ1 ¼ inffx : PðX nþ1 � xjF x;n Þ � 1 À ag; for α 2 (0, 1). Provided that the inverse of distribution function of X n+1 , given F x;n , exists, we may obtain VaR 1À a x;nþ1 ¼ F À 1 X nþ1 jF x;n ð1 À aÞ, according to Syuhada et al. [48]. This means that such VaR forecast is basically the (1 − α)-quantile of (conditional) distribution of the future return.
Since we use standard Student's t for innovation distribution and assume GARCH(1,1) for (negative) return model, the conditional distribution of X n+1 , given F x;n , is Student's t with degrees of freedom ν x and scale parameter s x;nþ1 ffi ffi ffi ffi ffi ffi ffi . This implies that its conditional mean and variance are respectively. Note that Student's t is a symmetrical distribution. Thus, the one-step-ahead VaR forecast for X n+1 , given F x;n , consists of such conditional mean and variance, i.e.
ffi ffi ffi ffi ffi ffi ffiffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi Parameter of model, (ω x , δ x , β x , ν x ), may be replaced by its estimate so that we have the estimative one-step-ahead VaR forecast given by d VaR 1À a x;nþ1 ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi In addition, the estimative ℓ(>1)-step-ahead VaR forecast is as follows x;nþ'À 1 Þn x À 2 whereŝ 2 x;nþ'À 1 ¼ô x þd xX 2 nþ'À 1 þb xŝ 2 x;nþ'À 2 . Note that forecasting VaR for Y t , given F y;tÀ 1 , and Z t , given F z;tÀ 1 , may be carried out through the similar procedure.
Aggregate VaR forecast. For the case of a portfolio or aggregate risk, we aim to forecast VaR for future aggregate risk of given previous information F x;y;z;n . We may employ vine copula to determine the joint distri- denotes its joint probability function having a certain decomposition. Then, the conditional distribution function of S n+1 , given F x;y;z;n , is determined as follows ðx; y; zÞ dz dy dx; s 2 R: By employing vine copula, forecasting the one-step-ahead individual VaRs is simultaneously carried out based on the joint distribution of ðX � nþ1 ; Y � nþ1 ; Z � nþ1 Þ. We collect all of them into a vector denoted by Then, we may forecast the one-step-ahead aggregate VaR, aggVaR 1À a s;nþ1 , by taking such individual VaR forecasts into account. We do this by first observing that where P ¼ (17), through simple algebraic operation we have This means that the aggVaR forecast above incorporates interactions between different returns by introducing their dependence measures. Note that the estimative aggVaR forecast may be obtained by using the individual estimative VaR forecasts and the estimate of dependence measures.
For perfect (positive) dependence, i.e. κ x,y = κ x,z = κ y,z = 1, it is obvious that the aggVaR forecast is equal to the simple sum of individual VaR forecasts written by Thus, simplesumVaR may simply be decided as the aggVaR forecast when the worse cases of all of our risks always occur simultaneously, see e.g. Li et at. [49] and Embrechts et al. [2]. In other words, the risks have dependence representation of the form M(u, v, w) = min(u, v, w) which may be called perfect-dependence copula. However, since the dependence measures satisfy |κ x,y |, |κ x,z |, |κ y,z |�1, we have aggVaR 1À a s;nþ1 � simplesumVaR 1À a s;nþ1 : This means that the aggVaR forecast in Eq (17) is bounded from above by simplesumVaR in Eq (18). The weaker the dependence among marginal risks, the lower the aggVaR forecast. This is inline with the fact that M(u, v, w) is the upper bound for all classes of copulas, see e.g. Joe [38], including pair-copula which determines the vine copula-based dependence among marginal risks for our portfolio risk.
In addition, it is interesting to investigate diversification benefits of using aggVaR and sim-plesumVaR. As in Li et at. [49], such diversification benefits may be measured by a so-called diversification coefficient (DC) defined as

PLOS ONE
From Eq (19), it may be understood that such coefficient is getting higher as the dependence among marginal risks is weaker. This makes the portfolio better diversified.
To find the ℓ-step-ahead aggregate VaR forecast, we consider the vector of ℓ-step-ahead individual VaR forecasts derived simultaneously. We denote it by each component is obtained similarly by Eq (15). Then, we have The calculation of DC for this aggVaR forecast is analogous to Eq (19) by involving the simple sum of individual VaR forecasts of Eq (20).

Backtesting VaR
After we find the (aggregate) VaR forecast, the assessment of the forecast accuracy is required. It may be carried out through backtesting. From a variety of the backtesting methods, we adopt the methods from Syuhada [34] and Jiménez et al. [20]. We discuss the following methods for the case of VaR 1À a x;nþ1 ; the discussion for the other individual VaR and the aggregate VaR forecasts are similar.
Probability-based backtesting. As stated before, the one-step-ahead VaR forecast is the forecasting limit of future risk, given information of risks in the past. To simply asses the accuracy of VaR forecast, we, therefore, may calculate the probability that the actual future risks do not violate the VaR forecast. It is called coverage probability (CP) defined by The closer the CP to the confidence level, 1 − α, the more accurate the VaR forecast. When the sequence of actual future risks, fx nþ1;k g N k¼1 , is obtained, we may define fI nþ1;k g N k¼1 , where I nþ1;k ¼ I fx nþ1;k �V aR 1À a x;nþ1 g , for k = 1, 2, . . ., N. We expect that fI nþ1;k g N k¼1 is the sequence of realizations of Bernoulli random variable with parameter 1 − α. Thus, the value of 1 N P N k¼1 I nþ1;k is required to estimate the CP.
From the sequence of fx nþ1;k g N k¼1 , we may also require the sequence of fI � nþ1;k g N k¼1 , where I � nþ1;k ¼ I fx nþ1;k >V aR 1À a x;nþ1 g , for k = 1, 2, . . ., N. Note that the term 1 refers to failure or violation of the VaR forecast; all terms are expected to be the realizations of Bernoulli random variable with failure proportion equal to α. Comparison of the actual number of failures, P N k¼1 I � nþ1;k , and the expected number of failures, Nα, may be considered to assess the accuracy of VaR forecast. Their ratio is called actual over expected (AE) ratio.
Conditional coverage test. The standard test for assessing the accuracy of VaR forecast is the test whether the failure proportion is equal to α. It is basically carried out through binomial test by considering normal approximation of binomial distribution. However, according to Christoffersen [50], the forecast accuracy may be assessed by determining whether the sequence of fI � nþ1;k g N k¼1 consists of independent and identical realizations with failure proportion equal to α. This means that we need to combine the test of proportion of failure (PoF) and the test whether this sequence satisfies the conditional coverage independence (CCI). Their combination is called conditional coverage (CC) test.

The PoF and CCI tests employ likelihood ratio (LR) test statistic. The LR statistic for PoF test is defined by
;k is the actual number of failures andâ ¼ N f N is its actual proportion. Meanwhile, the LR statistic for CCI test is where N ij is number of the term i followed by the term j and P ij is the corresponding transition probability, for i, j = 0, 1, from fI � nþ1;k g N k¼1 viewed as the sequence of realizations of a Markov chain model with binary states; π 1 is unconditional probability of the term 1 (failure). Both LR statistics are asymptotically chi-square distributed with 1 degree of freedom. As proposed by Christoffersen [50], the CC test employs the LR statistic defined by LR CC = LR PoF + LR CCI ; it is asymptotically chi-square distributed with 2 degrees of freedom. The null hypothesis of correct model specification fails to be rejected when p-value derived from this statistic is above significance level.
Backtesting through loss function. It is known that VaR is basically the quantile of (conditional) distribution of the future risk. In addition to the use of inverse of its distribution function, VaR as quantile may be formulated through different approach. Based on Kuan et al. [51], VaR 1À a x;nþ1 is the minimizer of the loss function defined by LðaÞ ¼ E½jð1 À aÞ À I fX nþ1 �ag j � jX nþ1 À aj; given F x;n �: In this case, the loss function evaluated at our VaR forecast may be relatively compared to that evaluated at VaR forecast obtained from other method as benchmark. They are called quantile losses. The value of their ratio below one means that our forecasting method outperforms the benchmark.

Daily returns, innovation distribution and marginal risks
The dynamic daily prices and returns for Bitcoin (BTC), Ethereum (ETH) and Litecoin (LTC) for a period of 730 days in 2017-2018 are given in Figs 3 and 4, respectively. The returns follow formula in Eq (1). Whilst the three cryptocurrencies seem to have similar prices behavior, Litecoin has a dramatic increase in day 350 (December 2017) before the decrease in day 400 (February 2018), in comparison to Bitcoin (less dramatic) and Ethereum (gradual increase). In addition, the returns of Litecoin have more high volatility compared to the low and medium volatility as shown in Bitcoin and Ethereum returns, respectively. The property of dynamic volatility significantly appears in all of these returns based on the result of hypothesis test for ARCH effect of Engle [52] in Table 3. This result leads us to assume conditional heteroscedasticity for each of our risk models. Furthermore, the stationarity assumption is also needed based on the result of ADF test in Table 3. In addition, the existence of (inverse) leverage effect may also be considered, especially for Bitcoin whose (negative) return is positively correlated with its squared volatility as visualized in Fig 5. In other words, we also need to assume asymmetric heteroscedastic model in order to capture the feature of asymmetrical volatility. This is in line with assumption in Bouri et al. [25] and Klein et al. [26].
We have assumed a GARCH(1,1), as in Eq (2), for each of our risk models. It consists of dynamic volatility and innovation. Table 4 gives the summary of statistics for the estimates of such innovation, defined in Eq (3), for the returns of Bitcoin, Ethereum and Litecoin. It is observed that the high empirical kurtosis leads us to employ heavy-tailed distribution for innovation.
Furthermore, we visualize the innovation estimates in histogram, Fig 6. Both standard normal and standard Student's t curves for probability function are fitted, where the estimated   degrees of freedom of standard Student's t are given in Table 5. In fact, it may be observed that standard normal distribution is not appropriate for each innovation. This is also described in Fig 7 for its distribution function. Based on AIC values also given in Table 5, standard Student's t distribution has lower value of AIC than standard normal distribution. This confirms the appropriateness of standard Student's t distribution for such innovation, as in Eq (4). Based on the assumption above, parameter estimates of GARCH(1,1) model are given in Table 6. By assuming symmetrical volatility, it is observed that the persistence parameters   Meanwhile, for asymmetric model with considering inverse leverage effect, we need an additional parameter γ to ensure that the positive-signed returns lead a rise in the volatility. In other words, such volatility is modeled by for (instance) Bitcoin returns. From Table 6, the stationarity condition is also satisfied for all asymmetric processes since the persistence parametersd þb þ 1 2ĝ ¼ 0:9247; 0:8207; 0:9239   are below one. For Litecoin, we obtainĝ z ¼ 0 which means that asymmetric GARCH(1,1) is equal to symmetric GARCH(1,1).

Dependence among innovations and the best copula selection
Modeling dependence of (ε x;t , ε y;t , ε z;t ) is used to compute and to find joint distribution of ; Y t jF y;tÀ 1 ; Z t jF z;tÀ 1 Þ. As stated before, suppose that ε x ¼ fε x;t g, ε y ¼ fε y;t g, and ε z ¼ fε z;t g. We calculate the dependence measure for each pair of innovation estimates and also calculate maximum of the sum of absolute dependence measures. Then, the appropriate trivariate copula density c ε x;t ;ε y;t ;ε z;t to determine joint distribution of (ε x;t , ε y;t , ε z;t ) is constructed. Fig 8 shows scatter plot (three dimension) of innovation estimates fðε x;t ;ε y;t ;ε z;t Þg whilst Fig 9 describes scatter plot of each two dimensional pair of innovation estimates: fðε x;t ;ε y;t Þg, fðε x;t ;ε z;t Þg and fðε y;t ;ε z;t Þg. Note that their dependence measures arek x;y ¼ 0:4501,k x;z ¼ 0:4910 andk y;z ¼ 0:5127, respectively. This means that innovation estimates with maximum of the sum of absolute dependence measures are fðε x;t ;ε z;t Þg and fðε y;t ;ε z;t Þg. Thus, the appropriate trivariate copula density c ε x;t ;ε y;t ;ε z;t to determine joint distribution of (ε x;t , ε y;t , ε z;t ) is Note that the first two bivariate copula densities, respectively, correspond to bivariate copula C ε x;t ;ε z;t and C ε y;t ;ε z;t for edges of the tree graph T 1 whilst the last one corresponds to conditional copula C ε x;t ;ε y;t jε z;t for edge of the next tree graph T 2 . The latter copula is defined through copula C ε x;t ;ε z;t and C ε y;t ;ε z;t .
According to Aas et al. [11] and Czado et al. [53], copula parameters are estimated through sequential maximum likelihood method as follows. First, the innovation estimatesε x;t ;ε y;t ;ε z;t are transformed through marginal distribution function F ε x;t ; F ε y;t ; F ε z;t that are The transformed data are then used to estimate parameter of copula C ε x;t ;ε z;t ð�; �; y x;z Þ and C ε y;t ;ε z;t ð�; �; y y;z Þ from several classes. The best copulas for C ε x;t ;ε z;t and C ε y;t ;ε z;t , based on several criteria, are used to find function C ε x;t jε z;t and C ε y;t jε z;t , respectively. Now, define u w;t ¼Ĉ ε x;t jε z;t ðu t jw t ;ŷ x;z Þ; v w;t ¼Ĉ ε y;t jε z;t ðv t jw t ;ŷ y;z Þ 2 ½0; 1�; so that we have the transformed data {(u w;t , v w;t )}. We then use these data to estimate parameter of copula C ε x;t ;ε y;t jε z;t ð�; �; y x;yjz Þ.
The scatter plot of the transformed data {(u t , w t )} and {(v t , w t )} is shown in Fig 10. It may be observed that these transformed data of the estimated innovations show asymmetrical dependence with strong dependence in upper tails. This means that an increase of extremely positive innovation for BTC (or ETH) returns is followed by an increase of that for LTC returns. This empirical fact indicates that Gumbel copula is appropriate for such data. This indication is inline with the result of copula selection based on AIC and goodness-of-fit criteria provided in Table 7. This is because Gumbel copula has the lowest value of AIC and the highest p-value (>0.05) of the goodness-of-fit test for both data {(u t , w t )} and {(v t , w t )}.
In addition, graphical approach visualized in Fig 11 explains that Kendall's distribution function (and the corresponding function of λ) derived from Gumbel copula seem to be close enough to the empirical Kendall's distribution function defined in Eq (7) (and the empirical function of λ n defined in Eq (8)). This means that Gumbel copula fits well to the upper-tail dependent data {(u t , w t )} and {(v t , w t )}.
From the result above, Gumbel copula is the best fitting copula for C ε x;t ;ε z;t and C ε y;t ;ε z;t . By employing Gumbel copula for such copulas, we obtain the transformed data {(u w;t , v w;t )} shown in Fig 12. We now observe that the new transformed data show symmetrical tail dependence. This empirical fact leads us to consider Student's t copula as the best fitting copula for respectively,n x ,n y ,n z . Fig 14 shows its graph structure. As a consequence, dependence model

VaR forecast
Now, the one-step-ahead VaR forecasts for Bitcoin, Ethereum and Litecoin returns based on (a)symmetric GARCH(1,1) model are provided in Table 9. Such VaR forecasts are calculated at 95%, 97% and 99% confidence levels (CLs) under two assumptions: (i) perfect dependence,  as benchmark, and (ii) vine copula-based dependence among different returns. For the former assumption, we simulate perfectly dependent innovations to calculate the previous returns and volatility and, hence, to find the individual VaR forecasts. This procedure is similarly applied for the latter assumption for which the innovations are simulated through vine copula we have   previously determined. Note that the individual VaR forecasts are found simultaneously by Eq (16) for one day only. The forecast accuracy under two assumptions above is assessed and compared through several backtesting methods adopted from Syuhada [34] and Jiménez et al. [20]. These methods include coverage probability (CP), actual over expected (AE) ratio, conditional coverage (CC) test and quantile loss (QL) ratio.
We find that the VaR forecasts for returns by assuming vine copula-based dependence have higher forecast accuracy. This is because the CP for these VaR forecasts are closer to the corresponding CL and their AE ratio are closer to one. Furthermore, this assumption successfully passes the CC test since the resulted p-value is higher than 5% significance level, except for Litecoin at the 97% CL with p-value 0.0322. The choice of significance level may be relaxed to 1% so that the CC test for Litecoin at such CL is also passed. Such results are obtained based on both symmetric and asymmetric GARCH(1,1) models. These show the positive impact of vine copula-based dependence on the VaR forecasts. These are confirmed from the low value of all QL ratios below one. As the CL increases, the QL ratio is getting lower. Meanwhile, the consideration of using asymmetric model do not give better impact on the VaR forecasts. This is because the VaR forecasts derived from both symmetric and asymmetric models do not obviously differ, but the asymmetric model provides higher value of QL ratio. As mentioned before, we have used the empirical data from 1-1-17 till 31-12-18. Now, we compare the VaR forecasts for several steps/days to the empirical data from 1-1-19 till 31-12-19 as visualized in Fig 16. The comparison is made based on (a)symmetric GARCH(1,1) model where the VaR forecasts are calculated simultaneously by Eq (20) through vine copula. According to the result in Table 9, we now provide, in Table 10, the one-step-ahead forecast of portfolio risk by using aggregate VaR (aggVaR) through vine copula in comparison to the simple sum of individual VaRs (simplesumVaR), as the benchmark, for one day only. Their forecast accuracy is compared in terms of CP, AE ratio, CC test and QL ratio. We may observe that the aggVaR forecast is lower in value with better accuracy at each CL. The CP as well as AE ratio are closer to the target and the CC test is successfully passed at 5% significance level. Meanwhile, in calculating simplesumVaR, the rejection of null hypothesis of the CC test is only at low level of significance, e.g. 1%, due to the low p-value. Furthermore, the QL ratio between aggVaR and simplesumVaR is below one and getting lower as the CL increases. These show that calculating simplesumVaR is too conservative and overestimates the portfolio risk. In other words, vine copula-based forecasting procedure outperforms perfect dependence assumption. The use of vine copula also makes the portfolio diversified, with diversification coefficient (DC) about 18%. Considering diversification typically provides a stable portfolio  since overall risk of the portfolio is properly reduced and, hence, risk allocation is allowed. However, as shown in Table 10, the DC is getting lower as the CL increases. This means that the higher the CL, the less stable the portfolio. The results under symmetric and asymmetric GARCH(1,1) models are in line. However, the asymmetric one performs worse with higher value of the QL ratio although the resulted DC is higher. In addition, we calculate the aggVaR forecasts for several steps/days. Such forecasts are visualized in Fig 17. This figure also provides their comparison to the simplesumVaR and the real data of aggregated returns from 1-1-19 till 31-12-19. Meanwhile, the corresponding DCs are displayed in Fig 18.   Fig 17. VaR forecasts of portfolio risk for several steps/days. The aggVaR forecasts are calculated at 95%, 97% and 99% confidence levels (CLs) based on symmetric (solid line) and asymmetric (dashed line) models through vine copula. They are compared to the simplesumVaR (purple) and to the real data (blue). https://doi.org/10.1371/journal.pone.0242102.g017

Conclusion
Model for dependent risks that form a portfolio risk has been constructed. Through vine copula, their marginal risks assumed to follow GARCH(1,1) model are coupled with complete representation through graph structure. This approach provides high flexible pair-copula model since it is able to capture dependence structure of all possible pairs of risks by using different bivariate copulas with the best criterion. As we hope, applying this dependence model to forecast VaR for Bitcoin, Ethereum and Litecoin returns provides good forecast accuracy. Furthermore, the resulted portfolio VaR forecast is low in value with high accuracy instead of simply summing the individual VaR forecasts under perfect dependence assumption. As a consequence, the portfolio is well diversified which means that its overall risk is well managed and reduced. The results under consideration of both symmetrical volatility and asymmetrical volatility in the marginal model are in line. However, the asymmetric model does not perform better although it makes the portfolio more diversified.
For further research, modeling dependence may be carried out for higher-dimensional risk data due to the increasing number of risks in the cryptocurrency market or other markets nowadays. Furthermore, the marginal risk model may be extended to another observable volatility models of GARCH. The use of latent volatility model of Stochastic Volatility Autoregressive (SVAR), see e.g. Han et al. [54] and Syuhada [34], perhaps gives more interesting results. This is due to volatility shock that appears in volatility process. As for risk measure forecast, it is important to consider (i) expected-based risk measure, namely Expected Shortfall or Tail-VaR, and (ii) expectile-based risk measure of EVaR [51]. Whilst the former (VaR) is more probability-based risk measure, the latter (ES/TVaR, EVaR) takes the magnitude of losses into account.