Graphical copula GARCH modeling with dynamic conditional dependence

Lupe S. H. Chan
The Hong Kong University of Science and Technology
and
Amanda M. Y. Chu
The Education University of Hong Kong
and
Mike K. P. So
The Hong Kong University of Science and Technology Corresponding author

Abstract

Modeling returns on large portfolios is a challenging problem as the number of parameters in the covariance matrix grows as the square of the size of the portfolio. Traditional correlation models, for example, the dynamic conditional correlation (DCC)-GARCH model, often ignore the nonlinear dependencies in the tail of the return distribution. In this paper, we aim to develop a framework to model the nonlinear dependencies dynamically, namely the graphical copula GARCH (GC-GARCH) model. Motivated from the capital asset pricing model, to allow modeling of large portfolios, the number of parameters can be greatly reduced by introducing conditional independence among stocks given some risk factors. The joint distribution of the risk factors is factorized using a directed acyclic graph (DAG) with pair-copula construction (PCC) to enhance the modeling of the tails of the return distribution while offering the flexibility of having complex dependent structures. The DAG induces topological orders to the risk factors, which can be regarded as a list of directions of the flow of information. The conditional distributions among stock returns are also modeled using PCC. Dynamic conditional dependence structures are incorporated to allow the parameters in the copulas to be time-varying. Three-stage estimation is used to estimate parameters in the marginal distributions, the risk factor copulas, and the stock copulas. The simulation study shows that the proposed estimation procedure can estimate the parameters and the underlying DAG structure accurately. In the investment experiment of the empirical study, we demonstrate that the GC-GARCH model produces more precise conditional value-at-risk prediction and considerably higher cumulative portfolio returns than the DCC-GARCH model.

1 Introduction

Portfolio selection for large portfolios is often a challenging problem. The mean-variance optimization introduced by \textcitemarkowitz provides a framework to decide the allocation of a set of assets by balancing the risk and return. The problem is also known as the minimum variance (MV) optimization if we do not specify the target return. To take the extreme scenarios into account (for example, catastrophic losses), \textciterockafellar2000optimization proposes an efficient formulation to minimize the Conditional Value-at-Risk (CVaR) of the portfolio. These optimization problems require us to estimate the covariance matrix and the CVaR of the portfolio. However, the number of parameters in the covariance matrix grows quadratically in the size of the portfolio, and the traditional models often omit the nonlinear dependence in the tails due to the assumption of multivariate normality or multivariate $t$ , which may give poor CVaR predictions.

The capital asset pricing model (CAPM) developed by \textcitesharpe1964capital and \textcitelintner1975valuation has often been adopted in the literature, and uses a smaller number of risk factors to capture the co-movements among a large portfolio of stocks, so as to reduce the number of parameters in the estimation. Several works attempt to apply the CAPM to reduce the dimension in the covariance estimation \parencitefan2008high,fan2011high,engle2019large. Nevertheless, these works are limited to multivariate normal distribution or multivariate $t$ distribution. These multivariate distributions fail to capture the nonlinear dependencies in the tails, which are crucial for the portfolio selection when we consider the minimum CVaR problem.

To effectively capture time series properties of high-dimensional returns, we also need to take the dynamic dependence features into account. \textciteDCC_model and \textciteDCC_GARCH propose the dynamic conditional correlation (DCC)-GARCH models to model the dynamic correlations. However, these models again fail to capture the non-linear dependencies in the tails. The literature provides evidence that assets are more likely to co-move in extreme scenarios \parencitejang2002asian,jiang2017financial. It motivates us to develop a model that can reduce the dimension of parameters, and in the meantime, also capture non-linear dependence in the tails dynamically.

\textcite

sklar1959fonctions introduces the framework of the copula, which promotes flexibility to model the marginal distributions of the variables and their multivariate dependencies separately. Several works introduce copula-GARCH models \parencitejondeau2006copula,huang2009estimating,aloui2013conditional,lee2009copula,mayer2023estimation, where the marginal time series are assumed to follow GARCH structures, and the multivariate dependencies among variables are modeled using multivariate copulas. Time-varying dependence structures can also be incorporated \parencitelu2014portfolio,wang2011dynamic,ozun2007portfolio. A drawback of the multivariate copula is that some copula parameters are applied to all variables to capture the non-linear dependence structure. For example, in the multivariate $t$ -copula, the degrees of freedom is shared among all variables to explain the tail dependencies. This could still be too restrictive as the tail dependencies may vary across different pairs of variables.

Pair-copula construction (PCC) has been used to further flexibilize the copula-GARCH models, where the joint distribution is decomposed into a product of bivariate copulas based on a hierarchical structure, allowing the copula parameters to vary for different pairs of variables. Vine decomposition proposed by \textcitejoe1997multivariate is one of the methods to define the hierarchy in the PCC. Vine-copula GARCH models are used in modeling financial returns in the literature \parencitebrechmann2013risk,vine_copula,min2010bayesian,acar2019flexible. \textcitemin2010bayesian,vine_copula use the Markov chain Monte Carlo (MCMC) methods for parameter estimation in vine-copula GARCH models. Besides, \textcitebauer2016pair,pairwise_copula_constructions attempt to use Bayesian networks (BN) to specify the hierarchical structures among variables in the PCC. Bayesian networks are graphical models to represent dependence structures among variables using directed acyclic graphs (DAG), where the nodes in the DAG represent variables, and the directed edges in the DAG represent the dependence structures. The advantage of using BN is that we can specify topological orders for the variables in the BN, where topological orders give a linear ordering of the variables in the BN. The topological orders of a BN of stock returns can be regarded as the flow of information or risk in the financial market \parencitechan2023moving,garvey2015analytical, which provide a natural way to set up a hierarchy in the variables for PCC. Furthermore, PCCs using BN are more parsimonious than that using vine copula models in general due to its focus on conditional independence \parencitebauer2016pair; the number of copulas to be estimated in PCCs using BN is smaller than or equal to that using vine decomposition.

In this paper, we aim to develop a Graphical Copula GARCH (GC-GARCH) framework to model the nonlinear dependencies in the tails among a large portfolio of stocks (say, 100 stocks). The GC-GARCH model consists of the following four components:

(1)

The conditional independence of stock returns given the risk factors motivated from the Capital Asset Pricing Model (CAPM).
(2)

The use of a directed acyclic graph to define the dependence structures of the risk factors.
(3)

To specify the conditional distributions of returns using pair-copula constructions.
(4)

To impose time-varying dependence structures through copula parameters.

Having these four components, we can greatly reduce the dimensionality in the parameter estimation; in the meantime, we can explain dynamically the nonlinear dependencies among stocks from the market indexes using the copulas, where the market indexes are hierarchized naturally using a DAG.

The rest of the paper is as follows. In Section 2, we provide details of the methodology of the GC-GARCH model and explain the logic of the above four components. In Section 3, we discuss the computational issues in the GC-GARCH model. In Section 4, we derive the likelihood functions and posterior distributions, and provide algorithms for parameter estimation and DAG estimation. In Section 5, we discuss how we use the GC-GARCH model in portfolio selection. In Section 6, we provide the results of a simulation study to illustrate that the estimation algorithms in Section 4 can recover the parameters and the underlying DAG accurately. In Section 7, we present findings of the empirical study by applying the GC-GARCH model to a portfolio of 92 stocks with 10 market indexes. We also compare the predictive performance between the GC-GARCH model and the DCC-GARCH model. In Section 8, we present our conclusions and a discussion of the paper.

2 Methodology

We denote the $m$ risk factor returns and $p$ stock returns on day $t$ by $r_{1,t},\ldots,r_{m,t}$ , and $r_{m+1,t},\ldots,r_{m+p,t}$ , respectively, for $t=1,2,\ldots,T$ . To model nonlinear dependence among risk factor returns and stock returns dynamically, we are interested in developing $F^{[t]}(r_{1,t},\ldots,r_{m+p,t})$ , the joint distribution function of $r_{1,t},\ldots,r_{m,t},r_{m+1,t},\ldots,r_{m+p,t}$ given $\mathcal{F}_{t-1}$ , the information set up to time $t-1$ . We further let $D=\{D_{1},\ldots,D_{T}\}$ be the set of data containing all returns from day $1$ to day $T$ , where $D_{t}=\{r_{1,t},\ldots,r_{m+p,t}\}$ contains all returns on day $t$ , $t=1,\ldots,T$ . This conditional distribution setting is inline with the usual GARCH modeling. In this paper, we develop the graphical copula GARCH (GC-GARCH) model. There are four main features in the GC-GARCH specification: (1) the conditional independence of stock returns, $r_{m+1,t},\ldots,r_{m+p,t}$ , given the risk factors, $r_{1,t},\ldots,r_{m,t}$ motivated from the Capital Asset Pricing Model (CAPM) theory; (2) the use of a directed acyclic graph to define the dependence structures of the risk factors; (3) the specification of the conditional distributions using pair-copula construction; and (4) the imposition of time-varying dependence structures in the modeling; more specifically, we allow the correlation parameters in the conditional bivariate copulas to be time varying in a similar way to the dynamic conditional correlation (DCC)-GARCH model \parenciteDCC_model,DCC_GARCH.

2.1 Conditional independence of stock returns

The first level of construction of the GC-GARCH model is through application of the conditional independence property motivated from the CAPM theory. The joint density function of stock returns and risk factor returns can be factorized into

	$\displaystyle f^{[t]}(r_{1,t},\ldots,r_{m+p,t})$	$\displaystyle=f^{[t]}(r_{m+1,t},\ldots,r_{m+p,t}\|r_{1,t},\ldots,r_{m,t})f^{[t]% }(r_{1,t},\ldots,r_{m,t})$		(1)
		$\displaystyle=\left[\prod_{j=m+1}^{m+p}f^{[t]}(r_{j,t}\|r_{1,t},\ldots,r_{m,t})% \right]f^{[t]}(r_{1,t},\ldots,r_{m,t}),$		(1)

where $f^{[t]}(r_{1,t},\ldots,r_{m,t})$ is the conditional joint density of the risk factor returns given $\mathcal{F}_{t-1}$ . It gives a crucial property that the stock return series $r_{i,t}$ is independent of other stock return series $r_{j,t}$ for $i\neq j$ and $i,j>m$ , conditional on the risk factor returns $r_{1,t},\ldots,r_{m,t}$ . A similar exploration under multivariate $t$ distribution assumption was found in \textciterisk_factor_mapping.

2.2 Directed acyclic graph models for risk factor returns

We formulate the conditional joint distribution of the $m$ risk factor returns, i.e., the term

f^{[t]}(r_{1,t},\ldots,r_{m,t})

in (1), using directed acyclic graphs (DAGs) \parencitekoller2009probabilistic. A principle of DAG is to represent possible causal or conditional independence relationships through a graph. The variables in the DAG are called nodes, and a relationship between any two variables is indicated by a directed edge. Figure 1 presents two DAGs with four nodes (i.e., $m=4$ ) and multiple directed edges. Whenever there is an edge from node $i$ to node $j$ , we say that node $i$ is a parent of node $j$ , and we notate $i\in\text{pa}(j)$ , where $\text{pa}(j)$ is the parent set of node $j$ . For example, in 1(a), $\text{pa}(1)=\text{pa}(2)=\varnothing$ (where $\varnothing$ denotes the empty set), $\text{pa}(3)=\{1\}$ , and $\text{pa}(4)=\{1,2,3\}$ . Variables in a DAG can be ordered topologically. A topological order, or simply an order, is valid if $i\in\text{pa}(j)$ , then node $i$ must be on the left of node $j$ in the order. We notate an order using the symbol $\prec$ . For example, an order of the DAG in 1(a) is $\prec=(1,2,3,4)$ . Another order $\prec=(2,1,3,4)$ is also valid for the DAG in 1(a); since nodes $1$ and $2$ are not connected by an edge, we can either order node $1$ first or node $2$ first. With the same order $\prec$ , the DAG, however, may not be unique. For example, the order $\prec=(1,2,3,4)$ is valid for both DAGs in 1(a) and 1(b).

By using the standard Markov assumptions in the DAG \parencitepairwise_copula_constructions, the joint density of $m$ risk factor returns can be decomposed into a product of univarate conditional densities

f^{[t]}(r_{1,t},\ldots,r_{m,t})=\prod_{i=1}^{m}f^{[t]}(r_{i,t}|\mathbf{r}_{% \text{pa}(i),t}),

(2)

where $\mathbf{r}_{\text{pa}(i),t}$ is the set of risk factor returns that are parents of node $i$ . For example, we have $f^{[t]}(r_{1,t},\ldots,r_{4,t})=f^{[t]}(r_{1,t})f^{[t]}(r_{2,t})f^{[t]}(r_{3,t% }|r_{1,t})f^{[t]}(r_{4,t}|r_{1,t},r_{2,t},r_{3,t})$ for the DAG in 1(a), and $f^{[t]}(r_{1,t},\ldots,r_{4,t})=f^{[t]}(r_{1,t})f^{[t]}(r_{2,t}|r_{1,t})f^{[t]% }(r_{3,t}|r_{1,t},r_{2,t})f^{[t]}(r_{4,t}|r_{3,t})$ for the DAG in 1(b).

Instead of adopting standard multivariate distribution such as multivariate normal or $t$ distributions, the decomposition in (2) enables us to flexibly specify the joint distribution of the $m$ risk factors through a graphical representation in the DAG.

(a) A DAG with order

\prec=(1,2,3,4)

\prec=(2,1,3,4)

is another valid order.

(b) Another DAG with order

\prec=(1,2,3,4)

Figure 1: Example of two DAGs.

2.3 Copulas and conditional copulas

By Sklar’s theorem \parencitesklar1959fonctions, $F(x_{1},\ldots,x_{d})$ , a $d$ -dimensional ( $d\in\{1,2,\ldots\}$ ) cumulative distribution function (CDF) for the variables $(x_{1},\ldots,x_{d})\in(-\infty,\infty)^{d}$ can be written as \parencitecopula_book

F(x_{1},\ldots,x_{d})=C(F_{1}(x_{1}),\ldots,F_{d}(x_{d})),

for some $d$ -dimensional copula $C:[0,1]^{d}\mapsto[0,1]$ and $F_{1},\ldots,F_{d}$ are the marginal CDFs of $x_{1},\ldots,x_{d}$ respectively. The joint density for $x_{1},\ldots,x_{d}$ is

f(x_{1},\ldots,x_{d})=c(F_{1}(x_{1}),\ldots,F_{d}(x_{d}))\prod_{i=1}^{d}f_{i}(% x_{i}),

(3)

where $f_{1},\ldots,f_{d}$ are the marginal density functions of $x_{1},\ldots,x_{d}$ respectively, and the copula density function $c$ is obtained by

c(u_{1},\ldots,u_{d})=\frac{\partial^{d}C(u_{1},\ldots,u_{d})}{\partial u_{1}% \ldots\partial u_{d}},

where $u_{i}=F_{i}(x_{i})$ for $i=1,\ldots,d$ .

The conditional distribution in (2) are often modeled using multivariate Gaussian distributions \parenciteGBN, which may performs poorly in extreme scenarios in financial applications. Instead, we factorize the conditional distributions in (2) using a similar logic of vine decomposition based on pair-copula construction (PCC) \parencitevine_copula. Using (3), we can express the conditional density of variables $x$ and $y$ given a random vector $z$ by $f(x,y|z)=c_{x,y|z}(F(x|z),F(y|z))f(x|z)f(y|z)$ , where $c_{x,y|z}$ is a conditional copula density function given $z$ . Therefore, we can express $f(x|y,z)$ as

f(x|y,z)=\frac{f(x,y|z)}{f(y|z)}=\frac{c_{x,y|z}(F(x|z),F(y|z))f(x|z)f(y|z)}{f% (y|z)}=c_{x,y|z}(F(x|z),F(y|z))f(x|z).

(4)

Using (4) and taking $x$ as $r_{j,t}$ , $y$ as $r_{m,t}$ , and $z$ as $r_{1,t},\ldots,r_{m-1,t}$ for $j=m+1,\ldots,m+p$ , the conditional density of return on the $j$ -th stock can be expressed as

		$\displaystyle f^{[t]}(r_{j,t}\|r_{1,t},\ldots,r_{m,t})$
	$\displaystyle=$	$\displaystyle c_{j,m\|1,\ldots,m-1}^{[t]}(F^{[t]}(r_{j,t}\|r_{1,t}\ldots,r_{m-1,% t}),F^{[t]}(r_{m,t}\|r_{1,t},\ldots,r_{m-1,t}))f^{[t]}(r_{j,t}\|r_{1,t},\ldots,r% _{m-1,t}),$

and recursively, for $i>1$ , we have

f^{[t]}(r_{j,t}|r_{1,t},\ldots,r_{m,t})=\left[\prod_{i=1}^{m}c_{j,i|1,\ldots,i% -1}^{[t]}(F^{[t]}(r_{j,t}|r_{1,t},\ldots,r_{i-1,t}),F^{[t]}(r_{i,t}|r_{1,t},% \ldots,r_{i-1,t}))\right]\cdot f_{j}^{[t]}(r_{j,t}),

(5)

where $c_{j,i|1,\ldots,i-1}^{[t]}(\cdot,\cdot)$ is the conditional bivariate copula density between the $j$ -th stock return $r_{j,t}$ and the $i$ -th risk factor return $r_{i,t}$ given the information on all the first $i-1$ risk factor returns $r_{1,t},\ldots,r_{i-1,t}$ and $\mathcal{F}_{t-1}$ , and $f_{j}^{[t]}(r_{j,t})$ is the conditional marginal distribution of $r_{j,t}$ given $\mathcal{F}_{t-1}$ , for $i=2,\ldots,m$ and $j=m+1,\ldots,m+p$ . For $i=1$ , we have

c_{j,i|1,\ldots,i-1}^{[t]}(F^{[t]}(r_{j,t}|r_{1,t},\ldots,r_{i-1,t}),F^{[t]}(r% _{i,t}|r_{1,t},\ldots,r_{i-1,t}))\equiv c_{j,1}(F^{[t]}(r_{j,t}),F^{[t]}(r_{1,% t})).

In the same way, to express the conditional distributions in the DAG for the risk factor returns in (2) into a product of conditional bivariate copula density, we label the $m$ risk factor returns to give the order $\prec=(1,2,\ldots,m)$ . Then, denote the parent set of $r_{i,t}$ as $\text{pa}(i)=\{i[1],\ldots,i[n(i)]\}$ , where $i[1]<\ldots<i[n(i)]$ , and $n(i)$ is the number of parents of $r_{i,t}$ . With the above setting, we can express the conditional density of the return of the $i$ -th risk factor to

f^{[t]}(r_{i,t}|\mathbf{r}_{\text{pa}(i),t})=\left[\prod_{k=1}^{n(i)}c_{i,i[k]% |i[1],\ldots,i[k-1]}^{[t]}(F^{[t]}(r_{i,t}|r_{i[1],t},\ldots,r_{i[k-1],t}),F^{% [t]}(r_{i[k],t}|r_{i[1],t},\ldots,r_{i[k-1],t}))\right]\cdot f_{i}^{[t]}(r_{i,% t}).

