Measuring systemic risk for financial time series: A dynamic bivariate Dvine model

1.   Introduction

Modeling multivariate time series is an important statistical application in different fields. There are two types of correlations embedded in multivariate time series: serial dependence, which refers to the dependence between different time points of a univariate time series, and cross-sectional dependence, which refers to the dependence between different dimensional time series at the same time point. Accurately measuring these two types of dependence is essential for evaluating the effectiveness of multivariate time series models.

Copula as a multivariate distribution with uniformly distributed marginals can flexibly model the dependence structure of multivariate variables, see Refs. [1–3] for details. The copula function can be used to model two types of dependence between multivariate variables. Currently, the application of the copula function to multivariate time series mainly focuses on cross-sectional dependence. When dealing with multivariate time series models, the popular method is to use the GARCH family models to analyze the serial dependence, and then employ the copula function to fit the correlation between the filtered data (see Refs. [4–7]). Moreover, Chen and Fan^[8] summarized the class of multivariate time series models as SCOMDY models and studied the theoretical properties of such models. In the framework of SCOMDY, Patton^[9] extended constant copula to time-varying copula by turning the copula parameters into functions of lagged variables and lagged exogenous variables. Giulio and Ergün^[10] proposed a dynamic GAS copula model whose iteration over time-varying parameters relies on a scaled score function of predictive model density rather than simply moments of observations. Time-varying copulas have succeeded in describing the dynamic cross-sectional dependence between multi-dimensional distributions in the past few years.

Modeling serial dependence based on copulas has also attracted attention. Joe^[11] proposed a class of stationary first-order Markov models based on bivariate copulas and marginal distributions. Chen and Fan^[12] and Beare^[13] generalized the class of models to semiparametric stationary Markov models and discussed asymptotic properties of the semiparametric Markov models. Zhao et al.^[14] extended the univariate copula-based time series model to higher order by using the D-vine copula. See Refs. [15–18] for more details of D-vine copula model. Bladt and McNeil^[19] constructed a new class of time series copula models using d-vines and v-transforms. Moreover, copula can also simultaneously capture the serial dependence and cross-sectional dependence in a class of vine copula models. Brechmann and Czado^[20], Smith^[21], and Beare et al.^[22] introduced different vine copula models that can account for both types of dependence in multivariate time series. Nagler et al.^[23] summarized a stationary vine copula model which can use a more general vine structure to fit multivariate time series including the above three models as special cases.

For financial time series, multivariate GARCH models focus on the linear dependence structure between variables; see, Refs. [24, 25]. However, copula functions can capture the nonlinear and asymmetric dependence that exists in multivariate time series. Thus, we design a dynamic bivariate Dvine model in which serial dependence and cross-sectional dependence are modeled based on copula functions. Specifically, inspired by Ref. [14], we consider a class of univariate copula-based stationary Markov models (Dvine), in which D-vine copulas are used to capture the serial dependence in univariate time series and marginal distributions are suitable parametric distributions. Subsequently, we propose the dynamic bivariate Dvine model in which the time-varying copula is used to capture the dynamic cross-sectional dependence by connecting the two-dimensional Dvine model. Based on the complicated and dynamic nature of the financial market, the dynamic bivariate Dvine model can provide a flexible and nonlinear copula-based method compared to the existing time series models.

In the financial field, the ultimate goal of modeling time series is to measure the future financial situation through past observations, such as the conditional quantile of interested financial companies. Therefore, forecasting the future behavior of financial assets has become the main method of risk measurement (e.g., Refs. [10, 26–28]). Given past observations, the corresponding estimator of VaR of financial assets can be constructed via the univariate Dvine model. For our dynamic bivariate Dvine model, we can also provide forecasting methods for CoVaR of the financial system when the financial company is in distress. Moreover, given the asymptotic properties of the parameters estimated by the univariate Dvine model, we can obtain the $\sqrt{T}$-consistency and asymptotic normality of the estimators of the conditional quantiles by the Delta method.

Our proposed model enriches the existing literature in three ways. First and foremost, we generalize the existing method of modeling the dependence of time series by proposing the dynamic bivariate Dvine model. Second, we provide an alternative method to predict the conditional quantile based on our proposed model. In the empirical analysis, the backtesting results show that our method can rival and outperform the traditional GARCH method. Third, by taking the two-stage likelihood method, we can obtain the estimators of our model parameters and study the statistical properties of the estimators of the Dvine model.

The rest of our paper is organized as follows. In Section 2, we introduce the dynamic bivariate Dvine model. Section 3 introduces the parameter estimation for our dynamic bivariate Dvine model and estimators of the conditional quantiles. In Section 4, we establish the asymptotic properties of the estimators of the univariate Dvine model. In Section 5, real data applications on financial stock prices show good performance compared to the traditional GARCH method. Finally, we offer some concluding remarks in Section 6. All the technical details are provided in Supporting information.

2.   Methodology

In this section, we propose a copula-based time series model which can provide a new tool for capturing joint behavior in the financial market. As we know, the dependence structure that exists in financial time series is always time-varying and nonlinear. We build the dynamic bivariate Dvine model to specify the dependence that exists in the bivariate time series $\left\{{\boldsymbol{X}}_{t} = \left(X_{t,1}, X_{t,2}\right)\right\}_{t = 1}^{T}$. Let $F_{t}\left(x_{t,1},x_{t,2}\mid {\cal{F}}_{t-1}\right)$ denote the joint distribution of ${\boldsymbol{X}}_{t}$ conditional on the past information ${\cal{F}}_{t-1} = \sigma\left({\boldsymbol{X}}_{t-1},\cdots,{\boldsymbol{X}}_{1}\right)$. The dynamic bivariate Dvine model for $\left\{{\boldsymbol{X}}_{t}\right\}$ is

where the marginal behavior of $X_{t,i}$ is captured by the univariate Dvine model with parameters $\left(\eta_{i},{\boldsymbol{\theta}}_{i}\right)$, and the joint behavior of ${\boldsymbol{X}}_{t}$ is characterized by the dynamic copula function $C_t$ with dependence parameter $\rho_{t}$ which means the cross-sectional dependence may change over time. We will introduce the specific forms in the following subsections.

