Abstract
Financial data usually have the features of complexity and interdependence structure, such as asymmetric, tail, and time-varying dependence. This study constructs a new multivariate skewed fat-tailed copula, namely, noncentral contaminated normal (NCCN) copula, to analyze the dependent structure of financial market data. The dynamic conditional correlation (DCC) model is also incorporated into constructing the time-varying NCCN copula model. This study comprehensively examines the effects of the DCC-NCCN copula and related models on fitting dependence structures of Hong Kong stock markets. The results show that the DCC-NCCN copula model can better depict the dependence structures of returns. Considering the flexibility and complexity, the DCC-NCCN copula model is a relatively ideal, time-varying, multivariate skewed fat-tailed copula model.
1. Introduction
After suffering from loss in the stock market, Isaac Newton ever said that “I can calculate the motions of the heavenly bodies, but not the madness of people.” This reflects the complexity of financial markets. In general, the financial asset return series have relatively complex interdependence structural features, such as asymmetric dependence, tail dependence, and time-varying dependence. According to whether it can depict asymmetric dependence and fat-tailed dependence, copula can be divided into four categories: symmetric thin-tailed copula, symmetric fat-tailed copula, skewed thin-tailed copula, and skewed fat-tailed copula. The examples above are normal copula, t-copula, skew-normal copula, and skew-t-copula. The multivariate skew-normal copula is the copula of the multivariate skew-normal distribution, such as Wei et al. [1] proposed the copula of the multivariate skew-normal distribution of Azzalini and Valle [2]. The multivariate skew-t-copula is the copula of the multivariate skew-t distribution, such as Demarta and McNeil [3] proposed the copula of the multivariate generalized hyperbolic skew-t (GHST) distribution of Barndorff-Nielsen [4]. Kollo and Pettere [5] propose the copula of the multivariate skew-t distribution of Azzalini and Capitanio [6]; Smith et al. [7] put forth the copula of the multivariate skew-t distribution of Sahu et al. [8]; and Liu et al. [9] advanced the copula of the multivariate extended skew-t (EST) distribution by Arellano-Valle and Genton [10].
Although these multivariate skew-t copulas are very flexible, they are also highly complex and challenging to apply. Considering both flexibility and complexity, these multivariate skew-t copulas may not be very ideal. This study constructs a new multivariate skewed fat-tailed distribution, namely, the multivariate noncentral contaminated normal (NCCN) distribution. The multivariate NCCN distribution can be interpreted as a multivariate noncentral normal scale mixture distribution, which is similar to the multivariate normal variance-mean mixture distribution and multivariate skew-normal scale mixture distribution. The multivariate NCCN distribution can also be interpreted as a simplified mixture of two multivariate normal distributions. Then, the copula of the multivariate NCCN distribution can be called the multivariate NCCN copula. Note that the multivariate NCCN copula cannot be interpreted as a mixture of two multivariate normal copulas. The NCCN copula may be relatively ideal. The advantages are shown as follows. First, the NCCN copula can flexibly describe positive and negative dependence. Second, according to Mardia’s kurtosis, the NCCN copula has stronger tail dependence than the normal copula. Third, the NCCN copula can flexibly describe asymmetric dependence. Fourth, the subclasses of the NCCN copula include normal copula and contaminated normal (CN) copula. Note that the flexibility of the CN copula is similar to the t-copula, but the complexity of the CN copula is significantly lower than that of the t-copula. Fifth, NCCN copula is suitable for two- and higher-dimensional dependence structure modeling. Sixth and the last, the flexibility of the NCCN copula is similar to that of skew-t copulas, but the complexity of the NCCN copula is significantly lower than that of skew-t copulas.
According to whether it can delineate the time-varying dependence, the copula can be divided into two classes: static copula and dynamic one. There are many options in modeling dynamic structures, including the time-varying parameter model [11], the dynamic conditional correlation (DCC) model [12] and the dynamic condition-related improvement (DCC-Student-t) model [13], the time-varying correlation (TVC) model [14], the asymmetric DCC (ADCC) model [15,16], and the generalized autoregressive score (GAS) model [17]. These dynamic models have pros and cons, and we compare them as follows. The advantage of the Patton model is that it does not limit the type of time-varying parameters, and the dynamic structure is relatively simple; the disadvantage is that the dimension is limited to two dimensions, and the meaning to interpret the dynamic structure is not very clear. The advantages of the DCC model and TVC model are as follows: dimension is unlimited, the dynamic structure is simple, and the interpretation meaning of the dynamic structure is clear; the disadvantage is to limit the type of time-varying parameters to the linear correlation matrix. The advantage of the ADCC model is that it further considers asymmetric dynamics based on the DCC model. For the GAS model, the advantages are unlimited dimension and time-varying parameter type. The interpretation of the dynamic structure is relatively definite. The disadvantage is that the dynamic structure is generally quite complex. There are only a few studies on the dynamic multivariate skewed fat-tailed copula, mainly including the dynamic asymmetric copula (DAC) model given by Christoffersen et al. [18], GAS-GHST copula model [19], and dynamic double asymmetric copula (DDAC) model [20]. The above dynamic structures provide ample potential options for building the time-varying NCCN copula.
The contributions of this paper are as follows. First, a new multivariate skewed fat-tailed distribution, namely, multivariate NCCN distribution, is constructed. Second, the copula of the multivariate NCCN distribution, namely, multivariate NCCN copula, is proposed. Third, we adopt the DCC model to construct a new time-varying copula model, namely, DCC-NCCN copula model. The last, employing the Hang Seng Index (HSI), Hang Seng China Enterprises Index (CEI), and Hang Seng China-Affiliated Corporations Index (CCI) as our sample data, we compare the fitting effects of the DCC-NCCN copula model with some other copula models and perform the visualized dependence analysis of Hong Kong stock markets.
2. Model Development
2.1. Fundamental Theory of the Copula
A copula is a multivariate cumulative distribution function (cdf) with uniform univariate margins and can be used to link univariate margins to a joint cdf. According to Sklar’s theorem, for a d-dimensional random vector with joint cdf and marginal cdfs , there exists a copula function such that
The copula is unique if the random vector is continuous. For a continuous random vector with joint cdf , joint probability density function (pdf) , marginal cdfs , marginal pdfs , and marginal quantile functions , the copula function and its pdf are, respectively, given bywhere . According to Sklar’s theorem, we can quickly get the copula of a given multivariate distribution. In particular, can be called the uniform scores of the random variables , and can be called the normal scores of the random variables , where is the quantile function of the univariate standard normal distribution. Clearly, the uniform score follows the univariate uniform distribution on , and the normal score follows the univariate standard normal distribution.
The copula function is closely related to many dependence measures, such as Kendall’s tau, quantile dependence (QD) coefficient, and Mardia’s skewness and kurtosis of normal scores. These dependence measures are briefly described in the following.
Kendall’s tau is also called Kendall’s rank correlation coefficient. It can be utilized to measure global dependence. Let and be independent and identically distributed random vectors. For the bivariate continuous random vector with joint cdf , copula function , and uniform scores , bivariate Kendall’s tau is given by
The value range of bivariate Kendall’s tau is . Bivariate Kendall’s tau can be interpreted as the probability of concordance minus the probability of discordance.
The quantile dependence coefficient can be used to measure local dependence. For a bivariate continuous random vector with copula function and uniform scores , the bivariate lower-lower (lower), upper-upper (upper), upper-lower, and lower-upper quantile dependence coefficients (LLQD, UUQD, ULQD, and LUQD) are, respectively, given bywhere is the quantile level. The value range of bivariate quantile dependence coefficients is . If , then . If , then and . If , we can get the bivariate tail dependence coefficients (LLTD, UUTD, ULTD, and LUTD).
Mardia [21] proposed Mardia’s skewness and kurtosis. For a multivariate continuous random vector with mean vector and covariance matrix , Mardia’s skewness and kurtosis are, respectively, given bywhere and are independent identically distributed random vectors. For a multivariate continuous random vector with the multivariate normal distribution, multivariate Mardia’s skewness and kurtosis are 0 and , respectively. Mardia’s skewness and kurtosis of normal scores can also be called Gaussian skewness and kurtosis. Gaussian skewness and kurtosis can be used to measure the asymmetric dependence and tail dependence, respectively. Note that the Gaussian skewness cannot measure the direction of asymmetric dependence. For a multivariate random vector with the multivariate normal copula, Mardia’s skewness and kurtosis are not clear, but the Gaussian skewness and kurtosis are 0 and , respectively. Note that the tail dependence coefficients cannot reasonably distinguish the strength of tail dependence. The two copulas with the same tail dependence coefficients may have different Gaussian kurtosis.
2.2. Nonlinear Asymmetric GARCH (NAGARCH) Model
Before modeling the dependence structure, we need to model the marginal distribution. This study adopts the NAGARCH model of Engle and Ng [22] to describe the dynamics of financial asset return series. The parameterization form of the NAGARCH model is not unique, and the distribution assumption is not unique. To easily explain the parameter of the NAGARCH model, this paper adopts a variance targeting (VT) form. To avoid the distribution specification error, this paper does not assume a specific distribution. We set the NAGARCH model:where is the return, is the unconditional mean, is the unconditional standard deviation, is the conditional standard deviation, is the residual with mean 0, and is the standardized residual with mean 0 and variance 1.
The conditional variance can be interpreted as the asymmetric information shock item plus the weighted average of unconditional variance and lagged conditional variance . This equation makes the parameter representation clearer. Parameter can control the dynamics of the conditional variance: the larger , the stronger the dynamics of the conditional variance. In particular, when = 0, the model is a constant volatility model. Parameter can control the clustering and mean reversion of the conditional variance: when is close to 1, the conditional variance shows stronger clustering and weaker mean reversion; when is close to 0, the conditional variance shows weaker clustering and stronger mean reversion. Parameter γ can control the asymmetric dynamics of the conditional variance: when > 0, the conditional variance has a positive asymmetry; when < 0, the conditional variance has a negative asymmetry, and the negative value impacts more on the conditional variance than the same degree of a positive one. In particular, when = 0, the model is the GARCH model. The conditional variance equation can also be expressed as
Clearly, all conditional variances can be insured to be greater than 0 under given parameter constraints.
In terms of parameter estimation, we adopt the quasi-maximum likelihood (QML) method to estimate the parameters of the NAGARCH model, that is, to apply the maximum likelihood (ML) method to conduct the estimation of the parameters of the NAGARCH-normal model.
The NAGARCH model is used to filter the return series to obtain the standardized residual series. Then, using the empirical cdf, we transform the standardized residual series into the uniform scores. For a standardized residual series , the empirical cdf iswhere is the indicator function. The uniform scores can be applied to further model the dependence structures.
2.3. Multivariate NCCN Distribution and Multivariate NCCN Copula
We firstly introduce the multivariate NCCN distribution. Because the location and scale parameters of the multivariate NCCN distribution cannot influence the multivariate NCCN copula [23], this paper does not take the location and scale parameters into account when defining the multivariate NCCN distribution.
Let and be the pdf and cdf of the univariate standard normal distribution , and be joint the pdf and joint cdf of the bivariate standard normal distribution with the linear correlation coefficient , and and be the joint pdf and joint cdf of the multivariate standard normal distribution with the linear correlation matrix R. , , and can be expressed as:
The multivariate NCCN distribution can be interpreted as a multivariate noncentral normal scale mixture distribution [24]. The stochastic representation of the multivariate NCCN distribution can be given bywhere the random vector follows the NCCN distribution. The random vector follows the multivariate noncentral normal distribution with correlation matrix and noncentral parameter vector . The random variable is a two-point distribution with probability parameter , shape parameter , and mean 1. and are independent of each other.
The multivariate NCCN distribution can also be interpreted as a simplified mixture of two multivariate normal distributions [25]. The stochastic representation of the multivariate NCCN distribution can also be given by
The parameters of the multivariate NCCN distribution can be divided into three parts: (1) a correlation matrix , (2) two tail parameters and , and (3) a skewness vector . According to the two stochastic representations, the joint cdf and pdf of the multivariate NCCN distribution can be easily given bywhere , , is a correlation matrix, , , and . The multivariate NCCN distribution family includes multivariate contaminated normal (CN) distribution when and multivariate normal distribution .
The mean vector, covariance matrix, and linear correlation matrix of the multivariate NCCN distribution are, respectively, given as equations (13)–(15):
Obviously, matrix R is not the linear correlation matrix of the multivariate NCCN distribution. For the multivariate CN distribution, , Mardia’s skewness is 0 and Mardia’s kurtosis is .
As to the two stochastic representations, any bivariate marginal distribution of the multivariate NCCN distribution is a bivariate NCCN distribution. The joint cdf and pdf of the bivariate NCCN distribution can be easily given bywhere , , , , , and . Note that this bivariate NCCN distribution can be regarded as the simplified bivariate mixed-normal distribution.
The linear correlation coefficient of the bivariate NCCN distribution is
Although the correlation parameter ρ is not equal to the linear correlation coefficient of the bivariate NCCN distribution, it changes in the same direction as the linear correlation coefficient. The two skewness parameters have a substantial influence on the linear correlation coefficient: when , ; when , ; and when , . The two tail parameters have a substantial influence on the linear correlation coefficient: when or , .
Each univariate marginal distribution of the multivariate NCCN distribution obviously follows the univariate NCCN distribution. The cdf and pdf of the univariate NCCN distribution are, respectively,where , , , , and . It is difficult to simplify the quantile function of the univariate NCCN distribution, . According to the cdf and pdf of the univariate NCCN distribution, the quantile function can be calculated by Newton’s method. Note that this univariate NCCN distribution can be regarded as the univariate mixed-normal distribution without location and scale parameters.
The mean, variance, skewness, and kurtosis of the univariate NCCN distribution are
For the univariate CN distribution, , , , and . Clearly, the two tail parameters and inversely change with the kurtosis.
According to the above description, the flexibility of the univariate, bivariate, and multivariate CN distributions are similar to that of the univariate, bivariate, and multivariate t-distributions, respectively. The flexibility of the univariate, bivariate, and multivariate NCCN distributions are similar to that of the univariate, bivariate, and multivariate skew-t distributions, respectively.
According to Sklar’s theorem, the multivariate NCCN copula function and its pdf can be easily expressed aswhere , is a correlation matrix, , , and . Similar to the multivariate NCCN distribution, subclasses of the multivariate NCCN copula include multivariate normal copula and multivariate CN copula. Clearly, the multivariate NCCN copula cannot be regarded as a mixture of two multivariate normal copulas.
The bivariate NCCN copula function and its pdf can be expressed aswhere , , , , and . The bivariate NCCN copula density function diagrams are given in Figures 1 and 2.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(a)