(6)

For example, for the DAG in 1(a), we have the parent set $\text{pa}(4)=\{1,2,3\}$ , the number of parents $n(4)=3$ , $4[1]=1$ , $4[2]=2$ and $4[3]=3$ . Therefore, (6) implies

	$\displaystyle f^{[t]}(r_{4,t}\|\mathbf{r}_{\text{pa}(4),t})$	$\displaystyle=c_{4,3\|1,2}^{[t]}(F^{[t]}(r_{4,t}\|r_{1,t},r_{2,t}),F^{[t]}(r_{3,% t}\|r_{1,t},r_{2,t}))c_{4,2\|1}^{[t]}(F^{[t]}(r_{4,t}\|r_{1,t}),F^{[t]}(r_{2,t}\|r% _{1,t}))$
		$\displaystyle c_{4,1}^{[t]}(F^{[t]}(r_{4,t}),F^{[t]}(r_{1,t}))\cdot f^{[t]}_{4% }(r_{4,t}).$

Similarly, we have $\text{pa}(3)=\{1\}$ , $n(3)=1$ and $3[1]=1$ . Applying (6) again implies

f^{[t]}(r_{3,t}|\mathbf{r}_{\text{pa}(3),t})=c_{3,1}^{[t]}(F^{[t]}(r_{3,t}),F^% {[t]}(r_{1,t}))\cdot f^{[t]}_{3}(r_{3,t}).

Note that $\text{pa}(r_{1,t})=\text{pa}(r_{2,t})=\varnothing$ and thus $f^{[t]}(r_{2,t}|\mathbf{r}_{\text{pa}(2),t})=f^{[t]}_{2}(r_{2,t})$ and $f^{[t]}(r_{1,t}|\mathbf{r}_{\text{pa}(1),t})=f^{[t]}_{1}(r_{1,t})$ . Combining these factorization results and (6), we can express the joint density of $m$ risk factor returns as

	$\displaystyle f^{[t]}(r_{1,t},r_{2,t},r_{3,t},r_{4,t})$	$\displaystyle=\Bigg{[}c_{4,3\|1,2}^{[t]}(F^{[t]}(r_{4,t}\|r_{1,t},r_{2,t}),F^{[t% ]}(r_{3,t}\|r_{1,t},r_{2,t}))c_{4,2\|1}^{[t]}(F^{[t]}(r_{4,t}\|r_{1,t}),F^{[t]}(r% _{2,t}\|r_{1,t}))$
		$\displaystyle c_{4,1}^{[t]}(F^{[t]}(r_{4,t}),F^{[t]}(r_{1,t}))\Bigg{]}\cdot c_% {3,1}^{[t]}(F^{[t]}(r_{3,t}),F^{[t]}(r_{1,t}))\cdot\prod_{i=1}^{4}f^{[t]}_{i}(% r_{i,t}).$

In general, combining the DAG decomposition of the $m$ risk factors in (2), the copula factorization of the stock returns in (5), and the copula factorization for the risk factor returns in (6), the joint density of all returns $r_{1,t}\ldots r_{m+p,t}$ in (1) can be decomposed into a product of conditional bivariate copulas and marginal densities:

$\displaystyle f^{[t]}(r_{1,t},\ldots,r_{m+p,t})$	$\displaystyle=\left[\prod_{j=m+1}^{m+p}\prod_{i=1}^{m}c_{j,i\|1,\ldots,i-1}^{[t% ]}(F^{[t]}(r_{j,t}\|r_{1,t},\ldots,r_{i-1,t}),F^{[t]}(r_{i,t}\|r_{1,t},\ldots,r_% {i-1,t}))\right]\cdot$	(7)
	$\displaystyle\left[\prod_{i=1}^{m}\prod_{k=1}^{n(i)}c_{i,i[k]\|i[1],\ldots,i[k-% 1]}^{[t]}(F^{[t]}(r_{i,t}\|r_{i[1],t},\ldots,r_{i[k-1],t}),F^{[t]}(r_{i[k],t}\|r% _{i[1],t},\ldots,r_{i[k-1],t}))\right]\cdot$
	$\displaystyle\prod_{j=1}^{m+p}f_{j}^{[t]}(r_{j,t}).$

In financial econometrics studies, the number of stocks, $p$ , is much larger than the number of risk factors, $m$ . Therefore, from the joint density in (7), the number of bivariate copulas is at most $mp+m(m-1)/2$ , equality holds if all pairs of nodes in the DAG are connected, in which case the number of edges in the DAG is $\binom{m}{2}=m(m-1)/2$ . The number of parameters is of $O(p)$ instead of $O(p^{2})$ in usual dynamic covariance modeling of $p$ stocks.

2.4 Dynamic conditional dependence and tail dependence

The decomposition in (7) also provides flexibility in the choice of conditional dependence, including linear correlation and the tail dependence parameters through conditional bivariate copulas. To incorporate time-varying dependence features in the GC-GARCH model, we follow the idea in \textcitevine_copula, which inspired by the DCC-GARCH models of \textciteDCC_model and \textciteDCC_GARCH, to allow the correlation parameters in the conditional copulas $c^{[t]}_{j,i|1,\ldots,i-1}$ and $c^{[t]}_{i,i[k]|i[1],\ldots,i[k-1]}$ to be time-varying. In this paper, we focus on the use of bivariate $t$ -copulas with correlation parameters $\varphi_{j,i|1,\ldots,i-1}^{[t]}$ and $\varphi^{[t]}_{i,i[k]|i[1],\ldots,i[k-1]}$ . In general, we denote the conditional copula and copula density at time $t$ for two arbitrary variables $r_{x,t}$ and $r_{y,t}$ conditional to the set $z$ as $C_{xy|z}^{[t]}$ and $c_{xy|z}^{[t]}$ respectively. In the simulation study and empirical study, we use the $t$ -copula to illustrate the GC-GARCH model. The conditional $t$ -copula at time $t$ with correlation parameter $\varphi=\varphi_{xy|z}^{[t]}$ and degrees of freedom $v=v_{xy|z}$ is given by \parencitedemarta2005t

\displaystyle C^{[t]}_{xy|z}(u_{x|z,t},u_{y|z,t})

\displaystyle=\int_{-\infty}^{t_{v}^{-1}(u_{x|z,t})}\int_{-\infty}^{t_{v}^{-1}% (u_{y|z,t})}f_{2;v,\varphi}(x,y)dydx,

where $u_{x|z,t}=F^{[t]}(r_{x,t}|r_{z,t})$ , $u_{y|z,t}=F^{[t]}(r_{y,t}|r_{z,t})$ , $t_{v}^{-1}(\cdot)$ is the inverse CDF of the univariate $t$ distribution with degrees of freedom $v$ , $f_{2;v,\varphi}(\cdot,\cdot)$ is the density function of the bivariate $t$ distribution with degrees of freedom $v$ and correlation parameter $-1<\varphi<1$ , defined as

f_{2;v,\varphi}(x,y)=\frac{1}{2\pi\sqrt{1-{\varphi}^{2}}}\left(1+\frac{x^{2}+y% ^{2}-2\varphi xy}{v(1-\varphi^{2})}\right)^{-\frac{v+2}{2}}.

We impose the finite covariance condition that $v>2$ \parencitedemarta2005t. Using (3), the conditional $t$ -copula density at time $t$ with correlation parameter $\varphi=\varphi_{xy|z}^{[t]}$ and degrees of freedom $v=v_{xy|z}$ is given by

\displaystyle c_{xy|z}^{[t]}(u_{x|z,t},u_{y|z,t})

\displaystyle=\frac{f_{2;v,\varphi}(t^{-1}_{v}(u_{x|z,t}),t^{-1}_{v}(u_{y|z,t}% ))}{f_{1;v}(t^{-1}_{v}(u_{x|z,t}))f_{1;v}(t^{-1}_{v}(u_{y|z,t}))},

where $f_{1;v}(\cdot)$ is the density function of the univariate $t$ -distribution with degrees of freedom $v$ , where

f_{1;v}(x)=\frac{\Gamma\left(\frac{v+1}{2}\right)}{\sqrt{v\pi}\Gamma\left(% \frac{v}{2}\right)}\left(1+\frac{x^{2}}{v}\right)^{-\frac{v+1}{2}}.

We assume that the conditional correlations change dynamically similar to the dynamics in \textcitevine_copula, inspired by the DCC-GARCH models of \textciteDCC_GARCH and \textciteDCC_model. The dynamic conditional correlation $\varphi^{[t]}_{xy|z}$ between returns $r_{x,t}$ and $r_{y,t}$ given $r_{z,t}$ and $\mathcal{F}_{t-1}$ is given by

\varphi^{[t]}_{xy|z}=(1-a_{xy|z}-b_{xy|z})\bar{\varphi}_{xy|z}+a_{xy|z}\xi_{xy% |z,t-1}+b_{xy|z}\varphi^{[t-1]}_{xy|z},

(8)

where $\bar{\varphi}_{xy|z}$ is the long-run correlation and $\xi_{xy|z,t-1}$ is the sample correlation at time $t-1$ . The stationary conditions are $0\leq a_{xy|z},b_{xy|z}<1$ , $a_{xy|z}+b_{xy|z}<1$ and $-1<\bar{\varphi}_{xy|z}<1$ . The sample correlation is obtained given the past $m_{sc}$ -period of data and $r_{z,t}$ , and is defined as

\xi_{xy|z,t-1}=\frac{\sum_{i=1}^{m_{sc}}\tilde{r}_{x|z,t-i}\tilde{r}_{y|z,t-i}% }{\sqrt{\sum_{i=1}^{m_{sc}}\tilde{r}_{x|z,t-i}^{2}\sum_{j=1}^{m_{sc}}\tilde{r}% _{y|z,t-j}^{2}}},

where $\tilde{r}_{x|z,t}=t_{v_{xy|z}}^{-1}(F^{[t]}(r_{x,t}|r_{z,t}))$ and $\tilde{r}_{y|z,t}=t_{v_{xy|z}}^{-1}(F^{[t]}(r_{y,t}|r_{z,t}))$ . We pick ${m_{sc}}=2$ in this paper for the simulation study and empirical study.

In financial application, we often want to estimate the unconditional correlations. In this case, computations are needed to convert the conditional correlation $\varphi_{xy|z}^{[t]}$ in the GC-GARCH model into unconditional correlation as in \textcitevine_copula. The unconditional correlation can be obtained iteratively using the formula in \textciterummel1976understanding,

\varphi_{xy|z_{-j}}^{[t]}=\varphi_{x,y|z}^{[t]}\sqrt{(1-{\varphi_{xz_{j}|z_{-j% }}^{[t]}}^{2})(1-{\varphi_{yz_{j}|z_{-j}}^{[t]}}^{2})}+\varphi_{xz_{j}|z_{-j}}% ^{[t]}\varphi_{yz_{j}|z_{-j}}^{[t]},

where $z_{j}$ is the $j$ th component of the vector $z$ , and $z_{-j}$ is the vector obtained by removing the $j$ th component in $z$ .

The $t$ -copula has often been used in financial return data modeling due to its ability to capture the phenomenon of dependent extreme values \parencitedemarta2005t,t_copula_rm_eg1,vine_copula. The tail dependence coefficient of the conditional $t$ -copula $C_{xy|z}^{[t]}$ with correlation parameter $\varphi_{xy|z}^{[t]}$ and degrees of freedom $v_{xy|z}$ is defined as

\lambda_{xy|z}^{[t]}=2t_{v+1}\left(-\sqrt{\frac{(v_{xy|z}+1)(1-\varphi_{xy|z}^% {[t]})}{(1+\varphi_{xy|z}^{[t]})}}\right),

where $t_{d}(\cdot)$ is the CDF of the standard $t$ -distribution with degrees of freedom $d>0$ . The tail dependence depends on two parameters, the correlation and the degrees of freedom. Note that, even when the correlation is zero, the $t$ -copula still gives dependence in the tails. The $t$ -copula captures the dependence of extreme values which is often observed in financial modeling. We report the estimate of the degrees of freedom $v_{xy|z}$ for each copula as a measure of tail dependence.

3 Computational issues

3.1 Marginal distributions

Similar to the approach in \textcitevine_copula, we assume that the each of the returns $r_{1,t},\ldots,r_{m+p,t}$ follows a GARCH(1,1)- $t$ model, i.e., a GARCH(1,1) model with $t$ distributed innovations,

	$\displaystyle r_{i,t}$	$\displaystyle=\sigma_{i,t}\varepsilon_{i,t},$		(9)
	$\displaystyle\sigma_{i,t}^{2}$	$\displaystyle=\omega_{i}+\alpha_{i}r_{i,t-1}^{2}+\beta_{i}\sigma_{i,t-1}^{2},$		(9)

where $\varepsilon_{i,t}$ ’s are independently and identically distributed standardized $t$ innovations with degrees of freedom $v_{i}$ , $i=1,\ldots,m+p$ . The constraints for positive variance and covariance stationary are $\omega_{i}>0$ , $\alpha_{i}\geq 0$ , $\beta_{i}\geq 0$ , and $\alpha_{i}+\beta_{i}<1$ for $i=1,\ldots,m+p$ . The conditional density function of $r_{i,t}$ under GARCH(1,1)- $t$ in (9) is given by

f_{i}^{[t]}(r_{i,t})=\frac{\Gamma\left(\frac{v_{i}+1}{2}\right)}{\sigma_{i,t}% \sqrt{(v_{i}-2)\pi}\Gamma\left(\frac{v_{i}}{2}\right)}\left(1+\frac{r_{i,t}^{2% }}{\sigma_{i,t}^{2}(v_{i}-2)}\right)^{-\frac{v_{i}+1}{2}},

(10)

for $t=1,\ldots,T$ . To ensure the variance is finite and (10) is well-defined, we require $v_{i}>2$ . (10) is derived using the transformation

r_{i,t}=\frac{\sigma_{i,t}}{\sqrt{\frac{v_{i}}{v_{i}-2}}}T_{v_{i}},

where $T_{v_{i}}$ follows a (unstandardized) $t$ -distribution with degrees of freedom $v_{i}$ .

3.2 Conditional distributions

Similar to \textcitevine_copula, the conditional distribution functions of the form $F^{[t]}(x_{t}|v_{t})$ in (7) have to be computed recursively using the $h$ -function associated with the conditional copulas $C_{x,v_{j}|v_{-j}}^{[t]}$ \parenciteJoe1996FamiliesO:

F^{[t]}(x_{t}|v_{t})=\frac{\partial C_{x,v_{j}|v_{-j}}^{[t]}(F^{[t]}(x_{t}|v_{% -j,t}),F^{[t]}(v_{j,t}|v_{-j,t}))}{\partial F^{[t]}(v_{j,t}|v_{-j,t})}:=h_{x,v% _{j}|v_{-j}}^{[t]}(F^{[t]}(x_{t}|v_{-j,t}),F^{[t]}(v_{j,t}|v_{-j,t})),

(11)

where $v_{t}$ is a random vector at time $t$ , $v_{j,t}$ is the $j$ th element of $v_{t}$ , and $v_{-j,t}$ is the vector obtained by removing $v_{j,t}$ from $v_{t}$ . Conversely, the inverse of the $h$ -function at time $t$ in (11) is also useful if we want to reduce the dimension of the conditioning set of $F^{[t]}(x_{t}|v_{t})$ . If both $F^{[t]}(x_{t}|v_{t})$ and $F^{[t]}(v_{j,t}|v_{-j,t})$ are available, then $F^{[t]}(x_{t}|v_{-j,t})$ can be evaluated using the inverse $h$ -function

F^{[t]}(x_{t}|v_{-j,t})={h^{[t]}}^{-1}_{x,v_{j}|v_{-j}}(F^{[t]}(x_{t}|v_{t}),F% ^{[t]}(v_{j,t}|v_{-j,t})).

(12)

In the simulation study and empirical study, we use the $t$ -copulas in the GC-GARCH model. The $h$ -function and inverse $h$ -function of the $t$ -copula with correlation parameter $\varphi=\varphi^{[t]}_{xy|z}$ and degrees of freedom $v=v_{xy|z}$ are respectively given by

h_{xy|z}^{[t]}(u_{xt},u_{yt})=t_{v+1}\left(\frac{t_{v}^{-1}(u_{xt})-\varphi t_% {v}^{-1}(u_{yt})}{\sqrt{\frac{({v}+[t_{v}^{-1}(u_{yt})]^{2})(1-\varphi^{2})}{{% v}+1}}}\right)

and

h^{{[t]}^{-1}}_{xy|z}(u_{xt},u_{yt})=t_{v}\left(t_{{v}+1}^{-1}(u_{xt})\sqrt{% \frac{({v}+[t_{v}^{-1}(u_{yt})]^{2})(1-\varphi^{2})}{{v}+1}}+\varphi t_{v}^{-1% }(u_{yt})\right),

where $u_{xt}=F_{x}^{[t]}(r_{xt})$ and $u_{yt}=F_{y}^{[t]}(r_{yt})$ are the cumulative distributions of two return variables $r_{xt}$ and $r_{yt}$ at time $t$ .

3.3 Computation of the likelihood function

The joint density in (7) contains three parts, which have to be computed sequentially: the marginal distributions, $f_{i}^{[t]}$ , the conditional copulas among the risk factors, $c_{i,i[k]|i[1],...,i[k-1]}^{[t]}$ , and the conditional copulas between the risk factors and stocks, $c_{j,i|1,...,i-1}^{[t]}$ . We define the following parameter sets:

1.

$\Theta_{1i}=\{\omega_{i},\alpha_{i},\beta_{i},v_{i}\}$ be the set of parameters in the marginal distribution of the $i$ th stock, where $i=1,\ldots,m+p$ . We also denote $\Theta_{1}=\{\Theta_{1i}\}_{i=1}^{m+p}$ .
2.

$\Theta_{2}=\{\bar{\varphi}_{i,i[k]|z_{ik}},a_{i,i[k]|z_{ik}},b_{i,i[k]|z_{ik}}% ,v_{i,i[k]|z_{ik}}:k=1,\ldots,n(i)\text{ and }i=2,\ldots,m\}$ be the set of parameters in the DAG copulas, where $z_{ik}=\{i[1],\ldots,i[k-1]\}$ , and
3.

$\Theta_{3j}=\{\bar{\varphi}_{j,i|1,\ldots,i-1},a_{j,i|1,\ldots,i-1},b_{j,i|1,% \ldots,i-1},v_{j,i|1,\ldots,i-1}:i=1,\ldots,m\}$ be the set of parameters in the copulas of the stock $j$ , where $j=m+1,\ldots,m+p$ . We also denote $\Theta_{3}=\{\Theta_{3j}\}_{j=m+1}^{m+p}$ .

The steps to compute of the joint density are depicted below. Starting at $t=1$ ,