The dynamic bivariate Dvine model is sensible and interpretable and links the univariate Dvine model by the dynamic copula function, which accurately reflects the nonlinear, time-varying dependence structure of financial time series. We emphasize the time-varying nature of the cross-sectional dependence in the dynamic bivariate Dvine model, which naturally extends the copula-based time series model in low-dimensional applications.

2.1   Univariate Dvine model

We develop a D-vine copula-based model for univariate time series $\left\{X_{t,i}\right\}_{t = 1}^{T}$ which satisfies the requirement of strict stationarity and Markov properties, $i = 1,2$. Inspired by the work of Ref. [14], a class of univariate Dvine models is employed to construct the conditional marginal distribution at time $t$ conditional on the past observations ${\cal{F}}_{t-1}^{i}$, $i = 1,2$.

A time series $\left\{X_{t,i}\right\}$ is called a Dvine($p_{i}$) copula process if, for all $T\geqslant 2$, the joint density of $\left\{X_{t,i}\right\}_{t = 1}^{T}$ can be written as

where $f(\cdot;\eta_{i})$ is the parametric marginal density with parameter $\eta_{i}$, $c_{k,i}$ is the bivariate copula in tree $k$ with parameter $\theta_{k,i}$ and ${\boldsymbol{\theta}}_{i} = \left\{\theta_{k,i}\right\}_{k = 1}^{p_{i}\vee (T-1)}$ are the set of all bivariate copulas parameters. The algorithm to obtain $F(x_{t,i}\mid {\cal{F}}_{t-1}^{i})$, the conditional distribution of $X_{t,i}$, is provided in Supporting information S.1.

The unconditional marginal distribution of our univariate time series model is estimated parametrically while Ref. [14] use a nonparametric approach. Utilizing the Dvine($p$) model to analyze the serial dependence present in univariate time series can provide a flexible way to analyze univariate time series in a nonlinear form without relying on the first or second moment of observations as GARCH-type models do.

2.2   Modeling the parameter of cross-sectional dependence

The modeling process for the dependence parameter of the cross-sectional copula is pivotal, and we introduce the GAS method^[29] to design time-varying mechanisms for the dependence parameter $\rho_{t}$. Specifically, we assume $C_{t}$ to be the Student’s t copula. To capture the dynamic cross-sectional dependence, we update the correlation parameter $\rho_{t}$ based on past observations. We transform it as $\rho_{t} = [1-\exp(-f_{t})]/ [1+\exp(-f_{t})]$ to ensure $\rho_{t}\in[-1,1]$, where $f_{t}$ is the time-varying parameter which can be updated by the following autoregressive process,

where $s_{t-1}$ denotes an appropriate function of the past information. Using the GAS framework, the driving mechanism $s_{t}$ is considered as the scaled score function of the observation density, and its specific form is given by

where ${\boldsymbol S}_{t}$ denotes the scaling matrix. In this paper, we define the scaling matrix as the square root matrix of the inverse of the Fisher information matrix with respect to the parameter $f_{t}$ and thus ${\boldsymbol S}_{t} = {\cal{J}}_{t \mid t-1}$, where ${\cal{J}}_{t \mid t-1} = \left({E}_{t-1}\left[\nabla_{t} \nabla_{t}^{\prime}\right]\right)^{-1/2}$. The detailed formula of the driving mechanism $s_{t}$ is presented in Supporting information S.2.

Together, the dynamic bivariate Dvine model for the time series $\left\{{\boldsymbol{X}}_{t}\right\}$ can be expressed as follows:

In this model, the conditional marginal distribution of $X_{t,i}$ is established by the univariate Dvine model with the parameters $\left(\eta_{i},{\boldsymbol{\theta}}_{i}\right)$, and the joint dynamic behavior of ${\boldsymbol{X}}_{t}$ is characterized by the bivariate copula $C_{t}$ with the time-varying coefficient parameter $\rho_{t}$. To our best knowledge, the dynamic bivariate Dvine model is the first copula-based time series model that can handle dynamic behavior. We emphasize the time-varying nature of the dynamic bivariate Dvine model, as it is built upon the assumptions of using a time-varying copula for the cross-sectional dependence. This makes it more suitable for financial time series modeling.

3.   Estimation and prediction

In this section, we provide the parameter estimation for our proposed bivariate Dvine model and introduce the prediction method via the proposed dynamic bivariate Dvine model.

3.1   Estimation of the proposed model parameters

Based on the model we introduced in the previous section, the parameters that need to be determined are ${\boldsymbol{\theta}} = \bigcup_{i = 1}^{2}\left(\eta_{i}, {\boldsymbol{\theta}}_{i}\right) \bigcup\left(w,\beta_{0}, \beta_{1},\sigma\right)$, where $\eta_{i}$ is the marginal distribution parameter of the Dvine($p_{i}$) model, ${\boldsymbol{\theta}}_{i}$ is the parameter set of the bivariate copulas in the Dvine($p_{i}$) model, $i = 1,2$, and ${\boldsymbol{\theta}}_{c} = \left(w,\beta_{0},\beta_{1},\sigma\right)$ denotes the parameter for the time-varying copula of the bivariate time series. The estimation of ${\boldsymbol{\theta}}$ is based on the likelihood approach. Denote $\left\{{\boldsymbol{x}}_{t} = (x_{t,1},x_{t,2})\right\}_{t = 1}^{T}$ as the $T$ observations of the bivariate time series, the complete log-likelihood function of $\left\{{\boldsymbol{x}}_{t}\right\}_{t = 1}^{T}$ is given by

where the conditional CDF (cumulative distribution function) $F_{i}\left(x_{t,i}\mid {\cal{F}}_{t-1}^{i};{\boldsymbol{\theta}}_{i}\right)$ can be derived from the Dvine($p_{i}$) model and $p = \max\left\{p_{1},p_{2}\right\}$. If we want to calculate the maximum likelihood estimates of all parameters of the model simultaneously, the calculation is complicated and inefficient. Joe^[30] and Tsukahara^[31] discussed the stepwise method in the copula-based models and proved its efficiency and theoretical properties. To improve our computational efficiency, we follow the stepwise maximum likelihood estimation method. The complete log-likelihood function in Eq. (6) is decomposed into several parts and optimized in successive steps.

