Approximate Factor Model with S-vine Copula Structure

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.

80,890 characters · 17 sections · 68 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

frontmatter\ead{[email removed]} \address[a]{School for Business and Society, University of York, UK} \cortext[cor1]{Corresponding author.} \ead{[email removed]} \begin{abstract} We propose a novel framework for approximate factor models that integrates an S‐vine copula structure to capture complex dependencies among common factors. Our estimation procedure proceeds in two steps: first, we apply principal component analysis (PCA) to extract the factors; second, we employ maximum likelihood estimation that combines kernel density estimation for the margins with an S‐vine copula to model the dependence structure. Jointly fitting the S‐vine copula with the margins yields an oblique factor rotation without resorting to ad hoc restrictions or traditional projection pursuit methods. Our theoretical contributions include establishing the consistency of the rotation and copula parameter estimators, developing asymptotic theory for the factor-projected empirical process under dependent data, and proving the uniform consistency of the projected entropy estimators. Simulation studies demonstrate convergence with respect to both the dimensionality and the sample size. We further assess model performance through Value-at-Risk (VaR) estimation via Monte Carlo methods and apply our methodology to the daily returns of S&P 500 Index constituents to forecast the VaR of S&P 500 index. \end{abstract} \begin{keyword} Distributional forecasting\sep Empirical process\sep Latent factor\sep Principal component analysis\sep Semiparametric copula model JEL Classification: C14\sep C32\sep C38\sep C55. \end{keyword}

Introduction

Background

Latent factor models and copula models are both essential tools in dependence modelling. Latent factor models are widely employed to uncover the underlying structure of high-dimensional data, offering a parsimonious representation of dependencies among large numbers of variables. Their ability to reduce dimensionality has enabled applications across diverse fields, including economics, social sciences, and environmental studies. In parallel, copula models have emerged as a powerful approach of capturing non-linearity for both serial and cross-sectional dependence among random variables, offering insights that traditional correlation measures may overlook. Vine copulas, in particular, have gained popularity due to their ability to combine multiple bivariate copulas into a cohesive framework, thereby accommodating complex dependence structures. By adding invariance conditions, the stationary vine (S-vine) copula can be obtained, which could facilitate inference in stationary time series.

This paper focuses on the dependence structure of the latent factors of an approximate factor model. Specifically, we consider a high-dimensional random vector, $\mathbf{X}_t$, generated by an approximate factor model:

equation[equation omitted — 138 chars of source]

where $\boldsymbol{\Lambda} = (\lambda_1, \cdots, \lambda_d)^\top$ is the matrix of factor loadings, $\mathbf{F}_t$ is an $K$-dimensional vector of common factors, and $\boldsymbol{\varepsilon}_t = (\varepsilon_{t1}, \cdots, \varepsilon_{td})^\top$ represents the idiosyncratic components uncorrelated with $\mathbf{F}_t$. By incorporating an S-vine copula structure within the common factors, we aim to extend the model's capacity to capture non-linear and asymmetric dependencies.This allows the model to better accommodate complex interactions between variables that may arise in high-dimensional data. The resulting model is capable of representing more intricate dependence patterns than either factor models or copulas alone.

Literature review

Copula models are widely used to model dependence structures in univariate or multivariate time series, capturing both cross-sectional and serial dependence. As the dimensionality of the time series increases, traditional copula models can become intractably complex. Vine copulas, introduced by joe1996families, offer a flexible framework by decomposing a multivariate copula into a sequence of bivariate copulas structured over a series of nested trees. The development of vine copula models has been advanced in a number of works, including cooke1997markov, joe2014dependence, kurowicka2006uncertainty, ACFB09,loaiza2018time and zhao2022modeling.

Various types of vine structure exist, such as D-vines, C-vines, and more general R-vines. In the context of time series modelling, CF06a investigate first-order Markov models as a special case of D-vines. ibragimov2009copula extends D-vine constructions to higher-order Markov models. For multivariate higher-order Markov chains, several alternative structures have been proposed, including D-vines smith2015copula, M-vines beare2015vine, and COPAR models brechmann2015copar. More recently, NKM22 establish general conditions under which a vine copula model yields a stationary multivariate time series. They show that D-vines and M-vines satisfy a translation invariance condition that ensures stationarity. While COPAR structures may be stationary in some cases, they do not satisfy this condition in general.

To address the challenge of modelling dependence in high-dimensional data, copula models have been extended to incorporate factor structures. joe2010tail introduces the factor copula model, in which the observed variables are assumed to be conditionally independent given a latent factor. This framework is generalized by krupskii2013factor to allow for multiple latent factors. A Gaussian factor model implies a Gaussian factor copula, but factor copulas allow for non-Gaussian dependence, enabling richer modelling of tail behaviour and asymmetry. These models have been widely applied in finance and insurance, as discussed in embrechts2009copulas, krupskii2013factor, chen2015multi and others. A special case is the linear factor copula model studied by OP17 and KG18, where observations are linear in factors and idiosyncratic errors. To handle dynamic dependence, OP18 propose copula-based time-varying factor models, and opschoor2021closed extend this to dynamic factor copulas with observation-driven time-varying loadings.

In parallel, the statistical literature on latent factor models offers a complementary approach to modeling high-dimensional data through dimension reduction. While factor copula models typically rely on conditional independence given latent variables, classical factor models assume that the errors are uncorrelated with the factors and with one another. A key advancement is the approximate factor model introduced by CR83, which allows for weakly correlated idiosyncratic errors. When both the number of variables and the sample size are large, approximate factor models can be efficiently estimated using principal component analysis (PCA) BN02. This approach is extensively employed in economics SW02 and enhances large covariance matrix estimation FLM13, which is crucial for effective risk management in finance.

Our approach adopts the approximate factor model framework but shifts focus toward understanding the dependence structure of the common latent factors. In contrast to factor copula models, which impose a conditional independence assumption between observed variables given the latent factors, the approximate factor model places less assumptions on the dependence structure between factors and idiosyncratic errors. Moreover, while the approximate factor model is primarily designed to capture conditional mean relationships, our work instead emphasizes the full distributional properties of the latent factors, including both their cross-sectional and serial dependence. This perspective allows us to capture temporal and cross-sectional dynamics in the underlying factors beyond what traditional mean-based models reveal.

The semiparametric estimation of copula models with non-parametric margins and parametric copula functions is studied by CF06a,CF06b, where the marginal distributions and the copula can be estimated separately. In our setting, however, the marginal distributions and the copula share a common set of parameters, so a full‐likelihood approach is required. Empirical process theory plays a central role in developing the asymptotic theory of maximum likelihood estimation for our model. KWXXY19 study the empirical process in approximate factor models under an i.i.d. factor assumption, whereas our work extends these results to accommodate dependent factors under a mixing condition. To resolve the factor rotation indeterminacy, we examine the projected empirical process of factors, a tool often used in model diagnostics ZL98, E06, X09. The consistency of our MLE hinges on the uniform convergence of this projected empirical process, specifically of the projections $\boldsymbol{v}^\top \boldsymbol{F}_t$, over all unit vectors $\boldsymbol{v}$.

Our estimation of the marginal log-likelihood is closely related to entropy estimation; see MP24 for a comprehensive review. Entropy estimators are widely used in projection pursuit to identify informative directions in high-dimensional data Huber85. In particular, for a one-dimensional projection, one may search for a direction $\boldsymbol{v}$ that maximizes the non-normality by maximizing the entropy $I(\boldsymbol{v}) = -\int g_{\boldsymbol{v}} \log g_{\boldsymbol{v}}$, where $g_{\boldsymbol{v}}$ denotes the density of the scalar projection $\boldsymbol{v}^\top \boldsymbol{F}_t$. Joe89 proposes the integral-form plug‐in entropy estimators based on kernel density estimation, with the leave-one-out variant discussed in IR81. HM93 further extends both estimators without resorting to numerical integration. In particular, HJWW20 establish minimax rates for entropy estimation and demonstrate that integral-form kernel-based plug-in estimators are nearly optimal, exhibiting only logarithmic suboptimality. However, most of the literature assumes i.i.d. samples and only a few studies consider dependent data Lim07,RAMM15. In this paper, we use entropy estimators to estimate the log-likelihood of the marginal distributions of rotated latent factors. Although the uniform consistency of the entropy estimator $I(\boldsymbol{v})$ over $\boldsymbol{v}$ has been established under i.i.d.\ assumptions ZWF95 and more recently in kernel-based likelihood estimation for independent component analysis by HHW25, we extend these results by establishing the uniform consistency of kernel-based plug-in entropy estimators under dependence and by accounting for the estimation error introduced through principal component analysis when estimating latent factors.