1.

(Marginal distribution) We first compute the cumulative distribution of the marginal distributions $F^{[t]}(r_{1,t}),\ldots,F^{[t]}(r_{m+p,t})$ and the log density functions $\log f^{[t]}_{1}(r_{1,t})$ , $\ldots,\log f^{[t]}_{m+p}(r_{m+p,t})$ .

(Conditional copulas of the DAG) For each $i=2,\ldots,m$ and if $r_{i,t}$ has at least one parent (i.e., $n(i)\geq 1$ ), then for $k=1,\ldots,n(i)$ ,

(i)

if $t=1$ , initialize $\varphi_{i,i[k]|i[1],\ldots,i[k-1]}^{[1]}$ as the long-run correlation $\bar{\varphi}_{i,i[k]|i[1],\ldots,i[k-1]}$ . Otherwise, $\varphi_{i,i[k]|i[1],\ldots,i[k-1]}^{[t]}$ is computed according to the dynamics in (8).

(ii)

Compute and store the log copula density with correlation parameter $\varphi_{i,i[k]|i[1],\ldots,i[k-1]}^{[t]}$ and degrees of freedom $v_{i,i[k]|i[1],\ldots,i[k-1]}$ :

\log c_{i,i[k]|i[1],\ldots,i[k-1]}^{[t]}(F^{[t]}(r_{i,t}|r_{i[1],t},\ldots,r_{% i[k-1],t}),F^{[t]}(r_{i[k],t}|r_{i[1],t},\ldots,r_{i[k-1],t})).

(iii)

If $k<n(i)$ , compute and store $F^{[t]}(r_{i,t}|r_{i[1],t},\ldots,r_{i[k],t})$ and $F^{[t]}(r_{i[k+1],t}|r_{i[1],t},\ldots,r_{i[k],t})$ using the $h$ -function in (11).

(Conditional copulas between the stocks and the risk factors) For each stock $j=m+1,\ldots,m+p$ , we compute as follows: For $i=1,\ldots,m$ ,

(i)

if $t=1$ , initialize $\varphi_{j,i|1,\ldots,i-1}^{[1]}$ as the long-run correlation $\bar{\varphi}_{j,i|1,\ldots,i-1}$ . Otherwise, $\varphi_{j,i|1,\ldots,i-1}^{[t]}$ is computed according to the dynamics in (8).

(ii)

Compute and store the log copula density with correlation parameter $\varphi_{j,i|1,\ldots,i-1}^{[t]}$ and degrees of freedom $v_{j,i|1,\ldots,i-1}$ :

\log c_{j,i|1,\ldots,i-1}^{[t]}(F^{[t]}(r_{j,t}|r_{1,t},\ldots,r_{i-1,t}),F^{[% t]}(r_{i,t}|r_{1,t},\ldots,r_{i-1,t})).

(iii)

If $i<m$ , compute and store $F^{[t]}(r_{j,t}|r_{1,t},\ldots,r_{i,t})$ and $F^{[t]}(r_{i+1,t}|r_{1,t},\ldots,r_{i,t})$ using the $h$ -function in (11).

4.

Set $t$ to $t+1$ and go back to step 1, until $t=T$ .

Finally, we can calculate the log likelihood function

		$\displaystyle\ell(\Theta\|D)=\sum_{t=1}^{T}\log f^{[t]}(r_{1,t},\ldots,r_{m+p,t})$
		$\displaystyle=\left(\sum_{i=1}^{m+p}\ell_{1i}(\Theta_{1i}\|D)\right)+\ell_{2}(% \Theta_{1},\Theta_{2}\|D)+\left(\sum_{j=m+1}^{m+p}\ell_{3j}(\Theta_{1},\Theta_{% 2},\Theta_{3j}\|D)\right),$

with

$\displaystyle\ell_{1i}(\Theta_{1i}\|D)$	$\displaystyle=\sum_{t=1}^{T}\log f_{i}^{[t]}(r_{i,t}),\text{for }i=1,\ldots,m+p,$	(13)
$\displaystyle\ell_{2}(\Theta_{1},\Theta_{2}\|D)$	$\displaystyle=\sum_{i=2}^{m}\sum_{k=1}^{n(i)}\sum_{t=1}^{T}\log c_{i,i[k]\|i[1]% ,\ldots,i[k-1]}^{[t]}(F^{[t]}(r_{i,t}\|r_{i[1],t},\ldots,r_{i[k-1],t}),F^{[t]}(% r_{i[k],t}\|r_{i[1],t},\ldots,r_{i[k-1],t})),$	(14)
$\displaystyle\ell_{3j}(\Theta_{1},\Theta_{2},\Theta_{3j}\|D)$	$\displaystyle=\sum_{i=1}^{m}\sum_{t=1}^{T}\log c_{j,i\|1,\ldots,i-1}^{[t]}(F^{[% t]}(r_{j,t}\|r_{1,t},\ldots,r_{i-1,t}),F^{[t]}(r_{i,t}\|r_{1,t},\ldots,r_{i-1,t}% )),$	(15)
	$\displaystyle\text{for }j=m+1,\ldots,m+p,$

and $\Theta=\{\Theta_{1},\Theta_{2},\Theta_{3}\}$ . In the estimation of parameters in DAG copulas, we need to compute starting from $i=2$ to $m$ sequentially since the conditional distributions depend on the parents. However, the parameters in each of the marginal distribution can be estimated separately. For the stock copulas, $r_{j,t}$ only depend on the risk factor returns $r_{1,t},\ldots,r_{m,t}$ but not other stock returns $r_{j^{\prime},t}$ , for all $j^{\prime}\in\{m+1,\ldots,m+p\}\setminus\{j\}$ . This implies that we can estimate the parameters in the stock copulas separately for each stock. Thus, we express the log likelihood function as a sum of marginal likelihood functions and a sum of stock copula likelihood functions because the parameters can be estimated separately for each term in the two summations. We discuss the detailed procedures for the estimation in Section 4.

3.4 The reduced DAG space

The conditional density of the $i$ th risk factor return, $i=1,\ldots,m$ , in (6) is dependent on the conditional distribution functions $F^{[t]}(r_{i,t}|r_{i[1],t},\ldots,r_{i[k-1],t})$ and $F^{[t]}(r_{i[k],t}|r_{i[1],t},\ldots,r_{i[k-1],t})$ , which are computed recursively using the $h$ -function in (11). However, the computation of the second conditional distribution function $F^{[t]}(r_{i[k],t}|r_{i[1],t},\ldots,r_{i[k-1],t})$ may be complicated for some DAGs. To illustrate this, Figure 2 contains a DAG with five variables as an example. Using the DAG factorization in (2), the joint density function of $r_{1,t},r_{2,t},r_{3,t},r_{4,t},r_{5,t}$ can be expressed as

f^{[t]}(r_{1,t},r_{2,t},r_{3,t},r_{4,t},r_{5,t})=f^{[t]}(r_{1,t})f^{[t]}(r_{2,% t})f^{[t]}(r_{3,t})f^{[t]}(r_{4,t}|r_{2,t},r_{3,t})f^{[t]}(r_{5,t}|r_{1,t},r_{% 3,t},r_{4,t}).

(16)

Choosing the topological order $\prec=(1,2,3,4,5)$ (as there are multiple orders that are compatible to the DAG in Figure 2), we further consider the pair-copula construction for the term $f^{[t]}(r_{5,t}|r_{1,t},r_{3,t},r_{4,t})$ in (16) using (6) :

	$\displaystyle f^{[t]}(r_{5,t}\|r_{1,t},r_{3,t},r_{4,t})=$	$\displaystyle c_{5,1}^{[t]}(F^{[t]}(r_{5,t}),F^{[t]}(r_{1,t}))c_{5,3\|1}^{[t]}(% F^{[t]}(r_{5,t}\|r_{1,t}),F^{[t]}(r_{3,t}\|r_{1,t}))\cdot$		(17)
		$\displaystyle c_{5,4\|3,1}^{[t]}(F^{[t]}(r_{5,t}\|r_{3,t},r_{1,t}),F^{[t]}(r_{4,% t}\|r_{3,t},r_{1,t})).$		(17)

On the other hand, the conditional density for $r_{4,t}$ in (16) is factorized as

f^{[t]}(r_{4t}|r_{2t},r_{3t})=c_{4,2}^{[t]}(F^{[t]}(r_{4t}),F^{[t]}(r_{2t}))c_% {4,3|2}^{[t]}(F^{[t]}(r_{4t}|r_{2t}),F^{[t]}(r_{3t}|r_{2t})),

(18)

which depends only on two copulas $c_{4,2}(\cdot,\cdot)$ and $c_{4,3|2}(\cdot,\cdot)$ . When we compute the conditional distribution function $F^{[t]}(r_{4t}|r_{3t},r_{1t})$ in the last term in (17), we apply the $h$ function in (11):

\displaystyle F^{[t]}(r_{4t}|r_{3t},r_{1t})=F^{[t]}(r_{4t}|r_{3t})=h_{4,3}^{[t% ]}(F^{[t]}(r_{4t}),F^{[t]}(r_{3t})),

(19)

where the first equality in (19) is due to the fact that $r_{4t}$ and $r_{1t}$ are independent in the Bayesian network in Figure 2 ( $r_{4t}$ and $r_{1t}$ are independent when $r_{5t}$ is not given). The second equality in (19) requires the $h$ function $h_{4,3}^{[t]}(u,v)=\partial C_{4,3}^{[t]}(u,v)/\partial v$ , where $u=F^{[t]}(r_{4t})$ and $v=F^{[t]}(r_{3t})$ , but the copula $C_{4,3}^{[t]}(\cdot,\cdot)$ is not defined in (18) under the DAG in Figure 2. Therefore, we cannot use the $h$ function to calculate $F^{[t]}(r_{4t}|r_{3t})$ . Alternatively, we can compute this conditional distribution function using the formula

	$\displaystyle F^{[t]}(r_{4t}\|r_{3t})$	$\displaystyle=\int_{-\infty}^{\infty}F^{[t]}(r_{4t}\|r_{3t},r_{2t})f^{[t]}(r_{2% t})dr_{2t}$		(20)
		$\displaystyle\equiv\int_{-\infty}^{\infty}h_{4,3\|2}^{[t]}(F^{[t]}(r_{4t}\|r_{2t% }),F^{[t]}(r_{3t}\|r_{2t}))f^{[t]}(r_{2t})dr_{2t},$		(20)

where the $h$ function $h_{4,3|2}^{[t]}(\cdot,\cdot)$ is directly computable from the copula $C_{4,3|2}^{[t]}(\cdot,\cdot)$ . However, the integration in (20) may not be easily computable and the use of numerical integration is required. To conclude, (19) is not easily computable because the condition in the conditional probability in (19) does not include the first parent of $r_{4t}$ , i.e., $r_{2t}$ . We need to marginalize $r_{2t}$ from the copula as in (20) to compute the conditional distribution $F^{[t]}(r_{4t}|r_{3t})$ .

To reduce computation burden, we consider the DAGs in the reduced space, such that each of the conditional probabilities in the equation

\displaystyle c_{i,i[k]|i[1],\ldots,i[k-1]}^{[t]}(F^{[t]}(r_{i,t}|r_{i[1],t},% \ldots,r_{i[k-1],t}),F^{[t]}(r_{i[k],t}|r_{i[1],t},\ldots,r_{i[k-1],t}))

(21)

contains a “cummulative” parent set. Appendix A contains an algorithm to subset the DAG space to reduce computation burden.

Figure 2: An example to illustrate the difficulty of calculating the conditional distributions.

3.5 Generating returns from the GC-GARCH model

In the simulation study, we simulate a hypothetical data set to evaluate the performance of the estimation procedures. To simulate $(m+p)$ -dimensional time series from the joint density function in (7), we use a method similar to the one in \textcitevine_copula. For each variable in the DAG, we generate a value from $U(0,1)$ , the uniform distribution over the interval $(0,1)$ , to represent a realization of the conditional CDF of $r_{i,t}$ given its parent set, i.e., $u_{i|\mathbf{r}_{\text{pa}(i)},t}=F^{[t]}(r_{i,t}|\mathbf{r}_{\text{pa}(i)})$ , for $i=1,\ldots,m$ . For each stock variable, we similarly generate a value from $U(0,1)$ to represent a realization of the conditional CDF of $r_{j,t}$ given the risk factor returns $r_{1,t}\ldots,r_{i-1,t}$ , notated by $u_{j|1,\ldots,i-1,t}=F^{[t]}(r_{j,t}|r_{1,t},\ldots,r_{i-1,t})$ . We then obtain $F^{[t]}(r_{1,t}),\ldots,F^{[t]}(r_{m+p,t})$ by applying the inverse $h$ -function in (12) to reduce the conditioning set by 1 each time until it is empty.

For example, suppose that we want to simulate returns from a GC-GARCH model with $m=4$ risk factors with the DAG structure in 1(a) and $p$ stocks. For $t=1,\ldots,T$ , we first generate $u_{1,t},u_{2,t},u_{3|1,t},u_{4|123,t}$ independently from $\text{Uniform}[0,1]$ for the risk factor returns. We then determine the unconditional cumulative distribution functions of the risk factor returns $F^{[t]}(r_{1,t}),\ldots,F^{[t]}(r_{4,t})$ :

	$\displaystyle F^{[t]}(r_{1,t})$	$\displaystyle=u_{1,t},$
	$\displaystyle F^{[t]}(r_{2,t})$	$\displaystyle=u_{2,t},$
	$\displaystyle F^{[t]}(r_{3,t})$	$\displaystyle=h_{3,1}^{{[t]}^{-1}}(u_{3\|1,t},u_{1,t}),$
	$\displaystyle F^{[t]}(r_{4,t})$	$\displaystyle=h_{4,1}^{{[t]}^{-1}}(h_{4,2\|1}^{{[t]}^{-1}}(h_{4,3\|1,2}^{{[t]}^{% -1}}(u_{4\|1,2,3,t},u_{3\|1,t}),u_{2,t}),u_{1,t}).$

For $j=m+1,\ldots,p$ , we generate $u_{j|1,2,3,4}$ from $\text{Uniform}[0,1]$ , which represents a realization of $F^{[t]}(r_{j,t}|r_{1,t},\ldots,r_{4,t})$ , the conditional distribution of the stock returns $r_{j,t}$ given the risk factor returns. Then, we apply the inverse $h$ -function in (12) recursively to obtain the unconditional cumulative distribution functions of the stock return $r_{j,t}$ :

F^{[t]}(r_{j,t})=h_{j,1}^{{[t]}^{-1}}h_{j,2}^{{[t]}^{-1}}(h_{j,3|1,2}^{{[t]}^{% -1}}(h_{j,4|1,2,3}^{{[t]}^{-1}}(u_{j|1,2,3,4},u_{4|1,2,3}),u_{3|1,t}),u_{2,t})% ,u_{1,t}).

Finally, as we assume the returns follow GARCH(1,1) models with $t$ distributed innovations as described in (9), by applying the inverse of the $t$ distribution for each $F^{[t]}(r_{1,t}),\ldots,F^{[t]}(r_{m+p,t})$ , we obtain a vector of simulated returns $r_{1t},\ldots,r_{m+p,t}$ .

4 Estimation

4.1 Parameter estimation

Recall that we divide the parameters in the GC-GARCH model into three groups in Section 3.3. In the literature, the estimation of the pair copulas is conducted individually for each copula \parencitebrechmann2013risk, referred to as sequential estimation. The sequential estimates can be obtained quickly due to the reduction in dimension. The sequential estimates are shown to be consistent and often perform well \parenciteHOBAEKHAFF201291,Haff2013ParameterEF.

For $\Theta_{1i}$ ’s, the parameters in the marginal distributions are first estimated using a numerical solver in R. This is the method of inference functions for margins (IFM) introduced by \textcitejoe1997multivariate,joe2005asymptotic. For $\Theta_{2}$ , the parameters in the DAG copulas, we use the sequential estimates of the parameters as starting values for a full estimation procedure using Markov chain Monte Carlo (MCMC) sampling \parenciteHaff2013ParameterEF. For $\Theta_{3j}$ ’s, the parameters in the stock copulas, however, we do not conduct MCMC sampling as the number of stocks in the portfolio, $p$ , is large in the simulation study and empirical study, therefore conducting MCMC sampling for each of the stocks is not feasible in practice.

Similar to the IFM, we first estimate $\Theta_{2}$ , the parameters in the DAG copulas. Then, we insert the estimates of $\Theta_{2}$ back into the likelihood function and estimate $\Theta_{3j}$ ’s, the parameters of the stock copulas. The rationale behind this is that the risk factors are used to explain the co-movements of the stock returns. Thus, the parameters in the DAG copulas should be estimated alone without including the stocks. For example, the parameters of the copulas in the DAG should be the same no matter whether we include $p=50$ or $p=100$ stocks in the GC-GARCH model.

We have already derived the likelihood functions in Section 3.3. We further derive the posterior distributions of the DAG copulas for MCMC sampling. The procedures of the parameter estimation in the GC-GARCH model are as follows.

1.

Each return series follows a GARCH(1,1) models with $t$ distributed innovations in (9). For $i=1,\ldots,m+p$ , we obtain the maximum likelihood estimates $\hat{\Theta}_{1i}$ from the log likelihood function of the $i$ -th return series in (13) using a numerical optimizer in R, with positive variances and stationary covariance constraints in (9).

Using the sequential estimates of $\Theta_{2}$ as starting values, we conduct MCMC sampling for the parameters $\Theta_{2}$ , whose posterior distribution is

		$\displaystyle\log f(\Theta_{2}\|\hat{\Theta}_{1},D)$		(22)
	$\displaystyle=$	$\displaystyle\ell_{2}(\hat{\Theta}_{1},\Theta_{2}\|D)+\log\pi(\Theta_{2})+\text% {constant},$		(22)

where $\ell_{2}(\hat{\Theta}_{1},\Theta_{2}|D)$ is given by (14) in regard to the DAG copulas, and $\pi(\Theta_{2})$ is the prior distribution for $\Theta_{2}$ . We choose a uniform prior for $\Theta_{2}$ , which is proportional to 1 if the parameters satisfy the stationary conditions and the finite covariance condition, i.e.,

		$\displaystyle\pi(\Theta_{2})$		(23)
	$\displaystyle\propto$	$\displaystyle\prod_{i=2}^{m}\prod_{k=1}^{n(i)}\mathbf{1}(0\leq a_{i,i[k]\|z_{ik% }},b_{i,i[k]\|z_{ik}}<1,a_{i,i[k]\|z_{ik}}+b_{i,i[k]\|z_{ik}}<1,-1<\bar{\varphi}_% {i,i[k]\|z_{ik}}<1,v_{i,i[k]\|z_{ik}}>2),$		(23)