In the first stage, for $i = 1,2$, the estimators of parameters $(\eta_{i},{\boldsymbol{\theta}}_{i})$ contained in the Dvine($p_{i}$) model are obtained using univariate time series $\left\{x_{t,i}\right\}_{t = 1}^{T}$. In this paper, we mainly adopt parametric methods to fit the marginal distributions and transform the observed data into pseudo-observations. In general, we define an appropriate parametric marginal density $f_{i}(\cdot;\eta_{i})$ for the $i$th univariate time series. The marginal parameters $\eta_{i}$ can be calculated by the following equation:

Then the pseudo-observations $\widehat{u}_{t,i}$ can be calculated by the parametric distribution function $\widehat{F}_{i}(\cdot;\eta_{i})$, where $\widehat{u}_{t,i} = \widehat{F}_{i}(x_{t,i};\eta_{i}),\; t = 1,\cdots,T, \;i = 1,2$. The vector $(\widehat{u}_{t,1},\widehat{u}_{t,2})$ can almost be seen as an i.i.d sample from the random vector $(U_{1},U_{2})^{\prime}$ and be used to estimate the parameters $\left({\boldsymbol{\theta}}_{1},{\boldsymbol{\theta}}_{2}\right)$ of bivariate copulas in Dvine models.

Given the pseudo-observations $\left\{\widehat{u}_{t,i}\right\}_{t = 1}^{T}$, $\widehat{{\boldsymbol{\theta}}}_{i}$ can be calculated by maximizing the joint log-likelihood of the Dvine($p_{i}$)model. Let $c_{k,i}(\cdot,\cdot;\theta_{k,i})$ denote the bivariate copula in tree $k$ with parameter $\theta_{k,i}$, the conditional distribution $F_{k,i}(x_{t+k,i}\mid x_{t+1,i},\cdots,x_{t+k-1,i})$ can be calculated by the bivariate copula function in tree $k$. Then, the log-likelihood function of the Dvine($p_{i}$) model for the pseudo-observations is given by

where $F_{k-1,i}\left(\widehat{u}_{t+k,i} \mid \widehat{\boldsymbol{u}}_{t+1,t+k-1}^{i};{\boldsymbol{\theta}}_{[k-1],i}\right)$ denotes the conditional distribution of $\widehat{u}_{t+k,i}$ given the intermediate variables $\widehat{\boldsymbol{u}}_{t+1,t+k-1}^{i} = (\widehat{u}_{t+1,i},\cdots,\widehat{u}_{t+k-1,i})$ and parameter set ${\boldsymbol{\theta}}_{[k-1],i} = (\theta_{1,i},\cdots,\theta_{k-1,i})$.

Through the tree-by-tree sequential selection methodproposed by Ref. [14], we determine the family and parameters of the copula function starting from the first tree, and obtain the estimated parameter of the copula in the second tree based on the result of the first tree. The process continuesuntil an independent copula is selected in a certain tree. Consequently, for $k = 1,\cdots,p_{i}$, we get

where the vector $\widehat{{\boldsymbol{\theta}}}_{[k-1],i}$ denotes the estimated parameters from the $k-1$ tree. Then the conditional CDF $F_{k,i}$ are calculated by the partial derivative or $h$-function of the copula $C_{k,i}$.

In the second stage, the parameter vector ${\boldsymbol{\theta}}_{c}$ of the time-varying copula function can be obtained by the following equation:

Compared to the full maximum likelihood estimator (MLE), the two-stage MLE provides a more efficient way to obtain the parameter estimates of our bivariate time series model.

3.2   Prediction

In the financial time series field, tail risk measures have been widely employed to measure financial asset risk. Given the past observations ${\cal{F}}_{t-1}$, the dynamic bivariate Dvine model can make predictions for the future tail behavior of ${\boldsymbol{X}}_{t}$. At time $t$, the conditional quantile or conditional VaR of $X_{t,i}$ at level $\alpha \in (0,1)$ given the past information ${\cal{F}}_{t-1}^{i}$ can be expressed as

For exposition, we first give the explicit equation for the $\alpha$th conditional quantile of $X_{t,i}$ when using the Dvine(2) model,

where $C_{k,i}(\cdot\mid \cdot;\theta_{k,i})$ denotes the conditional copula function in the tree $k$, and $C_{k,i}^{-1}(\cdot\mid \cdot;\theta_{k,i})$ denotes the inverse of the conditional copula function in tree $k$. We provide explicit expressions for conditional quantiles, but this is based on the fact that the copula functions included in the Dvine model all have the inverse of the conditional distribution. Equations for the $\alpha$th conditional quantiles of $X_{t, i}$ when the Dvine model is of order 1 or higher are in Supporting information S.3.

The CoVaR proposed by Adrian and Brunnermeier^[27] is a conditional quantile that measures the contribution of financial institutions to systemic risk. To facilitate CoVaR backtesting, Giulio and Ergün^[10] modified the A-B CoVaR to the G-E CoVaR, and measured the G-E CoVaR via the GARCH-copula method. Currently, the main methods used to estimate CoVaR include quantile regression and GARCH-copula models; see Refs. [27, 32]. In this section, we can also use our model to give an explicit equation for estimating the CoVaR. We denote the $\beta$th quantile of $X_{t,1}$ conditional on ${\rm{VaR}}_{t,2}^{\alpha}$ proposed by Ref. [10] as ${\rm{CoVaR}}_{t,1}^{\beta\mid \alpha}$, which can be seen in Eq. (11),

As mentioned earlier, we use the time-varying GAS copula function to measure ${\rm{CoVaR}}_{t,1}^{\beta\mid \alpha}$ which can be obtained by solving Eq. (12),

where $F_{x,1}\left({\rm{CoVaR}}_{t,1}^{\beta\mid \alpha} \mid {\cal{F}}_{t-1}^{1}\right)$ denotes the conditional distribution of $X_{t,1} = {\rm{CoVaR}}_{t,1}^{\beta\mid \alpha}$ given the past information ${\cal{F}}^{1}_{t-1}$. Then, we can obtain the corresponding ${\rm{CoVaR}}_{t,1}^{\beta\mid \alpha}$ by calculating the inverse of the conditional distribution using Eq. (17). In Section 5, we present the performance of the prediction method based on our model in an empirical study.

4.   The asymptotic properties of the proposed estimators

In this section, we discuss the consistency and asymptotic normality of our proposed estimators under sufficient conditions.