Contribution

In this paper, we propose a new semiparametric framework that integrates the approximate factor model with an S-vine copula structure to investigate the dependence among latent common factors. Our method extends principal component analysis by introducing a rotation of the initial factors that aligns them with a semiparametric vine copula model in the second step. This additional step enables the model to accommodate complex, non-Gaussian dependence structures that standard PCA fails to capture.

Unlike traditional approaches that rely on projection pursuit techniques such as quartimax or oblimax Huber85, or on imposed structural constraints like sparsity in factor loadings DD23, our procedure performs a likelihood-based estimation that simultaneously determines an optimal oblique rotation and fits a flexible S-vine copula model. As a result, our method captures rich dependence patterns across time and variables and facilitates distributional forecasting.

Our theoretical contributions are threefold. First, we establish consistency for both the factor rotation and the copula parameter estimates. Second, we develop the asymptotic theory for the empirical process of the latent factors in the presence of temporal dependence, extending existing results to mixing-dependent sequences. Third, we derive uniform convergence results for kernel-based entropy estimators, accounting for the compounded effects of projection, factor estimation error, and serial dependence. These developments collectively support the validity and robustness of our proposed estimation framework.

Organization of the Paper

The remainder of the paper is organized as follows. Section (ref) introduces the approximate factor model with S-vine copula-dependent common factors and outlines a two-step estimation procedure. Then we explores distributional forecasting based on the estimated model. Section (ref) establishes the asymptotic properties of the estimators. In particular, we study the convergence of the empirical process of the estimated factors under arbitrary rotations and prove the consistency of the proposed estimators. Section (ref) presents numerical results and Section (ref) illustrates an empirical application using real data. Section (ref) concludes.

(ref) contains the proofs of the main asymptotic theorems, while (ref) provides additional technical lemmas and their proofs. (ref) discusses the parameterization of the rotation matrix.

Throughout the paper, we use $\Vert\cdot\Vert$ to denote the Euclidean norm of a vector, the operator (or spectral) norm of a matrix, or the supremum norm of a function, depending on the context. Let $\lambda_{\max}(\cdot)$, $\lambda_{\min}(\cdot)$, $\operatorname{Tr}(\cdot)$, and $\operatorname{det}$ denote the maximum eigenvalue, minimum eigenvalue, trace, and determinant of a square matrix, respectively. For a $p\times p$ matrix ${\mathbf W}=(w_{ij})_{p\times p}$, we let $\Vert {\mathbf W}\Vert_F=\operatorname{Tr}^{1/2}\left({\mathbf W}^{^\intercal}{\mathbf W}\right)$ be its Frobenius norm, $\vert{\mathbf W}\vert_1=\sum_{i=1}^p\sum_{j=1}^p |w_{ij}|$, $\Vert{\mathbf W}\Vert_1=\max_{1\leq j\leq p}\sum_{i=1}^p |w_{ij}|$, $\Vert{\mathbf W}\Vert_{\infty}=\max_{1\leq i\leq p}\sum_{j=1}^p |w_{ij}|$, and $\Vert {\mathbf W}\Vert_{\max}=\max_{1\leq i\leq p}\max_{1\leq j\leq p} |w_{ij}|$. Denote $\nabla_{\boldsymbol{x}}{f}$ as the gradient of the function $f$ with respect to the vector $\boldsymbol{x}$. For a generic vector of functions $ \boldsymbol{G} = (G_1, \cdots, G_K) $, write $\boldsymbol{G}(\boldsymbol{x}) = (G_1(x_1), \cdots, G_K(x_K))$ with $\boldsymbol{x}=(x_1,\cdots,x_K)$.

Model and Estimation Procedure

In this section, we first introduce a factor model with S-vine copula dependence common factors and then we develop a two-step estimation procedure. Suppose we have time series observations $X_{it}$ for $t = 1, \cdots , T$ and $i = 1,\cdots, N$. Our goal is to estimate the copula structure of the latent factors.

Model setting: Factor model with S-vine copula dependence common factors

To accommodate strong cross-sectional dependence which is not uncommon for large-scale time series collected in practice, we assume that $X_{it}$ is generated by an approximate factor model:

equation[equation omitted — 99 chars of source]

where $\boldsymbol{\lambda}_i$ is a $K$-dimensional vector of the factor loading, ${\boldsymbol F}_t=(F_{1t},\cdots,F_{Kt})^\top$ is a $K$-dimensional vector of latent factors, and $\varepsilon_{it}$ is an idiosyncratic error. In matrix notation, (ref) can be written as

equation[equation omitted — 100 chars of source]

where ${\bf X}$ is a $T$-by-$N$ matrix, ${\bf F} = ({\boldsymbol F}_1,\cdots, {\boldsymbol F}_T)^\top$ is a $T$-by-$K$ matrix, $\boldsymbol{\Lambda} = (\boldsymbol{\lambda}_1,\cdots, \boldsymbol{\lambda}_N)^\top$ is a $N$-by-$K$ matrix , and $\boldsymbol{\varepsilon}$ is a $T$-by-$N$ error matrix.

In order to capture the dependence structure among factors, we assume that the factors of the approximate factor model follow a Markovian S-vine copula model NKM22.