(b)

(c)

(d)
According to the bivariate NCCN copula density function diagrams, the meaning of parameters can be easily understood. The correlation parameter ρ can affect global dependence: when ρ is larger, the negative dependence is weaker, and the positive dependence is stronger. The two tail parameters and can affect the tail dependence: when they are smaller, the tail dependence is stronger. The two skewness parameters can affect the asymmetric dependence: when and , the lower-lower tail dependence is stronger than the upper-upper tail dependence; when and , the upper-upper tail dependence is stronger than the lower-lower tail dependence; when and , the lower-upper tail dependence is stronger than the upper-lower tail dependence; and when and , the upper-lower tail dependence is stronger than the lower-upper tail dependence.
2.4. Time-Varying NCCN Copula
For the multivariate NCCN copula, this paper further considers that the correlation matrix may change over time and does not consider the dynamics of other parameters.
The DCC model is a basic dynamic correlation model. This study presents the DCC model aswhere T is the sample size, is the standardized residual vector, is the conditional correlation matrix of , is similar to the conditional covariance matrix, and is similar to the unconditional covariance matrix and can be estimated by using the sample mean of . Note that can be interpreted as an information shock term plus the weighted average of and .
Parameter α controls the dynamics of the conditional correlation matrix: when α is large, the dynamics of the conditional correlation matrix are strong. In particular, when α = 0, the model degenerates into the constant conditional correlation (CCC) model. Parameter β controls the clustering and mean reversion of the conditional correlation matrix: when β is close to 1, the conditional correlation matrix shows strong clustering and weak mean reversion; when β is close to 0, the conditional correlation matrix shows weak clustering and strong mean reversion. The renewal equation of can also be expressed as
In this renewal equation, can be interpreted as the weighted average of , and . Under the given parameter constraints, can be guaranteed to be a true correlation matrix.
The common DCC-normal copula model can be given bywhere follows the multivariate normal copula with time-varying parameter matrix and follows the standardized multivariate normal distribution with time-varying linear correlation matrix . For a sample of uniform scores , we can employ the ML method to estimate the parameter set .
Similar to the DCC-normal copula model, the DCC-NCCN copula model can be given bywhere is an identity matrix, follows the multivariate NCCN copula with time-varying parameter matrix , follows the multivariate NCCN distribution, and follows the standardized multivariate NCCN distribution with time-varying linear correlation matrix . Because the meaning of is not very clear, we first portray time-varying and then get time-varying .
Akaike information criterion (AIC) and Bayesian information criterion (BIC) can be used to compare the fitting effect of different models. The expressions for AIC and BIC arewhere LL is the log-likelihood, is the number of parameters, and is the sample size. The smaller the AIC and BIC are, the better the model is when we compare the models.
3. Empirical Results
3.1. Descriptive Statistics
This study employs the Hang Seng Index (HSI), Hang Seng China Enterprises Index (CEI), and Hang Seng China-Affiliated Corporations Index (CCI) from the Hong Kong stock market as our sample, abbreviated as HSI, CEI, and CCI thereafter, respectively. We obtained three daily closing price series from the period from Jan 1, 2005, to Dec 31, 2018, with 3451 data, respectively. The data source can be found at https://cn.investing.com, and we calculate the daily logarithmic return, , where is the daily closing price at time , and we have 3 return series, with 3450 observations each.
The relative price series are given in Figure 3.