Assume $\left\{{\boldsymbol{X}}_{t} = (X_{t,1},X_{t,2})\right\}_{t\in \mathbb{Z}}$ is a strictly stationary Markov process. In the first stage, for $i = 1,2$, $\widehat{\eta}_{i}$ and $\widehat{{\boldsymbol{\theta}}}_{i}$ are the estimators of the univariate time series model using the $i$-dimensional univariate time series. In the second stage, $\widehat{{\boldsymbol{\theta}}}_{c}$ is the estimator of the time-varying GAS copula given the conditional cumulative distribution function. In what follows, we provide the results on consistency and asymptotic normality of the proposed estimators $\widehat{\eta}_{i}$ and $\widehat{{\boldsymbol{\theta}}}_{i}$. In this paper, $\|\cdot\|$ denotes the Euclidean norm.

4.1   Asymptotic properties of the estimators in univariate time series model

Based on the parametric sequential maximum likelihood estimation method, we can estimate all parameters more simply and efficiently. Yet the theoretical results on the fully parametric MLE seem to be few, and we extend our results in this area. To establish the asymptotic properties of the estimators of the univariate time series model proposed by our paper, we need to introduce some notations.

where $F_{0}(\cdot)$ denotes the unconditional marginal distribution with parameter $\eta_{i}$, and ${\boldsymbol{\theta}}_{i} = \left(\theta_{1,i},\dots,\theta_{p,i}\right)$ denotes the parameter set of bivariate copulas in the Dvine($p$) model.

Under regularity conditions, the parametric sequential MLE $(\widehat{\eta}_{i},\widehat{{\boldsymbol{\theta}}}_{i})$ is defined as the solution of the following equation:

Under possible model mis-specification, $(\eta^{*}_{i},{\boldsymbol{\theta}}^{*}_{i})$ is defined as the pseudo-values via

where $(\eta^{*}_{i},{\boldsymbol{\theta}}^{*}_{i})$ can be seen as the true parameters when the model is correctly specified. Then we need to impose the following conditions:

$({\rm{C1}})$ (i) $(\eta^{*}_{i},{\boldsymbol{\theta}}^{*}_{i})\in ({\cal{H}},\varTheta)$, $({\cal{H}},\varTheta)$ is the parametric space of the Dvine($p$) model, $({\cal{H}},\varTheta)$ is the compact set of ${\cal{R}}^{b}$.

(ii) $E(\varPhi_{(\eta_{i}\in {\cal{H}},{\boldsymbol{\theta}}_{i}\in\varTheta)}(X_{1,i},\cdots,X_{1+p,i})) = 0$ if and only if $\eta_{i} = \eta^{*}_{i}, {\boldsymbol{\theta}}_{i} = {\boldsymbol{\theta}}^{*}_{i}$.

$({\rm{C2}})$ (i) $\varPhi_{(\eta_{i},{\boldsymbol{\theta}}_{i})}(X_{1,i},\cdots,X_{1+p,i})$ is well defined for $(X_{1,i},\cdots,X_{1+p,i};\eta_{i},{\boldsymbol{\theta}}_{i})\in {\cal{R}}^{1+p}\times {\cal{H}} \times \varTheta$; For all $\eta_{i}\in{\cal{H}},{\boldsymbol{\theta}}_{i}\in\varTheta$, the probability that $\varPhi_{(\eta_{i},{\boldsymbol{\theta}}_{i})}(X_{1,i},\cdots,X_{1+p,i})$ is Lipschitz continuous at $(\eta_{i},{\boldsymbol{\theta}}_{i})$ is 1.

(ii) $\nabla_{\left(X_{i}\right)}^{\prime}\varPhi_{(\eta_{i},{\boldsymbol{\theta}}_{i})}(X_{1,i},\cdots,X_{1+p,i})$ is well defined and continuous in $(X_{1,i},$ $\cdots,X_{1+p,i}; \eta_{i},{\boldsymbol{\theta}}_{i})\in {\cal{R}}^{1+p}\times {\cal{H}}\times \varTheta$.

$({\rm{C3}})$ The function $\varPhi_{\eta_{i}, {\boldsymbol{\theta}}_{i}}$ is continuously differentiable with respect to $(\eta_{i}, {\boldsymbol{\theta}}_{i})$ and satisfies

for any compact set ${\cal{H}} \times \varTheta$.

$({\rm{C4}})$ ${\boldsymbol{J}}_{\eta^{*}_{i}, {\boldsymbol{\theta}}^{*}_{i}} = E\left\{-\nabla_{\left(\eta^{*}_{i}, {\boldsymbol{\theta}}^{*}_{i}\right)}^{\prime} \varPhi_{\eta^{*}_{i}, {\boldsymbol{\theta}}^{*}_{i}}\left(X_{1,i}, \cdots, X_{1+p,i}\right)\right\}$ is an invertible matrix.

$({\rm{C5}})$ ${\boldsymbol{\varSigma}}_{\eta^{*}_{i}, {\boldsymbol{\theta}}^{*}_{i}} = \lim\limits_{T\rightarrow \infty}T\cdot {\rm var}\left( {\dfrac{1}{T-p}\displaystyle \sum\limits_{t = 1}^{T-p}\varPhi_{\eta^{*}_{i},{\boldsymbol{\theta}}^{*}_{i}}(X_{t,i},\cdots,X_{t+p,i})} \right)$ is positive definite.

$({\rm{C6}})$ $\nabla_{\left(\eta^{*}_{i}, {\boldsymbol{\theta}}^{*}_{i}\right)}^{\prime} \varPhi_{\eta^{*}_{i}, {\boldsymbol{\theta}}^{*}_{i}}\left(X_{1,i}, \cdots, X_{1+p,i}\right)$ is well defined and continuous in $(X_{1,i},$ $\cdots,X_{1+p,i}; \eta_{i},{\boldsymbol{\theta}}_{i})\in {\cal{R}}^{1+p}\times {\cal{H}}\times \varTheta$.

$({\rm{C7}})$ The $\beta$-mixing coefficients of $\left(X_{t,i}\right)_{t \in \mathbb{Z}}$ satisfy $\displaystyle \sum\limits_{t = 0}^{\infty} \int_{0}^{\beta(t)} Q^{2}(u) {\rm d} u < \infty$, where $Q$ is the inverse survival function of $\left\|\varPhi_{\eta^{*}_{i}, {\boldsymbol{\theta}}^{*}_{i}}\left(X_{1,i}, \cdots, X_{t+p,i}\right)\right\|$.