definition$({\boldsymbol F}_1,\cdots, {\boldsymbol F}_T)$ follows a Markovian S-vine copula model $({\cal V}^F,{\cal C}({\cal V}^F))$ of order $p$, where ${\cal V}^F$ is the vine structure consisting of a collection of trees $T_k$ with vertex sets $V_k^F$ and edge sets $E_k^F$, and ${\cal C}({\cal V}^F)$ is the collection of pair copulas associated with the edges in the vine, if the following conditions are satisfied: (i) The edges of the $(k-1)$-th tree will be the vertices of the $k$-th tree, where $V_k^F=E_{k-1}^F$ for $k=2,\cdots, Kp$. (ii) Assuming vertices $a,b \in V_k^F$ are connected by an edge $e\in E_k^F$ and the corresponding edges $a=\{a_1,a_2\},b=\{b_1,b_2\}\in E_{k-1}^F$, there is exactly one of the element in $a$ equals one of the element in $b$. (iii) There are $Kp-1$ conditioning set for for any edges in the $(Kp-1)$-th tree $T_{Kp-1}$. There is no conditioning set in the first tree. Every conditioned set appears exactly once. (iv) Assuming vertices $a_e,b_e \in V_k$ are connected by an edge $e\in E_k^F$, and edges $a_{e'}$ and $b_{e'}$ are the corresponding edges shift in time by $\tau$ steps, that is, $a_e = a_{e'} + (0, \tau )$ and $b_e = b_{e'} + (0, \tau )$, the pair copulas is translation invariant, that is, $c_{a_e, b_e | {\cal D}_e} = c_{a_{e'}, b_{e'} | {\cal D}_{e'}}$. Here, we adopt standard notation from the vine copula literature. For any edge $ e \in E_k $, the complete union of $ e $, denoted $ U_e $, is defined as \[ {\cal U}_e = \left\{ i \in V_1 \,\middle|\, i \in e_1 \in e_2 \in \cdots \in e \text{ for some } (e_1, \ldots, e_{k-1}) \in E_1 \times \cdots \times E_{k-1} \right\}. \] For a singleton $ i \in V_1 $, we define $ {\cal U}_i := \{i\} $. For an edge $ e \in E_k $ connecting two nodes $ v_1 $ and $ v_2 $, we define the conditioning set as $ {\cal D}_e := {\cal U}_{v_1} \cap {\cal U}_{v_2} $, and the conditioned set as the pair $ (a_e, b_e) := ({\cal U}_{v_1} \setminus {\cal D}_e, {\cal U}_{v_2} \setminus {\cal D}_e) $. We label such an edge by $ e = (a_e, b_e \mid {\cal D}_e) $. NKM22 investigate different S-vine structures, including M-vines, which are proposed by beare2015vine as a special case of S-vine. Within an M-vine framework, the cross-sectional dependence at each time point is represented by a D-vine configuration, wherein the D-vines corresponding to periods $t$ and $t+1$ are interconnected along a common boundary of their respective vine trees. Figure (ref) shows a three dimensional M-vine on three time points, where each column is a D-vine that captures the cross-sectional dependence. \begin{figure}[h!] \caption{Example of first tree level of three-dimensional M-vine on three time points, where $U_{it}$ denotes the probability integral transform of $F_{it}$.} \end{figure}

Estimation procedure: PCA and full maximum likelihood estimation

To estimate the copula structure of the common factors, we first estimate the unobserved factors ${\bf F}$ using principal component analysis. Specifically, we estimate the unobserved factors using $\widehat{\bf F}= (\widehat{{\boldsymbol F}}_1,\cdots, \widehat{{\boldsymbol F}}_T)^\top$ which is the matrix of $K$ eigenvectors (multiplied by $\sqrt{T}$) associated with the $K$ largest eigenvalues of the matrix ${\bf X}{\bf X}^\top/(TN)$ in decreasing order. Without loss of generality, we assume that the first non-zero element of each eigenvector is positive. The factor loading matrix is estimated by $\widehat{\boldsymbol{\Lambda}}={\bf X}^\top\widehat{\bf F}/T$. In practice, the number of factors, $K$, is usually unknown, and we determine it via an information criterion BN02, which minimizes the following function over different values of $k$,

equation*[equation* omitted — 208 chars of source]

Due to the identification issue inherent in factor models, the decomposition ${\bf X} = {\bf F} \Lambda + \boldsymbol{ \varepsilon}$ is not unique. Specifically, for any invertible matrix ${\bf H}$, the model can equivalently be written as ${\bf X} = ({\bf F}{\bf H})({\bf H}^{-1}\Lambda) + \boldsymbol{ \varepsilon}$. As a result, the estimated factors $\widehat{\bf F}$ approximate ${\bf F}$ only up to a rotation. Therefore, directly fitting $\widehat{\bf F}$ to an S-vine copula model may not yield a consistent estimate of the true vine structure.

To address this, we parametrize the rotation matrix as ${\bf H} = {\bf H}({\boldsymbol \theta}_H)$ and jointly estimate both the rotation parameters ${\boldsymbol \theta}_H$ and the copula parameters ${\boldsymbol \theta}_C$ via the maximum likelihood estimation. We assume that the likelihood function admits a unique global maximum (see Assumption (ref)). Specifically, the expected (qausi-)log-likelihood function of ${\bf F}_1,\cdots,{\bf F}_T$, $$ T\log(\operatorname{det}({\bf H}({\boldsymbol \theta}_H)))+T L_{M,HF}({\boldsymbol \theta}_H)+L_{C,HF}({\boldsymbol \theta}_C,{\boldsymbol \theta}_H)$$ attains a unique maximum at $({\boldsymbol \theta}_H^\circ, {\boldsymbol \theta}_C^*)$ with ${\bf H}({\boldsymbol \theta}_H^\circ)={\bf I}_K$, where $L_{M,HF}$ and $L_{C,HF}$ denotes the expected marginal log-likelihood and the expected copula log-likelihood function, respectively, based on the rotated factors ${\bf F}{\bf H}({\boldsymbol \theta}_H)$.

Noting that ${\bf H}$ is the Jacobian matrix Wi62 of the variable transformation from ${\boldsymbol F}_t$ to ${\bf H}{\boldsymbol F}_t$, ${\bf H}({\boldsymbol \theta}_H^\circ)={\bf I}_K$ implies that the log-likelihood function is uniquely maximized when no rotation is applied to the true factors. This ensures the identifiability of the true rotation and guarantees consistency in jointly estimating the rotation and copula parameters. Since the parameterization of the rotation matrix is not unique, the uniqueness of the maximum depends on the specific parameterization. We demonstrate that this identifiability is achievable through simulation in Section (ref). In practice, we normalize each column of the rotation matrix to have unit norm and carefully handle the potential sign-flip issue especially when the copula family is asymmetric. Further discussion on the parameterization of the rotation matrix is provided in (ref).

Assume that the PCA estimator $\widehat{\bf F}$ is an consistent estimator of ${\bf F}{\bf R}_0$, where ${\bf R}_0$ is a deterministic rotation matrix which is defined later in Proposition (ref). Define the c.d.f. and p.d.f. of ${\boldsymbol H}_j^\top{\bf R}_0{\boldsymbol F}_t$ as $G_j(x;{\boldsymbol \theta}_H)$ and $g_{jt}(x;{\boldsymbol \theta}_H)$, respectively, for $j=1,\cdots,K$, where ${\boldsymbol H}_j$ is the $j$-th column of ${\bf H}$. We can define the empirical distribution function of ${\boldsymbol H}_j^\top{\bf R}_0{\boldsymbol F}_t$ as follows,

equation[equation omitted — 166 chars of source]

for $j=1,\cdots,K$, and the pseudo-observations as $$\widehat {\boldsymbol U}_{t}({\boldsymbol \theta}_H)=\widehat {\boldsymbol G}({\bf H}\widehat{\boldsymbol F}_t; {\boldsymbol \theta}_H)=(\widehat G_1({\boldsymbol H}_j^\top\widehat{\boldsymbol F}_t; {\boldsymbol \theta}_H),\cdots,\widehat G_K({\boldsymbol H}_j^\top\widehat{\boldsymbol F}_t; {\boldsymbol \theta}_H))^\top.$$ To implement maximum likelihood estimation, we need also estimate the negative differential entropy of ${\boldsymbol H}_j^\top{\bf R}_0{\boldsymbol F}_t$, that is $\operatorname{E}[\log(g_{j}({\boldsymbol H}_j^\top{\bf R}_0{\boldsymbol F}_t;{\boldsymbol \theta}_H))]$. This can be estimated by using the kernel density estimator HMSW04, $$\widehat g_{j}(x;{\boldsymbol \theta}_H)= \frac{1}{Tb_j}\sum_{s=1}^T W(({\boldsymbol H}_j^\top\widehat{\boldsymbol F}_s-x)/b_j),$$ or the leave-one out kernel density estimator, $$\widehat g_{jt}(x;{\boldsymbol \theta}_H)= \frac{1}{(T-1)b_j}\sum_{s\neq t} W(({\boldsymbol H}_j^\top\widehat{\boldsymbol F}_s-x)/b_j),$$ which results in the leave-one out kernel entropy estimator MBND24,

equation[equation omitted — 186 chars of source]

for $j=1,\cdots,K$, where $W$ is a kernel function and $b_j$ is the bandwidth. Hall88 argues that for the nonnegative kernels we should use a window of size between $T^{- 1/3}$ and $T^{- 1/4}$ for entropy estimation. Thus we choose $b_j=\widehat\sigma_j T^{-1/4}$ as the bandwidth, where $\widehat\sigma_j$ is the standard deviation estimate of ${\boldsymbol H}_j^\top\widehat{\boldsymbol F}_t$. Then, the likelihood function of ${\bf F}$ can be defined as

equation[equation omitted — 336 chars of source]

with bandwidth ${\boldsymbol b}=(b_1,\cdots,b_K)$ and $\widehat L_{M}({\boldsymbol \theta}_H)=\sum_{j=1}^K\widehat L_{G_j}({\boldsymbol \theta}_H)$ being the estimate of $L_{M}({\boldsymbol \theta}_H):=\sum_{j=1}^K\operatorname{E}[\log(g_{j}({\boldsymbol H}_j^\top{\bf R}_0{\boldsymbol F}_t;{\boldsymbol \theta}_H))]$ and $\widehat L_C$ is the log-likelihood of the S-vine copula which is defined later in (ref).

By maximizing the objective function in (ref), we obtain the estimates $\widehat{\boldsymbol \theta}_C$ and $\widehat{\boldsymbol \theta}_H$. Let $\widehat{\bf H} = {\bf H}(\widehat{\boldsymbol \theta}_H)$. The factors are then recovered as $\widehat{\bf F}\widehat{\bf H}$ with the S-vine copula parameters $\widehat{\boldsymbol \theta}_C$. In practice, the vine structure is usually unknown, and we follow the procedure as in Appendix S3.1 of NKM22 to determine it, also see DBCK13.

Summary and further discussion

We summarize the estimation procedure:

description• Obtain $\widehat{\boldsymbol F}_t$, for $t=1,\cdots,T$ using PCA as discussed at the beginning of Section (ref). • For given ${\boldsymbol \theta}_H$, calculate the pseudo-observations, $\widehat {\boldsymbol U}_{t}({\boldsymbol \theta}_H)$ and estimate the parameter $\widehat{\boldsymbol \theta}_C({\boldsymbol \theta}_H)$ for the S-vine copula. Then the estimator of ${\boldsymbol \theta}_H$ is defined as \begin{equation} \widehat{\boldsymbol \theta}_H=\operatorname*{arg\,max}_{{\boldsymbol \theta}_H \in{\boldsymbol \Theta}_H} \widehat L_{T,{\bf b}}(\widehat{\bf F}{\bf H};{\boldsymbol \theta}_H,\widehat{\boldsymbol \theta}_C({\boldsymbol \theta}_H)), \end{equation} where $\widehat L_{T,{\bf b}}$ is defined in (ref). Then we let $\widehat{\boldsymbol \theta}_C=\widehat{\boldsymbol \theta}_C(\widehat{\boldsymbol \theta}_H)$. Thus, we obtain the estimates of the original factors and factor loadings (with S-vine dependence structure), respectively, $$\widetilde{\bf F}=\widehat{\bf F}{\bf H}(\widehat{\boldsymbol \theta}_H)\ \ \textit{ and }\ \ \widetilde{\boldsymbol{\Lambda}}:=\widehat{\boldsymbol{\Lambda}}({\bf H}(\widehat{\boldsymbol \theta}_H)^{-1})^\top.$$

Before concluding this section, we discuss the relationship between the stepwise estimation and Z-estimation. We say that two edges $ e $ and $ e' $ are translation invariance of each other if \[ a_e = a_{e'} + (\tau, 0), \quad b_e = b_{e'} + (\tau, 0), \quad {\cal D}_e = \{ v + (\tau, 0) : v \in {\cal D}_{e'} \}. \] We write $ e \sim e' $ to indicate that $ e $ and $ e' $ are translation invariance. This defines an equivalence relation on the set of edges. The corresponding equivalence class of an edge $ e $ is denoted by \[ [e] := \{ e' : e' \sim e \}. \]

