Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
155,508 characters · 63 sections · 86 citation commands
myblue Large Skew-t Copula Models and Asymmetric Dependence in Intraday Equity Returns
{0.15cm} {0.15cm} \pagestyle{empty}
\pagestyle{plain} \setcounter{equation}{0}
Copula models are widely employed for multivariate data because they separate the selection of marginals from the copula function. When modeling financial data, the copula should possess two features. First, it should allow for asymmetric dependence, where the level of dependence varies for different quantiles of the variables. Second, it should allow for high tail dependence, where extreme values of both variables are positively or negatively dependent. Skew-t copulas can capture both features in high dimensions. However, there exist multiple variants of skew-t copulas, with no consensus on the most effective choice, and their estimation in high dimensions is challenging. This paper addresses both issues by showing that the skew-t copula implicit in the Azzalini_Capitanio_2003 distribution, combined with Bayesian variational inference, offers an attractive solution. We employ the methodology in a study using 15 minute intraday returns on 93 U.S. equities over five years, and show that pairwise asymmetric dependence varies in a complex fashion over both equity pairs and time.
A skew-t copula is that implicit in a parametric skew-t distribution, and three types have been used previously. The first was proposed by demarta_mcneil_2005 and is based on the generalized hyperbolic (GH) skew-t distribution. It is the most popular of the three in financial studies, with applications by Christoffersen_Errunza_Jacobs_Langlois_2012,creal2015high,oh2017,lucas2017 and oh_patton_2021. The second was suggested by smith_gan_kohn_2010 and is based on the distribution of sahu_dey_branco_2003, while the third proposed by kollo2010 and Yoshiba_2018 is based on the distribution of Azzalini_Capitanio_2003, which we label the “AC skew-t copula”. All three skew-t distributions have latent generative representations for which factor models can be used, and their implicit copulas are called “factor copulas” by oh2017.\footnote{This type of factor copula is not to be confused with the similarly named vine-based factor copulas of krupskii2013 and others, which are not the implicit copulas of high-dimensional skew-t distributions studied here.} In this paper we show the AC skew-t factor copula allows for higher levels of asymmetry in pairwise dependence than the other two copulas, making it an attractive choice for financial data.
Yoshiba_2018 studies maximum likelihood estimation of the AC skew-t copula in low dimensions. However, the likelihood (and therefore also posterior) has a complex geometry and its direct optimization is hard in high dimensions, which has precluded its use with large financial panels. Here we propose a scalable Bayesian solution that uses a conditionally Gaussian generative representation of the AC skew-t distribution. Latent variables are introduced and augmented with the copula parameters to provide a tractable joint (or “augmented”) posterior. However, evaluation of this augmented posterior using Markov chain Monte Carlo (MCMC) methods is prohibitively slow in high dimensions. To overcome this we develop a new variational inference estimator for evaluating the proposed augmented posterior of the AC skew-t copula in high dimensions.
Variational inference (VI) approximates a target density with a tractable density called a variational approximation (VA). The VA is calibrated by solving an optimization problem where the Kullback-Leibler divergence between it and the target density is minimized. Well designed VI methods are effective for estimating models with large numbers of parameters and/or big datasets; see Blei_Kucukelbir_McAuliffe_2017 for an introduction. But here, the complex geometry of the posterior of the AC skew-t copula parameters is difficult to approximate directly using standard VI methods. Therefore, we instead propose using the augmented posterior as the target density, and then adopt a flexible VA of the form suggested by Loaiza-Maya_Smith_Nott_Danaher_2021. To solve the optimization problem stochastic gradient methods ranganath2014 are used with re-parameterization gradients Kingma_2014. To increase the speed of the method we derive these gradients analytically. The speed and accuracy of the resulting VI method is demonstrated using simulated data.
We use our new copula methodology in two studies using 15 minute intraday financial data. The first illustrates the flexibility of the AC skew-t copula by modeling returns on two equities and the VIX index during two 40 day periods: a low market volatility period, and the high market volatility period of the COVID-19 pandemic crash. We find strong asymmetries in the pairwise dependencies, which change in direction and scale from the low to high market volatility periods; asymmetries that are not identified using the two other variants of skew-t copula.
The second example is our main study, and uses an AC skew-t factor copula to capture dependence between returns on 93 equities. Using a rolling window of width 40 trading days from January 2017 to December 2021, we make three findings. First, we show that one-step-ahead forecast densities are more accurate than those from benchmark copula models. Forecast accuracy is maximized with 10 to 15 factors, suggesting a rich factor structure is worthwhile. This is consistent with oh2017,opschoor2021 and oh_patton_2021, who find that industry or other group-based factors increase accuracy. Second, a major finding is that pairwise dependence asymmetries can differ substantially over equity pairs in addition to overall correlations. Third, we show that (absent of trading costs) portfolio selection strategies based on pairwise quantile dependencies from the copula model can improve performance relative to index returns. Our findings are consistent with evidence of cross-sectional heterogeneity in asymmetric dependence of equity returns bollerslev2022realized,ando2022quantile and its importance in risk management patton2004out,harvey2010,giacometti2021.
Skew-t factor copulas have been used previously to capture dependence in high dimensional financial data. ohpatton2013 use the implicit copula of the GH distribution, and of a mixture of GH and t distributions. creal2015high and oh_patton_2021 consider the implicit copula of the GH distribution with a single global factor, which the latter authors augment with latent group specific factors. lucas2017 considers the same copula with a block equicorrelation structure. In comparison, ours is the first application of which we are aware of the AC skew-t copula to a large financial panel. We adopt an unrestricted specification for the skew parameters of the skew-t distribution, allowing for greater flexibility in its implicit copula. Previous studies consider dynamic latent factor specifications for fitting daily or weekly financial time series. In contrast, we employ intraday data at the 15 minute frequency in a rolling window study using a static factor copula model with up to 15 global factors. Our objective is to allow for a high degree of heterogeneity in the level of asymmetric dependence over equity pairs, rather than capturing time variation at a lower data frequency.
Our paper also contributes to the literature on Bayesian estimation of skew-t copulas. creal2015high use MCMC to estimate their skew-t factor copula, where a particle filter is used to evaluate the intractable likelihood. smith_gan_kohn_2010 use MCMC with data augmentation, but for the skew-t copula implicit in the distribution of sahu_dey_branco_2003. Both approaches can also be adopted to estimate the AC skew-t copula, but they incur a much higher computational burden than the VI approach suggested here. A few recent studies have used VI methods to estimate other high-dimensional copula models. loaiza2019variational use VI to estimate vine copulas with discrete margins, while nguyen2020variational use VI to estimate factor copulas based on pair-copula constructions. However, both studies employ simple mean field or other VAs that are ineffective approximations for the complex geometry of the posterior of the AC skew-t copula.
The rest of the paper is organized as follows. Section (ref) outlines the three skew-t copulas. Section (ref) develops the variational inference method for estimating the AC skew-t factor copula, and its efficacy is demonstrated using simulated data in Section (ref). Sections (ref) and (ref) contain the two empirical studies of intraday equity returns, while Section (ref) concludes.
A copula model Nelsen_1999,Joe_2014 expresses the joint distribution function of a continuous-valued random vector $\bm{Y}=(Y_1,\ldots,Y_d)^\top$ as
Here, $\text{\boldmath$y$}=(y_1,\ldots,y_d)^\top$, $F_{Y_j}$ is the marginal distribution function of $Y_j$, and $C:[0,1]^d \rightarrow \mathbb{R}$ is a copula function that captures the dependence structure. Differentiating (ref) gives the joint density
where $f_{Y_j}=\frac{\partial}{\partial y_j} F_{Y_j}(y_j)$, $c(\text{\boldmath$u$})=\frac{\partial }{\partial \text{\boldmath$u$}}C(\text{\boldmath$u$})$ is widely called the “copula density”, and $\text{\boldmath$u$}=(u_1,\ldots,u_d)^\top$.\footnote{The notations $C(u_1,\ldots,u_d)$ and $C(\text{\boldmath$u$})$ are used interchangeably, as are the notations $c(u_1,\ldots,u_d)$ and $c(\text{\boldmath$u$})$.}
The main advantage of using a copula model is that the marginals and the copula function can be modeled separately. When the elements of $\bm{Y}$ are asset returns, their distribution features asymmetric and high tail dependence; e.g. see Patton_2006 and Christoffersen_Errunza_Jacobs_Langlois_2012. A skew-t copula is one of only a few copulas that can account for these features in high dimensions.
Let the continuous random vector $\bm{Z}=(Z_1,\ldots,Z_d)^\top\sim F_Z$, with marginals $F_{Z_1},\ldots,F_{Z_d}$. Then if $U_j=F_{Z_j}(Z_j)$, the distribution function of $\bm{U}=(U_1,\ldots,U_d)^\top$ is $C(\text{\boldmath$u$}) = F_Z(F_{Z_1}^{-1}(u_1),...,F_{Z_d}^{-1}(u_d))$. This is called an implicit copula, and the elements of $\bm{Z}$ are called pseudo or auxiliary variables because their values are unobserved; see smith2021implicit for an introduction. The copula density is
where $\text{\boldmath$z$}=(z_1,\ldots,z_d)^\top$, $f_Z(\text{\boldmath$z$})=\frac{\partial}{\partial \text{\boldmath$z$}}F_Z(\text{\boldmath$z$})$ and $f_{Z_j}=\frac{\partial}{\partial z_j}F_{Z_j}(z_j)$. A skew-t copula is where $F_Z$ is a multivariate skew-t distribution. Because there are multiple skew-t distributions, there are also multiple copulas. Three types have been used previously, which we briefly outline below. The marginal moments of $\bm{Z}$ are unidentified in $C$, so that location parameters are unnecessary and the scale matrices are restricted to be a correlation matrix $\bar\Omega$ for all three skew-t copulas.
demarta_mcneil_2005 construct the implicit copula of the generalized hyperbolic (GH) skew-t distribution barndorff1977. This distribution has the generative representation
where $\bm{X} \sim N_d(\bm{0}, \bar{\Omega}), W \sim \textrm{Gamma}(\nu/2, \nu/2)$, and $\text{\boldmath$\delta$} \in \mathbb{R}^d$ is the skewness parameter. The joint density of $\bm{Z}_{\mbox{\tiny GH}}$ is available in closed form (see Part A of the Web Appendix). However, the marginal distributions and quantile functions are not and are difficult to compute using numerical methods, so that they are usually evaluated by Monte Carlo simulation from the generative representation. For high $d$ this is slow. oh_patton_2021 note that if $\text{\boldmath$\delta$}=\delta (1,1,\ldots,1)^\top$ with $\delta$ a scalar, then the $d$ marginals of the GH distribution are the same. This speeds computation greatly, but at the cost of restricting the flexibility of the pairwise asymmetric dependencies in its implicit copula.
smith_gan_kohn_2010 construct the implicit copula of the skew-t distribution proposed by sahu_dey_branco_2003. This distribution can be constructed from a t-distribution by hidden conditioning as follows. Let $\bm{X}$ and $\bm{L}$ be $d$-dimensional random variables jointly distributed
with $D = \operatorname{diag}(\text{\boldmath$\delta$}), \text{\boldmath$\delta$} \in \mathbb{R}^d$, $\bar\Omega$ a positive definite matrix, and $t_{2d}(\bm{0},\Omega,\nu)$ denoting a t-distribution with mean zero, scale matrix $\Omega$ and $\nu$ degrees of freedom. Then, $\bm{Z}_{\mbox{\tiny SDB}} \equiv (\bm{X}|\bm{L}>\bm{0})$ \footnote{The notation $\bm{X}|\bm{L}>\bm{0}$ corresponds to the conditional distribution of $\bm{X}$ given that all elements of $\bm{L}$ are positive. It does not denote the conditional distribution of $\bm{X}$ given a specific value for $\bm{L}$. This notational convention is used throughout the paper.} is a skew-$t$ distribution with skewness parameters $\text{\boldmath$\delta$}$ and joint density $f_{Z_{\mbox{\tiny SDB}}}(\text{\boldmath$z$} ; \bar\Omega, \text{\boldmath$\delta$}, \nu)= 2^d f_t(\text{\boldmath$z$};\bm{0},\bar\Omega+D^{2}, \nu) \mbox{Pr}(\bm{V}>\bm{0};\text{\boldmath$z$})$ where $f_t(\text{\boldmath$z$};\bm{0},\bar{\Omega},\nu)$ is the density of a $t_{d}(\bm{0},\bar{\Omega},\nu)$ distribution evaluated at $\text{\boldmath$z$}$, and $\bm{V}$ follows a $d$-dimensional t-distribution with parameters that are functions of $\text{\boldmath$z$}$. sahu_dey_branco_2003 show that the $j$th marginal has density $f_{Z_{\mbox{\tiny SDB}}}(z_j;1,\delta_j,\nu)$ with distribution and quantile functions that are computed numerically. Evaluation of the copula density at (ref) for large $d$ is difficult because computation of the term $\mbox{Pr}(\bm{V}>\bm{0};\text{\boldmath$z$})$ is also. However, smith_gan_kohn_2010 show how to estimate the distribution and its implicit copula using Bayesian data augmentation.
kollo2010 and Yoshiba_2018 construct the implicit copula of the skew-t distribution proposed by Azzalini_Capitanio_2003. This distribution is formed via hidden conditioning as follows. Let $\bm{X}$ be a $d$-dimensional random vector and $L$ a random variable with joint distribution
where $\text{\boldmath$\delta$} = (\delta_1,\dots, \delta_d)^\top$. Then, $\bm{Z}_{\mbox{\tiny AC}} \equiv (\bm{X}|L > 0)$ is a skew-{\em t} distribution with joint density
where ${\cal M}(\text{\boldmath$z$})=\text{\boldmath$z$}^{\top} \bar{\Omega}^{-1} \text{\boldmath$z$}$, and $T(x;\nu)$ is the distribution function of a univariate student-t with degrees of freedom $\nu$, evaluated at $x$. There is a one-to-one relationship between $\text{\boldmath$\alpha$}$ and $\text{\boldmath$\delta$}$, given by \[ \text{\boldmath$\alpha$} = (1 - \text{\boldmath$\delta$}^{\top} \bar{\Omega}^{-1} \text{\boldmath$\delta$})^{-1 / 2} \bar{\Omega}^{-1} \text{\boldmath$\delta$}\,, \;\mbox{ and }\;\; \text{\boldmath$\delta$} = (1+\text{\boldmath$\alpha$}^{\top} \bar{\Omega} \text{\boldmath$\alpha$})^{-1 / 2} \bar{\Omega} \text{\boldmath$\alpha$}\,. \] It is more convenient to employ $\text{\boldmath$\alpha$}\in \mathbb{R}^d$ as the skewness parameter, rather than $\text{\boldmath$\delta$}$, because it is unconstrained. The $j$th marginal of (ref) has density $f_{Z_{\mbox{\tiny AC}}}(z_j; 1,\delta_j, \nu)$, and the distribution function $F_{j}(z_j;\delta_j, \nu)=\int_{-\infty}^{z_j} f_{Z_{\mbox{\tiny AC}}}(\xi; 1,\delta_j, \nu)d\xi$ is computed using numerical integration.
A draw from the AC skew-t distribution can be obtained by first generating $W\sim \mbox{Gamma}(\nu/2,\nu/2)$, followed by $\tilde L\sim N(0,1)$ constrained so that $\tilde L>0$, and then from
where $\tilde L = W^{1/2}L$, and $L$ is defined at (ref). An alternative generative representation for the AC skew-t distribution is given in Appendix (ref). To convert a draw $\bm{Z}_{\mbox{\tiny AC}}=(Z_1,\ldots,Z_d)^\top$ from the skew-t distribution to a draw $\bm{U}\sim C$ from its implicit copula, set $U_j=F_j(Z_j;\delta_j, \nu)$ for all $j$.
The motivation for using any skew-t copula over a t-copula is to capture asymmetric dependence. To measure this here we compute the four pairwise quantile dependence metrics, and define measures of asymmetry in the major and minor diagonals as follows. Let $(Y_1,Y_2)$ follow a bivariate copula model\footnote{Because all three skew-t distributions are closed under marginalization, so are the skew-t copulas. Therefore, maximum asymmetric tail dependence in the bivariate case is equal to that for variable pairs in higher dimensions.} with copula function $C(u_1,u_2)=\mbox{Pr}(U_1\leq u_1,U_2 \leq u_2)$. Then, for $0<u<0.5$, we define the four quantile dependencies as: lower left $\lambda_{\mbox{\tiny LL}}(u) = P(U_2\leq u|U_1\leq u)$, upper right $\lambda_{\mbox{\tiny UR}}(u) = P(U_2>1-u|U_1>1-u)$, lower right $\lambda_{\mbox{\tiny LR}}(u) = P(U_2 \leq u| U_1 > 1-u)$, and upper left $\lambda_{\mbox{\tiny UL}}(u) = P(U_2 > 1 - u| U_1 \leq u)$. They are computed for the AC and SDB skew-t copulas using numerical integration, and by Monte Carlo simulation for the GH skew-t copula; see Appendix (ref) for more details. For a given quantile value $u$, the metrics \[ \Delta_{\mbox{\tiny Major}}(u)\equiv \lambda_{\mbox{\tiny UR}}(u)-\lambda_{\mbox{\tiny LL}}(u)\,,\;\;\; \mbox{ and }\;\;\; \Delta_{\mbox{\tiny Minor}}(u)\equiv \lambda_{\mbox{\tiny UL}}(u)-\lambda_{\mbox{\tiny LR}}(u)\,, \] measure asymmetry in the major and minor diagonals, respectively. ando2022quantile suggest similar quantile metrics, although based on a network model, rather than a copula model.
The maximum asymmetry for each of the three copulas is obtained by solving the optimization \[ \max_{\text{\boldmath$\delta$},\rho} \left\{\Delta(u) \right\} \] for given values of $0<u<0.5$ and $\nu$. The solution is the same for $\Delta_{\mbox{\tiny Major}}(u)$ and $\Delta_{\mbox{\tiny Minor}}(u)$, and optimization is with respect to $\text{\boldmath$\delta$}=(\delta_1,\delta_2)$ and the off-diagonal element $\rho$ of the $(2\times 2)$ correlation matrix $\bar \Omega$. Figure (ref) plots the maximums against $u$. At every quantile the GH skew-t copula has a higher maximum asymmetric dependence than the SDB skew-t copula. However, the AC skew-t copula has a higher maximum level of asymmetric dependence than both alternatives for all but the lowest values of $\nu$, at which $\Delta(u)$ is almost equal for the AC and GH skew-t copulas. Figure A1 in the Web Appendix plots contours of the bivariate densities for the three copulas with maximal $\Delta_{\mbox{\tiny Major}}(0.01)$, further illustrating the strong differences between the skew-t copulas
In summary, the AC skew-t copula allows for a higher level of asymmetric dependence, making it an attractive choice. In addition, the computational demands of evaluating the marginals of the GH skew-t distribution complicates the use of its implicit copula in higher dimensions. A comparison of the differing skew-t copulas is undertaken later in our study of intraday data.
For large $d$ an unrestricted correlation matrix $\bar \Omega$ is difficult to estimate and one solution is adopt a factor model for the auxiliary variables $\bm{Z}$. Following Murray_Dunson_Carin_Lucas_2013, a static factor copula is used here, which corresponds to adopting the factorization
Here, $G=\{g_{ij}\}$ is an $(d \times k)$ loadings matrix with $k<d$, zero upper triangle and positive leading diagonal elements (i.e. $g_{ii}>0$), and to identify the parameters $D=I_d$. The diagonal matrix $V_1=\mbox{diag}(V)^{-1/2}$ normalizes $V$ to a correlation matrix. If $\widetilde{G}$ equals $G$ with the leading diagonal elements replaced by their logarithmic values, then $\bar \Omega$ is parameterized by $\mbox{vech}(\widetilde{G})$. We adopt this transformation because our VI method is applied to unconstrained real-valued model parameters.
Consider observations $\text{\boldmath$y$}_i=(y_{i1},\ldots,y_{id})^\top$, for $i=1,\ldots,n$, drawn independently from the copula model at (ref) with copula parameters $\text{\boldmath$\theta$}$. The pseudo variables $\text{\boldmath$z$}_i=(z_{i1},\ldots,z_{id})^\top$ are given by
where $F_{Z_j}$ is a function of $\text{\boldmath$\theta$}_j \subseteq \text{\boldmath$\theta$}$. If $\text{\boldmath$y$}=\{\text{\boldmath$y$}_1,\ldots,\text{\boldmath$y$}_n\}$, then from (ref) and (ref), the likelihood is \[ p(\text{\boldmath$y$}|\text{\boldmath$\theta$})=\prod_{i=1}^n \left\{p(\text{\boldmath$z$}_i|\text{\boldmath$\theta$})\prod_{j=1}^d \frac{ f_{Y_j}(y_{ij})}{f_{Z_j}(z_{ij};\text{\boldmath$\theta$}_j)} \right\}\,. \] For the AC skew-t copula, this likelihood has a complex geometry, making direct maximization challenging for large $d$; see Part D.2 of the Web Appendix for an illustration. However, a more tractable extended likelihood of the AC skew-t copula model can be obtained from the generative representation for the AC skew-t distribution given in Section (ref) as follows.
Let $\text{\boldmath$w$}=(w_1,\ldots,w_n)^\top$ be a vector $n$ draws of $W$, and $\tilde\text{\boldmath$l$}=(\tilde l_1,\ldots,\tilde l_n)^\top$ be a vector of $n$ draws of $\tilde L$, then
is the extended likelihood. The product over $j$ in (ref) is the Jacobian of the transformation from $\text{\boldmath$z$}$ to $\text{\boldmath$y$}$ at (ref). From the generative representation of the AC skew-t distribution, the joint density
where $p(w_i|\nu)$ is the density of a Gamma($\nu/2,\nu/2)$ distribution, $p(\tilde l_i|w_i)= 2\phi_1(\tilde l_i;0,1)\mathds{1}(l_i>0)$ is the density of a constrained standard normal, and $p(\text{\boldmath$z$}_i|\tilde l_i,w_i,\text{\boldmath$\theta$})=\phi_d(\text{\boldmath$z$}_i;\text{\boldmath$\delta$} \tilde l_i w_i^{-1/2},w_i^{-1}(\bar\Omega - \text{\boldmath$\delta$}\text{\boldmath$\delta$}^\top))$. Here, $\phi_m(\cdot;\text{\boldmath$\mu$},\Sigma)$ denotes the density of a $N_m(\text{\boldmath$\mu$},\Sigma)$ distribution, and the indicator function $\mathds{1}(X)=1$ if $X$ is true, and zero otherwise. The copula model likelihood can be recovered by integrating over the two latent vectors; i.e. $p(\text{\boldmath$y$}|\text{\boldmath$\theta$})=\int p(\text{\boldmath$y$},\tilde \text{\boldmath$l$},\text{\boldmath$w$}|\text{\boldmath$\theta$})\mbox{d}(\tilde \text{\boldmath$l$},\text{\boldmath$w$})$.
This extended likelihood is tractable and fast to compute. Evaluation of $F_{Z_1}^{-1},\ldots,F_{Z_d}^{-1}$ at (ref) is undertaken using the spline interpolation method outlined by Smith_Maneesoonthorn_2018 and Yoshiba_2018 that uses the closed form densities $f_{Z_1},\ldots,f_{Z_d}$. These authors show this approach is highly accurate and is scalable with respect to $n$. An alternative extended likelihood can be specified using the other generative representation given in Appendix (ref). However, we found it to be less effective than that adopted here; see Parts C.3 and D.2 of the Web Appendix.
The copula parameters with a factor decomposition for $\bar \Omega$ are $\text{\boldmath$\theta$}=\{\mbox{vech}(\widetilde{G}),\text{\boldmath$\alpha$},\nu\}$, where $\text{\boldmath$\alpha$} = (\alpha_1,...,\alpha_d)^\top$ is a one-to-one function of $\text{\boldmath$\delta$}$. Bayesian estimation uses the posterior density $p(\text{\boldmath$\theta$}|\text{\boldmath$y$})$, which is the marginal in $\text{\boldmath$\theta$}$ of the augmented posterior
where $\text{\boldmath$\psi$}\equiv\{\text{\boldmath$\theta$},\tilde \text{\boldmath$l$},\text{\boldmath$w$}\}$. The augmented posterior is the product of the extended likelihood at (ref) and the prior $p(\text{\boldmath$\theta$})=p(\mbox{vech}(\widetilde{G}))p(\text{\boldmath$\alpha$})p(\nu)$. We adopt the generalized double Pareto distribution prior for each element of $\mbox{vech}(\widetilde{G})$ suggested by Murray_Dunson_Carin_Lucas_2013. The prior $\alpha_j \sim N(0, 5^2)$, which places 99% mass on $\delta_j\in(-0.997,0.997)$. The prior for $\nu$ is constrained so that $\nu>2$, with $(\nu-2) \sim\operatorname{Gamma}(3, 0.2)$. This places 99% mass on the range $\nu\in (3.69,48.47)$.
An MCMC scheme that produces draws from the augmented posterior of the AC skew-t copula is outlined in Appendix (ref). However, it is slow and for high $d$ (including when $d=93$ as in our equity application in Section (ref)) it does not evaluate the augmented posterior in reasonable computing time. Variational inference is a faster and scalable alternative, as we now discuss.
VI approximates the augmented posterior using a density $q(\text{\boldmath$\psi$})\in {\cal Q}$ called the “variational approximation” (VA) from a family of tractable approximations ${\cal Q}$. Typically, the VA is obtained by minimizing the Kullback-Leibler divergence (KLD) between the target density $p(\text{\boldmath$\psi$}|\text{\boldmath$y$})\propto p(\text{\boldmath$y$}|\text{\boldmath$\psi$})p(\text{\boldmath$\psi$})=h(\text{\boldmath$\psi$})$ and its approximation $q(\text{\boldmath$\psi$})$. It is easily shown that this is equivalent to maximizing the Evidence Lower Bound (ELBO)
over $q\in {\cal Q}$, where the expectation is with respect to $\text{\boldmath$\psi$}\sim q$.
The selection of families ${\cal Q}$ that balance accuracy of the VA and the time required to maximize the ELBO function, is the topic of much current research. For models with a large number of latent variables, such as the augmented posterior at (ref), Loaiza-Maya_Smith_Nott_Danaher_2021 point out that assuming simple fixed form approximations (such as the widely used mean field approximation) can be very inaccurate. Instead, they suggest using a VA family of the form
where $p(\tilde \text{\boldmath$l$},\text{\boldmath$w$}|\text{\boldmath$\theta$},\text{\boldmath$y$}) $ is the conditional posterior of the latent variables, and $q^0_{\lambda}(\text{\boldmath$\theta$})$ is the density of a fixed form VA with parameters $\text{\boldmath$\lambda$}$ which we discuss further below.\footnote{In the machine learning literature a parametric density $q_\lambda$ with parameters $\text{\boldmath$\lambda$}$ is often called a fixed form density.} These authors show that when using the VA at (ref), the ELBO of the target density $p(\text{\boldmath$\psi$}|\text{\boldmath$y$})$ is exactly equal to the ELBO for the target density $p(\text{\boldmath$\theta$}|\text{\boldmath$y$})$; that is, \[ \mathcal L = E_{q_\lambda}\left[\log h(\text{\boldmath$\psi$}) - \log q_\lambda(\text{\boldmath$\psi$})\right] = E_{q^0}\left[\log\left(p(\text{\boldmath$y$}|\text{\boldmath$\theta$})p(\text{\boldmath$\theta$})\right) - \log q^0_\lambda(\text{\boldmath$\theta$})\right]\,. \] Thus, maximizing $\cal L$ is equivalent to solving the variational optimization for the parameter posterior with the latent variables $\{\tilde \text{\boldmath$l$},\text{\boldmath$w$}\}$ integrated out exactly, which makes (ref) more accurate than other choices of VA for this target density.
To maximize ${\cal L}$ with respect to $\text{\boldmath$\lambda$}$, stochastic gradient ascent (SGA) is the most frequently used algorithm ranganath2014. From an initial value $\text{\boldmath$\lambda$}^{(0)}$, SGA recursively updates
until reaching convergence. Here, $\bm{\rho}^{(t)}$ is a vector of adaptive learning rates set using the momentum method of zeiler, `$\odot$' is the Hadamard product, and $\widehat{\nabla_\lambda\mathcal{L}}$ is an unbiased estimator of the gradient, which is evaluated at $\text{\boldmath$\lambda$}=\text{\boldmath$\lambda$}^{(t)}$. The key to fast convergence and computational efficiency of this SGA algorithm is that $\widehat{\nabla_\lambda\mathcal{L}}$ has low variability and is fast to evaluate. One of the most effective ways to achieve this is to use the re-parameterization trick of Kingma_2014. For the VA at (ref), the re-parameterization is a transformation from $\text{\boldmath$\theta$}$ to $\text{\boldmath$\varepsilon$}\sim f_\varepsilon$, where $\text{\boldmath$\theta$}=\tau(\text{\boldmath$\varepsilon$},\text{\boldmath$\lambda$})$ and $\tau$ is a deterministic function. In this case, Loaiza-Maya_Smith_Nott_Danaher_2021 show that
where the expectation is with respect to the density $f_{\varepsilon,\tilde l,w}(\text{\boldmath$\varepsilon$},\tilde \text{\boldmath$l$},\text{\boldmath$w$})=f_\varepsilon(\text{\boldmath$\varepsilon$}) p(\tilde \text{\boldmath$l$},\text{\boldmath$w$}|\text{\boldmath$\theta$},y)$. Typically only a single draw from $f_{\varepsilon,\tilde l,w}$ is required to obtain a low variance estimate of the expectation in (ref), resulting in a major computational advantage.
For $q^0_\lambda$ we use a Gaussian density with a factor model covariance matrix as suggested by miller2017 and Ong_Nott_Smith_2018 with $r$ factors. This is not to be confused with the factor model used to define the copula. Ong_Nott_Smith_2018 give the transformation $\tau$ required for the re-parametrization trick, fast to compute closed form expressions for the derivatives $\frac{\partial \text{\boldmath$\theta$}^\top}{\partial \text{\boldmath$\lambda$}}$ and $\nabla_\theta \log q_\lambda^0(\text{\boldmath$\theta$})$ in (ref), along with MATLAB routines for their evaluation. The derivative \[ \nabla_\theta \log h(\text{\boldmath$\psi$})=( \nabla_{{\rm vech}(\widetilde G)}^\top \log h(\text{\boldmath$\psi$}), \nabla_{\alpha}^\top \log h(\text{\boldmath$\psi$}), \nabla_\nu^\top \log h(\text{\boldmath$\psi$}) )^\top\,, \] is specific to the target density $p(\text{\boldmath$\psi$}|\text{\boldmath$y$})$, which is the augmented posterior at (ref) here. Table (ref) provides this gradient, which is evaluated recursively from the bottom of the columns upwards. Its derivation is given in Part B of the Web Appendix, and uses the trace operator to express the gradient in a computationally efficient directional derivative form. The gradient can also be computed using automatic differentiation, although this is much slower.
Loaiza-Maya_Smith_Nott_Danaher_2021 call an approach that combines the VA at (ref) with SGA optimization and the re-parameterized gradient at (ref), a “hybrid VI” method. This is because it nests an MCMC step within a well-defined stochastic optimization algorithm for solving the VI problem. Algorithm (ref) details our proposed hybrid VI algorithm for estimating the AC skew-t copula. At Step (b) a small number of Gibbs steps that first draw from $\tilde \text{\boldmath$l$}|\text{\boldmath$w$},\text{\boldmath$\theta$},\text{\boldmath$y$}$, and then from $\text{\boldmath$w$}|\tilde \text{\boldmath$l$},\text{\boldmath$\theta$},\text{\boldmath$y$}$, are used as outlined in Appendix (ref). We demonstrate that this works well here, although other methods, such as the Hamiltonian Monte Carlo sampler of hoffman2014no, can also be used at this step. The output of the algorithm is $q^0_{\widehat{\lambda}}(\text{\boldmath$\theta$})$, which is usually called the “variational posterior”, because it is the optimal VA to the parameter posterior $p(\text{\boldmath$\theta$}|\text{\boldmath$y$})$.
We first demonstrate the efficacy of the VI method using simulated data from lower dimensional examples where the exact posterior can be calculated using MCMC to measure accuracy.
A sample of size $n=1024$ is generated from each of two AC skew-t factor copulas. The first is a low-dimensional single factor model ($d=5,k=1$; Case 1), and the second is higher dimensional five factor model ($d=30,k=5$; Case 2). The copula parameters were obtained from fitting these models to returns data; see Part C of the Web Appendix for details of the data generating processes.
For both examples, the variational posterior is compared to the exact posterior $p(\text{\boldmath$\theta$}|\text{\boldmath$y$})$ computed using the (slower) MCMC scheme in Appendix (ref). MCMC evaluates the posterior up to an arbitrary error, which we made small using a large Monte Carlo sample size. A skew-t copula with very different parameter values can have similar densities $c(\text{\boldmath$u$})$ and thus dependence structures. Therefore, estimation accuracy is best measured using pairwise dependence metrics of the copula for all possible variable pairs, rather than parameter values. The metrics computed here include the quantile dependencies computed at the 1% and 5% quantiles, and the Spearman correlation. For the pair $(Y_i,Y_j)$, the latter is $\rho^S_{ij}= 12\int \int C_{ij}(u_i',u_j')du_i'du_j'-3$, with $C_{ij}$ the bivariate marginal copula for $(Y_i,Y_j)$ from the skew-t copula evaluated as in Appendix (ref).
The exact and approximate posterior means of the metrics are evaluated by averaging their values computed at Monte Carlo draws from $p(\text{\boldmath$\theta$}|\text{\boldmath$y$})$ and $q_{\widehat{\lambda}}^0(\text{\boldmath$\theta$})$, respectively. Figure (ref) plots the exact and variational posterior means of the pairwise Spearman correlations. Figures (ref) and (ref) do the same for the quantile dependence metrics at the 5% level. They show that VI produces a copula estimate with a dependence structure that is close to that of the exact posterior; further results are given in Part C of the Web Appendix. Variational inference is an approximate method and it typically gives estimates of the second moments of the Bayesian posterior that are less accurate than the posterior mean, which we also observe here, although this has limited impact on prediction; see frazier2023variational for a discussion.
In general, we implement Algorithm (ref) using 25 Gibbs draws at step (b). To illustrate that the algorithm is robust to the number of draws, Figure (ref) plots estimates for three implementations with 1, 10 and 25 Gibbs draws. The VI estimates are very similar, which is consistent with the finding of Loaiza-Maya_Smith_Nott_Danaher_2021 for some other models.
Table (ref) compares the computation time for MCMC and VI (with 25 Gibbs draws) using a standard laptop. The diagnostic of Geweke_1991 was used to determine the number of MCMC sweeps required. The number of steps in the SGA algorithm was 5000, which is a conservative number as judged by the stability of the optimal VA; for example, in Case 1 the VA based on only 2500 steps is very similar. The VI estimator is faster, and its computational advantage grows with $d$ and $k$. Finally, in the empirical work in Section (ref) where $d=93$, we found the VI estimator stable, whereas the MCMC algorithm in Appendix (ref) exhibited poor mixing and was unable to estimate the model in reasonable time.
Our approach uses the generative representation at (ref) (labelled `GR2') to construct a variational data augmentation method. One drawback is that a nested (albeit fast) MCMC sampler to draw from $p(\tilde{\text{\boldmath$l$}},\text{\boldmath$w$}|\text{\boldmath$\theta$},\text{\boldmath$y$})$ is required at step (b) of Algorithm (ref). An alternative generative representation `GR1' in Appendix (ref) with latent vector $\text{\boldmath$l$}=(l_1,\dots,l_n)^\top$ can also be used to construct another hybrid VI algorithm that targets $p(\text{\boldmath$l$},\text{\boldmath$\theta$}|\text{\boldmath$y$})$, where at step (b) direct generation from $p(\text{\boldmath$l$}|\text{\boldmath$\theta$},\text{\boldmath$y$})$ is possible without the use of a nested MCMC sampler. We implemented this alternative hybrid VI algorithm based on GR1, but found it to be less efficient than that based on GR2; see Part C.3 of the Web Appendix for more details. Additional discussion of why GR2 is an effective target distribution is also given in Part D.2 of the Web Appendix.
The first example applies the AC skew-t copula model, with $k=1$ factor and fitted using the proposed VI method, to two equity returns series and the VIX. The latter is an index of overall market volatility published by the Chicago Board of Exchange britten2000option. The objective is to illustrate the flexibility of the dependence structure of the AC skew-t copula.
Trade prices for Bank of America (BAC) and JP Morgan (JPM) were obtained from Revinitiv DataScope Select database. For each equity, the $N=26$ fifteen minute returns between 09:30 and 16:00 EDT of each trading day are constructed as follows. If $P_{\tau,t}$ denotes the last trade price of a given equity in the $\tau$th 15 minute interval on trading day $t$, then the return for that interval is $r_{\tau,t}=\log(P_{\tau,t})-\log(P_{\tau-1,t})$, with $P_{0,t}\equiv P_{N,t-1}$ the last trade price on the previous day. Days where trading is suspended are removed.
We follow Patton_2006 and many others by using the copula to capture cross-sectional dependence in the conditional distribution of the financial variables. The intraday GARCH model of Engle_Sokalska_2012 is used as a marginal model for the returns on each equity, which decomposes the return into components
Here, $h_t = \sum^N_{\tau=1}(r_{\tau,t})^2$ is the realized variance of day $t$, $s_\tau= \frac{1}{T}\sum_{t=1}^T(r_{\tau,t})^2 / h_t$ is an estimate of the diurnal pattern in volatility for the $\tau$th interval, and the conditional volatility movement
with $\epsilon_{\tau,t}=r_{\tau,t}/(h_t s_\tau)^{1/2}$, $\epsilon_{0,t}=\epsilon_{N,t-1}$ and $\sigma_{0,t}=\sigma_{N,t-1}$. In our analysis, we assume that the innovation $\varepsilon_{\tau,t} \sim t_1(0,1,\tilde{\nu})$, and estimate the parameters $(\beta_0,\beta_1,\beta_2,\tilde \nu)$ for each equity using maximum likelihood. The copula data is computed as $T(\varepsilon_{\tau,t};\tilde{\nu})$ at the model estimates.\footnote{To match this with the copula model notation in Section (ref), note that if $i=(N-1)t+\tau$ (so that $i=1,2,\ldots,NT$) and the equity return is the $j$th marginal, then $y_{ij}=r_{\tau,t}$ and $F_{Y_j}(y_{ij})=T(\varepsilon_{\tau,t};\tilde \nu)$ in (ref).}
Fifteen minute observations of the VIX were also obtained from the Revinitiv DataScope Select database. The VIX marginal model is a first order autoregression with a nonparametric disturbance estimated using a kernel density estimator.
The copula model was fit to data from two periods of 40 trading days: a low volatility period between 24 Oct. and 22 Dec. 2017, and a high volatility period between 11 Feb. and 15 Apr. 2020.\footnote{We select these because the first period contains the lowest VIX value from 2017 to 2021, whereas the second contains the highest VIX value that corresponds to the COVID-19 crash.} Table (ref) reports the estimates of the copula parameters, correlations and asymmetry measures. The two equity returns are positively correlated, whereas they are both negatively correlated with the VIX. There are strong asymmetries in the quantile dependencies, which change in direction and scale from the low to high volatility periods. For example, between the two equity returns, $\Delta_\text{Major}(0.05)$ changes from $0.169$ to $-0.279$. This suggests that during the period of high market volatility, negative returns are more dependent than positive returns. Whereas, between the equity returns and the VIX, $\Delta_\text{Minor}(0.05)$ changes from $-0.059, -0.061$ to $0.240, 0.211$. This suggests that during the period of high market volatility there is an increase in the dependence between high VIX values and negative equity returns.
Figures (ref) and (ref) further summarize the dependence structure of the estimated copula for the low and high market volatility periods, respectively. Following ohpatton2013, panels (d,g,h) give “quantile dependence plots” for the pairs BAC-VIX, JPM-VIX and JPM-BAC, respectively. These visualize asymmetric dependence along the major diagonal by plotting the quantile dependencies $\lambda_{\mbox{\tiny LL}}(u)$ & $\lambda_{\mbox{\tiny UR}}(1-u)$ against $u$ (blue & red lines), and along the minor diagonal by plotting $\lambda_{\mbox{\tiny LR}}(u)$ & $\lambda_{\mbox{\tiny UL}}(1-u)$ versus $u$ (yellow & purple lines). The dependence between the two equity returns in the major diagonal is positively asymmetric during the low volatility period in Fig. (ref)(h), and negatively asymmetric during the high volatility period in Fig. (ref)(h). In contrast, dependence in the minor diagonal between the returns on both equities and the VIX is close to symmetric in the low volatility period in Fig. (ref)(d,g), but positively asymmetric during the high volatility period in Fig. (ref)(d,g), as also indicated by the metrics in Table (ref).
To visualize the flexibility of the predictive distributions from the copula model, we compute these for the first 15-minute trading interval of the following day of each sample, which are 23 Dec. 2017 (low volatility period) and 16 Apr. 2020 (high volatility period). Figures (ref) and (ref) plot the marginal predictive densities in panels (a,e,i), and the bivariate slices of the joint predictive density in panels (b,c,f). In the high volatility period both the tails and asymmetric dependence of the distributions are greatly accentuated, compared to the low volatility period.
For comparison, we also fit the SDB and GH skew-t copula models with the same marginals. To fit the former we used the MCMC algorithm outlined by smith_gan_kohn_2010, and for the latter we used the code provided by oh_patton_2021. For SDB and AC we use $k=1$ factors for $\bar{\Omega}$, and for the GH we use 1 global factor plus 3 group-specific factors. Part D.1 of the Web Appendix reports the dependence metrics and parameter estimates for both copulas. We observe that all copulas have similar correlation estimates, but both SDB and GH exhibit much lower asymmetry in dependence compared to AC. Moreover, during the second period $\hat\nu =4.11$ for AC, whereas $\hat \nu=20.5$ and $21.2$ for SDB and GH, so that the AC skew-t copula also captures higher extremal tail dependence.
We now apply the AC skew-t copula model to the 15 minute returns on the constituents of the S&P100 index for the period 1 Jan. 2017 to 31 Dec. 2021, which is our main application. Only the $d=93$ equities that are listed over the entire period are used in our empirical analysis. The returns were constructed as in the previous example, and the same marginal models used. Rolling windows of width $T=40$ trading days are employed, and the copula model is estimated in each window using $n=26\times 40=1040$ intraday returns for each equity. The model is re-fit every 20 trading days, resulting in 60 overlapping windows. Our empirical study focuses on three research questions. The first is whether adopting the AC skew-t copula improves the accuracy of 15 minute ahead density forecasts of portfolio returns, relative to benchmarks. The second is what is the degree of asymmetric dependence over the five year period of our data. The third question is whether, or not, exploiting asymmetry in dependence when forming investment portfolios can improve returns.
We consider 15 minute (i.e. one-step-ahead) density forecasts of the return on a market value weighted portfolio of the equities. The weights are re-calculated at the same frequency as the model estimation; i.e. every 20 trading days. The forecasts are Bayesian posterior predictive distributions for every 15 minute period from 3 Mar. 2017 to 25 Jan. 2022, giving a total of $1200 \times 26= 31,200$ forecasts. Each of these is evaluated by drawing 10,000 iterates from the posterior predictive joint distribution of returns (see Part E.1 of the Web Appendix), from which draws of the portfolio return are obtained. Kernel density estimates of these are the density forecasts of the portfolio return.
Two density forecasting metrics are calculated: the log-score (LS) and the continuous ranked probability score (CRPS) of Gneiting_Balabdaoui_Raftery_2007. Higher values of the LS and lower values of the CRPS correspond to increased accuracy. Table (ref) reports the mean of these over all portfolio return density forecasts. Results are given for different numbers of factors $k$ in the copula model, and estimated using VI where the covariance matrix of the Gaussian VA $q^0_\lambda$ has different numbers of factors $r$. Note that while setting $r=0$ corresponds to a fully factorized VA for $q^0_\lambda$, we stress that the VA $q_\lambda$ at (ref) for the target density $p(\text{\boldmath$\psi$}|\text{\boldmath$y$})$ is not of a mean field type. We make two observations. First, using higher values of $k$ (i.e. $k=10$ and $k=15$) increases accuracy. This result is consistent with oh_patton_2021, who found latent group factors, over-and-above a single global factor, improve the dependence structure in their skew-t copula model. Second, the results are similar for $r=3,5,10$, which is consistent with Ong_Nott_Smith_2018 and subsequent authors who find lower values of $r$ typically work well. We focus on results for $k=10$ factors (which corresponds to a total of 979 copula parameters in $\text{\boldmath$\theta$}$) with $r=3$ for the VA. Estimation for a single window takes 29.7 hours using 20,000 steps of the SGA algorithm on a 2018 MacBook Pro.
The predictive accuracy of the AC skew-t copula model is compared with five benchmarks models. The first four are copula models with the same marginals. Three copulas are sub-types of the AC skew-t copula; namely, the AC skew-normal copula ($\nu\rightarrow \infty$), the t-copula ($\text{\boldmath$\delta$}=\bm{0}$), and the Gaussian copula ($\nu\rightarrow \infty$ and $\text{\boldmath$\delta$}=\bm{0}$). The fourth copula is the GH skew-t copula, implemented using the code and best fitting factor structure in oh_patton_2021; see Part E.3 of the Web Appendix for implementation details. The fifth benchmark model is an intraday Constant Conditional Correlation GARCH (CCC-GARCH) Bollerslev_1990 specified as in Section (ref), but with $\varepsilon_{\tau,t}\sim N(0,1)$ and cross-asset correlation estimated using the sample correlation. We treat the CCC-GARCH model as the baseline, and compute the cumulative difference of each metric for each copula model and this baseline model. Figure (ref) plots the differences over the validation period for (a) LS and (b) CRPS. The positive slope indicates that the five copula models all out-perform the baseline model throughout the period, and the AC skew-t factor copula dominates the other factor copula models on both metrics.
If $\Delta_{\text{Major},i,j}(u)$ denotes asymmetric quantile dependence for dimensions $i,j$, then we propose
as a measure of Total Asymmetric Dependence (TAD) in the major diagonal. We set $\epsilon=0.001$ and compute the integral over all quantiles numerically. Asymmetry along the major diagonal is measured, rather than the minor diagonal, because most equity pairs exhibit positive overall dependence. When dependence is symmetric $\mathcal T_{i,j}=0$, while the maximum value of the TAD for any pair of variables in the AC skew-t copula can be computed as $\max_{\rho,\text{\boldmath$\delta$}}\{\mathcal T_{i,j}\} = 0.193$ with $\nu = 2$.
This metric is computed for all $93\times 92/2=4278$ pairs of equities and for the AC skew-t copulas fitted at each of the 60 windows. Figure (ref) summarizes the evolution of the $\mathcal T_{i,j}$ values over the estimation windows as follows. The gray shaded area is the interval between $\min_{i,j}\{\mathcal T_{i,j}\}$ and $\max_{i,j}\{\mathcal T_{i,j}\}$, while the dashed black line depicts the mean value across all equity pairs. The windows with the highest level of TAD are 18 Dec. 2018 to 15 Feb. 2019, and 11 Feb. 2020 to 15 Apr. 2020. The former window corresponds to an escalating trade war between the U.S. and China, while the second window corresponds to the height of the COVID-19 equity market crash. The $\mathcal T_{i,j}$ values for the two pairs Apple--Microsoft and Google--Tesla are also plotted, illustrating how the level of TAD can vary substantially over different equity pairs.
Finally, to further highlight the heterogeneity in asymmetric quantile dependence, Figure (ref) depicts $\Delta_{\text{Major},i,j}(0.01)$ over the final ten non-overlapping windows for the ten equities with the largest market capitalization. Positive values of $\Delta_{\text{Major},i,j}(0.01)$ (green) indicate upside tail dependence, while negative values of $\Delta_{\text{Major},i,j}(0.01)$ (red) indicate downside tail dependence. In particular, during the early stages of the COVID-19 crash in panel (h), the copula captures extreme downside tail dependence across all equity pairs, with the exception of Apple.
Dependence between equity returns is one of the major attributes of asset selection strategies for investment portfolios. Traditionally, equities that are negatively correlated are preferred because they achieve a superior efficient investment frontier. More recently, consideration has been given to “tail risk” of extreme downside losses. For example, guidolin2008 and harvey2010 proposed portfolio selection based on skewness and/or kurtosis, giacometti2021 propose a portfolio allocation strategy with penalty terms based on tail risk measures, zhao2021 proposes a new portfolio optimization method using a bivariate peak-over-threshold approach; and bollerslev2022realized made use of semi-beta risk measures, separating market upside and downside movements, for portfolio construction. We consider using pairwise asymmetric quantile dependencies from the skew-t copula to construct portfolios.
The following two strategies are considered to select five equity pairs from the 4278 possible:
The equities that make up the top five pairs according to both criteria are used to construct portfolios using market-value weights.
Table (ref) reports backtesting results from our equity selection strategies implemented for the period from 3 Mar. 2017 to 25 Jan. 2022. For comparison, we also include the market-value weighted portfolio with all 93 equities (i.e. the S&P100 index corrected for sample selection), along with a strategy selecting the five equity pairs with the most negative Pearson's correlation, as our benchmarks. All strategies are re-balanced with the most current estimates and market values every 20 trading days. The high minor tail dependence strategy performs strongly, even against the benchmark S&P100 portfolio. It has the highest Sharpe ratio and lower downside risk as measured by the Value-at-Risk (VaR) at the 5% level.\footnote{The risk-free rate used to calculate the Sharpe Ratio is the Federal Reserve Effective Rate, while higher values of VaR(5%) indicate a decrease in downside risk at this quantile.} While this study does not incorporate trading costs, it suggests that trading strategies based on quantile dependencies from our AC skew-t copula model have the potential to improve the risk profile---particularly downside risk---of portfolios. Finally, we repeat the study using the GH skew-t copula, which performed poorly.
This paper makes three main contributions. First, we show that the AC skew-t copula can capture greater asymmetry in dependence than two other types of skew-t copulas. This is important because capturing asymmetry is the sole reason to adopt a skew-t copula over the t-copula. Second, we propose a VI method that can estimate the AC skew-t copula parameters for large panels of financial data. Third, our empirical work shows that the direction and extent of asymmetric dependence can vary over both equity pair and time, and that exploiting this in portfolio formation can improve investment performance. Below, we make some additional comments on the methodology.
While VI methods are increasingly employed in econometrics (see LMDN2022,chan2022fast and gefang2023 for examples) for more complex models the choice of both target and approximating densities is crucial. Standard VI methods cannot be used to approximate the posterior $p(\text{\boldmath$\theta$}|\text{\boldmath$y$})$ because it has a complex geometry. However, the augmented posterior $p(\text{\boldmath$\psi$}|\text{\boldmath$y$})$ is well approximated using (ref) and hybrid VI. While the SGA algorithm is widely used to solve variational optimization problems (for example, see ranganath2014) it is the efficient re-parameterization gradient at (ref) that makes our VI algorithm fast.
When estimating any implicit copula, it is necessary to compute the quantile functions $F_{Z_j}^{-1}$, for $j=1,\ldots,d$, at every observation. For the GH skew-t distribution this is difficult and large Monte Carlo samples (e.g. one million draws) are typically used to compute them to a necessary level of accuracy. One way to speed this computation is to set $\text{\boldmath$\delta$}=\delta(1,1,\ldots,1)^\top$, however this restricts variation in pairwise asymmetric dependence, which is the key focus of our study. This restriction, plus the difficulty in computing the MLE, is likely to contribute to its poor relative performance in Section (ref). In contrast, for the AC skew-t distribution $F_{Z_j}^{-1}$ can be computed with high accuracy numerically, and there is no need to restrict $\text{\boldmath$\delta$}$. Most previous factor copula studies use only one or two global factors, often enriched with an industry or latent group-specific factor. The computational efficiency of our VI method allows for a larger number of global factors to be used. This proves important in our empirical work, where we find 10 global factors improves performance, compared to a smaller number.
Finally, we conclude by mentioning two possible extensions of our work. In our study of intraday returns we use a static copula with a rolling window of width 40 trading days. But for studies of lower frequency returns, allowing for time variation in the AC skew-t copula can be achieved by following previous authors and adopting a dynamic specification. A second promising extension would be to consider generalizations of the AC skew-t distribution underlying the implicit copula to provide greater flexibility in capturing tail dependence. However, in both extensions identification of the copula parameters and an effective estimation methodology would need to be considered.
\oldappendix {\appendixname Part \Alph{section}\quad}
For each skew-t copula, the pairwise quantile dependence metrics in Section (ref) are evaluated from their bivariate copula function $C(u_1,u_2)=F_{Z_1,Z_2} \left( F_{Z_1}^{-1}(u_1), F_{Z_2}^{-1}(u_2) \right)$. Computation of $F_{Z_j}^{-1}(u_j)$ is discussed in the manuscript for each distribution. The function $F_{Z_1,Z_2}$ is obtained for the SDB skew-t distribution by numerical integration. For the AC skew-t distribution it is computed as the trivariate student t distribution function $ F_{Z_{\mbox{\tiny AC}}}\left( (z_1,z_2)^\top; \text{\boldmath$\delta$},\bar\Omega, \nu \right) = 2F_t\left((z_1,z_2,0)^\top; \Omega^*, \nu\right)\,, $ where
which can be derived from (ref). The GH skew-t distribution function cannot be computed in closed form, and numerical integration of its joint density is difficult. Therefore, a kernel-based approximation based on one million Monte Carlo draws is employed.
We list below two different generative representations for $\bm{Z}_{\mbox{\tiny AC}}$ that can be derived from (ref).
The extended likelihood in Section (ref) is based on generative representation GR2.
Algorithm (ref) is an MCMC sampling scheme to evaluate the augmented posterior of the AC skew-t factor copula parameters $\text{\boldmath$\theta$}$. Steps 1 and 2 are also used in the hybrid VI Algorithm (ref).
The conditional posteriors at Steps 3, 4 and 5 can be derived from (ref), but are unrecognizable, so we employ adaptive random walk Metropolis-Hastings schemes to generate each element.
To derive the posteriors at Steps 1 and 2, note that from (ref) and (ref),
with $\text{\boldmath$\mu$}_{z,i}=\text{\boldmath$\delta$} \tilde{l}_i w_i^{-1/2}$, $\Sigma_{z,i}=w_i^{-1}(\bar\Omega - \text{\boldmath$\delta$}\text{\boldmath$\delta$}^\top)$ and $p(w_i|\nu)$ is a Gamma$(\nu/2,\nu/2)$ density.
From the above, at Step 1 the conditional posterior $p(\text{\boldmath$\tilde l$}|\text{\boldmath$w$},\text{\boldmath$\theta$},\text{\boldmath$y$})= \prod_{i=1}^n p(\tilde{l}_i|w_i,\text{\boldmath$\theta$},\text{\boldmath$y$})$, where $\tilde{l}_i|\text{\boldmath$\theta$}, w_i, \text{\boldmath$y$}_i \sim N_{+}(A^{-1}B_i, A^{-1}) $, with $A = 1 + \text{\boldmath$\delta$}^\top (\bar\Omega - \text{\boldmath$\delta$}\text{\boldmath$\delta$}^\top)^{-1} \text{\boldmath$\delta$}, B_i = w_i^{1/2} \text{\boldmath$\delta$}^\top(\bar\Omega - \text{\boldmath$\delta$}\text{\boldmath$\delta$}^\top)^{-1}\text{\boldmath$z$}_i$ and $N_+$ denotes a univariate normal distribution constrained to positive values. Similarly, at Step 2 the conditional posterior $p(\text{\boldmath$w$}|\text{\boldmath$\tilde l$},\text{\boldmath$\theta$},\text{\boldmath$y$})=\prod_{i=1}^n p(w_i|\tilde{l}_i,\text{\boldmath$\theta$},\text{\boldmath$y$})$, with
An adaptive random walk Metropolis-Hastings step is used to draw from this posterior.
In Step (b) of Algorithm (ref), a draw of $\text{\boldmath$\tilde l$},\text{\boldmath$w$}$ is obtained by repeatedly drawing from the conditionals 25 times. We stress this is fast because all demanding computations do not involve $\text{\boldmath$l$}$ and $\text{\boldmath$w$}$, so that they only need to be computed once per SGA step. We found 25 draws to be adequate, which is consistent with the findings in Loaiza-Maya_Smith_Nott_Danaher_2021 for other models.
\singlespacing
\onehalfspacing \setcounter{page}{1}
\setcounter{figure}{0} \setcounter{table}{0} \setcounter{section}{0} \setcounter{equation}{0} \setcounter{algorithm}{1} This Online Appendix has five parts:
The joint density of $\bm{Z}_{\mbox{\tiny GH}}$ is
The function $K_\beta(s)$ is the modified Bessel function of the second kind with index $\beta$.
Figure (ref) plots the bivariate copula model densities with $N(0,1)$ marginals and copula parameters which give the maximal asymmetric dependence at the 1% quantile (i.e. maximum value of $\Delta(0,1)$.)
In this part of the Web Appendix we give further details on the VI methodology for the AC skew-t copula when using the generative representation GR2 with $\text{\boldmath$\psi$}\equiv\{\text{\boldmath$\theta$},\tilde \text{\boldmath$l$},\text{\boldmath$w$}\}$. We list the priors employed and then the log augmented posterior for this case. We derive the gradient $\nabla_\lambda \log h(\text{\boldmath$\psi$})$ for this case, the result of which is given in Table 2 of the manuscript. This derivation uses the trace operator to express gradients as directional derivatives, which provides a computationally efficient expression.
The parameters $\text{\boldmath$\theta$} = \{ \mbox{vech}(\widetilde{G}), \text{\boldmath$\alpha$}, \nu \}$ are transformed to unconstrained real values because the VI algorithm is applicable in this case. Thus we defined $\mbox{vech}(\widetilde{G}) = \{G_{p,k}, \widetilde{G}_{k,k}\}, p<k$, to represent off-diagonal and diagonal values of $\mbox{vech}(G)$ respectively, $\widetilde{G}_{k,k} = \log(\operatorname{Diag}(G))$ denotes the logarithms of leading diagonal values in $P \times K$ matrix $G$ with $K \ll P$. And $\tilde{\nu} = \log(\nu - 2)$ to promise VA inferenced $\nu > 2$.
The prior distribution is defined as $p(\text{\boldmath$\theta$}) = p(G_{P,K}) p(\widetilde{G}_{K,K}) p(\text{\boldmath$\alpha$}) p(\tilde{\nu})$ with
Below is the logarithm of the augmented posterior density for the generative representation GR2.
, where $\bar\Omega = V_1^{-1/2}VV_1^{-1/2}, V = GG^\top + D, V_1 = \operatorname{diag}(V)^{-1/2}$ and $\text{\boldmath$\delta$} = (1+\text{\boldmath$\alpha$}^{\top} \bar{\Omega} \text{\boldmath$\alpha$})^{-1 / 2} \bar{\Omega} \text{\boldmath$\alpha$} $
First we note that in expressing these gradients we adopt the following notation. If $f(\text{\boldmath$x$})$ is a scalar-valued function of a column vector $\text{\boldmath$x$}$, then $\nabla_x g(\text{\boldmath$x$})=\frac{\partial g}{\partial \text{\boldmath$x$}}^\top$ which is a column vector. If the vector-valued function $g(\text{\boldmath$x$})\in \mathbb{R}^d \times 1$ and $\text{\boldmath$x$} \in \mathbb{R}^n \times 1$, then $\frac{\partial g}{\partial \text{\boldmath$x$}}$ is a $d\times n$ matrix with $(i,j)$ element $\frac{\partial g_i}{\partial x_j}$.
with $\text{\boldmath$\delta$} = (1+\text{\boldmath$\alpha$}^{\top} \bar{\Omega} \text{\boldmath$\alpha$})^{-1 / 2} \bar{\Omega} \text{\boldmath$\alpha$}$ and $\bar\Omega = V_1VV_1, V = GG^\top + D, V_1 = \operatorname{diag}(V)^{-1/2} $
By Woodbury formula and Sherman–Morrison formula, we have
where $C_1 = 1+\text{\boldmath$\alpha$}^\top \bar{\Omega} \text{\boldmath$\alpha$}$, $C_2 = \text{\boldmath$\alpha$}\text{\boldmath$\alpha$}^\top \bar\Omega $, $C_3 = C_1^{-1} C_2 C_0 + C_1^{-1} C_0 C_2^\top - C_1^{-2} C_2C_0 C_2^\top $, $C_4 = C_0 - C_3$
, where $V_2 = \operatorname{diag}(V)^{-3/2}$.
From above we could derive the gradient of $T_{G1}$
where $C_{5i} = -w_i(\bar\Omega-\text{\boldmath$\delta$} \text{\boldmath$\delta$}^\top )^{-1} \text{\boldmath$y$}_i \text{\boldmath$y$}_i^\top (\bar\Omega-\text{\boldmath$\delta$} \text{\boldmath$\delta$}^\top )^{-1} = -w_i C_0 \text{\boldmath$y$}_i \text{\boldmath$y$}_i^\top C_0$
where $C_{6i} = C_1^{-1} C_2 C_{5i} + C_1^{-1} C_{5i} C_2^\top - C_1^{-2} C_2 C_{5i} C_2^\top $, $C_{7i} = C_{5i} - C_{6i} $
From above we could derive the gradient of $T_{G2i}$
\\
where $C_{8i} = - 2\tilde{l}_i w_i^{1/2} C_1 ^{-1/2} C_0 \text{\boldmath$y$}_i \text{\boldmath$\alpha$}^\top + \tilde{l}_i w_i^{1/2} C_1 ^{-3/2} C_2 C_0 \text{\boldmath$y$}_i \text{\boldmath$\alpha$}^\top $, $C_{9i} = 2\tilde{l}_i w_i^{1/2} C_0 \text{\boldmath$y$}_i \text{\boldmath$\delta$}^\top C_0$.
where $C_{10i} = C_1^{-1} C_2 C_{9i} + C_1^{-1} C_{9i} C_2^\top - C_1^{-2} C_2 C_{9i} C_2^\top $, $C_{11i} = C_{8i} + C_{9i} - C_{10i} $
\\
where $C_{12i} = 2\tilde{l}_i^2 C_1^{(-1/2)} \text{\boldmath$\alpha$} \text{\boldmath$\delta$}^\top C_0 - \tilde{l}_i^2 C_1^{-3/2} \text{\boldmath$\alpha$} \text{\boldmath$\delta$}^\top C_0 C_2^\top$, $C_{13i} = - \tilde{l}_i^2 C_0 \text{\boldmath$\delta$} \text{\boldmath$\delta$}^\top C_0$. \\
where $C_{14i} = C_1^{-1} C_2 C_{13i} + C_1^{-1} C_{13i}^\top C_2^\top - C_1^{-2} C_2 C_{13i} C_2^\top $, $C_{15i} = C_{12i} + C_{13i} - C_{14i}$
\\
\\
\\
where $C_{16i} = - 2\tilde{l}_i w_i^{1/2} C_1^{-1/2}\bar\Omega C_0 \text{\boldmath$y$}_i + 2\tilde{l}_i w_i^{1/2} C_1^{-3/2} \text{\boldmath$\alpha$}^\top \bar\Omega C_0 \text{\boldmath$y$}_i \bar\Omega \text{\boldmath$\alpha$}$, $C_{17i} = - 2\tilde{l}_i w_i^{1/2} C_0 \text{\boldmath$y$}_i \text{\boldmath$\delta$}^\top C_0$
\\
where $C_{18i} = \tilde{l}_i^2 C_1^{-1/2} \bar\Omega C_0 \text{\boldmath$\delta$} - \tilde{l}_i^2 C_1^{-3/2} \bar\Omega \text{\boldmath$\alpha$} \text{\boldmath$\alpha$}^\top \bar\Omega C_0 \text{\boldmath$\delta$} , C_{19i} = \tilde{l}_i^2 C_0 \text{\boldmath$\delta$} \text{\boldmath$\delta$}^\top C_0$
\\
\\
, where $\psi_0(t) = \diff{\log(\Gamma(t))}{(t)}$ is the digamma function
This web appendix provides additional details for the simulation in Section 4.
In Case 1 the dimension $d=5$ and number of factors $k=1$, with parameters obtained by fitting the skew-t copula with $\nu=10$ to $n=1040$ observations of 15 minute returns in the period between 2019-Sep-11 and 2019-Nov-6. The parameter values are:
In Case 2 the dimension $d=30$ and number of factors $k=5$, with parameters (reported below) obtained by fitting the skew-t copula with $\nu=10$ to $n=1040$ observations of 15 minute returns in the period between 2019-Sep-11 and 2019-Nov-6.
$G = $ $\left[
\right]$ \\
$\nu = 10$
Figure (ref) plots the marginal posterior means of the four quantile pairwise dependence metrics at the 1% level for Case 1. On the horizontal axis is the exact posterior mean computed using MCMC, and on the vertical axis is the approximate posterior mean computed using VI. Each point in a scatterplot corresponds to a specific pairwise dependence.
Figure (ref) plots the standard deviation of the marginal posterior (i.e. computed using MCMC) and its variational approximation (i.e. computed using VI) for the pairwise Spearman correlations. It is well-known that variational approximations are less well-calibrated to higher order moments, and this is apparent for some of the pairwise Spearman correlations. The accuracy of the variational posterior standard deviations can be further improved by considering more complex variational approximations, although the trade-off typically involves greater computation.
We find the generative representation GR2 to provide an augmented posterior that is most effective to estimate using our hybrid VI method, in comparison to GR1. To illustrate, Figure (ref) plots the logarithm of the posterior density against SGA step for both VI algorithms applied to Case 1. Results are given for the Gaussian factor approximation $q_\lambda^0$ using $r=0, 3$ and $5$ factors, and the optimization converges faster by this metric for GR2. Similar results (unreported) were found for Case 2 and in the analysis of the financial data.
\FloatBarrier
To show that the variational posterior mean is a good estimator of the true DGP, we reproduce Figures 3, 4 and 5 in the manuscript, but where we plot the variational posterior means against the true values. Figure (ref) does so for the Spearman correlations of both DGPs, while Figures (ref) and (ref) do so for quantile dependencies in both cases. This demonstrates the accuracy of the Bayesian posterior mean as an estimator.
\FloatBarrier
Below are estimates of the dependence metrics and parameters for the SDB and GH skew-t copulas fit to the $d=3$ dimensional example in Section 5.
The posterior of the AC skew-t copula parameters $\text{\boldmath$\theta$}$ has a complex geometry that we illustrate here using the fit to the high volatility period data in Section 5. In this example, $d=3$ and there is $k=1$ factor, so that the loading matrix
and the transformed degrees of freedom parameter $\tilde{\nu}=\log(\nu-2)$. Therefore, $\text{\boldmath$\theta$}=\{g_{11},g_{21},g_{31},\alpha_1,\alpha_2,\alpha_3,\tilde{\nu}\}$. With the vague proper priors in Section 3.3, the posterior is dominated by the likelihood. Figure (ref) plots the univariate and bivariate marginals of this posterior distribution, evaluated exactly using MCMC applied to the GR2 generative representation. Even in this low dimension, the complex geometry of the posterior is visible, so that it is a hard target distribution to approximate directly using SGA, or even traverse directly using MCMC. This motivates the use of the more tractable augmented posterior based on the generative representation.
The complex geometry of the posterior makes the variational optimization a difficult task. We solve this problem by proposing the use of the generative representation GR2 as the target posterior. To illustrate the improvement this provides, we apply our VI approach to three different target distributions for the high volatility period data. These target distributions are:
Figure (ref) plots values of the log posterior density against step number of the SGA optimization for all three target distributions. This allows us to monitor the effectiveness of the optimization. SGA struggles to converge for the GR0 and GR1 target distributions, whereas GR2 converges quickly. In our experience, the difficulties of applying SGA to the GR0 and GR1 target distributions only multiply for higher dimensional copulas. In contrast, we find that applying SGA to the target distribution GR2 is stable, which is demonstrated in our empirical work in Section 6.
Algorithm 5.1 outlines the steps necessary to generate draws from the posterior predictive return distribution for the density forecast evaluation study in Section (ref).
The $M$ draws of $(R_{\tau,t+1}^{[1]},...,R_{\tau,t+1}^{[M]})$ represents the draw from the posterior portfolio return distribution. These predictive distributions are evaluated using density forecast evaluation metrics in Section (ref).
Figure (ref) gives the quantile dependence plots between AAPL -- MSFT and GOOGL--TSLA from the fitted AC skew-t copula, calculated at the peak asymmetry period for each respective pair.
We repeat the prediction and investment studies in Sections 6 using the GH skew-t factor copula. The same marginals were used as for the AC skew-t copula. To estimate the GH skew-t copula the EM algorithm is used, as implemented in the routine “EM_static_Gcop.m” provided by oh_patton_2021 with one global factor plus heterogeneous loadings for 21 groups. This code determines group allocation and parameter estimation iteratively.
The optimal number of groups identified by oh_patton_2021 was 21 for their analysis of daily returns on the 110 stocks that were ever part of the S&P100 between 2010 and 2019. These authors found that their results were fairly robust to between 10 and 30 groups. The main difference with our study is that we employ intraday returns on (mostly) the same stocks, so that it is likely the optimum number of groups differs somewhat from 21. We tried a lower number of groups (10) and the prediction results were slightly less accurate, but (similar to these authors) we found the results relatively unchanged between 15 and 25 factors. Therefore, we simply present the results with same number of groups as the original study for simplicity.
Table 4 in the manuscript gives the equivalent investment performance using the quantile dependence based strategies in Section 6, but using the quantile dependence estimates from the GH skew-t copula. The relatively poor risk-adjusted performance for the GH skew-t copula is likely due to poor estimates of the quantile dependencies, and is consistent with the weaker one-step-ahead forecasts using this copula observed in Figure 7 in the manuscript.