Theorem 1. If the above conditions (C1)–(C3) hold for the univariate time series model, then we have $(\widehat{\eta}_{i},\widehat{{\boldsymbol{\theta}}}_{i})\xrightarrow{\rm p} (\eta^{*}_{i},{\boldsymbol{\theta}}^{*}_{i})$.

Theorem 2. If the above conditions (C1)–(C7) hold for the univariate time series model, then we have $\left\| (\widehat{\eta}_{i},\widehat{{\boldsymbol{\theta}}}_{i})-(\eta^{*}_{i},{\boldsymbol{\theta}}^{*}_{i})\right\| = O_{p}\left(T^{-1 / 2}\right)$ and

4.2   Asymptotic properties of the conditional quantile estimators

Similarly, one can show the asymptotic results of the estimator of the conditional quantile of $X_{t,i}$ under some appropriate conditions. The conditional quantile estimator of $X_{t,i}$ given ${\cal{F}}_{t-1}^{i}$ can be evaluated explicitly by the univariate Dvine model. More specifically, let $F\left(X_{t,i}\mid {\cal{F}}_{t-1}^{i} ;{\eta_{i},{\boldsymbol{\theta}}}_{i}\right)$ be the conditional distribution function parametrized by $(\eta_{i},{\boldsymbol{\theta}}_{i})$, and denote the value of our interest $\mu^{*} = \psi_{(\eta^{*}_{i},{\boldsymbol{\theta}}^{*}_{i})}\left(F\left(X_{t,i}\mid {\cal{F}}_{t-1}^{i} ;\eta^{*}_{i},{\boldsymbol{\theta}}^{*}_{i}\right)\right)$ for some functional $\psi$. Given the estimator $\left(\widehat{\eta}_{i},\widehat{{\boldsymbol{\theta}}}_{i}\right)$ of $\left(\eta_{i},{\boldsymbol{\theta}}_{i}\right)$, we define $\widehat{\mu} = \psi_{(\widehat{\eta}_{i},\widehat{{\boldsymbol{\theta}}}_{i})}\left(F\left(X_{t,i}\mid {\cal{F}}_{t-1}^{i};{\widehat{\eta}_{i},\widehat{{\boldsymbol{\theta}}}_{i}}\right)\right)$. In this paper, $\mu^{*}$ is defined as the conditional quantile.

Theorem 3. Assume the maps $\left(\eta_{i},{\boldsymbol{\theta}}_{i}\right)\rightarrow \psi_{(\eta_{i},{\boldsymbol{\theta}}_{i})}$ are Frechet differentiable at $\eta^{*}_{i},{\boldsymbol{\theta}}^{*}_{i}$ and $\left(\eta^{*}_{i},{\boldsymbol{\theta}}^{*}_{i}\right)$ with gradient $\varPsi_{\eta^{*}_{i},{\boldsymbol{\theta}}^{*}_{i}}$.

$({\rm{i}})$ If $T\rightarrow\infty$ and the above conditions (C1)–(C3) hold, then $\widehat{\mu} \xrightarrow{\rm p} \mu^{*}$.

$({\rm{ii}})$ If $T\rightarrow\infty$ and the above conditions (C1)–(C7) hold, then

The estimator $\widehat{\mu}$ of the conditional quantile converges to the pseudo-true $\mu^{*}$ at a rate $\sqrt{T}$. Theorem 3 exhibits the consistency and asymptotic normality of such predictions.

5.   Empirical analysis

5.1   Datasets and models

In this section, we illustrate the methods introduced in the previous sections with real financial datasets. The financial datasets are retrieved from Yahoo Finance and cover the period from January 1, 2015 to May 17, 2022, containing the turbulent period in financial markets caused by the outbreak of COVID-19. The first five datasets are mainly the daily closing prices of Amazon (AMAZ), Boeing (BA), Coca-Cola (COCA), General Motors (GM), and IBM companies stocks which belong to different industry branches. The last dataset is the daily closing price of the S$\&$P 500 index (SP500) which can be seen as a proxy for the financial system. There are 1856 trading days in each dataset and we take the log difference of the daily closing price as the empirical data in our analysis.

In the process of modeling the whole dataset, we first choose the Dvine($p$) model to analyze the serial dependence in the univariate time series. Table 1 shows the descriptive statistics for daily returns. The results indicate that all corresponding empirical distributions have asymmetric and heavy-tailed characteristics. By comparing different families of distributions, we employ the skewed-t distribution to fit the marginal distribution. After obtaining the pseudo-observations of individual time series, we select the order and the specific copula function for the parametric Dvine($p$) model according to the tree-by-tree selection method proposed in Ref. [14]. Table 2 exhibits the estimates for the Dvine($p$) models fitted to the financial datasets.

Table  1.  Descriptive statistics for daily returns.

	AMAZ	BA	COCA	GM	IBM	SP500
Mean	0.1062	0.0047	0.0368	0.0154	0.0099	0.0369
Skew	0.2037	−0.5538	−0.9523	−0.1728	−0.7491	−0.9664
Kurt	7.4601	19.8682	11.3986	9.5488	10.6625	18.5260
Std	1.9883	2.6450	1.1788	2.2385	1.5693	1.1578
Jarque-Bera	$<0.001$lt;0.001 $	$<0.001$lt;0.001 $	$<0.001$lt;0.001 $	$<0.001$lt;0.001 $	$<0.001$lt;0.001 $	$<0.001$lt;0.001 $
Mean refers to the average value; Skew refers to the skewness of the empirical distribution; Kurt refers to the kurtosis of the empirical distribution; Std refers to the standard deviation; Jarque-Bera refers to the p-value of the Jarque-Bera test.

 | Show Table

DownLoad: CSV

Table  2.  Summary of parameter estimates for the Dvine($p$) models.

	AMAZ	BA	COCA	GM	IBM	SP500