The graph shows that all price series have strong dynamics and an upward trend. Considering long-term investment, CCI is the best choice, and HSI is the worst choice. The return series are given in Figure 4.

(a)

(b)

(c)
Return series have significant dynamic characteristics, and the dynamic process shows significant mean reversion.
The univariate descriptive statistics of returns are presented in Table 1.
As expected, the minimum values are negative, and the maximum values are positive. The range of all returns is large. The median and mean values are close to zero, where the mean values are less than the median values. The standard deviation values are greater than 1. According to skewness, HSI is skewed to the left, and CEI and CCI are skewed to the right. The kurtosis values are significantly larger than 3, implying that returns have fatter tails than the normal distribution. The skewness and kurtosis tests show that the returns cannot follow the normal distribution.
Using the Ljung-Box Q(5) test method, the autocorrelation tests of the first four moments of returns are reported in Table 2.
As for the tests, the autocorrelation of returns is weak, but the autocorrelation of second-order, third-order, and fourth-order moments of returns is strong, indicating that the returns cannot have serial independence. The autocorrelation of squared returns is particularly prominent, indicating that the dynamic of volatility (variance or standard deviation) is the most important.
Using window length , the moving sample standard deviation series of returns are given in Figure 5.

(a)

(b)

(c)
For each return series, the time-varying volatility can be easily observed.
3.2. Fitting of the Marginal Distribution
We employ the NAGARCH model to describe the dynamics of the return series. Table 3 shows the estimation results of the NAGARCH model.
The values of parameter are greater than 0.05, indicating significant time-varying volatility. The values of parameter are close to 1, showing strong clustering and weak mean reversion. The values of parameter are less than 0, exhibiting volatility asymmetry. The standard deviation series of the NAGARCH model are given in Figure 6.