where $z_{ik}=\{i[1],\ldots,i[k-1]\}$ and $\mathbf{1}(\cdot)$ is an indicator function that takes 1 if the condition inside the indicator function is satisfied, and 0 otherwise. Note that we use the IFM method so that we insert $\hat{\Theta}_{1}$ , the maximum likelihood estimates of the parameters in the marginal distributions into $\log f(\Theta_{2}|\hat{\Theta}_{1},D)$ . We use $\overline{\Theta}_{2}$ , the medians of the samples drawn from the posterior distribution as the Bayes estimates of $\Theta_{2}$ .

3.

For each stock $j=m+1,\ldots,m+p$ , we conduct sequential estimation for $\Theta_{3j}$ by optimizing the function $\ell_{3j}(\hat{\Theta}_{1},\overline{\Theta}_{2},\Theta_{3j})$ in (15), with stationary condition constraints in (8). Note that we substitute $\hat{\Theta}_{1}$ into (15) due to the method of IFM, and $\overline{\Theta}_{2}$ into (15) as we estimate the DAG parameters separately.

4.2 MCMC sampling for parameters in DAG copulas

MCMC sampling is used for the estimation of $\Theta_{2}$ . We first obtain the sequential estimates of the parameters as the initial values for the MCMC sampling, denoted by $\Theta^{(0)}_{2}$ . We adopt the Robust adaptive Metropolis (RAM) algorithm \parenciteRAM for the MCMC sampling. Starting from $n=1$ , we do the following steps:

We set a proposal

\Theta_{2}^{*}=\Theta_{2}^{(n-1)}+S_{n-1}U_{n},

where $U_{n}$ is a $d$ -dimension vector drawn from the proposal distribution $q(\cdot)$ , $d=\absolutevalue{\Theta_{2}}$ is the number of parameters in $\Theta_{2}$ , and $S_{n-1}$ is a matrix that captures the correlation between the parameters. We use the multivariate normal distribution with mean equals a $d$ -dimensional vector of zeros and covariance matrix equals the $d$ -dimensional identity matrix as the proposal distribution. However, changing all parameters at once may lead to poor mixing in the MCMC sampling. We instead change a subset of parameters in $\Theta_{2}$ in each step. This can be done by setting the entries in $U_{n}$ corresponding to the variables that are not in the subset to 0.

2.

With a probability of $\alpha_{n}=\min\{1,f(\Theta_{2}^{*}|\hat{\Theta}_{1},D)/f(\Theta_{2}^{(n-1)}|% \hat{\Theta}_{1},D)\}$ , the proposal is accepted and we set $\Theta_{2}^{(n)}=\Theta_{2}^{*}$ . Otherwise, we set $\Theta_{2}^{(n)}=\Theta_{2}^{(n-1)}$ .

Update the lower-diagonal matrix with positive diagonal elements $S_{n}$ satisfying the equation

S_{n}S_{n}^{T}=S_{n-1}\left(I+\eta_{n}(\alpha_{n}-\alpha_{*})\frac{U_{n}U_{n}^% {T}}{\norm{U_{n}}^{2}}\right)S_{n-1}^{T},

where $I$ is the $d$ by $d$ identity matrix, $\alpha_{*}\in(0,1)$ is the target mean acceptance probability of the algorithm, and $\{\eta_{n}\}_{n\geq 1}$ is a step size sequence decaying to zero, where $\eta_{n}\in(0,1]$ . It can be done using Cholesky decomposition.

4.

Set $n\leftarrow n+1$ and go back to step 1.

The advantage of the RAM algorithm is that the matrix $S_{n}$ captures the correlation between parameters and allows the algorithm to attain a given mean acceptance rate $\alpha_{*}$ , which helps improve the convergence compared to the traditional MCMC sampling algorithms. $S_{0}$ is initialized with the $d$ by $d$ diagonal matrix. $\alpha_{*}$ is set to 23.4%, which is a typical choice in literature \parencite10.2307/2245134. The step size sequence $\eta_{n}$ is often defined as $\eta_{n}=n^{-\gamma}$ , where $\gamma\in(1/2,1]$ \parenciteRAM. In the simulation study and empirical study, we choose $\gamma=2/3$ . The sequence $\eta_{n}$ plays an important role to ensure the ergodicity of the MCMC algorithm.

4.3 Structural learning of the DAG

In practice, the underlying DAG of the risk factor returns is often unknown, and we need to estimate the DAG from the data. Score-based learning using MCMC sampling is adopted, which searches for DAGs that are around the posterior modes. In practice, there could be multiple DAGs that have similar scores. Specifically, to construct a score function, consider the term $\ell_{2}(\hat{\Theta}_{1},\Theta_{2}|D)$ in (14) with regard to the DAG parameters, in which we insert $\hat{\Theta}_{1}$ , the maximum likelihood estimates of the parameters in the marginal distribution due to the IFM method. Note that we do not include the likelihood function for stock copulas because we do not estimate the structure of the risk factors using the stock returns. To emphasize that the likelihood function is a function of a network $G$ , we rewrite the likelihood function as $\exp(\ell_{2}(\hat{\Theta}_{1},\Theta_{2},G|D))\equiv P(D|\hat{\Theta}_{1},% \Theta_{2},G)$ . In DAG structural learning, we need to calculate the marginal likelihood that marginalizes out the parameters:

P(D|\hat{\Theta}_{1},G)=\int P(D|\hat{\Theta}_{1},\Theta_{2},G)\pi(\Theta_{2}|% G)d\Theta_{2},

(24)

where $\pi(\Theta_{2}|G)$ is the prior of $\Theta_{2}$ given $G$ . Note that the prior distribution of $\Theta_{2}$ in (23) is dependent on the DAG $G$ as $G$ can be varying. Under this case, we rewrite the prior for $\Theta_{2}$ in (23) as $\pi(\Theta_{2}|G)$ in (24).

The integration in the marginal likelihood in (24) is, however, difficult to calculate in practice since the integrand involves a product of $t$ -copulas. An approximation is to use the Bayesian information criterion (BIC) derived by \textciteBIC:

\log P(D|\hat{\Theta}_{1},G)\approx\log P(D|\hat{\Theta}_{1},\hat{\Theta}_{2},% G)-\frac{\absolutevalue{\Theta_{2}}}{2}\log T,

where $\absolutevalue{\Theta_{2}}$ is the number of parameters in $\Theta_{2}$ , $\hat{\Theta}_{2}$ is the maximum likelihood estimates of $\Theta_{2}$ , and $T$ is the number of days of observations in the data set $D$ . However, calculating the maximum likelihood of $\Theta_{2}$ of the function $P(D|\hat{\Theta}_{1},\Theta_{2},G)$ is in general time consuming, since the dimension of $\Theta_{2}$ is moderately large. It is not feasible since it often involves at least thousands of iterations in the MCMC sampling. Thus, we can further approximate the BIC using the sequential estimates, which perform well in the literature \parenciteHOBAEKHAFF201291,Haff2013ParameterEF. Denote $\tilde{\Theta}_{2}$ to be the sequential estimates obtained by optimizing each copula in $P(D|\hat{\Theta}_{1},\Theta_{2},G)$ separately. The score function used in the DAG structural learning is

\log P(D|\hat{\Theta}_{1},G)\approx\log P(D|\hat{\Theta}_{1},\tilde{\Theta}_{2% },G)-\frac{\absolutevalue{\Theta_{2}}}{2}\log T.

(25)

Now we discuss the MCMC scheme for the DAG structural learning. We use the Structure MCMC proposed by \textcitestructure_mcmc. It samples graphs by conducting local edge movements. Let $\text{nbd}(G)$ be the set of neighborhood of the graph $G$ , defined as the set containing the graphs that are equal to $G$ except with a one-edge difference, and in the reduced DAG space introduced in Section 3.4. A neighbor graph $G^{\prime}$ is sampled randomly, and thus the transition kernel is given by

q(G^{\prime}|G)=\frac{1}{\absolutevalue{\text{nbd}(G)}},

where $\absolutevalue{\text{nbd}(G)}$ is the number of graphs in $\text{nbd}(G)$ . Under the Bayesian framework, we also need to derive the posterior distribution $P(G|\hat{\Theta}_{1},D)$ , which is given by

P(G|\hat{\Theta}_{1},D)\propto P(D|\hat{\Theta}_{1},G)\pi(G)\approx\exp(\log P% (D|\hat{\Theta}_{1},\tilde{\Theta}_{2},G)-\frac{\absolutevalue{\Theta_{2}}}{2}% \log T)\pi(G),

where we approximate the likelihood using (25) and $\pi(G)$ is the prior distribution of the graph $G$ . We set a flat prior $\pi(G)\propto 1$ if $G$ is a DAG and inside the reduced DAG space since we do not have any prior information with regard to the structure of the risk factor returns. We then apply the Metropolis-Hasting algorithm to sample graphs from the posterior distribution $P(G|\hat{\Theta}_{1},D)$ . We sample a graph $G^{\prime}$ from $\text{nbd}(G)$ . The graph $G^{\prime}$ is accepted with the probability

\alpha(G,G^{\prime})=\min\left\{1,\frac{P(G^{\prime}|\hat{\Theta}_{1},D)q(G|G^% {\prime})}{P(G|\hat{\Theta}_{1},D)q(G^{\prime}|G)}\right\}.

5 Portfolio selection

5.1 Minimum variance (MV) portfolio and minimum conditional value-at-risk (MCVaR) portfolio

The GC-GARCH model can be used for portfolio selection. Suppose that we take $G^{*}$ , the maximum a posteriori (MAP) network in the structural learning as the underlying network in the GC-GARCH model. Consider a portfolio of $p$ assets and $m$ risk factors. Suppose that we invest in a portfolio assigning weight $w_{j,t}$ to stock $j$ at the end of day $t$ , where $j=m+1,\ldots,m+p$ and $\sum_{j=m+1}^{m+p}w_{j,t}=1$ . The portfolio return on day $t+1$ (tomorrow) is given by

r_{PF,t+1}=\sum_{j=m+1}^{m+p}w_{j,t}r_{j,t+1}\equiv\mathbf{w}_{t}^{T}\mathbf{r% }_{t+1},

where $\mathbf{w}_{t}=(w_{m+1,t},\ldots,w_{m+p,t})^{T}$ is a vector of portfolio weights and $\mathbf{r}_{t+1}=(r_{m+1,t+1},\ldots,r_{m+p,t+1})^{T}$ is a vector of stock returns on date $t+1$ . Note that the returns are indexed by $t+1$ but the portfolio weights are indexed by $t$ , because the portfolio weight is determined before we observe $\mathbf{r}_{t+1}$ . We assume that the portfolio weights are determined 1 day before. The variance of the portfolio return on date $t+1$ is

\sigma_{PF,t}^{2}=\sum_{i=m+1}^{m+p}\sum_{j=m+1}^{m+p}w_{i,t}w_{j,t}\sigma_{ij% ,t+1}\equiv\mathbf{w}_{t}^{T}\Sigma_{t+1}\mathbf{w}_{t},

where $\sigma_{ij,t+1}$ is the conditional covariance between stocks $i$ and $j$ on date $t+1$ given $\mathcal{F}_{t}$ , and $\Sigma_{t+1}=[\sigma_{ij,t+1}]_{i,j=m+1}^{m+p}$ is the covariance matrix of the stock returns. We estimate $\Sigma_{t+1}$ by simulating returns from the fitted GC-GARCH model.

Using the minimum variance optimization proposed by \textcitemarkowitz, we determine the optimal portfolio weights $\mathbf{w_{t}}$ that minimize the 1-day ahead predicted variance of the portfolio return given $\mathcal{F}_{t}$ . Mathematically, we determine $\mathbf{w}_{t}$ for the problem

\min_{\mathbf{w_{t}}}\mathbf{w}_{t}^{T}\hat{\Sigma}_{t+1}\mathbf{w}_{t}~{}~{}% \text{ such that }\sum_{j=m+1}^{m+p}w_{j,t}=1,

(26)

where $\hat{\Sigma}_{t+1}$ is the 1-day ahead prediction of the conditional covariance of $\mathbf{r}_{t+1}$ matrix given $\mathcal{F}_{t}$ from the GC-GARCH model.

To estimate $\hat{\Sigma}_{t+1}$ for the problem (26), we simulate $K=20,000$ vectors of stock returns on day $t+1$ from the GC-GARCH model, $\mathbf{r}_{1,t+1},\ldots,\mathbf{r}_{K,t+1}$ , using the methodology in Section 3.5. We first estimated the 1-day ahead prediction of the conditional correlation matrix of $\mathbf{r}_{t+1}$ :

\hat{\Upsilon}_{t+1}=\text{diag}(\tilde{\Sigma}_{t+1})^{-1/2}\tilde{\Sigma}_{t% +1}\text{diag}(\tilde{\Sigma}_{t+1})^{-1/2},

where

\tilde{\Sigma}_{t+1}=\frac{1}{K}{\sum_{k=1}^{K}\mathbf{r}_{k,t+1}\mathbf{r}_{k% ,t+1}^{T}},

(27)

and $\text{diag}(\tilde{\Sigma}_{t+1})$ is the diagonal matrix containing the diagonal elements in $\tilde{\Sigma}_{t+1}$ . Then, we estimate the 1-day ahead conditional covariance matrix of $\mathbf{r}_{t+1}$ by

\hat{\Sigma}_{t+1}=\Lambda_{t}^{1/2}\hat{\Upsilon}_{t+1}\Lambda_{t}^{1/2},

(28)

where $\Lambda_{t}$ is the diagonal matrix, whose diagonals are the exact 1-day ahead predictions of the variances obtained from the marginal GARCH(1,1) models. Note that both (27) and (28) are estimators of the 1-day ahead covariance matrix, whereas (28) matches the variances to the exact 1-day ahead predictions of the variances. We prefer to use (28).

As discussed in \textcitebrechmann2013risk, since only the second moment is taken into account in the minimum variance portfolio problem, the information regarding to the tails is not considered. To also include the information in the tails, we solve the minimum conditional value-at-risk (MCVaR) portfolio problem:

\min_{\mathbf{w_{t}}}CVaR_{t+1}^{\alpha}(\mathbf{w}_{t}^{T}\mathbf{r}_{t+1})~{% }~{}\text{ such that }\sum_{j=m+1}^{m+p}w_{j,t}=1,

(29)

where

CVaR_{t+1}^{\alpha}(\mathbf{w}_{t}^{T}\mathbf{r}_{t+1})=-E(\mathbf{w}_{t}^{T}% \mathbf{r}_{t+1}|\mathbf{w}_{t}^{T}\mathbf{r}_{t+1}<q_{t+1}^{\alpha}(\mathbf{w% }_{t}^{T}\mathbf{r}_{t+1}),\mathcal{F}_{t})

and $q_{t+1}^{\alpha}(\mathbf{w}_{t}^{T}\mathbf{r}_{t+1})$ is the $\alpha$ th quantile of the 1-day ahead portfolio return $\mathbf{w}_{t}^{T}\mathbf{r}_{t+1}$ . \textciterockafellar2000optimization and \textciteVaR_optim_comp discuss an alternative method to approximate the MCVaR portfolio problem by generating returns from the model. We approximate the solution in (29) by sampling $K=20,000$ vectors of 1-day ahead stock returns from the GC-GARCH model, and solve

\min_{\mathbf{w_{t}},\mathbf{u},a}\left(-a+\frac{1}{K\alpha}\sum_{k=1}^{K}u_{k% }\right)

(30)

subject to the linear constraints

		$\displaystyle\sum_{j=m+1}^{m+p}w_{j,t}=1,$
		$\displaystyle u_{k}\geq 0,\text{ for }k=1,\ldots,K,\text{and}$
		$\displaystyle\mathbf{w}_{t}^{T}\mathbf{r}_{k,t+1}-a+u_{k}\geq 0,\text{ for }k=% 1,\ldots,K,$

where $\mathbf{u}=(u_{1},\ldots,u_{K})$ are auxiliary variables, $a$ is the VaR level, and $\mathbf{r}_{k,t+1}$ is the $k$ th sampled vector of 1-day ahead stock returns, $k=1,\ldots,K$ .

If we estimate the portfolio weights by solving the MCVaR portfolio problem (30), the minimized objective function $\left(-a+{(K\alpha)^{-1}}\sum_{k=1}^{K}u_{k}\right)$ is exactly the 1-day ahead prediction of CVaR. If we want to predict the 1-day ahead CVaR from the portfolio weights obtained in the MV portfolio problem in (26), this can be done by generating $K=20,000$ portfolio returns and estimate the 1-day ahead CVaR by

CVaR_{t+1}^{\alpha}=-\frac{1}{g_{t}}\sum_{k=1}^{K}\mathbf{w}_{t}^{T}\mathbf{r}% _{k,t+1}\mathbf{1}(\mathbf{w}_{t}^{T}\mathbf{r}_{k,t+1}<\hat{q}_{t+1}^{\alpha}% (\mathbf{w}_{t}^{T}\mathbf{r}_{t+1})),

(31)

where $\hat{q}_{t+1}^{\alpha}(\mathbf{w}_{t}^{T}\mathbf{r}_{t+1})$ is the $\alpha$ th sample quantile of the generated $K=20,000$ portfolio returns and $g_{t}=\sum_{k=1}^{K}\mathbf{1}(\mathbf{w}_{t}^{T}\mathbf{r}_{k,t+1}<\hat{q}_{t% +1}^{\alpha}(\mathbf{w}_{t}^{T}\mathbf{r}_{t+1}))$ is the number of generated portfolio returns that are smaller than $\hat{q}_{t+1}^{\alpha}(\mathbf{w}_{t}^{T}\mathbf{r}_{t+1})$ .

5.2 MV and MCVaR portfolios with model averaging

In Section 5.1, we solve the MV portfolio problem in (26) and the MCVaR portfolio problem in (30) using the MAP graph $G^{*}$ . However, in the structural learning, there could be multiple graphs that have similar network scores as $G^{*}$ . Then, we can conduct model averaging over the highest scoring $N_{g}$ networks for the portfolio selection instead of using only the MAP network. Model averaging methods are shown to outperform the model selection method in the literature \parencite10.5555/3023549.3023619,dash2004model. We denote $G^{(j)}$ to be the $j$ th highest scoring network (then, $G^{(1)}=G^{*}$ ). To estimate the 1-day ahead conditional covariance matrix by modeling averaging, we simulate $K=20,000$ vectors of stock returns from each estimated GC-GARCH model with network $G^{(j)}$ (we denote such estimated model as $\mathcal{M}^{(j)}$ ), $j=1,\ldots,N_{g}$ . We denote $r_{k,t+1}^{(j)}$ to be the $k$ th vector of stock returns generated from $\mathcal{M}^{(j)}$ . The 1-day ahead conditional covariance matrix under $\mathcal{M}^{(j)}$ , $\Sigma_{t+1}^{(j)}$ , is estimated similarly with the formula (28) using the vectors $r_{k,t+1}^{(j)}$ , $k=1,\ldots,K$ . The model-averaged 1-day ahead conditional covariance matrix is given by

\displaystyle\overline{\Sigma}_{t+1}