$\theta_{1,1}$	0.015	0.005	−0.046	0.036	−0.052	−0.065
$\theta_{1,2}$	4.654	2.755	3.411	3.914	4.121	2.535
$\theta_{2,1}$	−0.013	0.015	−0.044	−0.023	−0.021	−0.006
$\theta_{2,2}$	5.150	4.361	4.703	5.608	5.470	2.690
marginal
$\mu$	0.092	0.009	0.045	−0.007	0.023	0.0627
$\sigma$	2.199	3.568	1.271	2.491	1.760	1.5168
$\lambda$	0.966	0.987	0.948	0.981	0.963	0.9670
$v$	2.864	2.292	2.795	2.851	2.684	2.3651
For $i=1,2$, $\theta_{i,1}$ and $\theta_{i,2}$ refer to the estimators of the correlation coefficient and degrees of freedom of the t copula in tree $i$, respectively. The bottom four parameters denote estimators of the mean $\mu$, standard deviation $\sigma$, skewness $\lambda$, and degrees of freedom $v$ of the skewed t distribution.

 | Show Table

DownLoad: CSV

Certainly, the conditional CDF can be evaluated by the Dvine model. The conditional CDF $F(x_{t,i}\mid {\cal{F}}_{t-1}^{i})$ should be i.i.d. uniform if the fitted models are reasonable. Traditional modeling methods for univariate time series include GARCH-type models, we also fit the standard GARCH model proposed by Bollerslev^[4], the exponential GARCH model proposed by Nelson^[33], and the GJR-GARCH model proposed by Glosten et al^[34]. Hence, we can compare the performance of the parametric Dvine(p) model with GARCH-type models by testing the uniformity and independence of the conditional CDF. Table 3 shows the p-values of the Kolmogorov–Smirnov test and the Ljung-Box test using 10 lags. From the results of Table 3, we can see that the uniformity and independence of the conditional CDF are not rejected for those financial datasets, and the p-values of our model perform slightly better than the p-values of the GARCH-type models on both types of tests.

Table  3.  $p$-values for tests for uniformity and independence of u-PIT values.

Dataset	$n$	Dvine(p)			sGARCH			eGARCH			gjrGARCH
		$P_{\rm KS}$	$P_{ {\rm LB } }$		$P_{\rm KS}$	$P_{\rm LB}$		$P_{\rm KS}$	$P_{\rm LB}$		$P_{\rm KS}$	$P_{\rm LB}$
AMAZ	1853	$0.99$	$0.61$		$0.96$	$0.37$		$0.97$	$0.38$		$0.96$	$0.37$
BA	1853	$0.75$	$0.61$		$0.99$	$0.44$		$0.99$	$0.44$		$0.99$	$0.44$
COCA	1853	$0.89$	$0.09$		$0.77$	$0.00$		$0.77$	$0.00$		$0.77$	$0.00$
IBM	1853	$0.99$	$0.53$		$0.99$	$0.15$		$0.99$	$0.15$		$0.99$	$0.38$
SP500	1853	$0.21$	$0.71$		$0.11$	$0.25$		$0.12$	$0.02$		$0.10$	$0.28$
$P_{\rm KS}$ refers to the p-value of the Kolmogorov–Smirnov test and $P_{\rm LB}$ refer to the p-value of the Ljung-Box test using 10 lags.

 | Show Table

DownLoad: CSV

As for cross-sectional dependence, dynamic and nonlinear correlations between financial companies and the S$\&$P 500 index in financial markets can be inferred. Hence, we focus on the time-varying GAS copula functions including the Student’s t GAS copula and normal GAS copula. To show the strengths of the time-varying GAS copula functions over other types of copula functions, the time-varying DCC copula and the constant copula are fitted for the cross-sectional dependence. In Table 4, the corresponding statistics of the copula functions are listed. The minimum Akaike information criterion (AIC) values can be obtained by fitting the Student’s t GAS copula function. The results are reasonable because there exist certain upper and lower tail correlations between financial companies such as IBM and the financial system SP500. Moreover, the results show that the time-varying GAS functions are better than the constant copula and time-varying DCC copula functions. The Student’s t GAS copula function can effectively fit the dependence structure between our financial dataset. For demonstration, the time-varying average correlation among the five financial companies and the S&P 500 index estimated by the Student’s t GAS copula function is plotted in Fig. 1, which reflects a high positive correlation in financial markets.

Table  4.  Estimation results of time-varying copula function.

Dataset	Model	AIC	BIC	log-lik
AMAZ-SP500	t GAS	−1015.7	−993.58	511.84
	Normal GAS	−908.79	−886.69	458.40
	Constant	−988.24	−977.19	496.12
	t DCC	−988.00	−971.42	497.00
	Normal DCC	−922.29	−911.24	463.14
BA-SP500	t GAS	−817.73	−795.63	412.86
	Normal GAS	−773.01	−756.43	389.50
	Constant	−767.35	−756.30	385.67
	t DCC	−805.46	−788.89	405.73
	Normal DCC	−762.24	−751.19	383.12
COCA-SP500	t GAS	−573.69	−551.59	290.84
	Normal GAS	−559.88	−543.31	282.94
	Constant	−535.39	−524.35	269.69
	t DCC	−570.48	−553.90	288.24
	Normal DCC	−559.80	−548.75	281.90
GM-SP500	t GAS	−689.25	−667.15	348.62
	Normal GAS	−637.10	−653.67	329.83
	Constant	−661.44	−650.39	332.72
	t DCC	−685.37	−668.80	345.68
	Normal DCC	−656.79	−645.74	330.39
IBM-SP500	t GAS	−1075.00	−1050.90	541.52
	Normal GAS	−929.01	−945.58	475.79
	Constant	−1002.22	−991.17	503.11
	t DCC	−1038.41	−1021.83	522.20
	Normal DCC	−950.73	−939.68	477.36
The bold denotes the optimal AIC, BIC and log-likelihood.

 | Show Table

DownLoad: CSV

Figure  1.  The time-varying average correlation among the financial companies and system.

DownLoad: Full-Size Img PowerPoint

5.2   VaR and CoVaR

VaR and CoVaR as the mainstream approaches for measuring risks in financial markets are evaluated by our proposed models. In our paper, we aim to evaluate the one-day-ahead VaR and CoVaR at the $q = 5\%$ confidence level. We first calculate the $\text{VaR}_{q,t}^{i}$ for each financial dataset $i$ at time $t$. For comparison, we also consider GARCH-type models to calculate VaR for our univariate time series, including the standard GARCH model, the exponential GARCH model, and the GJR-GARCH model. Fig. 2 shows the estimation results of VaR for financial time series where a black line represents the original dataset, the VaR calculated by our model is marked with a red line, and lines in other colors represent the estimated VaR based on GARCH-type models. It shows that the estimated VaR calculated by our method can have excellent fitting results when the financial markets are in turbulence. Almost all financial datasets fluctuate sharply during the period from 2020 to 2022 due to the outbreak of COVID-19, yet our method accurately measures the risk value of financial data during periods of high volatility. As a result, the VaR building on our Dvine model is more sensitive to the occurrence of financial risk events. Thus, our method can be attractive to financial risk managers.