(a)

(b)

(c)
The standard deviation series of the NAGARCH model are consistent with the moving sample standard deviation series, indicating that the NAGARCH model can effectively describe the time-varying volatility. Based on the NAGARCH model, we can obtain the standardized residual series. The univariate descriptive statistics of standardized residuals are given in Table 4.
Compared with the return series, the range of standardized residual series is significantly reduced. The sample mean values of standardized residuals are almost equal to 0, and the sample standard deviation is almost equal to 1, which can meet the theoretical requirements. The skewness values of standardized residuals are quite different from returns. The kurtosis values of standardized residuals are smaller than returns. Based on the skewness and kurtosis tests, standardized residuals cannot follow a normal distribution.
Using the Ljung-Box Q(5) test method, the autocorrelation tests of the first four moments of standardized residuals are reported in Table 5.
Based on the tests, the autocorrelation of the first four moments of standardized residuals is not strong. The standardized residual series can basically meet the serial independence. In general, the NAGARCH model can effectively portray the dynamics of each return series.
Using the empirical cdf to transform standardized residuals into uniform scores, uniform scores satisfy the serial independence and follow a uniform distribution on [0,1].
3.3. Descriptive Analysis of Dependence Structures
To perform some sample analyses of the bivariate dependence structures, the sample bivariate dependence measures of uniform scores are reported in Table 6.
The bivariate global dependence measures are positive. The bivariate global dependence is the smallest for CEI-CCI and largest for HSI-CEI. The bivariate lower tail dependence measures are larger than the upper ones, implying that the bivariate dependence structures have stronger lower tail dependence. The upper-lower and lower-upper tail dependence measures are very close to zero. The Gaussian skewness tests show that the bivariate dependence structures are significantly asymmetric. The Gaussian skewness is the smallest for HSI-CCI and largest for HSI-CEI. The Gaussian kurtosis values are larger than 8, implying that the bivariate dependence structures have stronger tail dependence than the bivariate normal copula. The Gaussian kurtosis is the smallest for HSI-CCI and largest for HSI-CEI. The Gaussian skewness and kurtosis tests show that the bivariate dependence structures cannot follow the normal copula.
The bivariate scatter plots of uniform scores are given in Figure 7.