Since joint maximization of the full likelihood is often computationally intensive, following NKM22, we adopt the stepwise maximum likelihood estimator of ACFB09, which estimates the parameters of each pair-copula sequentially. Accordingly, we define

equation[equation omitted — 241 chars of source]

and for each representative edge $ e' \in E_k^F $, we define

eqnarray*[eqnarray* omitted — 500 chars of source]

where $ \boldsymbol{\theta}_C $ denotes the collection of all $ \boldsymbol{\theta}_{[e]} $, and ${\boldsymbol{\theta}}_{S_a([e])}$ and ${\boldsymbol{\theta}}_{S_b([e])}$ denote the parameter estimates from previous step.

Define

equation*[equation* omitted — 459 chars of source]

and

equation*[equation* omitted — 358 chars of source]

Let ${\boldsymbol s}_{{\boldsymbol \theta}_C, \boldsymbol {\nu}}$ be a vector of functions stacking $\boldsymbol{s}_{[e],{\boldsymbol \theta}_C}$ for $[e]\in E_i$ and $i=1,\cdots, (p+1)K$ and denote ${\boldsymbol s}_{{\boldsymbol \theta}_C, \boldsymbol {\nu}}\left(\boldsymbol{x}_t, \cdots, \boldsymbol{x}_{t+p}\right)={\boldsymbol s}_{{\boldsymbol \theta}_C}\left(\boldsymbol {\nu}(\boldsymbol{x}_t), \cdots, \boldsymbol {\nu}(\boldsymbol{x}_{t+p})\right)$. For given ${\boldsymbol \theta}_H$, we can rewrite the semiparametric step-wise MLE as the Z-estimation as NKM22, that is, $\widehat{\boldsymbol \theta}_C({\boldsymbol \theta}_H)$ solves the following estimation equation,

equation*[equation* omitted — 197 chars of source]

Furthermore, let

equation*[equation* omitted — 412 chars of source]

and ${\boldsymbol s}_{{\boldsymbol \theta}_H, \boldsymbol {\nu}}$ be a vector stacking $\boldsymbol{s}_{[e],{\boldsymbol \theta}_H,{\boldsymbol \nu}}$ for $[e]\in E_i$ and $i=1,\cdots, (p+1)K$, where

eqnarray[eqnarray omitted — 469 chars of source]

Then our proposed estimator can represent a Z-estimator of the following estimation equation,

empheq[left=\empheqlbrace]{align} &\frac{1}{T-p}\sum_{t=1}^{T-p}{\boldsymbol s}_{{\boldsymbol \theta}_C,\widehat{\boldsymbol{G}}}(\widehat{\boldsymbol F}_t,\cdots,\widehat {\boldsymbol F}_{t+p}) = {\boldsymbol 0}, \nonumber \\ &\boldsymbol{s}_{1,{\boldsymbol \theta}_H}+\widehat{\boldsymbol{s}}_{2,{\boldsymbol \theta}_H}+\frac{1}{T-p}\sum_{t=1}^{T-p} \boldsymbol{s}_{{\boldsymbol \theta}_H,\widehat{\boldsymbol{G}}}\left(\widehat{\boldsymbol F}_t, \cdots, \widehat{\boldsymbol F}_{t+p}\right)={\boldsymbol 0}. \nonumber

Forecasting

Following NKM22, we construct the predictive distribution via Monte Carlo simulation. Recall that $\widetilde{\boldsymbol{F}}_t$ and $\widetilde{\boldsymbol{\Lambda}}$ denote the estimated factors and the factor loadings matrix after proper rotation, respectively, and that $\widehat{\boldsymbol{\varepsilon}}_t$ denotes the residuals at time $t$. Our forecasting procedure is as follows:

description• Construct the empirical marginal distribution of $\widetilde{\boldsymbol{F}}_t$, and calculate the rank of $\widetilde{\boldsymbol{F}}_t$ for $t=1,\cdots, T$, and denoted as $\widetilde{\boldsymbol{U}}_t$. • Simulate $M$ paths $(\widetilde{\boldsymbol{U}}_{T+1}^{(m)}, \cdots,\widetilde{\boldsymbol{U}}_{T+h}^{(m)})$ for $m=1,\cdots, M$, from the S-vine copula model with starting values $\widetilde{\boldsymbol{U}}_{1}, \cdots,\widetilde{\boldsymbol{U}}_{T}$, where $h$ is the forecast horizon. • Transform $(\widetilde{\boldsymbol{U}}_{T+1}^{(m)}, \cdots,\widetilde{\boldsymbol{U}}_{T+h}^{(m)})$ into $(\widetilde{\boldsymbol{F}}_{T+1}^{(m)},\cdots, \widetilde{\boldsymbol{F}}_{T+h}^{(m)})$, for $m=1,\cdots, M$, in factor space by quantile mapping via the empirical marginal distribution of $\widetilde{\boldsymbol{F}}_t$. • Simulate $\widetilde{\boldsymbol{X}}_{T+k}^{(m)} = \widetilde{\boldsymbol{\Lambda}} \widetilde{\boldsymbol{F}}_{T+k}^{(m)} + \widetilde{\boldsymbol{\varepsilon}}_{T+k}^{(m)}$ for $k=1,\cdots,h$ and $m=1,\cdots,M$. where $\widetilde{\boldsymbol{\varepsilon}}_{T+k}^{(m)}$ is resampled from $\{\widehat{\boldsymbol{\epsilon}}_1,...,\widehat{\boldsymbol{\epsilon}}_T\}$ with replacement. • The predictive distribution at time $T+k$ is then given by the ensemble $\widetilde{\boldsymbol{X}}_{T+k}^{(m)}$ for $m=1,\cdots, M$. Quantiles of this distribution provide prediction intervals, while the ensemble mean serves as the point forecast.