\displaystyle=\sum_{j=1}^{N_{g}}\Sigma_{t+1}^{(j)}P(\mathcal{M}^{(j)})

where

	$\displaystyle P(\mathcal{M}^{(j)})$	$\displaystyle=\frac{P(G^{(j)}\|D,\hat{\Theta}_{1})}{\sum_{j=1}^{N_{g}}P(G^{(j)}% \|D,\hat{\Theta}_{1})}$
	$\displaystyle\approx$	$\displaystyle\frac{\exp(BIC^{(j)})}{\sum_{j=1}^{N_{g}}\exp(BIC^{(j)})}.$

Note that $\hat{\Theta}_{1}$ is the set of estimated parameters in the marginal distributions that we estimate separately due to the IFM method. As in the structural learning, the marginalized posterior probabilities are difficult to calculate, and we use the BICs in (25) to approximate them. $BIC^{(j)}$ is the score of $G^{(j)}$ in the structural learning in (25). The probability $P(\mathcal{M}^{(j)})$ is the marginalized posterior probability of $G^{(j)}$ normalized by the sum of the marginalized posterior probabilities of all candidates $G^{(1)},\ldots,G^{(N_{g})}$ (to ensure that $\sum_{j=1}^{N_{g}}P(\mathcal{M}^{(j)})=1$ ). We also assume a uniform prior for $G^{(j)}$ so that

P(G^{(j)}|D,\hat{\Theta}_{1})\propto P(D|\hat{\Theta}_{1},G^{(j)})\approx\exp(% BIC^{(j)}).

Then, we solve the MV problem in (26) using the covariance matrix $\overline{\Sigma}_{t+1}$ .

To solve the MCVaR problem, we draw $K=20,000$ vectors of returns from the GC-GARCH models $\mathcal{M}^{(1)},\ldots,\mathcal{M}^{(N_{g})}$ , where the probability of sampling from the model $\mathcal{M}^{(j)}$ is $P(\mathcal{M}^{(j)})$ . Then, we use the sampled $K$ vectors to solve the MCVaR portfolio problem in (30). The estimation methods of CVaR of the portfolio returns are discussed in Section 5.1. The sampled $K$ vectors can be viewed as the realizations from the mixture distribution of $\mathcal{M}^{(1)},\ldots,\mathcal{M}^{(N_{g})}$ with weights $P(\mathcal{M}^{(1)}),\ldots,P(\mathcal{M}^{(N_{g})})$ .

6 Simulation study

6.1 Simulation settings

In this simulation study, we aim to test our estimation method in Section 4. We conduct the following:

1.

When the underlying DAG is known, we show that the parameters in the copulas can be accurately estimated using the method in Section 4.1.
2.

We show that the structural learning algorithm in Section 4.3 can search for structures that are close to the underlying DAG.

We conduct the following two sets of simulation studies:

(S1)

$m=8$ risk factors, $p=100$ stocks and $T=1000$ days of observations. We pick a relative simple network.
(S2)

$m=10$ risk factors, $p=100$ stocks and $T=1000$ days of observations. We pick a more complicated graph in this case.

We repeat each experiment for 200 replications using the same DAG and the same set of parameters. We first simulated the common set of parameters:

(i)

The parameters in the marginal distributions in (9) are generated from the following distributions: $w_{i}\sim U(0.01,0.2)$ , $b_{i}\sim U(0.8,0.96)$ , $a_{i}\sim(1-b_{i})\times U(0.3,0.7)$ , and $v_{i}\sim U(5,10)$ . $U(k_{1},k_{2})$ denotes a uniform distribution over the interval $(k_{1},k_{2})$ for $k_{1}<k_{2}$ . The method ensures that $\omega_{i}>0$ , $a_{i},b_{i}\geq 0$ and $a_{i}+b_{i}<1$ for all $i=1,\ldots,m+p$ .
(ii)

The parameters in the DAG copulas and the stock copulas in (8) are generated from the following distributions: $\bar{\varphi}_{xy|z}\sim U(-1,1)$ , $b_{xy|z}\sim U(0.8,0.96)$ , $a_{xy|z}\sim(1-b_{xy|z})\times U(0.3,0.7)$ , and $v_{xy|z}\sim U(5,10)$ . The method ensures that $0\leq a_{xy|z},b_{xy|z}<1$ , $a_{xy|z}+b_{xy|z}<1$ , and $-1<\bar{\varphi}_{xy|z}<1$ .

In each replication, the data set is generated using the method in Section 3.5. We use a simpler setting under (S1) to illustrate the ideas in the proposed GC-GARCH model, whereas (S2) is a more complicated case that is closer to the cases in the empirical study. We present the simulation results for (S1) and (S2) below. We first present the performances of the structural learning. Then, we present the performances of the estimation of the parameters in the DAG copulas and stock copulas.

6.2 Structural learning evaluation

MCMC sampling is used to sample graphs from the posterior modes using the methods in Section 4.3. For each replication, we conduct the MCMC sampling for $200$ iterations. Let $G_{1},\ldots,G_{200}$ be the sampled graphs. We discard the first $B$ sampled graphs in the burn-in step, where $B$ is selected using the Geweke’s statistic (to be introduced later in this section). Using the graphs after burn-in, $G_{B+1},\ldots,G_{200}$ , we estimate the probability of the existence of an edge from $r_{i,t}$ to $r_{j,t}$ by

P(a_{ij}=1|D)=\frac{1}{200-B}\sum_{k=B+1}^{200}\mathbf{1}(a_{ij}(G_{k})=1),

(32)

where $a_{ij}$ denotes the $(i,j)$ th entry of the adjacency matrix of the true graph (which is assumed to be random in Bayesian analysis), and $\mathbf{1}(a_{ij}(G_{k})=1)$ equals 1 if the edge from $r_{i,t}$ to $r_{j,t}$ exists in the graph $G_{k}$ , and equals 0 otherwise. (32) counts the proportion of graphs that contain the edge from $r_{i,t}$ to $r_{j,t}$ , and is called an edge feature in the literature \parenciteorder_mcmc.

The edges features are used to predict the existences of edges. We predict the edge from $r_{it}$ to $r_{jt}$ to exist if $P(a_{ij}=1|D)>c$ for some threshold $c\in[0,1]$ . By altering the values of $c$ , we trace a receiver operating characteristic (ROC) curve. We use the completed partially directed acyclic graph (CPDAG) \parenciteCPDAG of the underlying network to evaluate the accuracy of the prediction, as the CPDAG represents an equivalence class of a set of graphs, where the graphs in the equivalence class have the same set of conditional independence. The CPDAG contains two types of edges. The first type is a directed edge, indicating that all graphs in the class agree with the direction. The second type is a bi-directed edge, indicating that some of the graphs disagree with the direction of an edge. The area under the ROC curve (AUROC) is used to measure the performance of the graph prediction.

On top of the AUROC, we also use the following metrics. Let $TP(c)$ , $TN(c)$ , $FP(c)$ and $FN(c)$ be respectively the numbers of true positive, true negative, false positive and false negative predictions when we use the threshold $c$ . We have the following five metrics

(1)

The accuracy

ACC(c)=\frac{TP(c)+TN(c)}{TP(c)+FN(c)+FP(c)+TN(c)},

(2)

the false discovery rate

FDR(c)=\frac{FP(c)}{TP(c)+FP(c)},

(3)

the false omission rate

FOR(c)=\frac{FN(c)}{FN(c)+TN(c)},

(4)

the sensitivity

SEN(c)=\frac{TP(c)}{TP(c)+FN(c)},\text{ and}

(5)

the specificity

SPE(c)=\frac{TN(c)}{TN(c)+FP(c)}.

The accuracy (ACC) measures the proportion of correct prediction, which gives an overall performance of the prediction. The false discovery rate (FDR) measures the proportion of incorrectly predicted edges out of all predicted edges. The false omission rate (FOR) measures the proportion of incorrectly omitted edges out of all omitted edges. The sensitivity (SEN) measures the accuracy of detecting edges that exist in the underlying network, and the specificity (SPE) measures the accuracy of omitting edges that do not exist in the underlying network. Yet the SEN and SPE have been used to calculate the AUROC by changing the threshold $c$ , and we also include the numerical results for reference.

To evaluate the convergence of $\{G_{k}\}_{k=B+1}^{200}$ , we calculate the Geweke’s statistic \parencitegeweke of a characteristic of the graph. We need to transform the network into a univariate time series $\{d(G_{k})\}_{k=B+1}^{200}$ , for some function $d:\{0,1\}^{m\times m}\mapsto\mathbb{R}$ , to calculate the Geweke’s statistic. We pick

d(G_{k})=\sum_{i=1}^{m}\sum_{j=1}^{m}\mathbf{1}(a_{ij}(G_{k})\neq a_{ij}),

where $a_{ij}(G_{k})$ is the $(i,j)$ th entry of the adjacency matrix of $G_{k}$ , and $a_{ij}$ is the $(i,j)$ th entry of the adjacency matrix of the true network, which is known in the simulation study. $\mathbf{1}(a_{ij}(G_{k})\neq a_{ij})$ equals 1 if the $(i,j)$ th entries in the adjacency matrices of $G_{k}$ and the true network, respectively, disagree with each other, and 0 zero otherwise. $d(G)$ is to keep track if the differences between the sampled graphs and the true graph are stabilized. The Geweke’s statistic tests if the average in the first 10% of $\{d(G_{k})\}_{k=B+1}^{200}$ is significantly different from the average in the last 50% of $\{d(G_{k})\}_{k=B+1}^{200}$ . If the p-value is greater than 0.01, this indicates that the difference of the two averages is not significant, and thus we can conclude that the samples are converged.

The underlying network in simulation (S1) is shown in Figure 3. We conduct 200 replications, where we generate different data sets from the same DAG and the same set of parameters. The burn-in $B$ in each replication is selected to minimize the Geweke’s statistic. The average AUROC out of 200 replications is 75.4%, showing that the structural learning performs well. The ACC, FDR, FOR, SEN and SPE in simulation (S1) using different thresholds ( $c=0.01,0.10,0.50,0.90,0.99$ ) is shown in Table 1. The ACC, SEN and SPE are moderatly high and the FDR and FOR are small, together with a good AUROC, indicating that many of the existing edges in the underlying graph are detected. The p-values of the Geweke’s statistics in 187 out of 200 replications are greater than 0.01, supporting that the burn-in sample is converged.

We do the same for the simulation (S2). The underlying network in simulation (S2) is shown in Figure 4. It contains more nodes and edges than the graph in simulation (S1), and we expect that the structural learning in simulation (S2) will be more challenging. The average AUROC out of 200 replications is 70.1%. The ACC, FDR, FOR, SEN and SPE in simulation (S1) is shown in Table 2. We observe that the false discovery rate becomes higher than the case in simulation (S1). This indicates that the structural learning algorithms tend to include more edges, especially in higher dimensions. Nevertheless, the algorithm has a sensitivity up to 0.642. The p-values of the Geweke’s statistics in 184 out of 200 replications are greater than 0.01. We still get a moderately good AUROC and convergence for a more challenging case.

However, we do not know the underlying network of a real data set. These two simulation results give some clues of the performance of the structural learning algorithm. The algorithm is able to cover a moderately high portion of existing edges in the underlying network. Then, in the empirical study, we also conduct 200 iterations for the structural learning.

Figure 3: The underlying network in the simulation (S1).

Threshold $c$	ACC	FDR	FOR	SEN	SPE
0.01	0.810	0.347	0.111	0.723	0.845
0.10	0.818	0.315	0.125	0.676	0.874
0.50	0.819	0.238	0.163	0.537	0.932
0.90	0.781	0.204	0.218	0.323	0.965
0.99	0.766	0.202	0.234	0.254	0.971

Table 1: The average ACC, FDR, FOR, SEN and SPE of the structural learning in simulation (S1) out of 200 replications at different thresholds.

Figure 4: The underlying network in the simulation (S2).

Threshold $c$	Accuracy	FDR	FOR	SEN	SPE
0.01	0.744	0.447	0.153	0.642	0.786
0.10	0.748	0.440	0.163	0.609	0.804
0.50	0.750	0.420	0.192	0.500	0.852
0.90	0.755	0.397	0.204	0.440	0.883
0.99	0.752	0.397	0.211	0.409	0.891

Table 2: The average ACC, FDR, FOR, SEN and SPE of the structural learning in simulation (S2) out of 200 replications at different thresholds.

6.3 Parameter estimation evaluation

The parameter estimations are separated into three phases: We first estimate $\Theta_{1}$ , the parameters in the marginal distributions. We then estimate $\Theta_{2}$ , the parameters in the DAG copulas using MCMC sampling. Finally, for each stock $j$ , we estimate $\Theta_{3j}$ , the parameters in the stock copulas $c_{j,i|1,\ldots,i-1}^{[t]}$ using sequential estimation starting from $i=1$ to $i=m$ . We do not report the estimation results for the parameters in the marginal distributions as our interest is the copula parameters.

6.3.1 Parameters in the DAG copulas

We conduct 20,000 iterations for learning the parameters in the MCMC sampling for the DAG copulas in each replication. The burn-in is chosen automatically by minimizing the average Geweke’s statistics of the chains of the parameters. We first present the estimation results for simulation (S1). We present the true values, the mean, the mean absolute error (MAE) of the mean, the median, the MAE of the median, the 5th quantile ( $q_{0.05}$ ), and the 95th quantile ( $q_{0.95}$ ) of the posterior median estimates out of the 200 replications in Table 3. A true value marked with an asterisk indicates that $(q_{0.05},q_{0.95})$ , the 90% credible interval contains the true value. All the true values are contained in the 95% credible intervals except for the parameter $\bar{\varphi}_{8,5}$ . The MAEs of the mean and median of $\bar{\varphi}_{xy|z}$ are very small, indicating that the estimation is precise. The MAEs for $a_{xy|z}$ $b_{xy|z}$ are, however, larger relative to the magnitudes of the parameters. This can be explained as follows. 5(b) shows the distributions of the posterior median estimates of (a) $a_{2,1}$ , (b) $b_{2,1}$ , (c) $\bar{\varphi}_{2,1}$ and (d) $v_{2,1}$ , across 200 replications, for the copula $c_{2,1}^{[t]}$ in simulation (S1) as an example. We observe that some of the estimates of $a_{xy|z}$ and $b_{xy|z}$ are close to zero. This is because the data sets contain noises, and the algorithm may not be able to detect the dynamic dependencies of some copulas from the data sets in a few replications. Nevertheless, the median of the estimates of $a_{xy|z}$ and $b_{xy|z}$ are close to the true values as seen in Table 3. The estimation results for simulation (S2) is in Appendix B.1 as the tables are too long to be included in the paper. The results are similar to those of simulation (S1), supporting that the MCMC estimation method in Section 4.2 is reliable.

Parameter	True value	Mean of estimates	Mean MAE	Median of estimates	Median MAE	$q_{0.05}$	$q_{0.95}$
$a_{2,1}$	$\mathbf{0.04^{*}}$	0.0298	0.0277	0.0269	0.0247	0	0.0829
$a_{3,2}$	$\mathbf{0.04^{*}}$	0.0331	0.025	0.0323	0.0198	0	0.0708
$a_{4,2}$	$\mathbf{0.04^{*}}$	0.0395	0.0214	0.0379	0.0168	0.0015	0.0745
$a_{5,2}$	$\mathbf{0.05^{*}}$	0.0466	0.0294	0.0458	0.0228	0.0004	0.0942
$a_{6,2}$	$\mathbf{0.02^{*}}$	0.017	0.0191	0.0111	0.0159	0	0.0581
$a_{7,2}$	$\mathbf{0.06^{*}}$	0.051	0.0361	0.0515	0.0301	0	0.1129
$a_{7,3\|2}$	$\mathbf{0.09^{*}}$	0.0898	0.0221	0.0896	0.0171	0.0577	0.1246
$a_{8,5}$	$\mathbf{0.12^{*}}$	0.0977	0.0307	0.096	0.0315	0.0497	0.1478
$b_{2,1}$	$\mathbf{0.89^{*}}$	0.6826	0.3549	0.8675	0.2568	0.0003	0.9947
$b_{3,2}$	$\mathbf{0.91^{*}}$	0.7973	0.2811	0.9066	0.1516	0.0165	0.9895
$b_{4,2}$	$\mathbf{0.88^{*}}$	0.7956	0.222	0.8664	0.1169	0.1002	0.9475
$b_{5,2}$	$\mathbf{0.84^{*}}$	0.6854	0.3141	0.8219	0.2102	0.0044	0.967
$b_{6,2}$	$\mathbf{0.96^{*}}$	0.6985	0.3853	0.9211	0.2747	0.0001	0.9962
$b_{7,2}$	$\mathbf{0.81^{*}}$	0.6308	0.3239	0.7719	0.2463	0.0002	0.9878
$b_{7,3\|2}$	$\mathbf{0.86^{*}}$	0.8412	0.0756	0.8495	0.0423	0.7593	0.9124
$b_{8,5}$	$\mathbf{0.8^{*}}$	0.7747	0.0851	0.7897	0.0661	0.6134	0.8878
$\bar{\varphi}_{2,1}$	$\mathbf{-0.1^{*}}$	$-$ 0.0945	0.0545	-0.0936	0.045	-0.1897	-0.009
$\bar{\varphi}_{3,2}$	$\mathbf{0.4^{*}}$	0.3774	0.0609	0.3758	0.0508	0.2867	0.481
$\bar{\varphi}_{4,2}$	$\mathbf{0.76^{*}}$	0.7346	0.0331	0.7365	0.0331	0.6806	0.7868
$\bar{\varphi}_{5,2}$	$\mathbf{0.59^{*}}$	0.5651	0.0462	0.562	0.0422	0.4874	0.6411
$\bar{\varphi}_{6,2}$	$\mathbf{0.34^{*}}$	0.3073	0.0657	0.3045	0.0572	0.2089	0.4085
$\bar{\varphi}_{7,2}$	$\mathbf{0.27^{*}}$	0.2496	0.0471	0.2497	0.041	0.171	0.3346
$\bar{\varphi}_{7,3\|2}$	$\mathbf{0.41^{*}}$	0.4011	0.0853	0.4042	0.0653	0.2541	0.5323
$\bar{\varphi}_{8,5}$	$-0.49$	$-$ 0.3225	0.0624	-0.3229	0.1678	-0.428	-0.2258
$v_{2,1}$	$\mathbf{5.29^{*}}$	5.5084	1.212	5.2921	0.8857	3.8757	7.7017
$v_{3,2}$	$\mathbf{5.22^{*}}$	5.765	1.1881	5.615	0.9539	4.1212	8.0312
$v_{4,2}$	$\mathbf{7.26^{*}}$	8.3123	2.359	7.6156	1.8309	5.568	12.7916
$v_{5,2}$	$\mathbf{9.71^{*}}$	10.3178	3.1879	9.7887	2.4571	6.3559	16.1122
$v_{6,2}$	$\mathbf{5.18^{*}}$	5.7413	1.5622	5.4904	1.1309	3.8794	8.7478
$v_{7,2}$	$\mathbf{7.07^{*}}$	7.9327	2.4533	7.3082	1.6834	5.3139	12.8052
$v_{7,3\|2}$	$\mathbf{6.59^{*}}$	7.342	2.3234	6.8102	1.7867	4.5986	11.9898
$v_{8,5}$	$\mathbf{7.27^{*}}$	10.6446	3.1184	10.0816	3.613	6.1601	16.0641