(a)

(b)

(c)
The points are mainly concentrated around the main diagonal. The points in the lower-lower and upper-upper tail regions are dense, but the points in the lower-upper and upper-lower tail regions are sparse.
To understand the bivariate local dependence of uniform scores, the sample quantile dependence curves of uniform scores are given in Figure 8. The horizontal axis shows quantile levels, and the vertical axis shows bivariate quantile dependence coefficients.

(a)

(b)

(c)
The sample quantile dependence curves show the following features: (1) LLQD and UUQD curves are significantly higher than ULQD and LUQD curves, indicating that all bivariate dependence structures have a strong positive dependence. (2) LLQD and UUQD curves are obviously not coincident, and the LLQD curve is significantly higher than the UUQD curve at the low quantile levels, indicating that all bivariate dependence structures are asymmetric, and the lower tail dependence is significantly higher than the upper tail dependence. (3) ULQD and LUQD curves are almost coincident.
Using window length , moving sample bivariate Kendall’s tau series of uniform scores are given in Figure 9.

(a)

(b)

(c)
The time-varying bivariate global dependence can be easily observed. From 2005 to 2007, the global dependence of HSI-CCI is the largest. However, since 2008, the global dependence of HSI-CEI is the largest.
3.4. Fitting of Dependence Structures
This paper considers CCC-N, CCC-CN, CCC-NCCN, DCC-N, DCC-CN, and DCC-NCCN copula models. Considering that the asymmetric dependence between uniform scores is mainly in the upper and lower tail, we can constrain all skewness parameters of the NCCN copula to be equal. In addition, we consider three common Archimedean copulas, namely, Clayton, Gumbel, and Frank copulas.(1)Clayton copula function: where Kendall’s tau of Clayton copula is .(2)Gumbel copula function: where . Kendall’s tau of Gumbel copula is .(3)Frank copula function:where . Kendall’s tau of Frank copula is .
For the Clayton, Gumbel, and Frank copulas, sample Kendall’s tau can be used to estimate their parameters. For other copulas, the maximum likelihood estimation (MLE) can be used.
To compare the fitting effect, the CCC-N copula model can be used as the benchmark model. Then, the parameter increment ∆P and LL increment ∆LL of each copula model can be calculated. Also, the models can be ranked by LL, AIC, and BIC, respectively. The fitting results of bivariate copula models are given in Tables 7–9.
The values of parameter α are greater than 0.02, and the values of parameter β are close to 1. The values of two tail parameters are not close to 1, implying that all bivariate dependence structures have stronger tail dependence than the normal copula. The values of the skewness parameter are negative, implying the lower tail dependence is stronger than the upper tail dependence for each bivariate dependence structure. The ∆LL values show that the fitting of Clayton, Gumbel, Frank, CCC-N, CCC-CN, CCC-NCCN, DCC-N, DCC-CN, and DCC-NCCN copula models is improved in turn. In terms of ranking, the DCC-NCCN copula model is the best.
To easily understand the fitting effect of the bivariate local dependence, the quantile dependence curves of bivariate copula models are given in Figures 10–15.