Alternative approaches for forecasting the idiosyncratic errors $\boldsymbol{\varepsilon}_t$ can also be considered. For instance, one may fit univariate time series models, such as ARMA-GARCH with specified innovation distributions, to each component of $\boldsymbol{\varepsilon}_t$. While such methods capture serial dependence within each series, they neglect potential cross-sectional dependence across components.

Asympotics

We next list some regularity conditions which are sufficient for deriving the uniform consistency results for $\widehat{\boldsymbol F}_t$, $\widehat{\boldsymbol\lambda}_i$, and $\widehat\varepsilon_{t,i}$.

assumption{ (i) The process $\{({\boldsymbol F}_{t}^\top, {\boldsymbol\varepsilon}_t^\top)^\top \}_{t=1}^\infty$ is stationary and $\beta$-mixing with the mixing coefficient $\beta_l$ satisfying $\beta_l=O\left(\rho^l\right)$, where $\rho$ is a constant that satisfies $ 0<\rho<1$.} { (ii) There exists a positive definite matrix ${\boldsymbol\Sigma}_\Lambda$ such that \[ \left\Vert \frac{1}{N}\sum_{i=1}^N {\boldsymbol\lambda}_{i}{\boldsymbol\lambda}_{i}^\top-{\boldsymbol\Sigma}_\Lambda\right\Vert=O\left(N^{-1/2}\right). \] and $\Vert{\boldsymbol\lambda}_{i}\Vert$ is uniformly bounded over $1\leq i\leq N$, where ${\boldsymbol\lambda}_{i}$ is the $i$-th column of ${\boldsymbol\Lambda}^\top$. } { (iii) The covariance matrix ${\boldsymbol\Sigma}_F:=\mathsf{Var}({\boldsymbol F}_t)$ is positive definite. In addition, there exists a constant $c_1>0$ such that $\operatorname{E}\left[\exp\left(c_1\Vert {\boldsymbol F}_t\Vert^2\right)\right]<\infty$.} { (iv) The idiosyncratic errors satisfy $$\operatorname{E}[{\boldsymbol\varepsilon}_t]={\bf 0},\ \ \ \operatorname{E}[\varepsilon_{t,i}{\boldsymbol F}_t]={\bf 0},\ \ \ {\rm and}\ \max_{1\leq i\leq N}\max_{1\leq t\leq T}\operatorname{E}\left[\exp\left(c_1 \varepsilon_{t,i}^2\right)\right]<\infty,$$ where $c_1$ is defined in (iii).} { (v) There exist $0<c_2<\infty$ and $\delta_{\varepsilon}>2$ such that} \[ \max_{1\leq t\leq T}\operatorname{E}\left[\left\Vert\sum_{i=1}^N {\boldsymbol\lambda}_{i}\varepsilon_{t,i}\right\Vert^{\delta_{\varepsilon}}\right]\leq c_2 N^{\delta_{\varepsilon}/2},\ \ \max_{1\leq s,t\leq T}\operatorname{E}\left[\left|\sum_{i=1}^N \left[\varepsilon_{s,i}\varepsilon_{t,i}-\operatorname{E}\left(\varepsilon_{s,i}\varepsilon_{t,i}\right)\right]\right|^{\delta_{\varepsilon}}\right]\leq c_2 N^{\delta_{\varepsilon}/2}. \] (vi) {\em The eigenvalues of ${\boldsymbol\Sigma}_\Lambda^{1/2}{\boldsymbol\Sigma}_F{\boldsymbol\Sigma}_\Lambda^{1/2}$ are bounded and distinct.} (vii) {\em For any unit vector $h \in \mathbb{R}^K$, the density function of the projection $h^\top {\boldsymbol F}_t$ exists and is denoted by $g(x)$. This density has bounded support $\operatorname{supp} g$, satisfies the following conditions for some $\delta > 0$: \[ \sup_{\| h \| = 1}\sup_{x\in \operatorname{supp}g} g(x)\geq c >0, \ \ \ \ \sup_{\| h \| = 1} \sup_{x \in \operatorname{supp}g} \left| \nabla_x g(x) \right| < \infty\] and \[ \sup_{\| h \| = 1} \operatorname{E} \left[ \left| \log\left( g(h^\top {\boldsymbol F}_t) \right) \right|^{1+\delta} \right] < \infty . \] }

Most of the above assumptions are standard in the approximate factor model estimation theory BN02, FLM13, CLLL25. The sub-Gaussian moment conditions in Assumption (ref)(iii) and (iv) are required to cover the ultra-high dimensional case where $N$ diverges at an exponential rate of $T$, and can be weakened if $N$ diverges at a polynomial rate of $T$. Assumption (ref)(vi) ensures the identification of the eigenvectors and Assumption (ref)(vii) are required to estimate the differential entropy of the rotated factors HHW25.

Define a $K\times K$ (stochastic) rotation matrix

equation[equation omitted — 237 chars of source]

where ${\boldsymbol V}_{NT}$ is a $K\times K$ diagonal matrix with the $K$ largest eigenvalues of ${\mathbf X}{\mathbf X}^\top/(NT)$ (arranged in a descending order) on its main diagonal. The following proposition gives the uniform consistency results for $\widehat{\boldsymbol F}_t$, $\widehat{\boldsymbol\lambda}_i$, and $\widehat\varepsilon_{t,i}$, which are comparable to those obtained in the existing literature BN02, FLM13,LLF23,CLLL25.

proposition{ Suppose that Assumption (ref)(i)--(v) is satisfied, $T^{4/\delta_{\varepsilon}} \ll N \ll \exp\{T^{1/5}\}$ with $\delta_{\varepsilon}$ defined in Assumption (ref)(v). Then, we have the following uniform consistency results: \begin{equation} \max_{1\leq t\leq T}\left\Vert\widehat{\bf{F}}_t-{\bf R}_{NT}{\bf{F}}_t\right\Vert=O_P\left(\frac{(\log T)^{1/2}}{T}+\frac{T^{2/\delta_{\varepsilon}}}{N^{1/2}}\right). \end{equation} Furthermore, if Assumption (ref)(vi) is satisfied, we have \begin{equation} \left\Vert{\bf R}_{NT}-{\bf R}_0\right\Vert =O_P\left(\frac{1}{T^{1/2}}+\frac{T^{2/\delta_{\varepsilon}}}{N^{1/2}}\right), \end{equation} and, as a consequence, \begin{equation} \max_{1\leq t\leq T}\left\Vert\widehat{\bf{F}}_t-{\bf R}_0{\bf{F}}_t\right\Vert=O_P\left((\log T)^{1/2}\left(\frac{1}{T^{1/2}}+\frac{T^{2/\delta_{\varepsilon}}}{N^{1/2}}\right)\right). \end{equation} where ${\bf R}_0={\boldsymbol V}_0^{-1/2}{\bf U}_0^\top{\boldsymbol \Sigma}_\Lambda^{1/2}$, ${\bf U}_0$ is a matrix consisting of the eigenvectors of ${\boldsymbol\Sigma}_\Lambda^{1/2}{\boldsymbol\Sigma}_F{\boldsymbol\Sigma}_\Lambda^{1/2}$, $v_{i,0}$ is the $i$-th largest eigenvalue of ${\boldsymbol\Sigma}_\Lambda^{1/2}{\boldsymbol\Sigma}_F{\boldsymbol\Sigma}_\Lambda^{1/2}$, and ${\boldsymbol U}_0=\operatorname{diag}\{v_{1,0},\cdots,v_{q,0}\}$. }

Proposition (ref) shows that $\widehat{\boldsymbol F}_t$ are consistent estimators of the rotated latent factors ${\bf R}_{0}{\bf{F}}_t$, rather than the factors themselves (unless ${\bf R}_{0}=\bf{I}$). In addition, if we assume that $N\gg T^{4/\delta_{\varepsilon}+1}$, the term $T^{2/\delta_{\varepsilon}}/N^{1/2}$ in ((ref)) would disappear.