Table 3: The true values, the mean, mean absolute error (MAE) of the mean, the median, the MAE of the median, the 5th quantile (

q_{0.05}

), and the 95th quantile (

q_{0.95}

) of the posterior median estimates out of the 200 replications for simulation (S1).

6.3.2 Parameters in the stock copulas

The parameters in the stock copulas are estimated using sequential estimation as discussed in Section 4.1. C contains a comparison with the sequential estimates and the MCMC estimates of the parameters in the stock copulas, and the results suggest that the sequential estimates work as well as the MCMC estimates. Since there are $m\times p$ estimates for, respectively, $a_{j,i|1,\ldots,i-1}$ , $b_{j,i|1,\ldots,i-1}$ , $\bar{\varphi}_{j,i|1,\ldots,i-1}$ , and $v_{j,i|1,\ldots,i-1}$ , it is more efficient to present the results using graphs. We also partition the parameters into $m$ levels, where the first level contains the parameters $a_{j,1}$ , $b_{j,1}$ $\bar{\varphi}_{j,1}$ and $v_{j,1}$ , and the $i$ th level contains $a_{j,i|1,\ldots,i-1}$ , $b_{j,i|1,\ldots,i-1}$ $\bar{\varphi}_{j,i|1,\ldots,i-1}$ and $v_{j,i|1,\ldots,i-1}$ , for $i=2,\ldots,m$ . We first present the sequential estimation results for the stock parameters in Simulation (S1). Figure 6, Figure 7, Figure 8, and Figure 9 show the graphs summarizing the parameter estimates for respectively $a_{xy|z}$ , $b_{xy|z}$ , $v_{xy|z}$ and $\bar{\varphi}_{xy|z}$ by level. The green curves indicate true values, the blue and orange curves indicate respectively the means and medians of the sequential estimates out of 200 replications, and the pair of dashed curves represent the 90% credible interval out of 200 replications. The x-axes are in ascending order according to the true values of the parameters in each level for a better visualization.

Both the means and medians of the sequential estimates of $a_{xy|z}$ are close to the true values for the first three levels as observed from Figure 6. The 90% credible intervals are also narrow, indicating that the estimates are precise. The estimates tend to be slightly smaller than the true values in deeper levels. This happens because the estimates in the deeper levels rely on the estimates in previous levels, which may lead to larger errors. We expect the errors in deeper levels would have mild effects on the model. This is because, as we will see later in the empirical study, usually the long-run correlations are small (and thus not useful to capture the linear dependence among stocks) and the degrees of freedom are large (and thus not useful to capture the tail dependence) in deeper levels. Then, as long as the estimates of the first several levels are accurate, the model can already capture the linear and non-linear dependencies among stocks accurately.

Both the means and medians of the estimates of $b_{xy|z}$ are close to the true values, even in the deeper levels as observed from Figure 7. However, the 90% credible intervals are wild, especially in deeper levels. This could be due to the “information” stored in the data set. Different data sets may contain different noises so that the estimation algorithm may not detect the temporal dependence for the correlation, and thus the algorithm may give a small value of $b_{xy|z}$ . Nevertheless, the algorithm can recover the true values of $b_{xy|z}$ on average as the means and medians of the sequential estimates out of 200 replications are close to their true values.

The means and medians of the estimates of $v_{xy|z}$ are close to the true values in all levels as observed in Figure 8. The 90% credible intervals are also narrow whenever the true degrees of freedom are small. This shows that the algorithm is able to detect strong tail dependence precisely, which is crucial in practice. However, the 90% credible intervals are wilder when the degrees of freedom are larger and in deeper levels. This is again due to the “information” stored in the data set. For larger degrees of freedom (and thus weaker tail dependence), the data may not contain much information about the co-moving tail cases, and thus the estimates will be less precise. But it is not an issue in practice, since we are able to capture the strong dependence from the data set precisely, which is more important in application.

The means and medians of the sequential estimates of $\bar{\varphi}_{xy|z}$ are close to their true values and the 90% credible intervals are very narrow as observed in Figure 9. This shows that the sequential estimation performs very well in estimating long-run correlations.

The sequential estimation results for the stock copulas in Simulation (S2) are contained in Appendix B.2. The estimation results are similar to those in this section, and the comments in this section are also applicable to the results in Appendix B.2. This shows that the algorithm also works well for a more complicated network with $m=10$ variables.

7 Empirical study

7.1 Data

In this section, we consider a portfolio of $p=92$ stocks in Table 1 in Appendix D and $m=10$ market indexes in Table 4. We conduct a full GC-GARCH estimation and also conduct an investment experiment using a moving-window approach, in which we determine the portfolio weights for rebalancing by refitting the GC-GARCH model weekly.

The daily closing prices of $m+p=102$ stocks/market indexes are available from 4 January 2017 to 5 May 2023 ( $T+1=1,561$ trading days of closing prices and $T=1,560$ returns can be calculated from the closing prices). Let $S_{i,t}$ be the closing price of the $i$ th stocks/market indexes on day $t$ for $i=1,\ldots,m+p$ , $t=0,\ldots,T$ . $t=0$ corresponds to 4 Jan 2017. The log returns on $T$ trading days of the stocks/market indexes, $r_{i,t}=100\times\left[\log(S_{i,t})-\log(S_{i,t-1})\right]$ for $i=1,\ldots,m+p$ and $t=1,\ldots,T$ , are used to fit the GC-GARCH model. Similar to the simulation study, we say that a copula is in level $k$ if the copula has $k$ conditioning variables.

Name	Symbol	Order	Number of constituents	Portfolio included
Hang Seng Composite LargeCap Index	HSLI	1	122	67
Hang Seng China Enterprises Index	HSCEI	2	50	33
HK-Listed Biotech Index	HSHKBIO	3	73	3
Hang Seng Index	HSI	4	76	58
Hang Seng ESG 50 Index	HSESG50	5	50	33
China H-Financials Index	H.FIN	6	30	20
Hang Seng Climate Change 1.5 $\deg$ C Target Index	HSC15TI	7	212	21
Hang Seng TECH Index	HSTECH	8	30	4
Hang Seng Composite MidCap Index	HSMI	9	52	4
Growth Enterprise Market Index	HKSPGEM	10	49	0

Table 4: The name, symbol, topological order, number of constituent and the number of stocks included in the portfolio of the

p=92

stocks for each risk factor. The topological orders of the stocks are based on the MAP network in Figure 11.

7.2 Full estimation

We first estimate the GC-GARCH model using all $T=1560$ trading days of returns. We conduct a structural learning with 200 iterations using the method in Section 4.3. Figure 10 shows the time series of the BIC in the structural learning of the full estimation. The burn-in is chosen automatically by minimizing the Geweke’s statistic of the BIC. The structure with the highest score is shown in Figure 11. The topological orders of the risk factors are shown in Table 4. The parameter estimates of the DAG copulas are shown in Table 5.

The Hang Seng Composite LargeCap index (HSLI) is at the top of the order in the network in Figure 11, and connected to seven risk factors, showing that HSLI can be regarded as the “source” of the market movements. The HSLI also contains 67 stocks in our portfolio as shown in Table 4. We expect that the HSLI can explain a high portion of the co-movements of the stocks. The HSCEI, HSI and HSHKBIO are also ranked in top orders in the network. Moreover, the risk factors in the network are highly interconnected, indicating that these risk factors are closely dependent.

The parameter estimates in the DAG copulas are shown in Table 5. Most of $a_{j,i|1,\ldots,i-1}+b_{j,i|1,\ldots,i-1}$ are close to 1 (as most of $b_{j,i|1,\ldots,i-1}$ themselves are close to 1), indicating that the correlations are highly persistent. The long-run correlations tend to be smaller in magnitude in deeper levels; for example, the long-run correlations of (HSMI,HSLI), (HSMI,HSCEI $|$ HSLI), (HSMI,HSHKBIO $|$ HSLI,HSCEI), and (HSMI,HSI $|$ HSLI,HSCEI,HSHKBIO) are respectively 0.9350, 0.0554, 0.4562, and $-$ 0.1924. This means that high portions of the correlations among risk factors are explained by the first level.

Some degrees of freedom of the risk factor pairs are small, indicating that the dependence in the tails for some pairs of risk factor are strong, captured by the GC-GARCH model. For example, the degrees of freedom of (HSCEI,HSLI), (HSI,HSCEI), (HSESG50,HSCEI) and (HSESG50,HSI $|$ HSCEI) are respectively 5.59, 3.81, 4.51 and 5.53. It indicates that there are high levels of tail dependence among the risk factors.

As there are $m\times p=920$ parameters for each $a_{j,i|1,\ldots,i-1}$ , $b_{j,i|1,\ldots,i-1}$ , $v_{j,i|1,\ldots,i-1}$ and $\bar{\varphi}_{j,i|1,\ldots,i-1}$ in the stock copulas, we use scatter plots in Figure 12, Figure 13 and Figure 14 to report the parameter estimates of $a_{j,i|1,\ldots,i-1}$ , $b_{j,i|1,\ldots,i-1}$ and $v_{j,i|1,\ldots,i-1}$ . $\bar{\varphi}_{j,i|1,\ldots,i-1}$ ’s are used as the x-axes. We report the estimates by level, where level $i$ contains the parameters in the copulas of the form $C^{[t]}_{j,i|1,\ldots,i-1}$ . By reporting by different levels, we can observe the effects of the risk factors to the stocks clearly. The y-axes correspond to the estimates of the above three parameters, and the x-axes correspond to the estimates of $\bar{\varphi}_{j,i|1,\ldots,i-1}$ . The reason for putting $\bar{\varphi}_{j,i|1,\ldots,i-1}$ in the x-axes is that we observe that the estimates of $a_{j,i|1,\ldots,i-1}$ , $b_{j,i|1,\ldots,i-1}$ and $v_{j,i|1,\ldots,i-1}$ are associated with the values of the long-run correlations $\bar{\varphi}_{j,i|1,\ldots,i-1}$ . $a_{j,i|1,\ldots,i-1}$ and $b_{j,i|1,\ldots,i-1}$ tend to be smaller, and $v_{j,i|1,\ldots,i-1}$ tends to be large whenever the long-run correlation $\bar{\varphi}_{j,i|1,\ldots,i-1}$ is close to zero. This could be due to the “information” contained in the data set; when the long-run correlation $\bar{\varphi}_{j,i|1,\ldots,i-1}$ is close to zero, the “information” stored in the data set would be relatively weak, and thus the algorithm may not be able to detect the temporal dependence, captured by $a_{j,i|1,\ldots,i-1}$ and $b_{j,i|1,\ldots,i-1}$ , and the tail dependence, captured by $v_{j,i|1,\ldots,i-1}$ (as a larger degrees of freedom indicates weaker tail dependence). From Figure 13, we observe that most of the estimates of $b_{j1}$ in level 1 (correspond to the risk factor HSLI) are close to 1, and most of the long-run correlations $\bar{\varphi}_{j1}$ are larger than 0.5, showing that the correlations of these stocks between HSLI are dynamic and most of these stocks are strongly correlated with HSLI. The estimates of $b_{ji|1,\ldots,i-1}$ tend to be smaller in deeper levels, and the long-run correlations in deeper levels also tend to move towards 0, showing that high portions of the correlations are explained by the risk factors in the first several levels. Figure 14 shows that the degrees of freedom in the first level are small, and tend to be larger in deeper levels, indicating that the tail dependence between the stocks and respectively the first several risk factors are stronger. High portions of tail dependence are explained by the first several levels. Recall we have mentioned, in Section 6.3.2, that the estimates could become slightly imprecise (i.e., the credible intervals are wider) in deeper levels. We claim that the effects would be mild as the correlations are getting smaller and the degrees of freedom are getting larger in deeper levels as observed in Figure 14. Nevertheless, we have shown that the estimates are accurate on average even in deeper levels in the simulation study.

Copula	$\bar{\varphi}_{xy\|z}$	$a_{xy\|z}$	$b_{xy\|z}$	$v_{xy\|z}$
HSCEI,HSLI	0.9911	0.0102	0.9776	5.59
HSHKBIO,HSLI	0.7601	0.0275	0.9347	12.65
HSI,HSCEI	0.9910	0.0122	0.9734	3.81
HSI,HSHKBIO $\|$ HSCEI	0.0560	0.0000	0.9996	24.34
HSESG50,HSCEI	0.9783	0.0132	0.9731	4.51
HSESG50,HSI $\|$ HSCEI	0.9361	0.0143	0.9732	5.53
H.FIN,HSLI	0.9666	0.0135	0.9847	9.90
H.FIN,HSCEI $\|$ HSLI	0.9114	0.0020	0.9975	17.72
H.FIN,HSHKBIO $\|$ HSLI,HSCEI	-0.0724	0.0022	0.9936	19.53
H.FIN,HSESG50 $\|$ HSLI,HSCEI,HSHKBIO	0.5499	0.0247	0.9643	18.36
HSC15TI,HSLI	0.9929	0.0062	0.9758	6.16
HSC15TI,HSCEI $\|$ HSLI	0.4995	0.0108	0.9679	9.57
HSC15TI,HSHKBIO $\|$ HSLI,HSCEI	0.5294	0.0149	0.9568	13.22
HSC15TI,HSESG50 $\|$ HSLI,HSCEI,HSHKBIO	0.4742	0.0191	0.9646	27.08
HSC15TI,H.FIN $\|$ HSLI,HSCEI,HSHKBIO,HSESG50	0.0680	0.0128	0.9624	18.11
HSMI,HSLI	0.9350	0.0153	0.9537	13.18
HSMI,HSCEI $\|$ HSLI	0.0554	0.0055	0.9818	19.91
HSMI,HSHKBIO $\|$ HSLI,HSCEI	0.4562	0.0243	0.8494	15.61
HSMI,HSI $\|$ HSLI,HSCEI,HSHKBIO	-0.1924	0.0264	0.9003	20.49
HSMI,HSESG50 $\|$ HSLI,HSCEI,HSHKBIO,HSI	0.0285	0.0112	0.9880	18.61
HSMI,H.FIN $\|$ HSLI,HSCEI,HSHKBIO,HSI,HSESG50	-0.2130	0.0087	0.9910	15.87
HSMI,HSC15TI $\|$ HSLI,HSCEI,HSHKBIO,HSI,HSESG50,H.FIN	0.7572	0.0016	0.9979	18.08
HSTECH,HSLI	0.9680	0.0226	0.9605	12.91
HSTECH,HSCEI $\|$ HSLI	-0.3211	0.0106	0.9783	19.93
HSTECH,HSHKBIO $\|$ HSLI,HSCEI	0.3041	0.0029	0.9960	23.84
HSTECH,HSESG50 $\|$ HSLI,HSCEI,HSHKBIO	-0.4858	0.0246	0.9692	15.94
HSTECH,HSC15TI $\|$ HSLI,HSCEI,HSHKBIO,HSESG50	0.3677	0.0078	0.9905	19.65
HKSPGEM,HSLI	0.5129	0.0810	0.5853	9.56
HKSPGEM,HSCEI $\|$ HSLI	-0.0213	0.0000	0.5636	21.31
HKSPGEM,HSMI $\|$ HSLI,HSCEI	0.1800	0.0667	0.0000	31.26
HKSPGEM,HSTECH $\|$ HSLI,HSCEI,HSMI	0.0563	0.0023	0.9968	26.78

Table 5: The parameter estimates of the DAG copulas in the MAP network in Figure 11.

7.3 Moving-window investment experiment

We conduct an investment experiment where we rebalance the portfolio on the first trading day of each week using the methodology in Section 5. We aim to compare the investment performances between the models (1) GC-GARCH model using the MAP network, (2) GC-GARCH model with model averaging using the $N_{g}=3$ highest scoring networks, and (3) DCC-GARCH model proposed by \textciteDCC_model and \textciteDCC_GARCH.

The experiment is conducted with a moving-window size of $w=750$ days. We rebalance the portfolio weights on the first trading day of a week. To determine the portfolio weights, we use the data on the previous $w=750$ trading days (up to the last trading days of the week, which is typically a Friday) to fit the model, and solve the MV and MCVaR problems in (26) and (30).

Suppose that $t_{1},t_{2},\ldots,t_{L}\in\{1,\ldots,T\}$ are the last trading days of all weeks in the data set, where $L=331$ is the number of weeks in the data set. The first model can be fitted using the first 749 trading days (in the week 159). The reason for using 749 and not 750 for the first window is that, the 749th trading day in the data set is a Friday, which fits the setting that we fit a model at the end of a week. We further reserve the first 32 weeks of model fits for the sake of an investment strategy, to be introduced later. Then, the first investment can be done on 31 August 2020, which is the first day in week 192, where the portfolio weights is determined at the end of week $159+32=191$ . Let $R_{PF,t}=\exp(r_{PF,t})-1$ be the simple rate of return of the portfolio on day $t$ . Note that the portfolio weights are determined on last days of each of week $t_{191},t_{192},\ldots,t_{330}$ , and we make investment on the first trading days of each week $t_{191}+1,t_{192}+1,\ldots,t_{330}+1$ . The simple rate of return of the portfolio in week $i$ ( $i=192,\ldots,331$ ) is given by $R_{PF,\text{week }i}=\prod_{\tau=t_{i-1}+1}^{t_{i}}(1+R_{PF,\tau})$ . Suppose that we invest $10,000$ dollars on day $t_{191}+1$ in week 192 (the first week we invest), the cumulative value of the portfolio as of week $i$ ( $i=192,\ldots,331$ ) is given by

P_{i}=10,000\prod_{\tau=t_{191}+1}^{t_{i}}(1+R_{PF,\tau}).

(33)

We also evaluate the quality of prediction, following \textciteES_estimation, we compute the cost function with coverage level $\alpha$ ,

C(\alpha)=\frac{1}{g(\alpha)}\sum_{i=192}^{331}\absolutevalue{-R_{PF,t_{i-1}+1% }-CVaR^{\alpha}_{t_{i-1}+1}}\mathbf{1}(-R_{PF,t_{i-1}+1}\geq CVaR^{\alpha}_{t_% {i-1}+1}),

(34)

for $\alpha=0.0005,0.001,0.005,0.01,0.05$ , where $CVaR^{\alpha}_{t_{i-1}+1}$ is the predicted conditional value-at-risk with coverage level $\alpha$ on day $t_{i-1}+1$ , estimated using the methods of the portfolio returns introduced in Section 5.1, and $g(\alpha)=\sum_{i=192}^{331}\mathbf{1}(-R_{PF,t_{i-1}+1}\geq CVaR^{\alpha}_{t_% {i-1}+1})$ is the number of days that losses exceed the predicted CVaRs from weeks $192$ to $331$ .