(a)

(b)

(c)

(a)

(b)

(c)

(a)

(b)

(c)

(a)

(b)

(c)

(a)

(b)

(c)

(a)

(b)

(c)
Compared with the sample QD curves, the bivariate Clayton copula model overestimates the degree of asymmetric dependence in the upper and lower tails. The bivariate Gumbel copula model has a wrong asymmetric direction on the upper and lower tail dependence, and it overestimates the strength of the asymmetric tail dependence. The bivariate Frank copula model cannot describe the asymmetric dependence and seriously underestimates the strength of the upper and lower tail dependence. The bivariate CCC-N copula model cannot describe the asymmetric dependence and significantly underestimates the degree of the lower tail dependence. The bivariate CCC-CN copula model cannot describe the asymmetric dependence. The bivariate CCC-NCCN copula model is basically correct.
To easily understand the fitting effect of the bivariate time-varying global dependence, Kendall’s tau series of the bivariate DCC-NCCN copula model are given in Figure 16.

(a)

(b)

(c)
Kendall’s tau series of the bivariate DCC-NCCN copula model is basically consistent with moving sample Kendall’s tau series. The results illustrate that the bivariate DCC-NCCN copula model can better depict the bivariate time-varying global dependence.
The 10% QD coefficient series of the bivariate DCC-NCCN copula model are given in Figure 17. Note that the 10% ULQD and 10% LUQD coefficient series are omitted because their values are very close to 0.