Recall that $\widehat G_j(x;{\boldsymbol \theta}_H)=T^{-1}\sum_{t=1}^T {\boldsymbol 1}({\boldsymbol H}_j^\top\widehat{\boldsymbol F}_t\leq x)$ is the marked empirical process of the common factors. Define $\widetilde G_j(x;{\boldsymbol \theta}_H)= T^{-1}\sum_{t=1}^T {\boldsymbol 1}({\boldsymbol H}_j^\top{\bf R}_0{\boldsymbol F}_t\leq x)$ as the oracle empirical process of the common factors. We have the following theorem.

theoremSuppose that Assumptions (ref) are satisfied and $T^{4/\delta_{\varepsilon}+1} \ll N \ll \exp\{T^{1/5}\}$ with $\delta_{\varepsilon}$ defined in Assumption (ref)(v). We have \begin{equation} \sqrt{T}\sup_{{\boldsymbol \theta}_H \in {\boldsymbol \Theta}_H}\sup _{x \in \mathbb{R}}\left|\widehat G_{j}(x;{\boldsymbol \theta}_H)- \widetilde G_j(x;{\boldsymbol \theta}_H)+\frac{1}{T}\sum_{t=1}^Tg_{j}(x;{\boldsymbol \theta}_H)z_{{\boldsymbol \theta}_H t}\right| \rightarrow_P 0. \end{equation} where $z_{{\boldsymbol \theta}_H t}={\boldsymbol H}_j^\top( \widehat{\boldsymbol F}_t-{\bf R}_0{\boldsymbol F}_t).$ In particular, we have \begin{equation} \sup_{{\boldsymbol \theta}_H \in {\boldsymbol \Theta}_H}\sup _{x \in \mathbb{R}}\left|\widehat G_{j}(x;{\boldsymbol \theta}_H)- G_j(x;{\boldsymbol \theta}_H)\right|\rightarrow_P 0. \end{equation}

remarkKWXXY19 establish the consistency of the marked empirical processes for the common components and idiosyncratic errors in an approximate factor model estimated via principal component analysis, under the assumption that the factors are independent and identically distributed (i.i.d.) random vectors, independent of the idiosyncratic component series. In comparison, our paper provides a novel consistency result for the marked empirical processes based on arbitrary rotations of the estimators of the common factors, relaxing these assumptions to accommodate a mixing condition.

Some further assumptions are needed for the maximum likelihood estimation.

assumption(i) The kernel $W(\cdot)$ is a twice continuously differentiable, symmetric probability density function with compact support on the interval $[-1, 1]$.\\ (ii) The bandwidth $b$ satisfies $b \to 0$ and $\ Tb^2/ \log^2(T) \to \infty$.

Define the pseudo-true values ${\boldsymbol \theta}^*=(({\boldsymbol \theta}^*_H)^\top,({\boldsymbol \theta}^*_C)^\top)^\top$ as the solution to

equation*[equation* omitted — 222 chars of source]

where

equation[equation omitted — 410 chars of source]

and

equation[equation omitted — 351 chars of source]

assumption(i) The pseudo-true values $(({\boldsymbol \theta}^*_H)^\top,({\boldsymbol \theta}^*_C)^\top)^\top$ lie in the interior of $\Theta_H\times \Theta_C$ and for every $\epsilon>0$, $$\inf_{\Vert {\boldsymbol\theta}-{\boldsymbol\theta}^*\Vert>\epsilon}\Vert \operatorname{E}\left[{\boldsymbol \phi}_{\boldsymbol \theta,\boldsymbol {G}}({\bf R}_0{\boldsymbol F}_1,\cdots,{\bf R}_0{\boldsymbol F}_{1+p})+{\boldsymbol \psi}_{{\boldsymbol \theta}_H}\right]\Vert>0.$$ (ii) The function ${\boldsymbol \phi}_{\boldsymbol{\theta},\boldsymbol{G}}$ is continuously differentiable with respect to $\boldsymbol\theta$ and its arguments and there is $\delta>0$ such that \begin{equation*}\operatorname{E}\left[\sup_{{\boldsymbol \theta}\in {\boldsymbol \Theta}} \left\Vert {\boldsymbol \phi}_{\boldsymbol{\theta},\boldsymbol{G}}({\bf R}_0{\boldsymbol F}_1,\cdots,{\bf R}_0{\boldsymbol F}_{1+p})\right\Vert<\infty,\quad \sup_{{\boldsymbol \theta}\in {\boldsymbol \Theta}}\left\Vert \nabla_{\boldsymbol{\theta}}{\boldsymbol \phi}_{\boldsymbol{\theta},\boldsymbol{G}}({\bf R}_0{\boldsymbol F}_1,\cdots,{\bf R}_0{\boldsymbol F}_{1+p})\right\Vert\right]<\infty, \end{equation*} and \begin{equation*}\operatorname{E}\left[\sup_{{\boldsymbol \theta}\in {\boldsymbol \Theta}}\sup_{{\boldsymbol{\nu}}\in \mathcal{F}_\delta} \left\Vert \frac{\partial}{\partial \{\nu_{j}({\boldsymbol \theta}_1,{\bf R}_0{\boldsymbol F}_s)\}}{\boldsymbol \phi}_{\boldsymbol{\theta},\boldsymbol{\nu}}({\bf R}_0{\boldsymbol F}_1,\cdots,{\bf R}_0{\boldsymbol F}_{1+p})\right\Vert\right]<\infty, \end{equation*} for $j=1,\cdots,K$, and $s=1,\cdots,1+p$, where $$\mathcal{F}_\delta=\{{\boldsymbol{\nu}}:\Vert{\boldsymbol{\nu}}-{\boldsymbol{G}}\Vert\leq\delta\}.$$ (iii) The matrix derivatives $\nabla_{\boldsymbol{\theta}}\nabla_{\boldsymbol{x}^\dag}{\boldsymbol \phi}_{\boldsymbol{\theta},\boldsymbol{G}}({\boldsymbol x}^\dag)$ are continuous in ${\boldsymbol x}^\dag$ where ${\boldsymbol x}^\dag=({\boldsymbol x}_1,\cdots,{\boldsymbol x}_{1+p})$. There is $\delta>0$ such that \begin{equation*} \sup_{{\boldsymbol{\nu}}\in \mathcal{F}_\delta} \Vert\nabla_{\boldsymbol{x}^\dag}{\boldsymbol \phi}_{\boldsymbol{\theta},\boldsymbol{\nu}}\Vert<\infty \end{equation*} and \begin{equation*}\operatorname{E}\left[\sup_{{\boldsymbol \theta}\in {\boldsymbol \Theta}}\sup_{{\boldsymbol{\nu}}\in \mathcal{F}_\delta} \left\Vert \frac{\partial}{\partial \{\nu_{j}({\boldsymbol \theta}_H,{\bf R}_0{\boldsymbol F}_s)\}}\nabla_{\boldsymbol{\theta}}{\boldsymbol \phi}_{\boldsymbol{\theta},\boldsymbol{\nu}}({\bf R}_0{\boldsymbol F}_1,\cdots,{\bf R}_0{\boldsymbol F}_{1+p})\right\Vert\right]<\infty, \end{equation*} for $j=1,\cdots,K$, and $s=1,\cdots,1+p$.

Assumption (ref)(i) ensures the identifiability of the model parameters and (ref)(iv) is a standard regularity condition. Assumption (ref)(ii)(iii) are standard in the semiparametric copulas model estimation theory NKM22.

theoremSuppose that Assumptions (ref)--(ref) are satisfied and $T^{4/\delta_{\varepsilon}+1} \ll N \ll \exp\{T^{1/5}\}$ with $\delta_{\varepsilon}$ defined in Assumption (ref)(v). We have $$\widehat{\boldsymbol \theta}_H {\rightarrow_p} {\boldsymbol \theta}^*_H,\ \ \widehat{\boldsymbol \theta}_C{\rightarrow_p} {\boldsymbol \theta}^*_C,$$ and \begin{equation} \max_{1\leq t\leq T}\left\Vert{\bf H}(\widehat{\boldsymbol \theta}_H)\widehat{\boldsymbol{F}}_t-{\boldsymbol{F}}_t\right\Vert=O_P\left((\log T)^{1/2}\left[\left(\frac{\log N}{T}\right)^{1/2}+\frac{T^{2/\delta_{\varepsilon}}}{N^{1/2}}\right]\right). \end{equation}