There are two investment strategies:

•

(Strategy 1) We rebalance the portfolio on the first day of every week. We determine the portfolio weights by solving (i) the MV portfolio problem in (26), and (ii) the MCVaR portfolio problems with $\alpha=0.0005,0.001,0.005,0.01,0.05$ .

•

(Strategy 2) We avoid investing in the stock market whenever we predict a high level of 1-day ahead portfolio CVaR. Mathematically, at the end of week $i$ , we first calculate the average 1-day ahead predicted CVaR with coverage level $\alpha$ in the past $w_{S}$ weeks (excluding week $i$ ):

\overline{CVaR}^{\alpha}_{(i-w_{S}):(i-1)}=\frac{1}{w_{S}}\sum_{i^{\prime}=i-w% _{S}}^{i-1}CVaR_{t_{i^{\prime}}+1}^{\alpha}.

(35)

We then predict $CVaR_{t_{i}+1}^{\alpha}$ . We decide to invest in the stock market only when $CVaR_{t_{i}+1}^{\alpha}$ $\leq$ $\overline{CVaR}^{\alpha}_{(i-w_{S}):(i-1)}$ , i.e., the potential loss is predicted to be lower than or equal to the $w_{S}$ -week average. We take $w_{S}=8,16,32$ to test if the strategy performs uniformly well with different choice of $w_{S}$ . We only solve the MCVaR portfolio problems for this strategy as this strategy is designed using CVaR.

To ensure that we are comparing the performance of strategies 1 and 2 on the same trading days, we reserve the first 32 weeks of model fits as we need $32$ weeks to calculate (35) when $w_{S}=32$ . The first investment is done on 31 August 2020, which is the first day in week 192, for both strategies 1 and 2 with $w_{S}=8,16,32$ and $\alpha=0.0005,0.001,0.005,0.01,0.05$ .

We compare the investment performances of between the models (1) GC-GARCH model using the MAP network, (2) GC-GARCH model with model averaging using the $N_{g}=3$ highest scoring networks (we call it MA for simplicity), and (3) DCC-GARCH model. We determine the portfolio weights using (i) MCVaR portfolio with different $\alpha$ values and (ii) MV portfolio for strategy 1. We also estimate the CVaR from the MV portfolios using (31) because we want to compare the performance of the MCVaR portfolio to the traditional MV portfolio. We expect the MV portfolio to work poorly as the MV portfolio problem did not consider the tail behavior. We only consider the MCVaR portfolio with different $\alpha$ values for strategy 2. We first compare the cost function and the number of days that losses exceed the CVaR. We then compare the investment performances by plotting the time series of the cumulative values of the portfolios in (33).

Table 6 shows the cost function in (34) and the number of days that the losses exceed the CVaR (or the number of exceedance) for models (1) to (3) where the portfolio weights are obtained by respectively solving the MCVaR portfolio and MV portfolio problems. First note that the cost function and the number of exceedance for the MV portfolios are all smaller than those for the MCVaR portfolios. This is not because the MV portfolios perform better, but rather due to the fact that the CVaRs are not minimized in the MV problem, and thus larger CVaRs give smaller cost function and number of exceedance. We should only compare the cost function and the number of exceedance between the GC-GARCH, MA and DCC for the same MCVaR or MV portfolios with the same $\alpha$ . We note that the cost functions of the portfolios with GC-GARCH and MA are all smaller than those of the DCC-GARCH model. The number of exceedance in the portfolios with GC-GARCH and MA are all smaller than those in DCC except for $\alpha=0.05$ . Furthermore, the MA has fewer days of exceedance than the GC-GARCH for MCVaR portfolio at $\alpha=0.01$ and MV at $\alpha=0.005,0.01$ . This suggests that the GC-GARCH model and MA work much better than the DCC-GARCH model in covering the losses using CVaR, and the model averaging using GC-GARCH can further improve the performance over the GC-GARCH using the MAP network.

We now compare the investment performance of strategy 1, where we rebalance the portfolio on the first day of every week. Figure 15 shows the time series of the cumulative values of the weekly rebalanced portfolios according to (33). The blue curves indicate the portfolios with the GC-GARCH model using the MAP network, the orange curves indicate the portfolios with the GC-GARCH model using model averaging, and the black curves indicate the portfolios obtained from the DCC-GARCH model. Figure 15 (a) shows the cumulative values of the minimum variance portfolios. The performances of GC-GARCH, MA and DCC are similar in this case. Figures 15 (b) to (f) show the minimum CVaR portfolios with $\alpha=0.0005,0.001,0.005,0.01,$ and $0.05$ . The values of the portfolios using GC-GARCH and MA are on average higher than those using DCC. The portfolios generate higher profits for smaller $\alpha$ . This supports that the GC-GARCH model outperforms the DCC-GARCH model when we take the tails into account. However, we observe that the cumulative values of the portfolios are falling from 2021 to 2022 in all portfolios in Figure 15. Even though the CVaRs have been minimized in the experiment, some of the losses could still be significant, and we need to modify the strategy to avoid huge potential losses.

Strategy 2 tries to avoid investing in the stock market whenever the predicted 1-day ahead CVaR is higher than the previous $w_{S}$ weeks average CVaR in (35). Figure 16, Figure 17, and Figure 18 show the cumulative values of the minimum CVaR portfolios with respectively $w_{S}=8,16,$ and $32$ . The cumulative values of the portfolios with GC-GARCH with the MAP network (blue curves) and GC-GARCH with model averaging (orange curves) are on average growing from 2020 to 2023, whereas the cumulative values of the DCC-GARCH portfolios (black curves) tend to be falling. To investigate in more detail, Table 7 shows the numbers of weeks invested in the stock market and the average returns on the weeks excluded from investing in the stock market for the investment experiment using strategy 2. The number of weeks is similar across different $\alpha$ and $w_{S}$ using different models. The average returns on the weeks excluded are all negative for GC-GARCH and MA, and are much smaller than those in DCC-GARCH. This suggests that the strategy with the GC-GARCH model successfully avoids days with high losses, and the portfolios are able to grow when we take the tails into account. The GC-GARCH model has good ability to capture the extreme scenarios from the tails. Strategy 2 performs well for different $w_{S}=8,16,$ and $32$ , especially when we take a small $w_{S}=8$ , in which case the average returns on the weeks excluded from investing in the stock markets are ranged from $-0.21$ % to $-0.16$ %, as shown in Table 7.

		Cost function			Number of CVaR exceedance
	$\alpha$	GC-GARCH	MA	DCC	GC-GARCH	MA	DCC
MCVaR portfolio	0.0005	0.0445	0.0454	0.0582	10	10	13
	0.001	0.0464	0.0451	0.0589	9	9	13
	0.005	0.0449	0.0502	0.0535	11	11	15
	0.01	0.0532	0.0533	0.0545	15	12	15
	0.05	0.0883	0.0765	0.0929	26	26	22
MV portfolio	0.0005	0.0026	0.0094	0.0158	2	2	6
	0.001	0.0109	0.0172	0.0273	3	3	10
	0.005	0.0319	0.0412	0.0434	9	8	14
	0.01	0.0450	0.0487	0.0493	13	11	14
	0.05	0.0864	0.0740	0.0913	26	25	22

Table 6: The cost function and the number of days that the losses exceed the CVaR for (1) GC-GARCH model using the MAP network, (2) GC-GARCH model with model averaging (marked MA in the table), and (3) DCC-GARCH model. The portfolio weights are obtained by respectively solving the MCVaR and MV problem.

		Number of weeks invested in the stock market			Average return excluded from investment (%)
	$\alpha$	GC-GARCH	MA	DCC	GC-GARCH	MA	DCC
$w_{S}=8$	0.0005	79	79	78	-0.21	-0.16	0.05
	0.001	79	79	78	-0.21	-0.15	0.05
	0.005	78	78	74	-0.16	-0.12	0.00
	0.01	79	79	74	-0.16	-0.12	0.02
	0.05	80	80	75	-0.17	-0.14	0.03
$w_{S}=16$	0.0005	80	80	72	-0.15	-0.14	0.00
	0.001	80	80	73	-0.14	-0.14	-0.02
	0.005	82	82	75	-0.12	-0.10	-0.05
	0.01	85	85	72	-0.11	-0.10	-0.01
	0.05	79	79	73	-0.12	-0.11	0.00
$w_{S}=32$	0.0005	89	89	76	-0.10	-0.09	0.08
	0.001	89	89	77	-0.09	-0.09	0.08
	0.005	87	87	78	-0.13	-0.08	0.05
	0.01	90	90	75	-0.10	-0.08	0.06
	0.05	87	87	74	-0.09	-0.07	0.07

Table 7: The number of weeks invested in the stock market and the average returns on the weeks excluded from investing in the stock market for the investment experiment using strategy 2.

8 Discussion and conclusion

We propose the GC-GARCH model to reduce the number of correlation parameters using the idea motivated from the CAPM theory, which uses the risk factors (market indexes) to explain the co-movements among a portfolio of stocks. The GC-GARCH model also captures the co-movements in the tails among risk factor and stock returns through the use of pair-copula construction. The risk factors are factorized using the Bayesian network, which allows us to “rank” the stocks from topological orders in BN. The edges in the BN indicate the flow of information. The time-varying specification of the parameters further flexibilizes the modeling. We adopt the three-stage estimation procedure, where we first estimate the parameters in the marginal distribution, then in the DAG copulas and finally the stock copulas. We have shown that the estimation procedures provide accurate estimates in Section 6. The moving-window investment experiment in Section 7 show that the GC-GARCH models using the highest scoring network and model averaging give better performances than the traditional DCC-GARCH model. The GC-GARCH has better ability to cover extreme losses, to avoid risk and to make profits from investing in the market.

There are several things we have not addressed in this paper. We assume that the underlying network in the GC-GARCH model is fixed. However, the underlying networks in financial markets are likely to be changing dynamically. We use the $t$ -copula, which is an elliptical copula, in the GC-GARCH model. However, the tail dependence of the returns could be asymmetric. We could employ the GC-GARCH model using other families of copulas, such as the Gumbel copula and the Clayton copula. We can address these points in future work.

Acknowledgement

This work was supported by the Hong Kong RGC General Research Fund (grant number 16507322).

\printbibliography

Appendix A An algorithm for subsetting the DAG space

We provide an algorithm to subset the DAG space to reduce computation burden. Let $\pi_{1}=\{r_{i[1]t},\ldots,r_{i[k]t}\}$ . We define $\pi_{s}^{+}$ ( $s=1,2,\ldots$ ) to be the last element in $\pi_{s}$ , and $\pi_{s}^{-}$ to be the set obtained by removing the last element in $\pi_{s}$ . Then, we have $\pi_{s}=\pi_{s}^{-}\cup\pi_{s}^{+}$ . We define a set of cumulative parents of a variable $r_{it}$ to be a set containing $\{r_{i[1],t},\ldots,r_{i[k],t}\}$ for some $k$ . For example, consider the variable $r_{5t}$ in the network in Figure 19, when we use the topological order $(1,2,3,4,5)$ for pair-copula construction, the sets $\{r_{1,t}\}$ , $\{r_{1,t},r_{3,t}\}$ and $\{r_{1,t},r_{3,t},r_{4,t}\}$ are considered as sets of cumulative parents, whereas the sets $\{r_{3,t}\}$ , $\{r_{1,t},r_{4,t}\}$ , $\{r_{3,t},r_{4,t}\}$ are examples that are not sets of cumulative parents.

To check whether the term $\pi_{s}^{-}$ contains a cumulative set of parents of $\pi_{s}^{+}$ , we do the following test, and we include the DAG in the reduced space when all conditional distributions in the network pass the test.

Step 1.

Take $\tilde{\pi}_{s}=\mathbf{r}_{\text{pa}(\pi_{s}^{+})}\cap\pi_{s}^{-}$ , where $\mathbf{r}_{\text{pa}(\pi_{s}^{+})}$ is the set containing the variables that are parents of $\pi_{s}^{+}$ .
Step 2.
- •
  
  If $\tilde{\pi}_{s}=\varnothing$ , the test is passed.
- •
  
  If $\tilde{\pi}_{s}\neq\varnothing$ , and if $\tilde{\pi}_{s}$ contains a set of cumulative parents of $\pi_{s}^{+}$ , then we take $\pi_{s+1}=\tilde{\pi}_{s}$ and update $s$ to $s+1$ . Go back to Step 1.
- •
  
  If $\tilde{\pi}_{s}\neq\varnothing$ , and $\tilde{\pi}_{s}$ does not contain a set of cumulative parents of $\pi_{s}^{+}$ , then we conclude that the calculation of $F^{[t]}(r_{i[k],t}|r_{i[1],t},\ldots,r_{i[k-1],t})$ requires numerical integration and the test is failed.

We reduce the DAG space by considering only those networks passing the above test. Even though there exist some cases where we still need to calculate conditional distributions using numerical integration, the possibility of involving numerical integration in the likelihood calculation has been greatly reduced.

For example, consider the DAG in Figure 19, which involves the copulas

c^{[t]}_{2,1}(\cdot,\cdot),c^{[t]}_{4,1}(\cdot,\cdot),c^{[t]}_{4,3|1}(\cdot,% \cdot),c^{[t]}_{5,1}(\cdot,\cdot),c^{[t]}_{5,3|1}(\cdot,\cdot),c^{[t]}_{5,4|3,% 1}(\cdot,\cdot),

we first try to check, directly without using the above algorithm, if we can compute a copula density directly without using integration. We consider the second conditional distribution in
$c_{5,4|3,1}^{[t]}(F^{[t]}(r_{5,t}|r_{3,t},r_{1,t}),F^{[t]}(r_{4,t}|r_{3,t},r_{% 1,t}))$ :

F^{[t]}(r_{4t}|r_{3,t},r_{1,t})=h_{4,3|1}^{[t]}(F^{[t]}(r_{4,t}|r_{1,t}),F^{[t% ]}(r_{3,t}|r_{1,t})).

Note that $F^{[t]}(r_{4,t}|r_{1,t})=h_{4,1}^{[t]}(F^{[t]}(r_{4,t}),F^{[t]}(r_{1,t}))$ can be directly computed using the $h$ function, and $F^{[t]}(r_{3,t}|r_{1,t})=F^{[t]}(r_{3,t})$ can be directly computed from the marginal distribution as $r_{3,t}$ and $r_{1,t}$ are independent.

If we use the algorithm above, the copula $c_{5,4|3,1}^{[t]}(F^{[t]}(r_{5,t}|r_{3,t},r_{1,t}),F^{[t]}(r_{4,t}|r_{3,t},r_{% 1,t}))$ corresponds to $i=5$ and $k=3$ . Then, $\pi_{1}=\{r_{5[1],t},r_{5[2],t},r_{5[3],t}\}=\{r_{1,t},r_{3,t},r_{4,t}\}$ , $\pi_{1}^{+}=\{r_{4,t}\}$ , and $\pi_{1}^{-}=\{r_{1,t},r_{3,t}\}$ . The second conditional distribution in the copula can be written as $F^{[t]}(\pi_{1}^{+}|\pi_{1}^{-})$ .

1.

For $s=1$ , we take $\tilde{\pi}_{1}=\mathbf{r}_{\text{pa}(\pi_{1}^{+})}\cap\pi_{1}^{-}=\mathbf{r}_% {\text{pa}(r_{4,t})}\cap\{r_{1,t},r_{3,t}\}=\{r_{1,t},r_{3,t}\}\cap\{r_{1,t},r% _{3,t}\}=\{r_{1,t},r_{3,t}\}\neq\varnothing$ . Since $\tilde{\pi}_{1}$ is cumulative parents of $\pi_{1}^{+}=r_{4,t}$ , we take $\pi_{2}=\tilde{\pi}_{1}$ and update $s=2$ .
2.

Now, $\pi_{2}=\{r_{1,t},r_{3,t}\}$ , $\pi_{2}^{+}=\{r_{3,t}\}$ and $\pi_{2}^{-}=\{r_{1,t}\}$ . Then, $\tilde{\pi}_{2}=\mathbf{r}_{\text{pa}(r_{3t})}\cap\{r_{1,t}\}=\varnothing\cap% \{r_{1,t}\}=\varnothing$ . Then, $F^{[t]}(r_{4,t}|r_{3,t},r_{1,t})$ passes the test.

We also give a counterexample. Consider the network in Figure 2. Consider the copula $c_{5,4|3,1}^{[t]}(F^{[t]}(r_{5t}|r_{3t},r_{1t}),$ $F^{[t]}(r_{4t}|r_{3t},r_{1t}))$ . We focus on the second argument $F^{[t]}(r_{4t}|r_{3t},r_{1t}))$ . Using the algorithm, the copula corresponds to $i=5$ and $k=3$ . Then, $\pi_{1}=\{r_{5[1],t},r_{5[2],t},r_{5[3],t}\}=\{r_{1,t},r_{3,t},r_{4,t}\}$ , $\pi_{1}^{+}=\{r_{4,t}\}$ , and $\pi_{1}^{-}=\{r_{1t},r_{3t}\}$ . For $s=1$ , we take $\tilde{\pi}_{1}=\mathbf{r}_{\text{pa}(\pi_{1}^{+})}\cap\pi_{1}^{-}=\{r_{2,t},r% _{3,t}\}\cap\{r_{1,t},r_{3,t}\}=\{r_{3,t}\}\neq\varnothing$ . However, the parents of $r_{4t}$ are $\{r_{2,t},r_{3,t}\}$ . Then, $\tilde{\pi}_{1}$ is not a set of cumulative parents of $\pi_{1}^{+}$ , and thus we exclude the DAG in the reduced space. This part corresponds to (20), where we need to compute the integration to marginalize out $r_{2t}$ . Without a set of cumulative parents, the computation must require the use of integration.

Figure 19: An example DAG that passes the test.

Appendix B Parameter estimation results in simulation (S2)

B.1 Parameters in the DAG copulas

The estimation results for simulation (S2) are presented in Table 1 for $a_{xy|z}$ and $b_{xy|z}$ , and Table 2 for $\bar{\varphi}_{xy|z}$ and $v_{xy|z}$ . Most of the true values are contained in the 95% credible intervals. The estimates are also close to the true values for most of the parameters. The performances are similar to the case in simulation (S1). It can be concluded that the MCMC estimation method in Section 4.2 is reliable.