(a)

(b)

(c)
The dynamic characteristics of 10% LLQD and 10% UUQD coefficient series are consistent with Kendall’s tau series. For each bivariate dependence structure, the time-varying lower and upper tail dependence can be easily observed.
In comparison with the fitting effects of the multivariate copula models, the fitting results of multivariate copula models are given in Table 10.
The fitting results of the multivariate copula models are consistent with the fitting results of the bivariate copula models. Based on LL values, the time-varying dependence, tail dependence, and asymmetric dependence all play an important role in improving the fitting effect of a multivariate dependence structure. The multivariate DCC-NCCN copula model is the best choice.
Some diagrams of the multivariate CCC-NCCN copula and DCC-NCCN copula models are given in Figures 18–20. Compared with bivariate CCC-NCCN and DCC-NCCN copula models, the effects of the multivariate CCC-NCCN copula and DCC-NCCN copula models have no significant differences.

(a)

(b)

(c)

(a)

(b)

(c)

(a)

(b)

(c)
4. Conclusion
This study examines the effects of the DCC-NCCN copula model and some other copula models on fitting dependence structures of Hong Kong stock markets. The main conclusions in this paper are as follows: First, according to descriptive statistics and fitting results of the marginal distribution, return series of HSI, CEI, and CCI all reveal significant time-varying volatility. NAGARCH model can well depict the dynamic characteristics of returns. Second, descriptive statistics and fitting results show that the bivariate dependence structures have strong positive dependence, asymmetric dependence, tail dependence, and time-varying dependence. For each bivariate dependence structure, the lower tail dependence is higher than the upper tail dependence. Third, through the comparison of the DCC-NCCN copula model and some other copula models, the DCC-NCCN copula model can well describe the bivariate dependence structures, but other copula models are not good. Considering the flexibility and complexity, the DCC-NCCN copula model is a relatively ideal copula model.
Data Availability
The data used to support the findings of this study are public and available on HK stock exchange or database.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this study.
Acknowledgments
This research was supported by the Key Program of National Social Science Foundation of China under Grant 18AGL001.