comment\begin{theorem} Suppose that Assumptions (ref)--(ref) are satisfied. We have $$T^{1/2}\left((\widehat{\boldsymbol \theta}_H, \widehat{\boldsymbol \theta}_C)^\top-({\boldsymbol \theta}^*_H, {\boldsymbol \theta}^*_C)^\top\right)\to_d {\cal N}(0,{\bf I}_{\boldsymbol{\theta}^*}{\bf J}_{\boldsymbol{\theta}^*}{\bf I}_{\boldsymbol{\theta}^*}),$$ where ${\bf I}_{\boldsymbol{\theta}^*}=\operatorname{E}[\boldsymbol{Z}_1\boldsymbol{Z}_1^\top]+2\sum_{t=2}^\infty\operatorname{E}[\boldsymbol{Z}_1\boldsymbol{Z}_t^\top]$ and $$\boldsymbol{Z}_t=$$ \end{theorem}

comment\section{Extension} We further assume $\{e_{1T},\cdots,e_{NT}\}$ follows a conditional vine copula model, $({\cal V}^e,{\cal C}({\cal V}^e))$, conditional on ${\boldsymbol F}_{t-1},\cdots,{\boldsymbol F}_{t-p}$, where ${\cal V}^e$ is the vine structure and ${\cal C}({\cal V}^e)$ is the sequence of pair copulas. With the conditional S-vine copula structure, we allow the serial dependence of $e_t$ and the dependence between $e_t$ and ${\boldsymbol F}_t$. \begin{assumption} ${\cal V}^e=(V_k^e, E_k^e)_{k=1}^{d}$ is a collection of trees where $V_k^e$ is the set of vertices and $E_k$ is the set edges of the $k$-th tree. Additionally, (i) $V_k^e=E_{k-1}^e$ for $k=2,\cdots$ (ii) if vertices $a,b \in V_k^e$ are connected by an edge $e\in E_k^e$, then the corresponding edges $a=\{a_1,a_2\},b=\{b_1,b_2\}\in E_{k-1}^e$ must share a common vertex. \end{assumption} \begin{assumption} (i) The conditional joint distribution of $\{e_{1t}, \cdots, e_{Nt}\}$ given history is the same as the conditional joint distribution of $\{e_{1t}, \cdots, e_{Nt}\}$ given ${\boldsymbol F}_{t-1},\cdots,{\boldsymbol F}_{t-p}$. (ii) The conditional joint distribution of $\{F_{1t}, \cdots, F_{Kt}\}$ given history is the same as the conditional joint distribution of $\{F_{1t}, \cdots, F_{Kt}\}$ given ${\boldsymbol F}_{t-1},\cdots,{\boldsymbol F}_{t-p}$. (iii) $\{e_{1t}, \cdots, e_{Nt}\}$ are independent with $\{F_{1t},\cdots,F_{Kt}\}$ given are conditional independent ${\boldsymbol F}_{t-1},\cdots,{\boldsymbol F}_{t-p}$. \end{assumption}

Simulation

In this section, we first discuss the practice issue in identification and then present two simulated examples to examine the finite-sample performance of our proposed methods. For each example, we assess the accuracy of the estimate when the number of common factors and the structure of the S-vine are correctly specified.

Simulation study on identification issues

The S-vine structure and the copula families remain invariant under scaling, permutation, and sign-flipping transformations (see Proposition (ref) in (ref)). While this invariance holds theoretically, practical implementations necessitate restricting the search to a finite set of copula families. We validate this proposition through simulation studies and explore whether additional transformations could introduce identification challenges when the set of admissible copula families is restricted.

To illustrate, we simulate two factors generated from an S-vine copula with 1,000 observations. The vine structure is demonstrated in Figure (ref). The copula families considered include Gaussian, Clayton, Frank, and Joe, while the marginal distributions are specified as either Gaussian or t4 distributions. The parameters of the copulas are shown in Table (ref).

figure[figure omitted — 282 chars of source]

table[table omitted — 1,940 chars of source]

The Figures (ref) and (ref) depict the log-likelihood (LLH) function contours after applying the affine transformation:

equation*[equation* omitted — 121 chars of source]

for $\theta_1,\theta_2\in[0,4\pi]$, in order to reveal the periodic pattern.

Figure (ref) shows that when the marginal distributions are $t_4$, all the LLH contour plots achieve $(0.5\pi, 0)$, $(0, 0.5\pi)$ and have a period of $\pi$ in both directions. This symmetry aligns with the estimation procedure’s allowance for copula rotations (e.g., 0$^{\circ}$, 90$^{\circ}$, 180$^{\circ}$, 270$^{\circ}$), as in NKM22, which corresponds to sign-flipping transformations of factors.

Regarding the Gaussian marginal distributions, Clayton, Frank, and Joe copulas retain the LLH maxima at $(0.5\pi, 0)$ and $(0, 0.5\pi)$, but not for the Gaussian copula. This discrepancy arises because rotating data in most case forces misspecification: the rotated joint distribution no longer resides within the original copula family, leading to inferior LLH values. The Gaussian copula with Gaussian margins, however, constitutes an exception. Affine transformations preserve the multivariate Gaussian structure, ensuring the copula family remains unchanged under rotation.

figure[figure omitted — 1,018 chars of source]

figure[figure omitted — 1,072 chars of source]

Although for Gaussian data, we may not be able to identify the specific rotation due to the invariance properties of the multivariate Gaussian distribution, this underspecification does not imply misspecification. As a result, the model can still achieve strong forecasting performance despite the lack of identifiability. We demonstrate this with an simulated example in Section (ref).

Performance of estimators

In order to quantify the assessment of the estimations, we compute the Root Mean Square Error (RMSE) for the parameters, common factors, and factor loadings as follows

eqnarray[eqnarray omitted — 624 chars of source]

where $ \dim({\boldsymbol \theta}_C) $ is the number of parameters in the vine copula. In practice, we sometimes obtain the flipping estimated factors, which has no influence on the same likelihood but leads to different estimates of the vine copulas and rotaion matrix parameters, ${\boldsymbol \theta}_C$ and ${\boldsymbol \theta}_H$. Therefore, we calculate the RMSEs based on the best sign-flipped estimates of the factors and the corresponding $\widehat{\boldsymbol \theta}_C$ and $\widehat{\boldsymbol \theta}_H$.

We now give a simulated example of an approximate factor model with two factors having a vine copula structure beare2015vine. The vine structure is presented in Figure (ref), which is an M-vine of second-order. All pair copulas in the vine are Frank copulas and Table (ref) shows the specification of the M-vine copula model with the vine structure and associated parameters for each edge.

figure[figure omitted — 283 chars of source]

table[table omitted — 823 chars of source]

The common factors in the approximate factor model are generated by $\Phi^{-1}(U_{it})$, for $i=1,2$ and $t=1,\cdots,T$, where $\Phi(\cdot)$ is the c.d.f. of the standard normal distribution and $U_{it}$ are generated from the M-vine model with parameters specified in Table (ref). Factor loadings are drawn independently from a normal distribution with mean and variance both equal to one. The idiosyncratic errors follow a zero-mean, unit-variance Gaussian AR(1) process with an autoregressive coefficient of 0.5. With the simulated factors, factor loadings, and idiosyncratic errors, we calculate the observations $X_{it}$ using the approximate factor model (ref), where the dimension $d$ is equal to 100, 200 or 500, and the number of observation $n$ is equal to 100, 250, 500, 750, 1000, or 2000. The simulations are repeated 200 times. The means of the RMSE for parameters, common factors, and factor loadings for the 200 simulations are presented in Table (ref)--(ref), respectively.

table[table omitted — 526 chars of source]

table[table omitted — 870 chars of source]

table[table omitted — 872 chars of source]

The three tables show that the RMSEs of the estimated parameters, common factors, and factor loadings converge toward zero as both the dimension and the number of observations increase. Our method performs well even in high-dimensional settings where $n$ is smaller than $d$. On the one hand, the accuracy of factor estimation improves as the dimension increases, due to more information being available for identifying the common structure. On the other hand, a larger number of observations leads to more accurate parameter estimation in the second step. Overall, the number of observations has a greater impact on estimation accuracy than the dimension.

Forecasting performance

