The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
81,422 characters
Spectral Dynamics and Regularization for High-Dimensional Copulas
\begin{frontmatter}
\title{Spectral Dynamics and Regularization for High-Dimensional Copulas}
\author[UvT]{Koos B. Gubbels}
\ead{[email removed]}
\author[VU]{Andre Lucas}
\ead{[email removed]}
\affiliation[UvT]{organization={Department of Econometrics and Operations Research, Tilburg University},
city={Tilburg},
country={The Netherlands}}
\affiliation[VU]{
organization={Department of Econometrics and Data Science and Tinbergen Institute, Vrije Universiteit Amsterdam},
city={Amsterdam},
country={The Netherlands}}
\begin{abstract}
We introduce a novel model for time-varying, asymmetric, tail-dependent copulas in high dimensions that incorporates both spectral dynamics and regularization.
The dynamics of the dependence matrix' eigenvalues are modeled in a score-driven way, while biases in the unconditional eigenvalue spectrum are resolved by non-linear shrinkage.
The dynamic parameterization of the copula dependence matrix ensures that it satisfies the appropriate restrictions at all times and for any dimension.
The model is parsimonious, computationally efficient, easily scalable to high dimensions, and performs well for both simulated and empirical data.
In an empirical application to financial market dynamics using 100 stocks from 10 different countries and 10 different industry sectors, we find that our copula model captures both geographic and industry related co-movements and outperforms recent computationally more intensive clustering-based factor copula alternatives.
Both the spectral dynamics and the regularization contribute to the new model's performance.
During periods of market stress, we find that the spectral dynamics reveal strong increases in international stock market dependence, which causes reductions in diversification potential and increases in systemic risk.\\[1ex]\end{abstract}
\begin{keyword}
Copulas
\sep principal components
\sep time-varying eigenvalues
\sep non-linear shrinkage
\sep high-dimensional dependence.
\end{keyword}
\end{frontmatter}
\section{Introduction}
\label{sec:intro}
High-dimensional models for dependence play an important role in quantitative risk management and finance; for an overview, see, e.g., \citet{McNeil2005}.
In such high-dimensional settings, standard multivariate densities are typically too tightly parameterized to describe the data well.
The standard solution to this problem, which is also adopted in this article, is to use the more flexible copula perspective.
Here, one splits the modeling process into two steps: a first stage where one builds univariate models for the marginal properties of each of the observed time series, and a second stage where one formulates a copula to capture the multivariate dependence structure between the different series \citep[see for instance][]{Joe2014}.
To capture financial market data well, there are three stylized facts that any good copula model should pick up, namely (i) time-variation of the dependence structure, (ii) asymmetry, and (iii) tail-dependence.
In addition, a good high-dimensional copula model should be able to deal with (iv) the increase in the number of parameters as the dimensionality grows, and (v) the biases arising in a high-dimensional context for typical dependence measures like covariance matrices or copula dependence matrices \citep[see, for instance,][]{LedoitWolf2004,LedoitWolf2012}.
Earlier literature typically deals with some, but not all of these challenges simultaneously, and a unified model dealing with all these challenges at the same time is currently lacking.
For instance, \citet{lucas2017modeling} use a dynamic, skewed and fat-tailed copula model for the time-varying dependence between 73 sovereigns and banks, but adopt a highly restrictive block-equicorrelation structure as in \citet{engle2012dynamic} to keep the number of parameters manageable in high dimensions.
\cite{Engle2019} on the other hand consider a vast-dimensional setting and use non-linear shrinkage to overcome the biases in the estimation of the long-term mean of the dynamic correlation matrix in the DCC specification of \citet{Engle2002}.
They do not take a copula perspective, however, and do not account for possible skewness in the dependence structure.
Moreover, if $d$ denotes the number of assets, their dynamics are imposed on the entire $d\timesd$ correlation matrix via the DCC specification, rather than on the dynamics of the much lower dimensional spectrum itself.
A different stream of literature use a factor-based copula approach with exogenous \citep{Creal2015,Oh2017,Opschoor2020} or endogenous \citep[via clustering;][]{Oh2023} assignment of cross-sectional units to groups with similar factor exposures.
These factor copulas are typically dynamic and fat-tailed (and sometimes skewed), whereas the dimensionality challenges are tackled via the allocation of assets to groups and the pooling of factor loadings.
Such group structures may, however, not be exact or may not even exist in particular applications, in which case a group-based factor copula model approach can become suboptimal or even biased.
Moreover, endogenously deciding on the number of clusters and factors in the context of a non-linear dynamic model can be quite challenging and become computationally prohibitively expensive in high dimensions.
In this article, we therefore propose a new high-dimensional copula model that addresses all the above challenges in one unified framework.
Our copula model includes skewness and tail dependence by starting from the generalized hyperbolic skewed $t$ copula as in, for instance, \citet{lucas2017modeling}.
We address the challenge of high-dimensionality by focusing on the spectral decomposition of the dependence matrix, modeling the dominant eigenvalues as dynamic, while keeping the remaining eigenvalues static and debiasing the eigenvalue spectrum using non-linear shrinkage techniques developed by \cite{Ledoit2022a, Ledoit2022b}.
\citet{hetland2023dynamic} study eigenvalue dynamics of a covariance matrix in a low-dimensional context and find that the score-driven dynamics of \cite{Creal2013} and \citet{Harvey2013} work well for eigenvalues.
We extend their framework (i) to a copula context with skewness and fat-tailedness, and (ii) to a high-dimensional context by including the regularization techniques required for debiasing the spectrum \citep{Ledoit2004,Ledoit2022a,Ledoit2022b}.
In contrast to \cite{hetland2023dynamic}, we explicitly distinguish between the marginal and the copula time series.
Therefore, our dynamic eigenvalues solely reflect the dependence structure, allowing us to show that market movements with increased dependence are followed by subsequent periods with high dependence.
We show in a simulation study that our copula model can recover cluster-based dependence structures with limited performance loss compared to cluster factor copulas of \citet{Oh2023} if the latter form the true data generating process (dgp).
However, when the dependence structure increasingly deviates from an exact group structure, our high-dimensional copulas pick up the dependence structure more accurately than the cluster-based factor copula alternative.
In an application to the dynamic dependence structure of global financial markets using 100 stocks from 10 different countries and 10 different industry sectors, we corroborate these results.
So far, high-dimensional copula studies have mainly focused on stocks from a single country, such as the US \citep{Oh2023}.
The large number of possible country and sector combinations in our current application complicates the detection of a proper factor structure using standard clustering approaches.
We find that high-dimensional copulas with regularized spectral dynamics perform well under these challenging circumstances.
For the empirical data considered in the application, both the dynamics of the spectrum and the regularization of the spectrum contribute to the new model's performance.
Moreover, the spectral copula dynamics reveal that in times of financial distress the stock market dependence structure stretches along the first spectral dimension, reducing diversification possibilities and increasing systemic risk.
The remainder of the paper is structured as follows.
In Section \ref{sec:Def} we introduce the modeling framework and discuss regularization of the dynamic spectra and estimation of the static parameters.
Section \ref{sec:Sim} presents simulation evidence.
Section \ref{sec:Emp} applies the new model to global financial markets.
Section \ref{sec:Concl} concludes.
\section{The model}
\label{sec:Def}
\subsection{Generalized hyperbolic skewed \textit{t} copula}
\label{subsec:GH}
Consider a vector-valued time series $\bm{y}_{t} = (y_{1,t},\ldots,y_{d,t})^\top \in \mathbb{R}^{d\times 1}$ observed for $t=1,\ldots,T$.
For each $y_{i,t}$, assume the availability of a conditional marginal model in the form of a cumulative distribution function (cdf) $F_{i}(\,\cdot\,|\mathcal{F}_{t-1})$ for a common information set $\mathcal{F}_{t} = \{\bm{y}_s\}_{s\let}$; see \citet{patton2006modelling}.
We use these marginal models to construct conditional probability integral transforms (PITs) $\bm{u}_{t} = (u_{1,t},\ldots,u_{d,t})^\top \in \mathbb{R}^{d \times 1}$ with $u_{i,t} = F_{i}(y_{i,t}\mid\mathcal{F}_{t-1})$, whose multivariate distribution is the main focus in the remainder of our analysis.
We assume the following dynamic conditional copula specification for the PITs:
\begin{align}
\label{eq:GH copula}
\bm{u}_{t} &\sim
c(\bm{u}_{t}\mid \mathcal{F}_{t-1}, \bm{\theta}_{c,t})
=
c(\bm{u}_{t}; \bm{\theta}_{c,t})
=
\frac{
g\left(
G_1^{-1}(u_{1,t}; \bm{\theta}_{c,t}),
\ldots,
G_{d}^{-1}(u_{d,t}; \bm{\theta}_{c,t})
\right)
}{
g_1\left( G_1^{-1}(u_{1,t}; \bm{\theta}_{c,t}) \right)
\cdots
g_{d}\left( G_{d}^{-1}(u_{d,t}; \bm{\theta}_{c,t}) \right)
}
,
\end{align}
where the dynamic copula parameter $\bm{\theta}_{c,t}$ summarizes all the dependence information from the information set $\mathcal{F}_{t-1}$ as needed for the copula, and where $g(\,\cdot\,)$ is an appropriate multivariate density with marginal pdfs $g_{i}(\,\cdot\,;\,\cdot\,)$ and cdfs $G_{i}(\,\cdot\,;\,\cdot\,)$ for $i=1,\ldots,d$.
As an appropriate choice for $g(\,\cdot\,;\,\cdot\,)$, we consider the multivariate generalized hyperbolic skewed $t$ copula with degrees-of-freedom parameter $\nu$, skewness parameter $\bm{\gamma}$, and scale matrix $\bm{R}_{t} = \bm{R}(\bm{\theta}_{c,t})$, as generated by
\begin{align}
\label{eq:GH pdf}
g\left( \bm{y^{\star}}_{t}; \bm{\theta}_{c,t} \right)
&=
\frac{
2^{-\nu/2+1}\nu^{\nu/2}\alpha_{t}^{(d+\nu)/2}
}{
(2\pi)^{d/2}\Gamma(\nu/2)|\bm{R}_{t}|^{1/2}
}
\frac{
e^{\bm{y^{\star}}_{t}{}^\top\bm{\beta}_{t}}
K_{(\nu+d)/2} \left(
\alpha_{t} \sqrt{\nu+\bm{y^{\star}}_{t}{}^\top\bm{R}_{t}^{-1}\bm{y^{\star}}_{t}}
\right)
}{
\left( \nu+\bm{y^{\star}}_{t}{}^\top\bm{R}_{t}^{-1}\bm{y^{\star}}_{t} \right)^{(d+\nu)/4}
}
,
\end{align}
where $\bm{y^{\star}}_{t} = (y^{\star}_{1,t},\ldots,y^{\star}_{d,t})^\top$, with $y^{\star}_{i,t} = G_{i}^{-1}\left( u_{i,t}; \bm{\theta}_{c,t} \right)=G_{i}^{-1}\left( u_{i,t}; \nu,\gamma_{i} \right)$ for $i=1,\ldots,d$, and where $\bm{\beta}_{t}=\bm{R}_{t}^{-1}\bm{\gamma}$ and $\alpha_{t} = (\bm{\gamma}^\top\bm{R}_{t}^{-1}\bm{\gamma})^{1/2}$; see also \citet{Demarta2005}.
As the marginal scales of the copula are not identified, we restrict the diagonal of $\bm{R}_{t}$ to be equal to one, such that $\bm{R}_{t}$ has a correlation matrix form.
The multivariate generalized hyperbolic skewed $t$ distribution in \eqref{eq:GH pdf} can be constructed as a mean-variance mixture
\begin{align}
\label{eq:GH mv mixture}
\bm{y^{\star}}_{t} &=
\varv_{t}\,\bm{\gamma} + \sqrt{\varv_{t}}\, \bm{R}_{t}^{1/2}\bm{z}_{t},
\end{align}
where $\bm{z}_{t} \sim {\rm N}(\bm{0}_{d}, \mathbf{I}_{d})$ is a multivariate standard normally distributed random variable, and $\varv_{t} \sim {\rm IG}(\nu/2,\nu/2)$ is Inverse Gamma distributed and independent of $\bm{z}_{t}$.
Here $\bm{R}_{t}^{1/2}$ denotes the symmetric root of $\bm{R}_{t}$.
The above distribution has been a popular choice to study dynamic dependence structures in financial markets \citep{Lucas2014,Opschoor2020, Oh2023}.
The distribution is both flexible and analytically tractable and accommodates skewness as well as fat tails.
It is immediately clear from \eqref{eq:GH mv mixture} that the marginal distribution of $y^{\star}_{i,t}$ is univariate skewed $t$ with shape parameter $\nu$, skewness parameter $\gamma_i$, and scale parameter 1.
The marginal pdfs $g_i(\,\cdot\,;\,\cdot\,)$ in \eqref{eq:GH copula} are thus known analytically.
\subsection{Spectral dynamics}
\label{subsec:dynamics and shrinkage}
The core challenge for a high-dimensional time-varying copula model is the curse of dimensionality and the explosion of the number of free parameters in $\bm{R}_{t}$.
Different solutions have been proposed.
\citet{Oh2018} and \citet{Opschoor2020} study factor copulas where $\bm{R}_{t}$ is decomposed into a low-rank matrix plus a diagonal matrix, both of which can vary over time.
Further parsimony can be imposed by pooling the dynamic parameters across pre-specified groups, for example, industries, or by choosing the groups in a data-driven way, e.g., using clustering techniques as in \citet{Oh2023}.
Clustering in the context of a non-linear model can quickly become time-consuming if either the sample size or the number of assets becomes large.
Moreover, clustering techniques may face challenges if the cluster structure is only approximate; see Section~\ref{sec:Sim} for the effects of such deviations in a controlled setting.
To overcome these issues, we do not impose a group structure, but instead use a normalized spectral decomposition of the copula dependence matrix $\bm{R}_{t}$ by imposing
\begin{align}
\label{eq:R parameterization}
\bm{R}_{t}
&=
\operatorname{diag}\left(
\bm{W} \bm{\Lambda}_{t} \bm{W}^\top
\right)^{-1/2}\
\bm{W} \bm{\Lambda}_{t} \bm{W}^\top
\operatorname{diag}\left(
\bm{W} \bm{\Lambda}_{t} \bm{W}^\top
\right)^{-1/2}
,
\end{align}
where $\bm{W}$ is orthogonal and $\bm{\Lambda}_{t}$ is time-varying and diagonal with strictly positive entries.
Note that the parameterization of $\bm{R}_{t}$ automatically ensures that $\bm{R}_{t}$ has the format of a correlation matrix as long as the diagonal elements of $\bm{\Lambda}_{t}$ are strictly positive.
The latter can be ensured by using an exponential link function for the diagonal elements of $\bm{\Lambda}_{t}$.
Alternatives for parameterizing dynamic correlation matrices are, for example, the hypersphere parameterization of \citet{Jaeckel2000} as used in \citet{Creal2011} and \citet{buccheri2021score}, or the log-correlation matrix parameterization of \citet{ArchakovHansen2021} as used in \citet{HafnerWang2023}.
By concentrating on the eigenvalues, the number of free dynamic parameters is considerably less than in these alternative parameterizations, which helps for the model's tractability in high-dimensional settings.
Moreover, the spectral parameterization of the correlation matrix in \eqref{eq:R parameterization} still allows for explicit derivative expressions with respect to all nonzero elements in $\bm{\Lambda}_{t}$.
We describe the dynamics of $\bm{\Lambda}_{t}$ using the score-driven approach of \citet{Creal2013} and \citet{Harvey2013},
\begin{align}
\label{eq:first score eq}
\bm{f}_{t+1} &= \bm{\omega} + \bm{B}\,\bm{f}_{t} + \bm{A}\,\bm{\nabla}_{t},
\end{align}
where $\bm{\nabla}_{t} = \partialg(\bm{y^{\star}}_{t};\bm{\theta}_{c,t})/\partial\bm{f}_{t}$ is the score of the copula density, $\bm{f}_{t} = \log\bm{\lambda}_{t}$, and where we use unit scaling as defined by \citet{Creal2013}.
The relevant score equations for the new parameterization from \eqref{eq:R parameterization} in the copula setting are given by the following proposition, the proof of which can be found in the appendix.
\begin{proposition}[\textbf{score equations for spectral dynamics}]
\label{prop:score dynamics equations}
Let $c\left( \bm{u}_{t}; \bm{\theta}_{c,t} \right)=c\left( \bm{u}_{t}; \bm{R}_{t}, \nu, \bm{\gamma} \right)$ be the skewed $t$ copula density from \eqref{eq:GH copula} and \eqref{eq:GH pdf}, let $f_{i,t} = \log\lambda_{i,t}$, $\lambda_{i,t} = \Lambda_{i,i,t}$ and $\tilde\alpha_{t} = \alpha_{t}\sqrt{\nu + \bm{y^{\star}}_{t}{}^\top\bm{R}_{t}^{-1}\bm{y^{\star}}_{t}}$ with $\alpha_{t} = \sqrt{\bm{\gamma}^\top\bm{R}_{t}^{-1}\bm{\gamma}}$.
Then, using $\bm{\nabla}_{t} = (\nabla_{1,t},\ldots,\nabla_{d,t})^\top = \partial \log c\left( \bm{u}_{t}; \bm{\theta}_{c,t} \right)/\partial \bm{f}_{t}$, we have that the unit-scaled score dynamics are given by Eq.~\eqref{eq:first score eq}, with
\begin{align}
\label{eq:GH score expression}
\nabla_{i,t} &=
\frac{\partial \log c(\bm{u}_{t}; \bm{R}_{t},\nu,\bm{\gamma})}{\partial f_{i,t}}
=
-\tfrac12\
\frac{\partial \log|\bm{R}_{t}|}{\partial f_{i,t}}
-\tfrac12\
\bm{y^{\star}}_{t}{}^\top
\frac{\partial\bm{R}_{t}^{-1}}{\partialf_{i,t}}
\left( \tilde\nu_{t}\ \bm{y^{\star}}_{t} - 2\bm{\gamma} \right)
\\
\nonumber
&\qquad\qquad
+
\tfrac12\,
\tilde\alpha_{t}\cdot
k'_{(\nu+d)/2}\left( \tilde\alpha_{t} \right) \times
\Bigg(
\frac{\bm{\gamma}^\top
\left( \partial\bm{R}_{t}^{-1}/\partialf_{i,t} \right)\bm{\gamma}}{\bm{\gamma}^\top \bm{R}_{t}^{-1}\bm{\gamma}}
+
\frac{\bm{y^{\star}}_{t}{}^\top \left( \partial\bm{R}_{t}^{-1}/\partialf_{i,t} \right)\bm{y^{\star}}_{t}}{\nu+\bm{y^{\star}}_{t}{}^\top \bm{R}_{t}^{-1}\bm{y^{\star}}_{t}}
\Bigg)
,
\end{align}
where $\tilde\nu_{t} = (d + \nu)/(\nu + \bm{y^{\star}}_{t}{}^\top\bm{R}_{t}^{-1}\bm{y^{\star}}_{t})$, and
\begin{align}
\nonumber
&\bm{R}_{t} =
\operatorname{diag}(\bm{\Sigma}_{t})^{-1/2}\ \bm{W} \bm{\Lambda}_{t} \bm{W}^\top\ \operatorname{diag}(\bm{\Sigma}_{t})^{-1/2},
\qquad
\bm{\Sigma}_{t} = \bm{W} \bm{\Lambda}_{t} \bm{W}^\top,
\qquad
\\
\nonumber
&\frac{\partial \log|\bm{R}_{t}|}{\partial f_{i,t}}
= 1-\lambda_{i,t} \sum_{j=1}^{d}
\frac{w^2_{j,i}}{\Sigma_{j,j,t}},
\qquad
\dot\bm{\Sigma}_{i,t} = \tfrac12\lambda_{i,t}
\operatorname{diag}\left(
\frac{w_{1,i}^2}{\Sigma_{1,1,t}}
, \ldots,
\frac{w_{d,i}^2}{\Sigma_{d,d,t}}
\right),
\\
\nonumber
&\frac{\partial\bm{R}_{t}^{-1}}{\partialf_{i,t}}
=
-\operatorname{diag}(\bm{\Sigma}_{t})^{1/2}\, \frac{\bm{w}_{i}\bm{w}_{i}^\top}{\lambda_{i,t}}\, \operatorname{diag}(\bm{\Sigma}_{t})^{1/2}
+
\dot\bm{\Sigma}_{i,t}\,\bm{R}_{t}^{-1} +
\bm{R}_{t}^{-1}\,\dot\bm{\Sigma}_{i,t},
\end{align}
and $k_{\nu}'(x) = \partial \log \left( x^{\nu}\cdot K_{\nu}(x) \right)/\partial x$ for given $\nu > 0$ and $x\in\mathbb{R}^+$, where $\bm{w}_{i}$ denotes the $i$th column of $\bm{W}$ and $w_{j,i}$ its $(j,i)$th element, and $\Sigma_{j,j,t}$ denotes the $(j,j)$th element of $\bm{\Sigma}_{t}$.
\end{proposition}
The eigenvalue dynamics in Proposition~\ref{prop:score dynamics equations} are substantially different from the Gaussian and Student's $t$ based eigenvalue dynamics in \citet{hetland2023dynamic} in at least three respects: the dynamics in Proposition~\ref{prop:score dynamics equations} relate to (i) a high-dimensional copula setting, (ii) a correlation matrix rather than a covariance matrix parameterization, and (iii) a more general class of distributions.
Still, the dynamics of $\bm{\lambda}_{t}$ given in Proposition~\ref{prop:score dynamics equations} have an intuitive form.
To see this, we introduce a scaled and rotated version of $\bm{y^{\star}}_{t}$, namely $\tilde{\bm{y}}^{\star}_{t} = \bm{W}^\top\operatorname{diag}(\bm{\Sigma}_{t})^{1/2}\bm{y^{\star}}_{t}$.
We also define the corresponding rescaled and rotated version of $\bm{\gamma}$, namely $\tilde{\bm{\gamma}}_{t} = \bm{W}^\top\operatorname{diag}(\bm{\Sigma}_{t})^{1/2}\bm{\gamma}$.
We can now rewrite \eqref{eq:GH score expression} as
\begin{align}
\nabla_{i,t}
=&
\tfrac12\,\left(
\tilde\nu_{t}
\frac{\tilde y^{\star}_{i,t}{}^2}{\lambda_{i,t}}
- 1
\right)
-
\frac{\tilde{\gamma}_{i,t}\tilde y^{\star}_{i,t}}{\lambda_{i,t}}
-\tfrac12\,
\tilde\alpha_{t}\cdot
k'_{(\nu+d)/2}\left( \tilde\alpha_{t} \right) \cdot
\left(
\frac{\tilde{\gamma}_{i,t}^2/\lambda_{i,t}}{\tilde{\bm{\gamma}}_{t}^\top\bm{\Lambda}_{t}^{-1}\tilde{\bm{\gamma}}_{t}}
+
\frac{\tilde y^{\star}_{i,t}{}^2/\lambda_{i,t}}{\nu+\tilde{\bm{y}}^{\star}_{t}{}^\top\bm{\Lambda}_{t}^{-1}\tilde{\bm{y}}^{\star}_{t}}
\right)
\nonumber \\
&-\,\left(
\tilde\nu_{t}\,
\tilde{\bm{y}}^{\star}_{t}{}^\top\bm{\Lambda}_{t}^{-1}\bar{\bm{y}}^{\star}_{i,t}
-
\operatorname{trace}(\dot\bm{\Sigma}_{i,t})
\right)
+
\left(
\tilde{\bm{\gamma}}_{t}{}^\top\bm{\Lambda}_{t}^{-1}\bar{\bm{y}}^{\star}_{i,t}
+
\bar{\bm{\gamma}}_{i,t}^\top\bm{\Lambda}_{t}^{-1}\tilde{\bm{y}}^{\star}_{t}
\right) \nonumber \\
&+
\tilde\alpha_{t}\cdot
k'_{(\nu+d)/2}\left( \tilde\alpha_{t} \right) \cdot
\left(
\tfrac{\tilde{\bm{\gamma}}_{t}^\top\bm{\Lambda}_{t}^{-1}\bar{\bm{\gamma}}_{i,t}}{\tilde{\bm{\gamma}}_{t}^\top\bm{\Lambda}_{t}^{-1}\tilde{\bm{\gamma}}_{t}}
+
\tfrac{\tilde{\bm{y}}^{\star}_{t}{}^\top\bm{\Lambda}_{t}^{-1}\bar{\bm{y}}^{\star}_{i,t}}{\nu+\tilde{\bm{y}}^{\star}_{t}{}^\top\bm{\Lambda}_{t}^{-1}\tilde{\bm{y}}^{\star}_{t}}
\right)\label{eq:GH score expression rewrite}
,
\end{align}
where $\bar{\bm{y}}^{\star}_{i,t}$ and $\bar{\bm{\gamma}}_{i,t}$ are similar rescaled and rotated versions of $\bm{y^{\star}}_{t}$ and $\bm{\gamma}$, respectively, as defined in the proof of \eqref{eq:GH score expression rewrite} in \ref{app:proofs}.
The score in \eqref{eq:GH score expression rewrite} consists of two main parts: terms 1 to 3, and terms 4 to 6.
The first term holds the familiar weighted EGARCH-like volatility dynamics $\tilde\nu_{t}\/\tilde y^{\star}_{i,t}{}^2/\lambda_{i,t} - 1$ that is familiar from the literature on score-driven models for a time-varying log-variance \citep[see][]{Creal2013,Harvey2013}.
This is intuitive, as the $i$th eigenvalue is the variance of the the $i$th spectral projection $\tilde y^{\star}_{i,t}$.
The weighting factor $\tilde\nu_{t}$ causes extreme observations to have less impact on the eigenvalue dynamics for finite $\nu$, and it collapses to unity if $\nu\to\infty$.
The second and third term of \eqref{eq:GH score expression rewrite} are due to the skewness of the copula specification.
They vanish if $\tilde{\gamma}_{i,t}$ equals zero.
If $\tilde{\gamma}_{i,t} > 0$, then a large positive $\tilde y^{\star}_{i,t}$ has a smaller impact on the next eigenvalue $\lambda_{i,t+1}$ than a large negative $\tilde y^{\star}_{i,t}$.
This is because large positive outcomes are more likely under positive skewness and are, therefore, not attributed to volatility increases along the corresponding spectral dimension.
For negative skewness, the opposite holds.
The second main part of the score in \eqref{eq:GH score expression rewrite}, terms 4--6, stem from our parameterization of $\bm{R}_{t}$ as a correlation matrix.
These three terms mimic the first three terms, but with an opposite sign.
We first see a weighted EGARCH type term $(\tilde\nu_{t}\tilde{\bm{y}}^{\star}_{t}{}^\top\bm{\Lambda}^{-1}\bar{\bm{y}}^{\star}_{i,t} - \operatorname{trace}(\dot\bm{\Sigma}_{i,t})$, followed by a term related to the skewness, and a term related to the combination of skewness and kurtosis.
These additional, more complex terms adjust the dynamics of the eigenvalues to account for the fact that $\bm{R}_{t}$ has unit diagonal elements by construction.
They therefore involve all the spectral projections in $\tilde{\bm{y}}^{\star}_{t}$ simultaneously, combined with a related spectral projection $\bar{\bm{y}}^{\star}_{i,t}$ that includes an additional scaling by $\dot\bm{\Sigma}_{i,t}$.
Despite its more complex expression, the score is still easy to compute and available in analytical form and decomposable into different terms that are attributable to the features of the skewed $t$ copula.
Due to the explicit normalization of $\bm{R}_{t}$ in Eq.~\eqref{eq:R parameterization} and Proposition~\ref{prop:score dynamics equations}, $\bm{R}_{t}$ is a proper correlation matrix for any update of $\bm{f}_{t}$.
The number of time-varying parameters in $\bm{R}_{t}$, however, is considerably less than in a full dynamic hypersphere or log-correlation matrix parameterization as in \citet{Creal2011}, \citet{buccheri2021score}, \citet{ArchakovHansen2021}, or \citet{HafnerWang2023}.
Also note that not all parameters in the parameterization of $\bm{R}_{t}$ in \eqref{eq:R parameterization} can be identified simultaneously.
For instance, both $\bm{\lambda}_{t}$ and $k\cdot\bm{\lambda}_{t}$ for some constant $k>0$ give the same copula dependence matrix $\bm{R}_{t}$ in \eqref{eq:R parameterization}.
We solve this by restricting the elements of $\bm{\omega}$ in the recurrence relation for $\bm{f}_{t+1}$ later on using a targeting approach.
The spectral copula approach of Proposition~\ref{prop:score dynamics equations} comes with several advantages.
First, the approach is purely data-driven, as opposed to using a pre-defined group classification based on, for instance, industries as in \citet{Oh2018} and \citet{Opschoor2020}.
Second, the spectral decomposition is computationally much less demanding compared to a full clustering-based approach as in, for example, \citet{Oh2023}.
Third, the approach imposes no restrictions on the dependence matrix.
Fourth, the spectral approach facilitates an exploratory phase of the modeling process.
For instance, using an initial guess of $\nu$ and $\bm{\gamma}$ (such as the standard normal $\bm{\gamma} = \bm{0}$ and $\nu\to\infty$), the inverse cdf transforms $\bm{y^{\star}}_{t}$ of the PITs $\bm{u}_{t}$ immediately lead to an initial estimate of $\hat{\bm{W}}$ based on the unconditional correlation matrix of $\bm{y^{\star}}_{t}$.
This, in turn, allows us to construct preliminary estimates of the spectral projections of the data in $\tilde y^{\star}_{i,t}$, which one can inspect visually to get an impression of the extent of time-variation in the volatility of $\tilde y^{\star}_{i,t}$, i.e., in $\lambda_{i,t}$.
We can use such information to impose further parsimony on the model, e.g., by restricting some of the spectral dimensions to have constant rather than time-varying volatility.
We come back to this in Section \ref{subsec:data}.
\subsection{Non-linear shrinkage and model selection}
\label{subsec:nonlinear shrinkage}
In high-dimensions it becomes challenging to estimate the copula correlation matrix.
It is known that for large dimension $d$ the sample eigenvalues $\{\hat\bm{\lambda}_{t}\}_{t=1}^{d}$ become a biased estimator of the true spectrum $\{\bm{\lambda}_{t}\}_{t=1}^{d}$ when $d,T\to\infty$ and $d/T$ converges to some positive, nonzero constant \citep{marcenkopastur1967,Johnstone2001Spiked,LedoitWolf2012,Ledoit2022a,Ledoit2022b}.
\citet{Engle2019} use shrinkage techniques to correct these biases by adjusting the (targeted) high-dimensional intercept in their DCC transition equation.
They show that such a shrinkage procedure for the $d\timesd$ volatility intercept provides better results than a targeting procedure without shrinkage.
We follow a similar procedure, using the de-biased spectrum based on the more recent quadratic shrinkage techniques of \citet{Ledoit2022b} as our target.
If $\hat\bm{\lambda}_{t}$ denotes the sample spectrum at time $t$, the quadratic shrinkage formula $ f_{\rm QS}$ is given by
\begin{align}
\check{\lambda}_{i,t}^{-1}
&=
f_{\rm QS}(\hat\bm{\lambda}_{t}^{-1})
=
\left( 1-q \right)^2\hat\lambda_{i,t}^{-1}
+
2q\left( 1-q \right) \hat\lambda_{i,t}^{-1}\frac1{d}
\sum_{j=1}^{d}
\hat\lambda_{j,t}^{-1}
\frac{
\hat\lambda_{j,t}^{-1} - \hat\lambda_{i,t}^{-1}
}{
(\hat\lambda_{j,t}^{-1} - \hat\lambda_{i,t}^{-1})^2
+
h^2\hat\lambda_{j,t}^{-2}
}
\nonumber
\\
\label{eq:Shrink}
& \qquad
+
q^2\hat\lambda_{i,t}^{-1}
\left(
\left[
\frac1{d} \sum_{j=1}^{d}
\hat\lambda_{j,t}^{-1}
\frac{
\hat\lambda_{j,t}^{-1} - \hat\lambda_{i,t}^{-1}
}{
(\hat\lambda_{j,t}^{-1} - \hat\lambda_{i,t}^{-1})^2
+
h^2\hat\lambda_{j,t}^{-2}
}
\right]^2
+
\left[
\frac1{d} \sum_{j=1}^{d}
\lambda_{j,t}^{-1}
\frac{
h\,\hat\lambda_{j,t}^{-1}
}{
(\hat\lambda_{j,t}^{-1} - \hat\lambda_{i,t}^{-1})^2
+
h^2\hat\lambda_{j,t}^{-2}
}
\right]^2
\right)
,
\end{align}
where $\check{\bm{\lambda}}_{t} = \big(\check{\lambda}_{1,t},\ldots,\check{\lambda}_{d,t}\big)^\top$ denotes the de-biased spectrum, $q=d/T$, and $h$ is a smoothness parameter.
\cite{Ledoit2022b} show that this type of shrinkage mapping is asymptotically optimal under various loss metrics, such as Frobenius loss.
They advise to take $h = \min(q^2, q^{-2})^{0.35}\cdot d^{-0.35}$ and demonstrate that their non-linear shrinkage estimator has excellent in-sample and out-of-sample performance in a wide range of settings.
In combination with the above shrinkage approach to overcome eigenvalue biases in high dimensions, we also impose parsimony on the model by limiting the number of dynamic eigenvalues.
We do so using standard model-selection criteria.
In applications to stock return data as in Section~\ref{sec:Emp}, one typically finds that the first few eigenvalues are much larger than the remaining ones.
The notion of a few dominant eigenvalues that drive the dependence structure is conceptually in line with the perspective of a `spiked' eigenvalue dependence model \citep[see, e.g.,][]{Johnstone2001Spiked,FanLiaoMincheva2013POET,donoho2018optimal} and corresponds closely to the typical covariance structure present in financial markets \citep[see, for instance,][and many more]{FamaFrench1993,FamaFrench1998International}.
Capturing the dynamics of the largest eigenvalues is therefore most important and contributes most to the model's fit.
For the lowest eigenvalues, it is more important to prevent them from becoming too close to zero due to the high-dimensional biases, which would result in poor out-of-sample performance.
To select the number of dynamic eigenvalues, we therefore implement the following model selection strategy.
Starting from a model with only static eigenvalues, we increase the number of dynamic eigenvalues one-by-one, starting from the largest eigenvalue.
We then compute the change in BIC relative to the static model by accounting for the change in the in-sample log-likelihood and the change in the number of parameters.
We select the dynamic model with the lowest BIC.
The procedure typically results in a limited number of dynamic eigenvalues and a large number of static ones.
Combined with the high-dimensional shrinkage procedure for targeting the intercepts, this results in a considerable reduction of the parameter space.
We show later using simulated and empirical data that this approach efficiently captures the salient features of the data, including their dynamics.
\subsection{Parameter estimation}
\label{subsec:parameter estimation}
We split the estimation procedure in two standard steps.
We first estimate the marginal behavior of each original univariate series $r_{i,t}$ by an AR(1)-GARCH(1,1) model using standard quasi maximum likelihood (QMLE).
From the devolatilized residuals $y_{i,t}=\epsilon_{i,t}/\sigma_{i,t}$, we use the non-parametric rank transformation $u_{i,t} = \text{rank}(y_{i,t})/(T+1/2)$ to obtain the PITs.
This makes the construction of the PITs more robust to any potential mis-specification of the marginal distributions, in line with the use of QMLE to estimate the AR(1)-GARCH(1,1) marginal models.
Next, we estimate the dynamic dependence structure by maximizing the copula density \eqref{eq:GH copula}.
For a given tail shape parameter $\nu$ and skewness parameter $\bm{\gamma}$, we can compute the inverse marginal cdf projections $\bm{y^{\star}}_{t} = \bm{y^{\star}}_{t}(\nu,\bm{\gamma}) = G^{-1}(u_{i,t};\bm{\gamma},\nu)$.
Let $\bm{S}^{\star} = \mathbb{E}[\bm{y^{\star}}_{t} \bm{y^{\star}}_{t}{}^\top]$ and $\bm{R} = \mathbb{E}[\bm{R}_{t}]$, then for the skewed $t$ distribution we have that
\begin{align}
\label{eq:target omega}
\bm{S}^{\star} &=
\frac{\nu}{\nu-2} \bm{R} +
\frac{2\nu^2}{(\nu-2)^2(\nu-4)}
\bm{\gamma}\bm{\gamma}^\top
\quad \Leftrightarrow \quad
\bm{R} =
\frac{\nu-2}{\nu} \bm{S}^{\star} -
\frac{2\nu}{(\nu-2)(\nu-4)}
\bm{\gamma}\bm{\gamma}^\top
.
\end{align}
Replacing $\bm{S}^{\star}$ in \eqref{eq:target omega} by the sample estimator $\hat{\bm{S}}^{\star} = T^{-1}\sum_{t=1}^{T} \bm{y^{\star}}_{t}\bm{y^{\star}}_{t}{}^\top$, we then define the estimators
\begin{equation}
\label{eq:target omega2}
\begin{split}
\hat\bm{\Sigma} &=
\hat\bm{\Sigma}(\bm{\gamma},\nu) =
\hat\bm{W}(\bm{\gamma},\nu)\ \hat\bm{\Lambda}(\bm{\gamma},\nu)\ \hat\bm{W}(\bm{\gamma},\nu)^\top
=
\frac{\nu-2}{\nu} \hat{\bm{S}}^{\star} -
\frac{2\nu}{(\nu-2)(\nu-4)}
\bm{\gamma}\bm{\gamma}^\top
,
\\
\hat\bm{R} &= \operatorname{diag}(\hat\bm{\Sigma})^{-1/2}\ \hat\bm{\Sigma}\ \operatorname{diag}(\hat\bm{\Sigma})^{-1/2}.
\end{split}
\end{equation}
The intercepts of the score-driven transition Eq.~(\ref{eq:first score eq}) are targeted as explained in Section~\ref{subsec:nonlinear shrinkage} using the quadratic shrinkage procedure of \citet{Ledoit2022b} from Eq.~(\ref{eq:Shrink}), yielding the final (shrunken) estimators
\begin{align}
\label{eq:target omega3}
\check{\bm{\Sigma}} &=
\check{\bm{\Sigma}}(\bm{\gamma},\nu) =
\hat\bm{W}(\bm{\gamma},\nu)\ \check{\bm{\Lambda}}(\bm{\gamma},\nu)\ \hat\bm{W}(\bm{\gamma},\nu)^\top
,
&
\check{\bm{R}} &=
\operatorname{diag}(\check{\bm{\Sigma}})^{-1/2}\ \check{\bm{\Sigma}}\ \operatorname{diag}(\check{\bm{\Sigma}})^{-1/2}
,
\end{align}
where $\check{\bm{\Lambda}} = \check{\bm{\Lambda}}(\bm{\gamma},\nu)$ holds the quadratically shrunken spectrum from Eq.~(\ref{eq:Shrink}).
Note that all these quantities are still functions of $\bm{\gamma}$ and $\nu$.
Finally, we compute the copula log-likelihood function as
\begin{align}
\label{eq:loglik iteration k}
\mathcal{L}(\bm{\psi}) &=
\sum_{t=1}^T
\left(
\log g\left( \bm{y^{\star}}_{t} ; \check{\bm{R}}_{t}, \bm{\gamma}, \nu \right)
-
\sum_{i=1}^d \log g_{i}\left( y^{\star}_{i,t} ; \gamma_{i}, \nu \right)
\right)
,
\\
\nonumber
\check{\bm{R}}_{t} &= \operatorname{diag}(\check{\bm{\Sigma}}_{t})^{-1/2}\,
\check{\bm{\Sigma}}_{t}\, \operatorname{diag}(\check{\bm{\Sigma}}_{t})^{-1/2},
\qquad
\check{\bm{\Sigma}}_{t} = \hat\bm{W}\, \check{\bm{\Lambda}}_{t}\, \hat\bm{W}^\top,
\\
\nonumber
\log \check{\lambda}_{i,t+1}
&= \left\{ \begin{array}{ll}
(1-b_{i}) \log \check{\lambda}_{i} + b_{i}\, \log \check{\lambda}_{i,t} + a_{i} \nabla_{i,t},
& \text{for } i=1,\ldots,d_0,
\\
\log \check{\lambda}_{i}, & \text{for } i=d_0+1,\ldots,d,
\end{array} \right.
\end{align}
where $\bm{\psi}$ gathers all the static parameters $\bm{A}$, $\bm{B}$, $\bm{\gamma}$, and $\nu$ of the model, and where we use the BIC to select the number of time-varying eigenvalues $d_0$, as explained in Section~\ref{subsec:nonlinear shrinkage}.
The static parameter vector $\bm{\psi}$ is then estimated using Maximum Likelihood.
\section{Simulation study}
\label{sec:Sim}
To assess the performance of the high-dimensional dynamic copula model with spectral regularization in a controlled environment, we perform two simulation experiments.
In the first experiment, we focus on the performance of the spectral copula structure vis-\`a-vis a grouped factor copula structure \citep[as used in, e.g.,][]{Oh2017,Oh2018,Oh2023,Opschoor2020}. We investigate how the different copula structures and the shrinkage method behave under controlled deviations from a grouped factor structure.
In our second experiment, we focus on the parameter estimation performance of our new model.
\subsection{Experiment 1: spectral versus grouped copula structures}
\label{subsec:SimStatic}
In both experiments, we specify the unconditional correlation matrix structure $\bm{R}$ in our data generating process (dgp) as
\begin{equation}
\label{eq:sim dependence structure}
R_{i,j} =
\frac{
\beta_{M}^2 +
\beta_{G_{i}}^2 \delta_{G_{i},G_{j}} +
\beta_{I}^2 \delta_{i,j}+
\beta_{C}^2 e^{-|C_{i} - C_{j}|/2}
}{
\sqrt{
\beta_{M}^2 +
\beta_{G_{i}}^2 +
\beta_{C}^2 +
\beta_{I}^2
}\ \
\sqrt{
\beta_{M}^2 +
\beta_{G_{j}}^2 +
\beta_{C}^2 +
\beta_{I}^2
}\
}
.
\end{equation}
Here, $\delta_{i,j}$ is the Kronecker delta with $\delta_{i,j} = 1$ if $i=j$, and zero else.
We can interpret this correlation structure as corresponding to a factor model of the form:
$ \beta_{M}\,F^{M}_{t} + \beta_{G_{i}}\,F^{G_{i}}_{t} + \beta_{C_{i}}\,F^{C_{i}}_{t} + \beta_{I}\,\eta_{i,t}$,
where $\beta_{M}$ is the loading coefficient of a common market factor $F^{M}_{t}$ to which all assets are exposed;
$\beta_{G_{i}}$ is the factor loading for a group factor $F^{G_{i}}_{t}$ corresponding to group $G_{i}$ of asset $i$, e.g., an industry factor;
$\beta_{C}$ is the factor loading for the factor $F^{C_{i}}_{t}$, which defines, for instance, a country factor, i.e., a second grouping structure over and above the first (industry) grouping effect;
and finally $\beta_{I}$ represents the size of the idiosyncratic noise component.
Higher values of $\beta_{I}$ decrease the signal-to-noise ratio.
In our example, the (industry) indicators $G_{i}$ and (country) indicators $C_{i}$ both range from 1 to 10, and the factors have zero means and are uncorrelated with each other.
The only exception is the country factor $F^{C_{i}}_{t}$, which is correlated between countries $C_{i}$ with a correlation that varies with the distance $|C_i-C_j|$.
If $\beta_{C}=0$, the dgp collapses to a pure group structure as in \citet{Oh2017,Oh2018,Oh2023} and \citet{Opschoor2020}.
For $\beta_{C} > 0$, however, the group structure is only approximate, which allows us to study the effect of such deviations on both our new copula structure as well as alternatives from the literature.
We consider $d=100$ assets, where each asset corresponds to a unique pair $(G_{i},C_{i})$.
If $\beta_{C} > 0$, the correlation structure can no longer be easily restored by a simple extension of the number of factors, while retaining the block-diagonality of the factor loading matrix and the orthogonality of the group factors; see Supplementary \ref{app:oh methodology}.
\begin{table} [t]
\caption{\label{tab:SimExpDynamic} Correlation properties of the empirical data set and the stylized dependence model. We show the properties of Gaussian rank correlations, such as the cross-sectional mean, standard deviation and the range over all correlations. We show also the mean, standard deviation and the range over two specific sectors, namely the industry sector with highest average correlation and lowest average correlation. }
\begin{center}
\begin{tabular}{ l ccc c ccc}
\hline
& \multicolumn{3}{c}{Empirical} & & \multicolumn{3}{c}{Simulation}\\
\cmidrule{2-4} \cmidrule{6-8}
& mean & std. dev. & range & & mean & std. dev. & range \\
\hline
All correlations & 0.29 & 0.12 & $[0.00-0.78]$ & & 0.29 & 0.15 & $[0.05-0.73]$ \\
Lowest correlated sector & 0.30 & 0.12 & $[0.11-0.56]$ & & 0.31 & 0.13 & $[0.14-0.58]$ \\
Highest correlated sector & 0.54 & 0.08 & $[0.41-0.71]$ & & 0.59 & 0.07 & $[0.49-0.73]$ \\
\hline
\end{tabular}
\end{center}
\end{table}
We select the values of the parameters in the dgp to resemble the empirical complexity of the international stock market data as used in Section~\ref{sec:Emp}.
We set $\beta_{I}=1$, $\beta_{M} = 0.75$ and let $\beta_{G_{i}}=1.75-0.15i$ to have groups with different intragroup (industry) copula correlations.
For $\beta_{C}$, we consider a parameter range $\beta_{C} \in \{0,0.75,1.5\}$.
This means that country-effects are either fully absent, or are smaller than industry-effects, or have a comparable size.
For $\beta_{C} =0$, the grouped factor structure defined by the groups $G_{i}$ is exact and all within-group correlations have the same value.
For $\beta_{C}>0$, the within-group correlations start to differ and the group structure based only on $G_{i}$ no longer captures the full heterogeneity in the dependence structure.
Having considerable heterogeneity is in line with our empirical data study in Section~\ref{sec:Emp}.
Table \ref{tab:SimExpDynamic} shows various properties of unconditional Gaussian rank correlations for the empirical data from Section~\ref{sec:Emp}, and for simulated Gaussian copula data based on the stylized correlation matrix with $\beta_{C} = 1.5$.
The two sets of correlations and their heterogeneity show similar properties between the simulated dgp and the empirical data.
The stylized parameters result in a cross-sectional average correlation of 0.29, which includes both intrasector correlations and cross-sector correlations.
The correlations averaged over each sector in the simulation vary from 0.31 to 0.59.
The intrasector correlations themselves vary from 0.14 to 0.58 within the lowest correlated sector and from 0.49 to 0.73 within the highest correlated sector.
Similar values and ranges are found for the data set used in Section~\ref{sec:Emp}.
The simulation setting thus mirrors the stylized facts in the empirical data quite well.
As our benchmark model in the simulations and in the empirical application later on, we use the dynamic factor copula approach of \citet{Opschoor2020} and \citet{Oh2023}.
We select the optimal number of clusters endogenously using the clustering methodology of \cite{Oh2023} based on their computer code.
See Supplementary \ref{app:additional simulations} for more details on the dependence structure imposed by their methodology.
We consider a $d = 100$ dimensional time series and set the number of simulated observations to $T \in \{250, 1000\}$.
We compute the models' performance metrics both in-sample and out-of-sample to illustrate the effect of over-fitting and of shrinkage.
On top of the $T$ in-sample observations, we therefore generate an additional 1000 observations from the same dgp and re-calculate the log-likelihood without re-estimating the model's parameters to obtain the out-of-sample log-likelihood.
To simulate the data, we use a skewed $t$ with $\nu=25$ and $\bm{\gamma} = \gamma\,\bm{\iota}$, where $\gamma = -0.25$ and $\bm{\iota}$ is a vector of ones, as in \citet{Oh2023}.
These parameter values are similar to those obtained from empirical data.
\begin{table}[t]
\begin{center}
\caption{
\label{tab:SimExpStatic}
Comparison of in-sample and out-of-sample log-likelihoods ($\ell_{in}$ and $\ell_{out}$) for different correlation structures, estimation methods and sample sizes $T$ in $d=100$ dimensions.
The true copula dependence structure $\bm{R}$ is given in \eqref{eq:sim dependence structure} and either has a single clear grouped correlation structure (such as industry-only, $\beta_{C}=0$), or a stylized, but realistic additional intra and inter-group correlation structure (e.g. country, $\beta_{C}=0.75$ or 1.5).
The $T$ in-sample observations are used to estimate the models.
An additional $T$ simulated out-of-sample observations are used to compute $\ell_{out}$.
In the dgp, data are generated using the skew $t$ copula with $\nu=25$ and $\bm{\gamma} = -0.25\,\bm{\iota}$.
Best performing models per combination of dependence structure and sample size (row-wise) are bolded.
}
\begin{tabular}{cc cc c cc c cc c cc}
\hline
$\beta_{C}$ & $T$ & \multicolumn{2}{c}{True } & & \multicolumn{2}{c}{Regularized} & & \multicolumn{2}{c}{Sample} & & \multicolumn{2}{c}{Factor} \\
\cmidrule{3-4} \cmidrule{6-7} \cmidrule{9-10} \cmidrule{12-13}
& & $\ell_{in}$ & $\ell_{out}$ & & $\ell_{in}$ & $\ell_{out}$ & & $\ell_{in}$ & $\ell_{out}$ & & $\ell_{in}$ & $\ell_{out}$\\
\midrule
0 & $1000$ & 35,168 & 36,206 & & 36,899 & 34,798 & & 37,563 & 33,366 & & 34,811 & $\bm{35,791}$ \\
0 & $250$ & 8,562 & 8,577 & & 9,925 & 7,694 & & 11,496 & 4,380 & & 8,555 & $\bm{8,507}$ \\
0.75 & $1000$ & 39,392 & 39,665 & & 41,062 & $\bm{38,135}$ & & 41,800 & 36,805 & & 30,922 & 31,304 \\
0.75 & 250 & 9,422 & 9,414 & & 10,903 & $\bm{8,283}$ & & 12,399 & 5,168 & & 7,402 & 7,313 \\
1.5 & $1000$ & 55,693 & 55,748 & & 57,185 & $\bm{54,459}$ & & 58,132 & 52,906 & & 33,943 & 33,439 \\
1.5 & 250 & 13,611 & 13,529 & & 15,051 & $\bm{12,464}$ & & 16,601 & 9,295 & & 8,343 & 8,120 \\
\hline
\end{tabular}
\end{center}
\end{table}
Table~\ref{tab:SimExpStatic} presents the results.
There are three main takeaways.
First, if there is a true factor structure ($\beta_{C}=0$), then the factor cluster copula model of \citet{Oh2023}, labeled Factor, performs best based on the out-of-sample log-likelihood $\ell_{out}$, with the regularized spectral model in second place.
This is encouraging, as the out-of-sample performance loss for the spectral method vis-\`a-vis a correctly specified more parsimonious model is limited.
To appreciate this, we note that the spectral dynamic copula model is agnostic about the true model structure and therefore still requires the estimation of 4950 different correlation parameters.
Second, we see that MLE combined with the sample correlation matrix (in the column Sample) suffers from the curse of dimensionality for large correlation matrices, as expected \citep{Ledoit2022a}.
The out-of-sample performance of the MLE is significantly worse than its in-sample performance.
The biases in the spectrum become more pronounced if the number of in-sample data points decreases to, e.g., $T=250$.
The problem of overfitting is resolved either when a parsimonious model structure is imposed through a cluster factor structure (in the column Factor) or when we use the quadratic shrinkage approach of \citet{Ledoit2022b} (in the column Regularized).
Third, when we depart from the true group structure, the results for $\beta_{C}=0.75$ and 1.5 show that the spectral copula approach with regularization performs best by a wide margin.
Again, when shrinkage is not applied, we see similar biases as before.
Also, the effect of shrinkage becomes larger as the ratio $d/T$ increases.
We also see that the differences in out-of-sample log-likelihood performance between the clustering and the regularized approach is large for $\beta_{C} = 0.75$ and 1.5.
This holds even if one allows for an endogenous choice of the number of groups and asset group assignments as in \citet{Oh2023}.
In the stylized correlation matrix from Eq.~\eqref{eq:R parameterization}, there are variations in correlations both within groups and across groups, while the cluster factor approach imposes that all correlations for stocks belonging to a group are the same.
Moreover, the cluster factor structure does not span the full space of correlation matrices (see Supplementary \ref{app:additional simulations}).
As a result, the cluster factor structure can miss out on important forms of heterogeneity in the data, both in the simulation and the empirical study.
\subsection{Experiment 2: quality of point and path estimates}
\label{subsec:SimDynamic}
\begin{table} [t]\centering
\caption{\label{tabSimExpDynamic} Simulation results for the copula parameter estimates based on the dynamic skew $t$ copula as dgp with the dependence structure from \eqref{eq:sim dependence structure}. The sample size is $T=1000$ in $d=100$ dimensions.
Sample ML uses the sample spectrum of $\hat\bm{R}$ in Eq.~\eqref{eq:target omega2}, while Regularized ML uses the quadratic shrinkage formula of \citet{Ledoit2022b}. The simulations are repeated 100 times. The mean and the standard deviation of the estimates are reported.
}
\begin{tabular}{ l r c r r c r r }
\hline
& & & \multicolumn{2}{c}{Regularized ML} & & \multicolumn{2}{c}{Sample ML} \\
\cmidrule{4-5} \cmidrule{7-8}
Parameter & True & & Mean & Std. dev. & & Mean & Std. dev. \\
\midrule
$\lambda_1$ & 29.8 & & 29.4 & (1.8) & & 29.7 & (1.7) \\
$\lambda_2$ & 10.7 & & 10.7 & (0.8) & & 10.8 & (0.7) \\
$\lambda_{99}$ & 0.160 & & 0.171 & (0.008) & & 0.101 & (0.003) \\
$\lambda_{100}$ & 0.160 & & 0.164 & (0.009) & & 0.097 & (0.003) \\
$a_1$ & 0.10 & & 0.09 & (0.01) & & 0.09 & (0.01) \\
$b_1$ & 0.90 & & 0.89 & (0.03) & & 0.89 & (0.03) \\
$a_2$ & 0.10 & & 0.09 & (0.02) & & 0.08 & (0.02) \\
$b_2$ & 0.90 & & 0.88 & (0.06) & & 0.88 & (0.07) \\
$\nu$ & 25.0 & & 24.6 & (1.8) & & 29.2 & (2.2) \\
$\gamma$ & -0.25 & & -0.28 & (0.07) & & -0.31 & (0.08) \\
\hline
\end{tabular}
\end{table}
So far, the simulation study concentrated on the unconditional dependence structure of the copula and the performance of the shrinkage procedure.
Next, we consider the quality of the point estimates of the copula parameters, the selection of the number of dynamic components, and the path estimates of the dynamic eigenvalues.
We take the score-driven dynamic spectral copula model from Section~\ref{sec:Def} as our dgp.
The stylized correlation matrix from \eqref{eq:sim dependence structure} is again used with $\beta_{C}=1.5$.
The other parameter values (see Table~\ref{tab:SimExpDynamic}) are chosen to be similar to the empirical results from Section \ref{sec:Emp}.
The first 2 eigenvalues in the dgp are dynamic, while the remaining 98 eigenvalues in the dgp are static, also in line with the empirical results.
We estimate the static parameters using the Sample ML as well as the Regularized ML procedure as described in Section~\ref{subsec:nonlinear shrinkage}.
The estimation results can be found in Table~\ref{tabSimExpDynamic}.
The results in Table~\ref{tabSimExpDynamic} confirm that the Sample ML approach suffers from biases in the estimation of the eigenvalues at the lower end of the spectrum, i.e., for small eigenvalues.
At the upper end of the spectrum, i.e., for the large eigenvalues, the effect of shrinkage is less pronounced.
This means that the Regularized ML estimator successfully removes the biases at the low end of the spectrum, while leaving `spiked' eigenvalues intact.
We note that also all other parameters, such as $a_{i}$, $b_{i}$, $\nu$, and $\bm{\gamma}=\gamma\,\bm{\iota}$, are estimated accurately.
\begin{figure}[t]
\centering
\includegraphics[width=1.0\columnwidth]{fig_DynSim.eps}
\caption{Performance of the skewed $t$ copula with regularized spectral dynamics. The left plot shows changes in the (in-sample) BIC of a model with static eigenvalues versus a model with eigenvalues $1,\ldots,i$ being dynamic, as a function of the spectral index $i$.
The BIC is minimal for $i=2$, which is also the true number of dynamic eigenvalues in the dgp. The average result of 10 Monte Carlo simulations is shown. The right-hand plot shows the dynamics of the first eigenvalue for a single MC simulation, where the true dynamics of the dgp is compared with the predicted dynamics from the estimated model. The period after 4 years is an out-of-sample forecast. }\label{sim:figEigValDyn}
\end{figure}
Figure~\ref{sim:figEigValDyn} presents the results for the model selection procedure to determine the number of dynamic eigenvalues.
In each Monte Carlo simulation and for $i=1,\ldots,7$, we successively estimate models with eigenvalues $1,\ldots,i$ as dynamic, and eigenvalues $i+1,\ldots,d$ as static.
We compute the in-sample BIC decrease with respect to a fully static model ($i=0$), as well as its out-of-sample log-likelihood increase.
The left-hand panel in Figure~\ref{sim:figEigValDyn} shows that the average $\Delta$BIC curve has its minimum at the correct number of $d_0=2$ dynamic components.
This analysis validates our choice to use BIC for selecting the number of dynamic eigenvalues in the empirical study of Section~\ref{sec:Emp}.
From the out-of-sample log-likelihood increases (right-hand axis in Figure~\ref{sim:figEigValDyn}), we see that it is particularly important to include the dynamics of the first eigenvalue, since it has the largest impact.
For spectral indices above $i=2$, there are no further increases in log-likelihood if we make these eigenvalues dynamic, and the BIC grows linearly in the number of parameters.
The stable out-of-sample log-likelihood suggests that the risk in overfitting the number of dynamic components is relatively low, since overstating the number of dynamic eigenvalues does not have a negative effect on the out-of-sample log-likelihood.
The right-hand panel in Figure~\ref{sim:figEigValDyn} shows a simulated path of the first dynamic eigenvalue with the dgp and the estimated path from the dynamic copula model.
This shows that the eigenvalue dynamics are captured accurately.
An additional simulation experiment in Supplementary~\ref{app:sim misspec} confirms that the dynamic score-driven model can still recover the true, unobserved eigenvalue dynamics even if it is mis-specified, in line with theoretical consistency results in \citet{beutnerlinlucas2023}.
\section{Empirical study}
\label{sec:Emp}
\subsection{Data}
\label{subsec:data}
In our empirical study, we consider daily log return data for 100 stocks from 10 different industry sectors and 10 different European countries over a period of 10 years.
The data are obtained from LSEG (formerly known as Refinitiv), so that our industry sectors are based on The Refinitiv Business Classification (TRBC), which distinguishes between Financials (FI), Industrials (IN), Technology (TE), Basic Materials (BM), Consumer Cyclicals (CC), Consumer Non-Cyclicals (CN), Utilities (UT), Health Care (HC), Energy (EN) and Real Estate (RE).
We select stocks that are issued and traded in 10 different European countries: Germany (DE), United Kingdom (UK), France (FR), Spain (ES), Italy (IT), Sweden (SE), Norway (NO), The Netherlands (NL), Belgium (BE) and Switzerland (CH). For each sector and country combination, we select stocks with the largest market capitalizations after the end of the observation period, while also requiring that each stock is fully observed over the 10 year period from January 2015 to 31 December 2024.
We have a few (8) missing country-industry combinations that do not satisfy our criteria, e.g., because they are not observed over the full sample period.
In such cases, we take the stock from another country in that same industry (with second highest market capitalization).
We are able to do so in such a way that we end up with each of the 10 countries and the 10 sectors occurring precisely 10 times, where some country-sector combinations will be missing, while other combinations occur twice.
Table \ref{tabIndices100} of Supplementary~\ref{app:additional empirics} lists the full set of all 100 stocks across countries and sectors.
Thus far, high-dimensional copula studies have mainly focused on stocks from a single country, such as the US \citep[see, e.g.,][]{Oh2023}.
The large number of possible country and sector combinations increases the complexity of the dependence structure in the current data set and poses challenges to standard factor clustering algorithms.
As a result, copulas that do not impose strong restrictions on the correlation structure can achieve substantial performance gains.
\begin{table} [t]\centering
\caption{\label{tabEmpiricalMarginal} Summary statistics of univariate times series. Panel A presents the unconditional mean, standard deviation, skewness and kurtosis of the daily log-returns. Panel B shows the estimated parameters of the marginal AR(1)-GARCH(1,1) models for the univariate time series,
$
r_{i,t}= \delta_{i} + \phi_{i} r_{i,t-1}+ \epsilon_{i,t}$ and $\sigma^2_{i,t} = \omega_{i} + \alpha_{i} \epsilon^2_{i,t-1}+ \beta_{i} \sigma^2_{i,t-1}$, estimated by QMLE.
Panel C shows the correlations of the devolatilized returns $y_{i,t} = \epsilon_{i,t}/\sigma_{i,t}$. In all cases the mean and the range over the cross-sectional dimension of 100 stocks is given.}
\begin{tabular}{ l p{0.1cm} cc c l cc}
\hline
Cross section & & mean & range & & & mean & range \\
\hline
\multicolumn{8}{l}{Panel A: log-returns $r_{i,t}$}\\
\hline
Mean & & 0.000 & $[-0.001,\,0.001]$ & & Skewness & $-0.406$ & $[-1.82,\,1.25]$ \\
Std. dev. & & 0.018 & $[0.011,\,0.030]$ & & Kurtosis & 12.95 & $[5.72,\,35.25]$ \\
\hline
\multicolumn{8}{l}{Panel B: univariate AR(1)-GARCH(1,1) parameter estimation results}\\
\hline
$\sqrt{\omega_{i}}$ & & 0.004 & $[0.001,\,0.011]$ & & $\delta_{i}$ & 0.000 & $[-0.001,\,0.001]$ \\
$\alpha_{i}$ & & 0.090 & $[0.013,\,0.317]$ & & $\phi_{i}$ &-0.023 & $[-0.195,\, 0.061]$ \\
$\beta_{i}$ & & 0.846 & $[0.548,\,0.979]$ & & & & \\
\hline
\multicolumn{8}{l}{Panel C: Average cross-sectional correlations devolatilized residuals $y_{i,t} = \epsilon_{i,t}/\sigma_{i,t}$}\\
\hline
$\rho_{i,j}$ & & 0.275 & $[0.004,\,0.786]$ & & & & \\
\hline
\end{tabular}
\end{table}
We merge the time series of the log-returns from different countries based on common trading days.
This results in a data set with $T = 2,425$ daily returns $r_{i,t}$ for each asset $i=1,\ldots,d$.
Summary statistics of the data are presented in Panel A of Table \ref{tabEmpiricalMarginal}, confirming standard stylized facts such as unconditional left-skewness and substantial excess kurtosis of the raw log-returns.
We filter each time series $r_{i,t}$ using a standard AR(1)-GARCH(1,1) filter
and proceed the analysis with the devolatilized residuals $y_{i,t}=\epsilon_{i,t}/\sigma_{i,t}$.
Panel B of Table \ref{tabEmpiricalMarginal} summarizes the univariate estimation results.
The first two columns of panels in Fig. \ref{figAcf} show the autocorrelation functions for absolute log-returns $|r_{i,t}|$ and for absolute devolatilized returns $|y_{i,t}|$ for 3 (arbitrarily chosen) stocks.
As expected, we observe strong serial correlation for the $|r_{i,t}|$, which indicates clear volatility clustering effects.
After applying the univariate AR(1)-GARCH(1,1) filters, the standardized absolute residuals $|y_{i,t}|$ no longer indicate any substantial volatility clustering.
We take these standardized residuals $y_{i,t}$ as input for the Probability Integral Transformations (PITs) using the rank-based transforms $u_{i,t} = \text{rank}(y_{i,t})/(T+1/2)$ as motivated in Section~\ref{subsec:parameter estimation}.
Our final sample consists of $T = 2,425$ pseudo-copula observations in $d=100$ dimensions.
We use the first half of the data for estimation, and the second half for our out-of-sample performance evaluation.
\begin{figure}[!t]
\centering
\includegraphics[width=1.0\columnwidth]{fig_acf.eps}
\caption{Autocorrelation functions for log-returns $|r_{i,t}|$, devolatilized residuals $|y_{i,t}|$, $|y^{\star}_{i,t}|= |G^{-1}(u_{i,t})|$ and spectral projections $|\tilde y^{\star}_{i,t}|$ for $i=1,2,10$. The Gaussian copula ($\gamma=\bm{0}_{d}$ and $\nu^{-1}=0$) is used to determine $y^{\star}_{i,t}$ and $\tilde y^{\star}_{i,t}$.}\label{figAcf}
\end{figure}
Interestingly, we can use the pseudo-copula observations $y^{\star}_{i,t} = G_{i}(u_{i,t})$ to explore whether there will be time-variation in the eigenvalues of the copula dependence matrix.
The third column in Figure~\ref{figAcf} provides the autocorrelation function of the absolute pseudo-copula observations $|y^{\star}_{i,t}|$.
Again, we see no substantial signal of volatility clustering in the values of $y^{\star}_{i,t}$.
However, if we compute the eigenvalue-eigenvector decomposition of the sample correlation matrix of $\bm{y}_{t}$ and use it to compute initial estimates of the spectral projections $\tilde y^{\star}_{i,t}$, the last column of Figure~\ref{figAcf} clearly shows volatility clustering effects for the first spectral projection, which is indicative of time-variation in $\lambda_{1,t}$.
As another example, also the second spectral projection is shown, for which the time-varying volatility, i.e., a time-varying $\lambda_{2,t}$, is less strong.
The further down in the spectrum we consider the spectral projections, the less evidence we find for a time-varying $\lambda_{i,t}$.
For instance, for $i=10$ the lower-right autocorrelation function largely remains within the confidence band.
This indicates that a modeling approach with a limited number of time-varying $\lambda_{i,t}$s for the copula dependence parameters is both parsimonious and congruent with the financial data at hand.
Before presenting the estimation results for the dynamic spectral regularized copula specification, we note that our analysis differs conceptually from existing (spectral) GARCH models in the literature, such as the orthogonal GARCH or the $\lambda$ GARCH models \citep[e.g.,][]{hetland2023dynamic}.
The latter do not explicitly distinguish between the marginal and the copula time series.
As a result, the conditional volatility effects in spectral directions in those models could stem from marginal volatility increases as well as from increased correlation effects.
In our framework, by contrast, a larger conditional variance $\lambda_{1,t}$ solely reflects the dynamics of the copula, allowing us to answer whether strongly dependent market movements (as measured by the first spectral projections) are followed by periods of higher dependence.
\subsection{Estimation results}
Given that the time-variation in $\lambda_{i,t}$ is concentrated in the first few spectral projections according to the autocorrelograms in Figure~\ref{figAcf}, we use a BIC guided selection procedure to determine the final number of dynamic components.
We start from the static skew $t$ copula, after which we make the first eigenvalue dynamic.
We report the change in (unregularized) BIC relative to the static model.
After that, we also make the second value dynamic and report the change in BIC relative to the static model, and so on.
We do not apply shrinkage at this stage.
Figure~\ref{figEigValDyn} shows the changes in the BIC vis-\`a-vis the static model when adding the dynamic eigenvalues for the skewed $t$ copula one by one for $i=1,...,7$.
Clearly, the largest gain is obtained by allowing the first $\lambda_{i,t}$ to be dynamic, followed by a smaller improvement for $i=2$.
The BIC is at its minimum for $i=2$ dynamic components, which therefore is the value we use in the remainder of the analysis.
We use the same number of 2 dynamic eigenvalues for the three different copula densities, endowing them with the score-driven dynamics from Section~\ref{subsec:dynamics and shrinkage}.
\begin{table}[t]
\caption{\label{tabLoglCase}
Performance of various static and dynamic copula specifications with either a factor structure or a regularized spectral structure for the correlation matrix. We consider 100 European stocks from 10 different countries and industry sectors ($d=100$) with 10 years of daily data. The first 5 years are used for estimation and the last 5 years for out-of-sample performance.
Panel A shows the in-and out-of-sample log-likelihood for dynamic copula specifications.
Panel B shows the same results based on static copula specifications.
The best performing copula specification is indicated in bold, which is the dynamic skew $t$ copula with non-linear shrinkage.
}
\begin{tabular}{l c c c c c c c c c c c}
\hline
& \multicolumn{3}{c}{Regularized} & & \multicolumn{3}{c}{Sample} & & \multicolumn{3}{c}{Factor} \\
\cmidrule{2-4} \cmidrule{6-8} \cmidrule{10-12}
& \it{G} & $t_d$ & skew $t_d$ & & \it{G} & $t_d$ & skew $t_d$ & & \it{G} & $t_d$ & skew $t_d$ \\
\hline
\multicolumn{10}{l}{Panel A: dynamic copula likelihoods}\\
\hline
$\ell_{in}$ & 38,068 & 38,537 & 38,574 & & 38,220 & 39,012 & 39,045 & & 28,376 & 29,641 & 29,686 \\
$\ell_{out}$ & 27,882 & 29,725 & $\bm{29,752}$ & & 27,139 & 28,705 & 28,733 & & 24,882 & 26,531 & 26,555\\
\hline
\multicolumn{10}{l}{Panel B: static copula likelihoods}\\
\hline
$\ell_{in}$ & 37,397 & 38,058 & 38,104 & & 37,664 & 38,650 & 38,691 & & 28,088 & 29,481 & 29,529 \\
$\ell_{out}$ & 26,769 & 29,529 & 29,549 & & 25,513 & 28,230 & 28,253 & & 23,605 & 25,696 & 25,727 \\
\hline
\end{tabular}
\end{table}
Table \ref{tabLoglCase} shows the performance for the three different approaches (Regularized spectral copula, Sample-based spectral copula, and the Factor copula with clustering) and three different choices for the copula density (Gaussian, Student's $t$, and the skewed $t$).
We estimate the factor model with optimal cluster assignment in 100 dimensions using the algorithm of \cite{Oh2023}, where 17 clusters turn out to be optimal in terms of BIC.
Panels A and B of Table \ref{tabLoglCase} show the in-sample and out-of-sample log-likelihood values for the different dynamic and static models, respectively.
First, we find that the spectral dynamic copula models (panel A) perform significantly better than the static copula models (panel B), both in-sample and out-of-sample.
Modeling the dynamics of the dependence structure thus significantly improves the fit and predictive power of the models, even in the current parsimonious setting where only 2 out of the 100 eigenvalues are dynamic.
Second, the dynamic factor copula specifications with endogenous clustering perform significantly worse compared to the regularized dynamic spectral copula specifications.
This is mainly due to their more restrictive dependence structure compared to the spectral copula specification.
Third, as expected and as mentioned earlier, the $t$ copula performs substantially better than its Gaussian counterpart.
Performance differences between the symmetric and the skew $t$ copulas, by contrast, are modest, though there is a slight increase in the log-likelihood, both in-sample and out-of-sample.
Finally, non-linear shrinkage of the unconditional eigenvalues leads to significant improvements of $+1000$ points in the out-of-sample log-likelihood compared to the sample eigenvalues.
Spectral dynamics and regularization thus emerge as the most important aspects in modeling the copula dependence spectrum for large $d$ and $T$.
In Table \ref{tabEstCase}, we show the parameter estimates for the two dynamic spectral indices of the best performing dynamic skew $t$ copula with regularization.
To obtain confidence intervals for the estimates, the block bootstrap method was used on pseudo-copula observations $u_{i,t}$ using 20 trading days as block length and 200 bootstrap samples.
For both eigenvalues, the dynamics are highly persistent with $b_1$ and $b_2$ equal to 0.9 or higher.
The value of $a_1$ can be accurately estimated and appears clearly significant, signaling a time-varying dependence structure.
The value of $a_2$ can be less accurately estimated and the corresponding improvement in log-likelihood is smaller, but it is still yields a lower BIC and also the out-of-sample log-likelihood improves.
Such observations are also in line with Figure~\ref{figAcf}.
The estimated value for the skewness parameter $\gamma$ is significantly negative, while the estimated tail parameter $\nu$ is fairly large with values around 45. Still, the improvement in log-likelihood is substantial when allowing for tail-dependence, both in sample and out-of-sample, as can be seen from Table \ref{tabEstCase} when comparing the (skewed) $t$ copula with its Gaussian counterpart.
\begin{table}[t]\centering
\caption{\label{tabEstCase}
Estimation results for the skew $t$ copula with spectral dynamics based on 100 European stocks and 5 years of historic data.
The 90\% confidence intervals are determined with the block bootstrap method on the pseudo-copula observations $u_{i,t}$ using 20 trading days as block length.
}
\begin{tabular}{lcc ccc ccc}
\hline
Par. & Est. & 90\% CI & Par. & Est. & 90\% CI & Par. & Est. & 90\% CI \\
\hline
$\check{\lambda}_1$\rule{0ex}{2.5ex} & 30.6 & [26.2,\,\,35.5] & $a_{1}$ & 0.06 & [0.05,\,\,0.07] & $b_1$ & 0.90 & [0.86,\,\,0.94] \\
$\check{\lambda}_2$ & 5.8 & [ 5.3,\,\, 6.6] & $a_2$ & 0.06 & [0.03,\,\,0.21] & $b_2$ & 0.97 & [0.26,\,\,0.99] \\
$\nu$ & 44.1 & [41.7,\,\,54.8] & $\gamma$ & -0.37 & [-0.47,\,\,-0.21] \\
\hline
\end{tabular}
\end{table}
\subsection{Eigenvalues and eigenvectors}
\begin{figure}[t]
\centering
\includegraphics[width=1.0\columnwidth]{fig_DynEmp.eps}
\caption{Spectral dynamics of the first spectral eigenvalue, where the ratio of the first eigenvalue and the sum of the higher eigenvalues is plotted over time. Under normal market circumstances, the first eigenvalue explains about 30\% of the spectral variance, leading to a ratio below 1/2. When the ratio increases, it signals enhanced spectral variance in the parallel direction (first eigenvector) compared to the total spectral variance in the orthogonal directions, indicating less diversification and higher systemic risk or return. The three most pronounced events are given an economic interpretation. The period from January 2020 to December 2024 represents an out-of-sample forecast. }\label{figEigValDyn}
\end{figure}
To conclude the empirical analysis, we zoom in on several aspects of the estimation results, such as the eigenvalue dynamics and its relation to financial crises, as well as the interpretation of the first eigenvectors.
We use the best performing copula model, namely the skew $t_d$ copula with regularized spectral dynamics.
In the right plot of Fig.~\ref{figEigValDyn}, we show the dynamics of the first eigenvalue $\lambda_{1,t}$, where we plot the ratio of the first eigenvalue to the sum of the smaller eigenvalues ($i\leq 2$) over time.
Under normal market circumstances, the first eigenvalue explains about 30\% of the spectral variance, leading to a ratio around 30/70, which is below 0.5.
There is strong serial dependence in the first eigenvalue, in line with the earlier result in the right-hand panels in Figure~\ref{figAcf} and Panel A in Table~\ref{tabEstCase}.
As is seen in Figure~\ref{figEigVec}, the first eigenvector takes a more or less equally weighted position in each of the assets, thus reflecting a parallel movement in the market.
Increases in the first eigenvalue compared to the remaining ones therefore result in correlation matrices with less diversification potential, leading to periods with higher systemic risk.
The right-hand panel in Figure~\ref{figEigValDyn} reveals that such parallel-to-orthogonal variance ratio increases can reach up to a factor two or higher during crisis periods, implying the first eigenvalue explains up to 70\% of the spectral variance during such times.
During these crisis events, the cross-sectional average correlation $\bar{\bm{R}}_t$ peaks at 0.65, while the unconditional average is only 0.30.
The three most pronounced events in the time period of Figure~\ref{figEigValDyn} are the flash crash in 2015, the Brexit in 2016 and the Covid pandemic in 2020, which corroborate the lack of diversification potential during such periods.
\begin{figure}[tb]
\centering
\includegraphics[width=1.0\columnwidth]{fig_ev.eps}
\includegraphics[width=1.0\columnwidth]{ev2.png}
\caption{Heatmap of estimated eigenvectors weights $\hat{w}_{i,j}$ of the unconditional copula correlation matrix $\hat\bm{R}$ from Eq.~\eqref{eq:target omega2}. The first two eigenvectors are shown ($j=1,2$). The stocks $i$ are ordered in terms of countries along the $x$-axis and in terms of sectors along the $y$-axis. The white color indicates that the country-sector combination is absent in the data set. In case of two stocks per country-sector combination, the result for the first stock from Table \ref{tabIndices100} is shown. The first eigenvector has positive weights for all stocks, corresponding to a parallel market movement. We compare the second eigenvector with optimal cluster assigments using \citet{Oh2023} giving 17 clusters, indicated by the numbers in each sector-country cell. The colors relate to the value of the corresponding element of the eigenvectors in the spectral copula specification.
The bottom panel contains the same information as the right-hand panel, but the stocks are ordered per cluster (in order of their average value of the second eigenvector elements).
}\label{figEigVec}
\end{figure}
In Figure~\ref{figEigVec}, we analyze the eigenvectors corresponding to the largest two dynamic eigenvalues in more detail and compare them to the cluster assignments based on \citet{Oh2023}.
This is an extention to \citet{Gubbels2025}, who perform the eigenvector analysis solely for country-effects in a static spectral copula context.
The left panel shows the first eigenvector as a heatmap.
All elements are positive and of roughly equal magnitude, so that the first eigenvector represents the collective movement of the European stock market as a whole.
This is in line with the market factor interpretation of the first factor by \citet{Oh2023}: the estimated intercepts of the factor loadings for the optimal number (using the BIC) of 17 clusters are all positive, see Table \ref{tabEstClus100D} of Supplementary~\ref{app:additional empirics}.
Moreover, we also see from Figure~\ref{figEigVec} that the stocks with lowest weight in the first eigenvector (green color in the left panel) all correspond to cluster 11, which has the lowest market co-movement in Table~\ref{tabEstClus100D}.
More generally, the inner product of the first eigenvectors of the unconditional correlation matrix from both approaches is 0.994.
To further compare the spectral eigenvectors with the optimal cluster assignments based on the methodology of \cite{Oh2023}, we label each of the assets by their cluster number 1,\ldots,17, where 17 denotes the optimal number of clusters in terms of the BIC.
The right-hand panel in Figure~\ref{figEigVec} shows the heatmap of the second eigenvector, which is orthogonal to the first eigenvector.
It represents the most important cross-sectional direction for diversification in the European stock market.
The eigenvector mainly captures the diversification between different sectors, where the financial, energy and basic material sectors have opposite sign to the real estate, health, utility and consumer non-cyclical sectors.
From Table~\ref{tabEstClus100D} and Figure~\ref{fig:omega OhPatton vs eigvec 2}, we see that clusters 3, 7, 12, 14 have the largest absolute cluster loadings $\omega^{C}_{i}$) and represent the sectors real estate, financials, utilities and energy.
For these clusters, the corresponding eigenvector weights are most negative for sectors 7 and 14, and most positive for sectors 3 and 12.
This means that there is alignment between the two models in identifying co-moving stocks.
Although weaker than sector effects, we also observe country effects in the second eigenvector.
The colors that correspond to Norway are more aligned with Sweden than with the UK, for example.
This shows that in our heterogeneous data set both industry sectors and countries have intertwined effects.
To further illustrate the alignment between cluster assignments and the second eigenvector, the bottom panel in Figure~\ref{figEigVec} shows the heatmap of the weights grouped per cluster and sorted by their (cluster) average weight.
The cluster factor selects clusters with similar second eigenvalue weights in the spectral decomposition.
The main difference between the two approaches, however, lies in the handling of the between-cluster correlations: in the cluster factor approach with its block-diagonal cluster loading matrix (see Supplementary~\ref{app:oh methodology}) the off-diagonal factor loadings are restricted to zero, whereas no such restriction applies in the spectral approach.
This is the primary cause why the spectral approach obtains a better fit to the data.
\begin{figure}[t]
\centering
\includegraphics[width=1.0\columnwidth]{fig_shrinkage.eps}
\caption{Estimated eigenvalues excluding and including shrinkage, $\hat{\lambda}_{i}$ and $\hat{\mu}_{i}$, of the unconditional copula correlation matrix for the dynamic skew $t$ copula. The sample eigenvalues are ordered from high to low. The left plot shows the first 10 eigenvalues and the middle plot the other 90 eigenvalues. The sample eigenvalues are shown in blue, while the eigenvalues after applying shrinkage are shown in red. The right plot shows the out-of-sample log-likelihood improvement when the first $i$ sample eigenvalues are replaced with shrinkage eigenvalues. The largest improvements stem from the highest spectral indices, whose sample eigenvalues are most severely biased. }\label{figShrinkage}
\end{figure}
As final analysis, we zoom in on the effect of shrinkage on the empirical results.
When $d = 100$, the concentration ratio $d/T$ is 8\%, since we use $T=1,213$ in-sample observations for estimation.
Though modest, the biasing effect on the spectrum is clearly present, as we have seen in the earlier results.
To understand how the quadratic shrinkage procedure contributes to the fit of the model to the empirical data, Figure~\ref{figShrinkage} shows the largest 10 and the remaining 90 shrunken and unshrunken targeted intercepts $\hat\omega_{i}$ from the score-driven transition dynamics in \eqref{eq:first score eq}.
For the largest eigenvalues we observe that shrinkage leads to a (slight) downward adjustment, since the sample eigenvalues have an upward bias.
For the lowest eigenvalues, the converse happens.
The \textit{relative} impact of shrinkage is much larger at the lower end of the spectrum.
It is precisely this correction at the lower end of the spectrum that causes the substantial increases in out-of-sample log-likelihood.
To see this, Figure~\ref{figShrinkage} plots the improvement in \textit{out-of-sample} log-likelihood for the dynamic skew $t$ copula due to non-linear shrinkage as a function of the spectral index.
We calculate this change in log-likelihood by only replacing the sample eigenvalues with regularized eigenvalues up to spectral index $i$, for $i=1,\ldots,100$, such that $i=100$ corresponds to the fully regularized model.
The figure shows that the bulk of the out-of-sample likelihood improvement comes from the highest spectral indices ($i>70$).
It underlines the importance of the regularization step for in the new model's set-up for the small eigenvalues jointly with introducing dynamics for the largest eigenvalues.
\section{Conclusion}
\label{sec:Concl}
In this article, we proposed a new dynamic copula model for high-dimensional time-dependent financial applications based on the skewed $t$ copula.
The model does not require variables to be grouped in clusters, but rather uses a spectral decomposition of the copula dependence matrix.
Asymptotic biases in the spectrum are avoided by regularization of the correlation matrix, while the dynamics are captured via time-varying volatilities in the principal spectral dimensions, where the number of dynamic components is kept to a minimum via a model selection procedure based on the BIC.
The combination of all three elements results in a parsimonious, yet flexible time-varying dependence model with limited complexity and increased interpretability.
A simulation study confirmed that the regularized dynamic copula performs well out-of-sample in high-dimensional settings with general dependence structures.
In an empirical study, we showed that regularized score-driven spectral dynamics can be used to study contractions in the dependence structure of the international financial market during crisis times, which reduces diversification potential and enhances systemic risk.
We also showed that regularization outperforms cluster assignments in terms of in-sample and out-of-sample performance for our heterogeneous data set containing a broad range of countries and sectors.
In particular, regularization turned out to be most important at the low end of the spectrum, while the introduction of dynamics was most important for the highest eigenvalues.
\bibliography{references}
\newpage