Parameter	True value	Mean of estimates	Mean MAE	Median of estimate	Median MAE	$q_{0.05}$	$q_{0.95}$
$a_{10,3}$	$\mathbf{0.05^{*}}$	0.0439	0.0221	0.043	0.018	0.0004	0.0797
$a_{10,9\|3}$	$\mathbf{0.03^{*}}$	0.0123	0.0244	0.0011	0.0276	0	0.0743
$a_{2,1}$	$\mathbf{0.04^{*}}$	0.0417	0.0217	0.0406	0.017	0.0052	0.0791
$a_{3,1}$	$\mathbf{0.06^{*}}$	0.0609	0.0202	0.0587	0.0154	0.0295	0.0933
$a_{3,2\|1}$	$\mathbf{0.02^{*}}$	0.0245	0.029	0.0161	0.0215	0	0.0805
$a_{4,1}$	$\mathbf{0.03^{*}}$	0.0275	0.0283	0.0207	0.0227	0	0.0732
$a_{4,2\|1}$	$\mathbf{0.03^{*}}$	0.032	0.0343	0.0265	0.0229	0	0.1115
$a_{5,1}$	$\mathbf{0.06^{*}}$	0.0573	0.0284	0.0567	0.0222	0.0123	0.1043
$a_{5,3\|1}$	$\mathbf{0.03^{*}}$	0.026	0.028	0.0185	0.0224	0.0001	0.0827
$a_{5,4\|1,3}$	$\mathbf{0.09^{*}}$	0.0817	0.0432	0.0826	0.0344	0.0138	0.1538
$a_{6,5}$	$\mathbf{0.08^{*}}$	0.0596	0.0466	0.0541	0.0415	0.0002	0.1512
$a_{7,1}$	$\mathbf{0.09^{*}}$	0.0845	0.0209	0.0831	0.0169	0.0538	0.1188
$a_{7,2\|1}$	$\mathbf{0.05^{*}}$	0.051	0.0382	0.0461	0.0282	0.0007	0.1096
$a_{7,5\|1,2}$	$\mathbf{0.09^{*}}$	0.0845	0.0376	0.0836	0.0278	0.0314	0.1393
$a_{8,1}$	$\mathbf{0.07^{*}}$	0.072	0.0297	0.069	0.0223	0.0345	0.1177
$a_{8,3\|1}$	$\mathbf{0.06^{*}}$	0.0582	0.0338	0.0577	0.0253	0.0015	0.1178
$a_{9,1}$	$\mathbf{0.07^{*}}$	0.0685	0.0225	0.0686	0.0173	0.0339	0.1062
$a_{9,4\|1}$	$\mathbf{0.02^{*}}$	0.0143	0.018	0.0086	0.0155	0	0.0465
$a_{9,7\|1,4}$	$\mathbf{0.02^{*}}$	0.0199	0.025	0.0128	0.0172	0	0.074
$b_{10,3}$	$\mathbf{0.92^{*}}$	0.8629	0.1962	0.9196	0.0812	0.4181	0.9614
$b_{10,9\|3}$	$\mathbf{0.9^{*}}$	0.5729	0.4059	0.7276	0.3801	0.0001	0.9993
$b_{2,1}$	$\mathbf{0.86^{*}}$	0.7981	0.185	0.842	0.1066	0.3114	0.9503
$b_{3,1}$	$\mathbf{0.91^{*}}$	0.8974	0.0681	0.9123	0.0339	0.8237	0.9508
$b_{3,2\|1}$	$\mathbf{0.95^{*}}$	0.7064	0.3385	0.8519	0.2598	0.001	0.9946
$b_{4,1}$	$\mathbf{0.92^{*}}$	0.7247	0.3232	0.8879	0.2292	0.0243	0.9953
$b_{4,2\|1}$	$\mathbf{0.95^{*}}$	0.7817	0.3202	0.9366	0.1857	0.0208	0.9933
$b_{5,1}$	$\mathbf{0.86^{*}}$	0.7725	0.2147	0.8378	0.1238	0.1073	0.9369
$b_{5,3\|1}$	$\mathbf{0.94^{*}}$	0.7385	0.324	0.9026	0.224	0.0028	0.9947
$b_{5,4\|1,3}$	$\mathbf{0.81^{*}}$	0.7245	0.2275	0.7883	0.1429	0.0438	0.9443
$b_{6,5}$	$\mathbf{0.84^{*}}$	0.6638	0.3074	0.7869	0.2269	0.004	0.9621
$b_{7,1}$	$\mathbf{0.86^{*}}$	0.8457	0.0742	0.8549	0.0376	0.7567	0.9115
$b_{7,2\|1}$	$\mathbf{0.86^{*}}$	0.7247	0.2877	0.8401	0.1923	0.0195	0.9834
$b_{7,5\|1,2}$	$\mathbf{0.84^{*}}$	0.786	0.1899	0.8399	0.1043	0.3452	0.9407
$b_{8,1}$	$\mathbf{0.86^{*}}$	0.8085	0.175	0.857	0.0863	0.4541	0.9301
$b_{8,3\|1}$	$\mathbf{0.86^{*}}$	0.7497	0.2502	0.8415	0.1544	0.0419	0.9521
$b_{9,1}$	$\mathbf{0.84^{*}}$	0.8152	0.1225	0.8359	0.0649	0.679	0.924
$b_{9,4\|1}$	$\mathbf{0.93^{*}}$	0.7019	0.3498	0.8918	0.2601	0.0124	0.9952
$b_{9,7\|1,4}$	$\mathbf{0.96^{*}}$	0.7551	0.3261	0.9283	0.2201	0.0124	0.996

Table 1: The true values, the mean, mean absolute error (MAE) of the mean, the median, the MAE of the median, the 5th quantile (

q_{0.05}

), and the 95th quantile (

q_{0.95}

) of the posterior median estimates of

a_{xy|z}

and

b_{xy|z}

out of the 200 replications for simulation (S2).

Parameter	True value	Mean of estimates	Mean MAE	Median of estimate	Median MAE	$q_{0.05}$	$q_{0.95}$
$\bar{\varphi}_{10,3}$	$\mathbf{0.24^{*}}$	0.2831	0.1249	0.2781	0.0965	0.0987	0.4595
$\bar{\varphi}_{10,9\|3}$	$-0.41$	$-$ 0.2534	0.0469	$-$ 0.256	0.1573	$-$ 0.327	$-$ 0.186
$\bar{\varphi}_{2,1}$	$\mathbf{-0.78^{*}}$	$-$ 0.7702	0.0283	$-$ 0.7711	0.0232	$-$ 0.8198	$-$ 0.7228
$\bar{\varphi}_{3,1}$	$\mathbf{-0.36^{*}}$	$-$ 0.355	0.1126	$-$ 0.3551	0.0847	$-$ 0.5275	$-$ 0.1755
$\bar{\varphi}_{3,2\|1}$	$\mathbf{-0.02^{*}}$	$-$ 0.0202	0.0807	$-$ 0.0187	0.06	$-$ 0.1422	0.0986
$\bar{\varphi}_{4,1}$	$\mathbf{0.26^{*}}$	0.2593	0.0645	0.2602	0.0495	0.1531	0.3621
$\bar{\varphi}_{4,2\|1}$	$\mathbf{0.02^{*}}$	0.0185	0.1042	0.0203	0.074	$-$ 0.1257	0.1474
$\bar{\varphi}_{5,1}$	$\mathbf{0.54^{*}}$	0.5194	0.0548	0.5168	0.0473	0.4325	0.6079
$\bar{\varphi}_{5,3\|1}$	$\mathbf{0.12^{*}}$	0.1144	0.0722	0.1122	0.0559	0.0094	0.2255
$\bar{\varphi}_{5,4\|1,3}$	$\mathbf{-0.26^{*}}$	$-$ 0.2428	0.0687	$-$ 0.2419	0.0558	$-$ 0.3642	$-$ 0.1406
$\bar{\varphi}_{6,5}$	$\mathbf{0.29^{*}}$	0.195	0.0811	0.1994	0.1055	0.0771	0.3028
$\bar{\varphi}_{7,1}$	$\mathbf{-0.47^{*}}$	$-$ 0.4649	0.0742	$-$ 0.468	0.0598	$-$ 0.5728	$-$ 0.3435
$\bar{\varphi}_{7,2\|1}$	$\mathbf{-0.42^{*}}$	$-$ 0.4121	0.0632	$-$ 0.4126	0.0489	$-$ 0.5067	$-$ 0.3083
$\bar{\varphi}_{7,5\|1,2}$	$\mathbf{0.03^{*}}$	0.0127	0.1103	0.0185	0.0777	$-$ 0.1545	0.1786
$\bar{\varphi}_{8,1}$	$\mathbf{0.43^{*}}$	0.4018	0.0776	0.4074	0.0638	0.2633	0.532
$\bar{\varphi}_{8,3\|1}$	$\mathbf{0.5^{*}}$	0.489	0.0651	0.4921	0.0527	0.3823	0.6006
$\bar{\varphi}_{9,1}$	$\mathbf{-0.69^{*}}$	$-$ 0.6862	0.0413	$-$ 0.6894	0.0327	$-$ 0.7534	$-$ 0.6112
$\bar{\varphi}_{9,4\|1}$	$\mathbf{0.2^{*}}$	0.1871	0.0506	0.1869	0.04	0.1082	0.2713
$\bar{\varphi}_{9,7\|1,4}$	$0.74$	0.5693	0.0704	0.5564	0.1725	0.4759	0.7001
$v_{10,3}$	$\mathbf{7.9^{*}}$	9.9883	3.1891	9.6558	2.8627	5.7493	15.934
$v_{10,9\|3}$	$6.67$	14.649	3.764	14.3274	8.0097	8.9711	21.157
$v_{2,1}$	$\mathbf{7.97^{*}}$	7.9942	2.2888	7.7657	1.8153	4.9135	12.4083
$v_{3,1}$	$\mathbf{7.91^{*}}$	8.6289	2.7028	8.0335	1.9919	5.6216	14.5503
$v_{3,2\|1}$	$\mathbf{6.19^{*}}$	6.6826	2.0876	6.2403	1.5549	4.1974	10.8303
$v_{4,1}$	$\mathbf{7.27^{*}}$	8.0482	2.4939	7.5606	1.8723	5.0499	12.1542
$v_{4,2\|1}$	$\mathbf{9.98^{*}}$	9.9089	3.3581	9.082	2.6316	5.7849	16.4214
$v_{5,1}$	$\mathbf{5.41^{*}}$	5.8013	1.5465	5.4467	1.1149	3.8427	8.8646
$v_{5,3\|1}$	$\mathbf{5.26^{*}}$	6.0324	1.6404	5.6626	1.3035	4.0945	8.8652
$v_{5,4\|1,3}$	$\mathbf{7.47^{*}}$	8.14	2.2743	7.8704	1.7779	5.0633	13.045
$v_{6,5}$	$\mathbf{5.27^{*}}$	8.34	2.4175	7.6211	3.1237	5.2655	12.8108
$v_{7,1}$	$\mathbf{9.75^{*}}$	10.0117	3.058	9.1846	2.4324	6.2343	15.6522
$v_{7,2\|1}$	$\mathbf{8.17^{*}}$	8.7747	2.8075	8.2079	2.0742	5.1824	13.6875
$v_{7,5\|1,2}$	$\mathbf{6.35^{*}}$	8.0564	2.6115	7.6028	2.271	4.7868	12.7598
$v_{8,1}$	$\mathbf{9.98^{*}}$	10.0132	3.3093	9.4096	2.5236	6.014	16.2632
$v_{8,3\|1}$	$\mathbf{6.79^{*}}$	7.6349	2.6396	7.0186	1.9867	4.4484	12.7467
$v_{9,1}$	$\mathbf{6.86^{*}}$	7.1971	1.9426	6.7725	1.5081	4.8446	10.8707
$v_{9,4\|1}$	$\mathbf{6.65^{*}}$	7.7157	2.4767	7.2372	1.9109	4.5809	12.0884
$v_{9,7\|1,4}$	$\mathbf{5.52^{*}}$	8.6611	2.6835	8.2426	3.2541	4.9682	13.7083

Table 2: The true values, the mean, mean absolute error (MAE) of the mean, the median, the MAE of the median, the 5th quantile (

q_{0.05}

), and the 95th quantile (

q_{0.95}

) of the posterior median estimates of

v_{xy|z}

and

\bar{\varphi}_{xy|z}

out of the 200 replications for simulation (S2).

B.2 Parameters in the stock copula

This section gives the sequential estimation results for simulation (S2). Figure 1, Figure 2, Figure 3 and Figure 4 are four graphs summarizing the parameter estimation for $a_{xy|z}$ , $b_{xy|z}$ , $v_{xy|}z$ , and $\bar{\varphi}_{xy|z}$ . These graphs are similar to the graphs in Figure 6, Figure 7, Figure 8 and Figure 9 presented in Section 6.3.2, and the observations regarding the estimation results are similar to those in Section 6.3.2.

Appendix C Comparing the sequential estimator and the posterior mean/median estimators

In our empirical study, we conduct a moving-window study with a total of 173 windows, in each of which we include 92 stocks. Estimating stock copulas using MCMC sampling is infeasible in this case. For $m=10$ , $p=92$ and $T=750$ , for each stock, it took 6 hours to run 20,000 iterations in the MCMC sampling. The sequential estimation, however, only takes 30 seconds.

In this section, we compare the sequential estimates to the posterior mean/median estimates. We conduct 20 replications with different data sets generated from the GC-GARCH model with the same underlying network in Figure 4 as in simulation (S2). The true parameters are also fixed across replications. We calculate the MAEs of the estimates of all 100 stocks, by level. Figure 1 shows the time series plots of the average MAE, within the same level, of the sequential estimates and the posterior mean/mode/median estimates of (a) $\bar{\varphi}_{xy|z}$ , (b) $a_{xy|z}$ , (c) $b_{xy|z}$ , and (d) $v_{xy|z}$ . The x-axes are the levels.

•

In (a), we observe that the MAEs of the sequential estimates of $\bar{\varphi}_{xy|z}$ (dashed blue line) are similar to the MAEs of the mean and median MCMC estimates (respectively the solid red and dotted green lines).
•

In (b), the MAEs for $a_{xy|z}$ using sequential estimation are smaller than the MAEs of the MCMC estimates in most of the levels.
•

In (c), we observe that the sequential estimates of $b_{xy|z}$ have slightly larger MAEs for the first four levels, and the sequential estimates actually work better for the subsequent levels.
•

In (d), the MAEs of the sequential estimates and the MCMC median estimates of $v_{xy|z}$ are similar.

Summarizing the above four points, we observe that sequential estimates work as well as the MCMC estimates.

Appendix D List of stocks

Name	Stock symbol	Name	Stock symbol
CKH HOLDINGS	0001.HK	CSPC PHARMA	1093.HK
CLP HOLDINGS	0002.HK	CHINA RES LAND	1109.HK
HK & CHINA GAS	0003.HK	CK ASSET	1113.HK
HSBC HOLDINGS	0005.HK	YANKUANG ENERGY	1171.HK
POWER ASSETS	0006.HK	SINO BIOPHARM	1177.HK
HANG SENG BANK	0011.HK	CHINA RAIL CONS	1186.HK
HENDERSON LAND	0012.HK	BYD COMPANY	1211.HK
SHK PPT	0016.HK	ABC	1288.HK
NEW WORLD DEV	0017.HK	AIA	1299.HK
GALAXY ENT	0027.HK	PICC GROUP	1339.HK
MTR CORPORATION	0066.HK	ICBC	1398.HK
HANG LUNG PPT	0101.HK	PSBC	1658.HK
TSINGTAO BREW	0168.HK	CRRC	1766.HK
GEELY AUTO	0175.HK	GF SEC	1776.HK
ALI HEALTH	0241.HK	CHINA COMM CONS	1800.HK
CITIC	0267.HK	PRADA	1913.HK
WH GROUP	0288.HK	COSCO SHIP HOLD	1919.HK
CHINA RES BEER	0291.HK	SANDS CHINA LTD	1928.HK
TINGYI	0322.HK	CHOW TAI FOOK	1929.HK
SINOPEC CORP	0386.HK	MINSHENG BANK	1988.HK
HKEX	0388.HK	COUNTRY GARDEN	2007.HK
CHINA RAILWAY	0390.HK	ANTA SPORTS	2020.HK
TECHTRONIC IND	0669.HK	CHINA VANKE	2202.HK
CHINA OVERSEAS	0688.HK	SHENZHOU INTL	2313.HK
TENCENT	0700.HK	PING AN	2318.HK
CHINA TELECOM	0728.HK	MENGNIU DAIRY	2319.HK
AIR CHINA	0753.HK	PICC P&C	2328.HK
CHINA UNICOM	0762.HK	LI NING	2331.HK
LINK REIT	0823.HK	GREATWALL MOTOR	2333.HK
PETROCHINA	0857.HK	PRU	2378.HK
XINYI GLASS	0868.HK	SUNNY OPTICAL	2382.HK
ZHONGSHENG HLDG	0881.HK	BOC HONG KONG	2388.HK
CNOOC	0883.HK	CPIC	2601.HK
HUANENG POWER	0902.HK	CHINA LIFE	2628.HK
CONCH CEMENT	0914.HK	ENN ENERGY	2688.HK
CHINA LONGYUAN	0916.HK	STANCHART	2888.HK
CCB	0939.HK	ZIJIN MINING	2899.HK
CHINA MOBILE	0941.HK	BANKCOMM	3328.HK
LONGFOR GROUP	0960.HK	CICC	3908.HK
XINYI SOLAR	0968.HK	CM BANK	3968.HK
SMIC	0981.HK	BANK OF CHINA	3988.HK
CITIC BANK	0998.HK	CMOC	3993.HK
CKI HOLDINGS	1038.HK	CITIC SEC	6030.HK
HENGAN INT’L	1044.HK	CMSC	6099.HK
CHINA SOUTH AIR	1055.HK	CEB BANK	6818.HK
CHINA SHENHUA	1088.HK	HTSC	6886.HK

Table 1: The

p=92

stocks included in the portfolio.

	$\displaystyle f^{[t]}(r_{5,t}\|r_{1,t},r_{3,t},r_{4,t})=$	$\displaystyle c_{5,1}^{[t]}(F^{[t]}(r_{5,t}),F^{[t]}(r_{1,t}))c_{5,3\|1}^{[t]}(% F^{[t]}(r_{5,t}\|r_{1,t}),F^{[t]}(r_{3,t}\|r_{1,t}))\cdot$		(17)
		$\displaystyle c_{5,4\|3,1}^{[t]}(F^{[t]}(r_{5,t}\|r_{3,t},r_{1,t}),F^{[t]}(r_{4,% t}\|r_{3,t},r_{1,t})).$		(17)

	$\displaystyle F^{[t]}(r_{4t}\|r_{3t})$	$\displaystyle=\int_{-\infty}^{\infty}F^{[t]}(r_{4t}\|r_{3t},r_{2t})f^{[t]}(r_{2% t})dr_{2t}$		(20)
		$\displaystyle\equiv\int_{-\infty}^{\infty}h_{4,3\|2}^{[t]}(F^{[t]}(r_{4t}\|r_{2t% }),F^{[t]}(r_{3t}\|r_{2t}))f^{[t]}(r_{2t})dr_{2t},$		(20)