Figure  2.  The VaR of AMAZ, BA, COCA, GM, IBM, and SP500 at level $p = 0.05$, respectively.

DownLoad: Full-Size Img PowerPoint

Moreover, we evaluate our estimated VaR using the two-sided binomial score test, Kupiec test^[35] and Christoffersen test^[36]. Table 5 shows the results of the backtesting of VaR and compares it with other models. There are 1853 samples for backtesting and the actual number of violations where the true loss exceeds the estimated VaR are listed in column 3 (labeled $n_{\rm e}$). The p-values of the two-sided binomial score test, Kupiec test and Christoffersen test are given in columns 4, 6 and 8 (labeled $P_{\rm B}$, $P_{\rm uc}$ and $P_{\rm cc}$), respectively, and the largest ones are listed in bold. From the results of Table 5, the unconditional coverage property and conditional coverage property are not rejected at the $\alpha = 0.05$ level of significance for all datasets. The number of violations is in line with the expected number of violations under the null hypothesis that the Dvine($p$) model is the correct model. In addition, but the p-values obtained by our method are also slightly better than the results calculated by the GARCH method.

Table  5.  Backtesting of the calculated VaR at level $\alpha = 0.05$.

Dataset	Model	$n_{\rm e}$	$P_{\rm B}$	$S_{\rm uc}$	$P_{uc}$	$S_{cc}$	$P_{cc}$
AMAZ	Dvine(p)	92	0.99	0.00	0.94	0.45	0.79
	sGARCH	101	0.37	0.75	0.38	2.88	0.24
	eGARCH	91	0.92	0.03	0.86	0.57	0.77
	gjrGARCH	97	0.63	0.20	0.65	0.38	0.82
BA	Dvine(p)	93	0.96	0.00	0.97	0.11	0.95
	sGARCH	102	0.31	0.94	0.33	1.96	0.38
	eGARCH	96	0.71	0.12	0.72	0.97	0.62
	gjrGARCH	103	0.26	1.15	0.29	1.45	0.48
COCA	Dvine(p)	90	0.83	0.08	0.78	0.17	0.92
	sGARCH	89	0.75	0.16	0.69	0.80	0.41
	eGARCH	88	0.67	0.26	0.60	5.01	0.09
	gjrGARCH	88	0.67	0.26	0.60	2.04	0.36
GM	Dvine(p)	87	0.59	0.36	0.54	0.58	0.75
	sGARCH	108	0.11	2.51	0.11	5.84	0.05
	eGARCH	102	0.31	0.94	0.33	4.13	0.12
	gjrGARCH	101	0.36	0.75	0.38	4.14	0.12
IBM	Dvine(p)	90	0.83	0.08	0.77	0.12	0.94
	sGARCH	79	0.15	2.25	0.13	4.11	0.13
	eGARCH	78	0.12	3.84	0.10	6.15	0.05
	gjrGARCH	79	0.15	2.25	0.13	4.11	0.13
SP500	Dvine(p)	95	0.71	0.06	0.80	0.06	0.97
	sGARCH	107	0.24	2.49	0.14	2.98	0.18
	eGARCH	96	0.71	0.11	0.73	0.33	0.84
	gjrGARCH	98	0.56	0.30	0.58	0.44	0.80
$n_{\rm e}$ is the number of violations where the true loss exceeds the estimated VaR, $P_{\rm B}$ is p-value of two-sided binomial score test, $P_{\rm uc}$ and $S_{\rm uc}$ denote the p-values and statistics of the LR test of unconditional coverage (uc), and the subscript $\rm cc$ denotes the LR test of conditional coverage and independence (cc).

 | Show Table

DownLoad: CSV

To demonstrate the performance of our model in predicting VaR, the period of the training dataset is set from January 1, 2015 to February 1, 2022, while the period of the testing dataset is set from February 2, 2022 to May 17, 2022. Table 6 shows the results of the predicted VaR obtained by our model. The p-values of the two-sided binomial score test, Kupiec test, and Christoffersen test demonstrate the accuracy of the predicted VaR. When we calculate the mean square error (MSE) and correlation between the predicted VaR and the calculated VaR, the results show a small difference and high correlation between the predicted VaR and calculated VaR. The backtesting results prove that our predicted VaR is reliable.

Table  6.  The performance of the predicted VaR at level $\alpha = 0.05$.

Dataset	$P_{\rm B}$	$S_{\rm uc}$	$P_{\rm uc}$	$S_{\rm cc}$	$P_{\rm cc}$	MSE	Corr
AMAZ	0.41	0.51	0.47	1.27	0.53	0.0007	0.9994
BA	0.78	0.05	0.79	0.63	0.72	0.0008	0.9995
COCA	0.41	0.51	0.47	1.50	0.47	0.0006	0.9996
GM	0.41	0.36	0.54	1.58	0.43	0.0067	0.9987
IBM	0.59	0.88	0.35	1.00	0.61	0.0005	0.9997
SP500	0.78	0.04	0.83	0.52	0.77	0.0014	0.9991
$P_{\rm B}$ is p-value of two-sided binomial score test, $P_{\rm uc}$ and $S_{\rm uc}$ denote the p-values and statistics of the LR test of unconditional coverage (uc), and the subscript $\rm cc$ denotes the LR test of conditional coverage and independence (cc), MSE denotes the mean square error between predicted VaR and calculated VaR, and Corr denotes the correlation between the predicted VaR and the calculated VaR.

 | Show Table

DownLoad: CSV

Subsequently, we estimate the $\text{CoVaR}_{q,t}^{s\mid j}$ of the financial system at time $t$ conditional on the financial company $j$ being in financial distress (when the company’s return falls below its $5\%$ VaR). To compare the performance of our model, we also fit the eGARCH-t GAS copula model to evaluate CoVaR. The results of the estimated CoVaR of SP500 conditional on the financial company being in financial distress are shown in Fig. 3.