In this section, we compare the one-step prediction performance of our S-vine factor models (SF) with Frank copula sequences to the SF model with mispecified pair copulas (Gaussian copulas), GARCH processes with skewed student t innovation distribution, and the dynamic factor model (DFM), where for forecasting using DFM see HT23. Factors are generated from the two-variable second-order M-vine with pair copulas and parameters presented in Table (ref). We forecast with the proposed SF model (SF_1) and compare it with the SF model without factor rotation (SF_2) and the SF model with a misspecified pair copula (SF_wrong) uses Gaussian copulas as the selected copula family. With the kernel function, the rotation matrix, the simulated factor loadings, and the idiosyncratic component generated from a normal distribution, we simulate $X_{it}$. The dimension is set to 100 or 200 and the number of training data is 100, 250, 500 and 750, respectively. We set the number of test data to 200 for all these cases.

In order to compare the accuracy of the distributional predictions, we calculate the quantile scores for each model. The score function proposed by gneiting2011comparing, to assess the one-step prediction performance, which can be written as

equation[equation omitted — 88 chars of source]

where $y$ is the forecast value at $\alpha$-quantile and $x$ is the realized value. Suppose that the one-step predictive distribution of $X_{it}$ is $\widehat{G}_{i,t+1}$ in S-vine factor models. We compute the one-step-ahead forecast of the Value-at-Risk $\widehat{VaR}_\alpha^t = \widehat{G}^{-1}_{i,t+1}(\alpha|\mathbf{X}_{it})$ at $\alpha$-quantile, for $t = 1,....,n$, and calculate the mean of quantile scores mcneil2015quantitative,

equation[equation omitted — 107 chars of source]

In our simulation, we only present the results for $i=1$, as the results across dimensions are similar.

The predictions are repeated 100 times. We calculate the average value of the mean of quantile scores for the 100 repetitions in each case and demonstrate it in Table (ref). According to Table (ref), the SF models have smaller average quantile scores than DFM and GARCH, which is reasonable, since the data are simulated from the SF models. SF_1 has smaller average quantile scores compares to SF_2, which means the SF model with rotation matrix perform better in prediction. Therefore, the rotation is an effective step for SF models in the estimation and prediction process. The difference of the average quantile scores between the SF models with the "true" and "wrong" copulas is not pronounced, which implies that the type of copulas may not have significant influence on the SF models.

table[table omitted — 1,339 chars of source]

The box plots of the mean quantile scores across 100 simulations are shown in Figure (ref). We present results only for the case with $ d = 200 $. The SF_2 model is omitted from the plots due to its poor forecasting performance without rotation, as reported in Table (ref). The inadequacy of the incorrect SF model can be observed in the box plots, where the green box (representing SF_1) is slightly higher than the blue box (representing the true model). Similarly, both the DFM and GARCH models exhibit higher quantile scores and greater variability compared to the SF models.

figure[figure omitted — 524 chars of source]

Application

In this section, we apply the proposed method to estimate the dependence structure of the common factors driving the daily volatilities of the S&P 500 Index constituents, and to predict the Value-at-Risk (VaR) of the daily return of the S&P 500 Index. The return data are obtained from the London Stock Exchange Group database and cover the period from 1 January 2021 to 31 December 2024. We begin the sample in 2021 to avoid structural breaks in the S&P 500 Index constituents, which occurred on 6 March 2020 and were detected by LLF23. Missing values are filled by interpolation, and each stock’s returns are standardised to have zero mean and unit variance. Following barigozzi2016generalized, we adopt an approximate factor model for volatility, using the absolute values of returns as a proxy. To highlight the role of cross-sectional dependence and latent factors, we compare the predictive performance of our proposed factor copula model with that of a univariate benchmark. Specifically, we consider a GARCH(1,1) model with skewed Student-t innovations bollerslev1986generalized, which is a widely used and effective method for modelling financial return series.

The training set consists of data from the first three years (2021–2023), with the 2024 data reserved for one-step-ahead out-of-sample prediction. A total of 484 firms were listed on the S&P 500 Index during the entire period. Additionally, we include the S&P 500 Index itself as the 485th variable, with a particular focus on forecasting its Value-at-Risk (VaR). Specifically, we forecast the distribution of the absolute value of returns using the volatility factor model and assume equal probabilities for negative and positive returns to construct the VaR forecast.

Let $ X_{it} $ denote the absolute stock return of asset $ i $ at time $ t $, where $ i = 1, \ldots, 485 $ and $ t = 1, \ldots, 753 $ corresponds to the training period. The remaining 252 observations are reserved for the test period. We employ a factor model with a fifth-order M-vine structure to fit the absolute returns, where the vine order is selected based on the autocorrelation function (ACF) plots of the estimated factors. The number of factors is determined using the information criterion proposed by BN02, which identifies six latent factors. The vine structure is shown in Figure (ref).

figure[figure omitted — 286 chars of source]

We fit the model using our proposed method under different pair copula families. The candidate pair copulas include the Gaussian, Frank, Clayton and "all" copulas, where the pair copulas in different edges are selected according to the AIC values. Table (ref) shows the log-likelihood of the selected SF models. Using the estimated factor model with M-vine structure, we forecast the one-step-ahead distribution and calculate the Value-at-Risk (VaR) at 0.05, 0.10, 0.90, 0.95-quantile via an expanding window. The parameters and rotation matrix are fixed and without refit in the prediction process.

table[table omitted — 1,204 chars of source]

We calculate the violations of VaR of the one-step prediction in all these models and present it in Table (ref). The SF models with Frank and Clayton copulas have similar violations as GARCH process at 0.05 quantiles, while the SF with Gaussian and "all" copulas overestimate the violations. The GARCH process does not perform well at 0.95-quantile. On the contrary, the SF models with Clayton or "all" copulas show the advantages at 0.95-quantile. The violations at the lower and upper quantiles reveal the asymmetry of the left and right tail of the predictive distribution. For example, the SF models with Gaussian copulas overestimate the violations at lower quantile and underestimate the ones at upper quantile, which lead to understate of risk of loss and expectation of profit.The performance of SF models can be improved by selecting suitable combinations of pair copulas. The "all" copula model serves as a representative example, showing more precise violations across quantiles.

figure[figure omitted — 375 chars of source]

Figure (ref) illustrates the forecasted quantiles at levels 0.05, 0.10, 0.90, and 0.95 for the GARCH model and SF models with the "all", Gaussian, Frank, and Clayton copulas, based on test data. The VaR predicted by both the GARCH process and the SF models is influenced by extreme historical values. However, the GARCH model exhibits a smoother prediction, while the SF model displays a less persistent pattern and recovers to normal levels more quickly. Overall, the quantile plots on the test data show satisfactory performance.

Conclusion

This paper introduced a novel framework by integrating an approximate factor model with an S-vine copula structure to provide extra flexibility and enhanced the precision of dependency modelling in high-dimensional data. We proposed a two-step estimation procedure and developed the one-step-ahead distributional prediction using the S-vine factor models. A factor rotation matrix was applied to the PCA factor estimates in order to accommodate their dependence and align them with a semiparametric S-vine copula model. The proposed rotation and the parameters of the S-vine copula are estimated via maximum likelihood.

We establish the consistency of rotation and parameter estimates in S-vine copulas. To achieve this, we develop asymptotic theory for the empirical process of latent factors under mixing conditions, thereby extending existing results from i.i.d. sequences to dependent settings. Furthermore, we derive uniform convergence results for kernel-based entropy estimators, explicitly accounting for the combined effects of projection, factor estimation error, and serial dependence. These theoretical justifications provide strong support for the validity and robustness of our proposed estimation framework.

The numerical experiments embodied the fact that the estimated parameters converge to the true parameters with the increase in dimensionality and sample size. We examined the forecasting performance using VaR estimation via Monte-Carlo methods in simulation studies and empirical applications. The practical identification issue is discussed and illustrated via simulation in Section (ref), which sheds light on that the model can attain prominent and robust performance in one-step-ahead forecast process even when identification issues exists. We apply the daily returns of S&P 500 Index constituents as an example to assess the forecasting performance of our model. The approximate factor model with S-vine structure demonstrated a pronounced advantage in VaR prediction, where it outperformed the GARCH process at various quantiles.