Figure  3.  The CoVaR of SP500 conditional on the financial companies AMAZ, BA, COCA, GM, and IBM being in financial distress.

DownLoad: Full-Size Img PowerPoint

By analyzing the CoVaR in Fig. 3, we find that the CoVaR sharply plunge during the COVID-19 pandemic period. The values of VaR and CoVaR obtained by our model are sensitive to the high volatility period. Its risk contagion is more significant in crisis periods using our method compared to the GARCH method. The finding reflects that companies in different industries are closely linked in the international market and that different companies have a risk contagion effect on the financial system. The magnitude is different but the trend is generally similar. To better interpret the risk spillover effect on the financial system conditional on financial company $j$ being financial distress, we define $\Delta\text{CoVaR}_{q,t}^{s\mid j}$ as Eq. (15).

where $\text{CoVaR}_{q\mid 0.5,t}^{s\mid j}$ denotes the $q$th quantile of the financial system conditional on financial company $j$ being its median state.

By evaluating the $\Delta\text{CoVaR}_{q,t}^{s\mid j}$ measure, we can more intuitively determine the risk spillover effect of the financial company $j$ on the financial system and its magnitude. The estimation results of the value of $\Delta \text{CoVaR}$ are shown in Fig. 4. The estimated $\Delta \text{CoVaR}$ can reflect the impact of the financial company on the S$\&$P500 index in a similar way to the previous volatility. Meanwhile, the results suggest that IBM and GM companies can have a greater influence on the financial system while other companies have a relatively small influence on the US financial market. The estimation results demonstrate that some industrial companies may be more prone to contribute to financial risk contagion among financial markets than Internet companies in crisis periods.

Figure  4.  The $\Delta \text{CoVaR}$ of SP500 conditional on the financial companies AMAZ, BA, COCA, GM, and IBM.

DownLoad: Full-Size Img PowerPoint

To ensure the prediction accuracy under our model, we also backtest our CoVaR according to the method proposed in Ref. [10]. For each financial company $j$, we obtain its $T$ observations when the actual loss at time $t$ is greater than its estimated $\text{VaR}_{q,t}^{j}$. Then, in those $T$ observations, we can construct a hit sequence again which compares the estimated $\text{CoVaR}_{q,t}^{s\mid j}$ at time $t$ with the actual loss of the financial system $R_{t}^{s}$,

The hit sequence returns 1 if the actual loss of the financial system is greater than its estimated $\text{CoVaR}_{q,t}^{s\mid j}$ on day $t$ in which financial company $j$ is in financial distress. Then we can backtest our CoVaR measure by using Kupiec and Christoffersen tests. Table 7 reports the test results for CoVaR. As show in Table 7, our CoVaR measure can satisfy the unconditional and conditional coverage properties at the $5\%$ confidence level. The backtest results of our method are also quite good compared to traditional eGARCH methods. This shows that the calculation of our model is accurate and efficient. Therefore, using our model enables accurate measurement of risk spillovers in extreme risk situations.

Table  7.  Backtesting of CoVaR at level $\alpha = 0.05$.

Dataset	Model	$n_{\rm e}$	$P_{\rm B}$	$S_{\rm uc}$	$P_{\rm uc}$	$S_{\rm cc}$	$P_{\rm cc}$
AMAZ-SP500	Dvine	5	0.81	0.03	0.87	1.41	0.49
	eGARCH	9	0.81	3.61	0.05	3.62	0.16
BA-SP500	Dvine	5	0.81	0.02	0.87	0.60	0.74
	eGARCH	9	0.81	3.84	0.08	4.99	0.08
COCA-SP500	Dvine	7	0.22	1.25	0.26	1.62	0.44
	eGARCH	6	0.81	0.55	0.45	1.44	0.48
GM-SP500	Dvine	4	0.99	0.03	0.86	2.14	0.35
	eGARCH	9	0.81	2.51	0.11	2.56	0.27
IBM-SP500	Dvine	5	0.81	0.05	0.81	0.65	0.72
	eGARCH	6	0.81	1.02	0.31	1.59	0.45
The note is the same as that of Table 5.

 | Show Table

DownLoad: CSV

Table 8 shows the results of the predicted CoVaR obtained by our model. As we can see from the table, our predicted CoVaR can satisfy the unconditional coverage property and conditional coverage property at the $5\%$ confidence level. The results of MSE and correlation also show the similarity between the predicted CoVaR and calculated CoVaR. Therefore, we can also obtain reliable results by using our model to estimate the conditional quantile.

Table  8.  The performance of the predicted CoVaR at level $\alpha = 0.05$.

Dataset	$P_{\rm B}$	$S_{\rm uc}$	$P_{\rm uc}$	$S_{\rm cc}$	$P_{\rm cc}$	MSE	Corr
AMAZ-SP500	0.22	0.61	0.37	1.41	0.49	0.0039	0.9889
BA-SP500	0.18	0.34	0.51	1.39	0.49	0.0180	0.9817
COCA-SP500	0.22	0.36	0.54	0.62	0.71	0.0071	0.9883
GM-SP500	0.99	0.05	0.81	1.42	0.49	0.0007	0.9824
IBM-SP500	0.10	2.56	0.09	4.11	0.13	0.0159	0.9761
The note is the same as that of Table 6.

 | Show Table

DownLoad: CSV

6.   Conclusions

Our paper proposes a dynamic bivariate Dvine model to deal with the dependence that exists in bivariate time series. We apply a class of stationary Dvine($p$) models to fit serial dependence in univariate time series and then model cross-sectional dependence via a time-varying copula model. To improve the computational efficiency, we use a stepwise estimation method to obtain the parameter estimates we need and provide the theoretical proof of the asymptotic properties of the parameter estimates of the Dvine($p$) model under sufficient conditions. Our method is a new extension to the modeling of time series, and we consider the copula-based models to forecast the VaR and CoVaR of financial time series.

Using applications to the financial stock prices, our VaR and CoVaR can provide predictions of tail risk for financial time series, especially during periods of high volatility caused by the outbreak of the novel coronavirus. Backtesting the conditional quantiles calculated by our method proves the reliability of the results and the improvement over the previous methods. Further research is needed to produce better algorithms to improve the efficiency of models on higher dimensional data.

Supporting information

The supporting information for this article can be found online at https://doi.org/10.52396/JUSTC-2023-0014.