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.
77,799 characters · 22 sections · 46 citation commands
Latent community paths in VAR-type models via dynamic directed spectral co-clustering
\affil[1]{Cornell University} \affil[2]{Sungkyunkwan University}
\baselineskip17pt
The analysis of high-dimensional vector autoregressive (VAR) models is often constrained by the curse of dimensionality. Even when regularization techniques yield sparse estimates basu2015regularized,brownlees2018realized,barigozzi2023fnets, direct inspection of individual coefficients---each representing predictive influence in the sense of granger1969investigating---rarely provides a coherent account of system-wide dependence. In many macroeconomic and financial applications, the main question is not simply which coefficients are nonzero, but how the system organizes itself into interpretable groups and how that organization changes across seasons or dependence horizons.
To address this, we study a dynamic network framework for latent community paths: ordered trajectories of directional groups that capture how clusters of “senders” and “receivers” of shocks split, merge, persist, and recompose over time. By shifting attention from isolated coefficients to the dynamic evolution of group structure, we turn a high-dimensional estimation problem into one of recovering interpretable trajectories of directional influence.
The core of our methodology embeds a degree-corrected stochastic co-blockmodel (ScBM; qin2013regularized) into the transition matrices of VAR-type systems. While recent literature has advanced block recovery in VAR models via spectral methods gudhmundsson2021detecting, brownlees2022community, these approaches largely rely on symmetrization, which inherently suppresses directional roles, or treat network snapshots as static entities. More recently, gudhmundsson2025detecting moved toward directional groups through an ScBM representation, yet primarily focused on lag aggregation. We overcome these limitations by integrating directed spectral co-clustering, which preserves the asymmetry between sending and receiving roles rohe2016co, wang2020spectral, with eigenvector smoothing across ordered stages in the spirit of liu2018global. This yields a dynamic community-recovery scheme tailored to track the continuous evolution of dependent multivariate time series.
We develop this framework in two distinct high-dimensional settings. First, in the periodic VAR (PVAR) model ursu2009modelling, baek2018periodic, we uncover how directed community structures evolve cyclically across seasons, moving beyond static partitions. Furthermore, in the generalized VHAR model, we replace the conventional daily--weekly--monthly aggregation corsi2009simple, baek2021sparse with arbitrary fixed lengths to track how global dependencies transition across ordered short-, medium-, and long-run horizons LeeBaek2023Crypto.
The primary theoretical contribution of this paper is establishing rigorous non-asymptotic misclassification bounds for directed community recovery within the ScBM--PVAR and generalized ScBM--VHAR frameworks. The approach is modular by design. We first derive estimation error bounds in operator norm for regularized lasso estimators in high-dimensional time series wong2020lasso. These bounds are then integrated with random-graph perturbation arguments for degree-corrected co-blockmodels qin2013regularized to control the distance between the estimated and population transition matrices. In a second step, we extend the eigenvector smoothing perturbation argument of Liu et al. (2018) to accommodate ordered chains of left and right singular-vector projectors across evolving seasons or horizons. The resulting bound on the Frobenius norm of row-normalized singular-vector perturbations directly translates into a misclassification rate, completing a clean theoretical pipeline from first-stage estimation error through spectral perturbation to clustering error. The modular structure also allows the co-clustering theory to be updated as tighter first-stage estimators become available. Finally, we provide verifiable sufficient conditions under which the core spectral assumptions hold, connecting the abstract conditions to primitive properties of the ScBM network parameters.
The empirical relevance of our framework is demonstrated in two applications. In the U.S. nonfarm payroll analysis (PVAR), the estimated paths reveal a recurrent business-centered core together with broader seasonal reallocation among more mobile sectors. In the global realized-volatility application (VHAR), the estimated paths reveal a compact U.S.-centered long-horizon block, a Europe-heavy developed core, and pronounced cross-horizon reorganization at shorter horizons. These results show that the proposed framework provides an interpretable account of how directional dependence is reorganized across seasons and horizons.
The rest of the paper is organized as follows. Section (ref) introduces the models and population networks. Section (ref) presents the estimators and the dynamic co-clustering procedure. Section (ref) establishes the theoretical misclassification bounds. Sections (ref) and (ref) report simulation results and empirical applications, respectively. Section (ref) concludes.
\noindentNotation. For a vector $v \in \mathbb{R}^d$, we denote its Euclidean norm by $\|v\| = \big(\sum_{i=1}^d v_i^2\big)^{1/2}$. For a matrix $M$, we write $\|M\|$ for the spectral norm, $\|M\|_F$ for the Frobenius norm, and $\|M\|_1$ and $\|M\|_{\infty}$ for the matrix $\ell_1$ and $\ell_{\infty}$ norms, that is, the maximum absolute column sum and maximum absolute row sum, respectively. Let $\rho(M)$ denote the spectral radius, write $M^\prime$ for the transpose of $M$, and use $\mathrm{Tr}(M)$ for its trace. We denote by $\mathbb{Z}$ the set of integers, by $\mathbb{N}$ the set of positive integers, and by $\mathbb{N}_0 := \mathbb{N}\cup\{0\}$ the set of nonnegative integers. Finally, big-$\mathcal{O}$ and big-$\Omega$ denote deterministic asymptotic upper and lower bounds, respectively, while $\mathcal{O}_{\mathbb{P}}$ denotes stochastic boundedness in probability. We write $a_n \asymp b_n$ if there exist constants $0 < c < C < \infty$ such that $c\, b_n \le a_n \le C\, b_n$ for all sufficiently large $n$.
We introduce three ScBM-based specifications for multivariate autoregressions. ScBM-VAR is the basic building block; ScBM-PVAR and generalized ScBM-VHAR extend it to seasonal and ordered-horizon settings.
Let $\mathcal{G}=(\mathcal{V},\mathcal{E},\mathcal{W})$ be a directed weighted graph with node set $\mathcal{V}=\{1,\ldots,q\}$, edge set $\mathcal{E}\subset \mathcal{V}\times\mathcal{V}$, and weights $\mathcal{W}=\{[w]_{ij}\in\mathbb{R}\setminus\{0\}:(i,j)\in\mathcal{E}\}$. The adjacency matrix $A \in \mathbb{R}^{q \times q}$ satisfies $[A]_{ij} = [w]_{ij}$ if $(i,j) \in \mathcal{E}$ and $0$ otherwise, where the nonzero weights satisfy $|[w]_{ij}| \le w_{\max}$ almost surely. We allow $K_y$ sending communities and $K_z$ receiving communities, with membership matrices $Y\in\{0,1\}^{q\times K_y}$ and $Z\in\{0,1\}^{q\times K_z}$. Under the ScBM, edge probability between nodes $i$ and $j$ satisfies
where $\Theta^y$ and $\Theta^z$ are diagonal propensity matrices and $B\in[0,1]^{K_y\times K_z}$ is the community link matrix. Let $\{\mathcal{V}_k\}$ be a partition of the node set. For an identifiability purpose,
The corresponding population adjacency matrix is
for some constant $\mu>0$.
Let $\{Y_t\}_{t=1}^T$ be a $q$-dimensional mean-zero time series. The VAR$(p)$ model is
At each lag $h$, the ScBM is embedded in the transition matrix as
where $\phi_h$ is a scalar to ensure the stability of the VAR($p$) model and
with $\tau_h := q^{-1}\sum_{i,j}[A_h']_{ij}$ denoting the average out-degree regularizer chaudhuri2012spectral,qin2013regularized, which is known to improve clustering performance under high heterogeneity in degrees. The transpose in (ref) ensures that columns represent lagged senders and rows current receivers. We assume that the latent network is common across lags $h=1,\ldots,p$, and that $\{\phi_h\}$ is chosen to satisfy
so that the VAR$(p)$ model is stable.
For each season $m = 1,\ldots,s$, let $G_m = (\mathcal{V}_m,\mathcal{E}_m,\mathcal{W}_m)$ be a directed weighted graph with adjacency matrix $A_m \in \mathbb{R}^{q \times q}$, where $[A_m]_{ij} = [w_m]_{ij}$ if $(i,j) \in \mathcal{E}_m$ and $0$ otherwise. We assume that there exists a finite constant $w_{\max,m} > 0$ such that $|[w_m]_{ij}| \le w_{\max,m}$ almost surely for all $(i,j) \in \mathcal{E}_m$. Let $Y_m\in\{0,1\}^{q\times K_{y_m}}$ and $Z_m\in\{0,1\}^{q\times K_{z_m}}$ denote the sending and receiving membership matrices, respectively. With propensity matrices $\Theta_m^y,\Theta_m^z$ and community link matrix $B_m$, the population adjacency matrix is
In addition to the setup of each season, we assume a cyclic evolution of communities across seasons, namely,
Note that “sending” and “receiving” denote directional connectivity “within” each season, following the convention of rohe2016co,wang2020spectral. The restriction (ref) implies that the receiving membership in one season is carried over as the sending membership in the next season. Due to this circular structure, the community evolution remains identical with the reversed flow; $Z_{m}=Y_{m-1}$, $m=2,\ldots,s$, and $Z_1=Y_s$.
The periodic VAR model is
where $p_m$ may depend on the season. Writing $t = m + ns$ yields the equivalent representation
The seasonal transition matrices are parameterized as
with a scalar $\phi_{h,m}$ and
Within each season, the latent network is taken to be common across lags, and $\phi_{h,m}$ is chosen so that the process is stable.
For both estimation and theory, it is convenient to stack one full seasonal cycle. Assume $T=Ns$ so that the sample contains $N$ complete cycles, and define
Then the ScBM-PVAR model can be written in the stacked form
where $p^*=\lceil (\max_{m=1,\ldots,s} p_m)/s \rceil$ and $\Psi_h^* \in \mathbb{R}^{qs\times qs}$ is obtained from the original seasonal coefficient matrices $\{\Phi_{h,m}\}_{m=1}^s$; see, for example, Section 2 of ursu2009modelling. Thus, the ScBM-PVAR model is represented as a VAR$(p^*)$ model on the stacked process $\{Y_n^*\}$. When $s=1$, this reduces to the ScBM-VAR model in (ref)--(ref).
Let $1<b_M<b_L$. The generalized VHAR model is
where
The effective sample size is $N=T-b_L$. When $(b_M,b_L)=(5,22)$, the generalized VHAR model reduces to the classical daily-weekly-monthly specification corsi2009simple. Equation (ref) is equivalent to a constrained VAR($b_L$). If $\Phi_1,\ldots,\Phi_{b_L}$ denote the implied VAR($b_L$) transition matrices, then
Hence, the number of free coefficients decreases from $b_Lq^2$ in an unrestricted VAR($b_L$) model to $3q^2$. The ScBM restriction is then constructed as, similar to (ref),
with a scalar $\phi_{(h)}$ and
For each $h\in\{S,M,L\}$, the population adjacency matrix $A_{(h)}$ is specified by an ScBM with horizon-specific parameters $(\mu_{(h)},\Theta^y_{(h)},\Theta^z_{(h)},Y_{(h)},Z_{(h)},B_{(h)})$ associated with the underlying graph $\mathcal{G}_{(h)}$, similar to (ref). We also assume a bounded-support condition: there exists a finite constant $w_{\max,(h)}>0$ such that $|[w_{(h)}]_{ij}|\le w_{\max,(h)}$ almost surely for all $(i,j)\in\mathcal{E}_{(h)}$. Unlike the cyclic PVAR case, these horizons are transient and are interpreted in the natural long-to-short order. We impose the restrictions $Y_{(L)}=Z_{(L)}=Y_{(M)}$ and $Z_{(M)}=Y_{(S)}$, so that the long-horizon community persists within the long-horizon block and is then inherited by the medium-horizon sending block, while the medium-horizon receiving block is inherited by the short-horizon sending block. This yields a hierarchical flow from broader horizons to shorter ones, where long-horizon structure constrains medium- and short-horizon interactions muller1997volatilities,corsi2009simple.
This section presents the estimation-clustering pipeline for the ScBM-VAR, ScBM-PVAR, and generalized ScBM-VHAR models, separating transition-matrix estimation from community recovery.
As the first step, we estimate the transition matrices via ordinary least squares (OLS) as the baseline, particularly in low-dimensional settings. In higher dimensions, we use lasso estimation.
\noindentScBM-VAR and ScBM-PVAR. For the ScBM-VAR model, either OLS or lasso estimation basu2015regularized is applied directly to the regression in (ref). For the ScBM-PVAR model, we use the stacked representation in (ref). Define
and set
The lasso estimator for the ScBM-PVAR model is
where $\lambda_{N,P}>0$ is a regularization parameter. In the low-dimensional case, the OLS estimator is obtained by removing the $\ell_1$ penalty and solving the corresponding least-squares problem.
\noindentGeneralized ScBM-VHAR. For the generalized ScBM--VHAR model, define $$ Z = (Y_{b_L+1}', \ldots, Y_T')' \in \mathbb{R}^{N\times q}, \qquad B =
\in \mathbb{R}^{3q\times q}, \qquad \beta_V^* = \operatorname{vec}(B), $$ with $N=T-b_L$. Let $X \in \mathbb{R}^{N\times qb_L}$ be the unrestricted VAR$(b_L)$ lag-design matrix, and set $$ X_e = X R_{b_M,b_L} \in \mathbb{R}^{N\times 3q}, $$ where $R_{b_M,b_L} \in \mathbb{R}^{qb_L\times 3q}$ is the deterministic restriction matrix implied by \eqref{e:VHAR_coefficients}. Then $$ Z = X_e B + E, \qquad \operatorname{vec}(Z) = (I_q \otimes X_e)\beta_V^* + \operatorname{vec}(E). $$ The restricted lasso estimator is
where $\lambda_{N,V}>0$ is a regularization parameter. The restricted OLS estimator is obtained by dropping the $\ell_1$ penalty.
The co-clustering step has the same backbone in all models. Given estimated autoregressive matrices $\hat\Phi_m$, compute their leading left and right singular vectors, row-normalize them, optionally smooth the associated projection matrices across ordered seasons or horizons, and apply $k$-means to the matched left-right pairs.
\noindentScBM-VAR. We begin with the ScBM-VAR model, combining the directed co-clustering idea of rohe2016co with the VAR Blockbuster approach of gudhmundsson2021detecting. Let
and compute the leading $K$ left and right singular vectors $\hat X_L,\hat X_R\in\mathbb{R}^{q\times K}$. We then row-normalize them as
If $K_y=K_z=K$, we run $k$-means on $(\hat X_L^*\ \hat X_R^*)\in\mathbb{R}^{q\times 2K}$ and use the resulting labels for both sending and receiving communities. Otherwise, we cluster $\hat X_L^*$ and $\hat X_R^*$ separately using $K_y$ and $K_z$ clusters, respectively.
\noindentScBM-PVAR. We now extend the procedure to the ScBM-PVAR model and incorporate smoothing across seasons via PisCES liu2018global. For simplicity, we take $p_m = p$ for all $m=1,\ldots,s$. For each season $m$, form
and compute its top $K_{y_m}$ left and $K_{z_m}$ right singular vectors $\hat X_{m,L}$ and $\hat X_{m,R}$. We impose the cyclic rank constraints $K_{z_{m-1}}=K_{y_m}$ for $m=2,\ldots,s$ and $K_{z_s}=K_{y_1}$. Define projector matrices
and initialize $\bar U_{m,L}^{(0)}=\hat U_{m,L}$ and $\bar U_{m,R}^{(0)}=\hat U_{m,R}$. For a smoothing parameter $\alpha_N>0$, PisCES updates the left sequence by
and the right sequence analogously:
where $\Pi(M;K)=\sum_{k=1}^K \nu_k\nu_k'$ is the projector onto the span of the top-$K$ eigenvectors of $M$. After convergence, write $\bar U_{m,L}=\bar X_{m,L}\bar X_{m,L}'$ and $\bar U_{m,R}=\bar X_{m,R}\bar X_{m,R}'$, row-normalize $\bar X_{m,L}$ and $\bar X_{m,R}$ as in (ref), and link adjacent seasons by $k$-means on $(\bar X_{m-1,R}^*\ \bar X_{m,L}^*)$ for $m=2,\ldots,s$ and on $(\bar X_{s,R}^*\ \bar X_{1,L}^*)$ for the cyclic pair. When $\alpha_N=0$, this reduces to season-wise spectral co-clustering.
\noindentScBM-VHAR. For the generalized ScBM-VHAR model, we follow the same procedure as for ScBM-PVAR, except that the three “seasons” are now replaced by short-, medium-, and long-horizon components with aggregation lengths $1$, $b_M$, and $b_L$. Define
We then apply the same singular-vector extraction and PisCES smoothing routine with $s=3$, imposing the rank constraints $K_{y_{(L)}}=K_{z_{(L)}}=K_{y_{(M)}}$ and $K_{z_{(M)}}=K_{y_{(S)}}$. Under the long-to-short convention, the linked stages are
and terminal $Z_{(S)}$. After smoothing, we jointly cluster $(\bar X_{(L),R}^*\ \bar X_{(L),L}^*\ \bar X_{(M),R}^*)$ to recover the first stage, jointly cluster $(\bar X_{(M),L}^*\ \bar X_{(S),R}^*)$ to recover the second stage, and cluster $\bar X_{(S),L}^*$ separately. When $\alpha_N=0$, the procedure reduces to horizon-wise spectral co-clustering.
We select the PisCES smoothing parameter $\alpha_N$ by cross-validation, adapting the procedure of liu2018global. Specifically, we apply the method to either seasonal or horizon-specific autoregressive matrices $\{\hat\Phi_m\}_{m=1}^s$ as follows:
This criterion may be interpreted as a quadratic approximation to the von Neumann entropy; see, for example, equations (10)–(11) in ye2014approximate.
This section develops the theory for the ScBM-PVAR and ScBM-VHAR procedures. The key idea is modularity: once an operator-norm bound for $\hat{\Phi}-\Phi$ is available, the downstream co-clustering theory follows by substitution. All proofs are deferred to Appendix (ref).
To study the stability of the ScBM-PVAR model, we work with the stacked representation in (ref). Under this representation, the ScBM-PVAR model is a VAR$(p^*)$ model on the stacked process $\{Y_n^*\}$, where the coefficient matrices $\Psi_h^*$ are inherited from the original seasonal system. This allows us to analyze stability through the corresponding stacked companion form.
Condition (i) is the standard VAR stability condition, and condition (ii) is a norm bound that controls the overall scale of the stacked transition matrices, in line with assumptions used in the related models gudhmundsson2021detecting,yin2023general.
The generalized ScBM-VHAR model in (ref) can be written as a VAR($b_L$) with constrained coefficients as in (ref). Hence, by treating $\{\Phi_h\}_{h=1}^{b_L}$ as the VAR($b_L$) transition matrices, stability is ensured by the usual characteristic-root condition
Since $b_L$ is fixed, we additionally impose a deterministic norm envelope on the constrained lag matrices; this will be used later in the estimation bounds.
Condition (i) is the standard stability condition for a VAR($b_L$) companion system. Condition (ii) is a convenient envelope used below to control the norms of the constrained lag matrices uniformly.
In this section, we derive the convergence rates for both models. The ScBM-PVAR and generalized ScBM-VHAR models can be analyzed analogously by stacking their respective seasons or horizons. For the ScBM-PVAR model, we define the stacked autoregressive matrix as
Let $\hat{\Phi}$ and $\tilde{\Phi}$ denote the sample estimator and its population counterpart, respectively, where $\tilde{\Phi}$ is the expectation taken over the network randomness. The total error then decomposes as
This section focuses on bounding the estimation error, $\|\hat{\Phi} - \Phi\|$, while the network randomness term in (ref) is addressed in the subsequent section.
First, consider when the OLS estimator $\hat{\Phi}^{\mathrm{ols}}$ is used. The idea is to view $\{Y_n^*\}$ as a weakly dependent VAR process. A similar approach is also used in the standard PVAR process ursu2009modelling,boubacar2023estimating. Let $\Sigma_Y^* := {\rm Var}(Y_n^*)$ and $\Sigma_\varepsilon^* := {\rm Var}(\varepsilon_n^*)$. Define the normalized process $(\Sigma_Y^*)^{-1/2} Y_n^*$ and its strong mixing coefficients $\alpha_{Y^*}(\ell)$ in the usual way.
For the generalized ScBM-VHAR model, we use the restricted VAR($b_L$) representation induced by the horizon lengths $1<b_M<b_L$, which yields the following rate.
Next we replace OLS by lasso estimation $\hat{\Phi}^{\mathrm{lasso}}$. For the ScBM-PVAR model, this is the stacked VAR($p^*$) regression, whereas for the generalized ScBM-VHAR model, it is the restricted VAR($b_L$) regression induced by $R_{b_M,b_L}$. In both cases, the downstream co-clustering theory depends on the first stage only through an operator-norm bound for $\hat{\Phi}-\Phi$. For the ScBM-PVAR model, let $\alpha_P^*$ denote the stacked coefficient vector; since $p^*$ is fixed, we suppress the distinction between $N$ and $N-p^*$ in the rates.
These are precisely the deterministic inputs required by Proposition 4.1 of basu2015regularized. Note that the analysis in basu2015regularized is conducted under Gaussian assumptions, whereas Assumption (ref) can be interpreted as imposing conditions on sub-Gaussian random vectors. However, subsequent studies (e.g., wong2020lasso,basu2024high) show that the suggested convergence rates remain the same without additional conditions. Hence, the theoretical guarantees derived under Gaussian assumptions extend. When the stacked ScBM-PVAR is Gaussian and stable, Proposition 4.2 and Proposition 4.3 of basu2015regularized verify Assumption (ref)(ii)--(iii) with process dimension $p=qs$ and lag order $d=p^*$. Since $p^*$ is fixed, the rate simplifies to $\lambda_{N,P}\asymp\sqrt{\log(qs)/N}$.
For generalized ScBM-VHAR model, let $\beta_V^*$ denote the restricted coefficient vector in (ref).
For the fixed-horizon VHAR design, Proposition 1 of baek2021sparse verifies the analogue of the Gram-matrix construction in their equation (2.10) and the RE and DB in their equation (2.11) by embedding VHAR into a stable high-order VAR model. The same argument carries over to arbitrarily fixed $(b_M,b_L)$ after replacing the specific aggregation matrix in their equation (2.12) by $R_{b_M,b_L}$. Because $b_M$ and $b_L$ are fixed, this changes only constants, so Proposition 4.1 of basu2015regularized yields the same rate order.
As a continuation of Section (ref), we now bound $\|\tilde{\Phi} - \Phi\|$, i.e., the deviation arising from random graphs under the ScBM assumptions.
This assumption corresponds to the sparse assumption of networks due to the regularization with expected degrees of order at least $\log(sq)$, which allows exact recovery in ScBM-type models rohe2016co,abbe2018community.
We now translate the operator-norm bounds from Sections (ref) and (ref) into perturbation bounds for the row-normalized singular vectors and the resulting misclassification rates. For season $m=1,\ldots,s$, let $\hat X^*_{m,L}$ and $\hat X^*_{m,R}$ denote the estimated row-normalized left and right singular vectors, let $X^*_{m,L}$ and $X^*_{m,R}$ denote their population counterparts, and let $\bar X^*_{m,L}$ and $\bar X^*_{m,R}$ denote the corresponding PisCES-smoothed versions. Likewise, for horizon $h\in\{S,M,L\}$, let $\hat X^*_{(h),L}$ and $\hat X^*_{(h),R}$ denote the estimated row-normalized singular vectors, let $X^*_{(h),L}$ and $X^*_{(h),R}$ denote their population counterparts, and let $\bar X^*_{(h),L}$ and $\bar X^*_{(h),R}$ denote the corresponding PisCES-smoothed versions. Since singular vectors are identified only up to orthogonal transformations, all bounds below are understood after suitable rotations.
The following conditions are often implicit in prior analyses of regularized spectral clustering and directed co-clustering; cf. qin2013regularized,rohe2016co,gudhmundsson2021detecting. In the present asymmetric co-clustering framework, we state them explicitly because stable recovery of both the sending and receiving singular spaces is needed for the singular-vector perturbation and misclassification arguments below.
We report finite-sample results for the proposed procedures using sparsely generated ScBM-PVAR and generalized ScBM-VHAR models. In each setup, the networks embedded in the transition matrices are generated once and held fixed across replications, ensuring that the innovations are the sole source of randomness. Furthermore, we do not impose degree correction in this simulation study; this allows us to isolate the effects of latent path evolution from the confounding effects of node-specific degree heterogeneity. We present the results based on the PisCES-smoothed lasso estimator, while the corresponding OLS estimation and adjusted Rand index vinh2009information comparisons are deferred to Appendix (ref).
For each season $m=1,\ldots,4$, we first generate a binary support matrix $A_m$. If node $i$ belongs to sending block $k$ and node $j$ belongs to receiving block $r$, then for $i\neq j$, \[ [A_{m}]_{ij}\sim \mathrm{Bernoulli}(\pi_{kr,m}), \qquad \pi_{kr,m}=\min\!\left\{\frac{\kappa_{kr}}{n^{z}_{m,r}-\mathbf{1}(k=r)},\,0.95\right\}, \] where $n^{z}_{m,r}$ is the size of receiving block $r$ in season $m$. Hence the expected off-diagonal support per row is $O(1)$, and the matrix density is $O(q^{-1})$. Conditional on $A_m$, the raw seasonal coefficient matrix $\widetilde\Phi_m$ is generated entrywise. The own-lag diagonal is always retained, \[ [\widetilde\Phi_{m}]_{ii}=a_{\mathrm{self}}u_{ii}, \qquad u_{ii}\sim \mathrm{Unif}(0.95,1.05), \] and for $i\neq j$ with $k=y_{m,i}$ and $r=z_{m,j}$, \[ [\widetilde\Phi_{m}]_{ij}=
\qquad u_{ij}\sim \mathrm{Unif}(0.90,1.10). \] Thus $A_m$ determines the support, while the block pair determines the coefficient magnitude. For Type 1 we use $(a_{\mathrm{self}},a_{\mathrm{diag}},a_{\mathrm{upper}},a_{\mathrm{lower}})=(0.30,0.14,0.04,0.06)$, $(\kappa_{\mathrm{diag}},\kappa_{\mathrm{upper}},\kappa_{\mathrm{lower}})=(2.6,0.6,0.9)$, and for Type 2 we use $(a_{\mathrm{self}},a_{\mathrm{diag}},a_{\mathrm{upper}},a_{\mathrm{lower}})=(0.28,0.12,0.05,0.07)$, $(\kappa_{\mathrm{diag}},\kappa_{\mathrm{upper}},\kappa_{\mathrm{lower}})=(2.2,0.7,1.3)$. Type 2 is therefore intentionally more difficult, because it allows relatively stronger between-block connectivity. After constructing $\widetilde\Phi_1,\ldots,\widetilde\Phi_4$, we apply a common rescaling so that $\sigma_{\max}(\Phi_4\Phi_3\Phi_2\Phi_1)=0.90$. This preserves the relative seasonal pattern while ensuring a stable but non-trivial signal.
Figure (ref) displays the three path designs. Path 1 is static, with four communities throughout. Path 1 also serves as a static benchmark: when the block structure does not change over seasons, the proposed co-spectral clustering may be viewed as a directed extension of the symmetrized spectral clustering of gudhmundsson2021detecting. Including this case allows us to verify that the method works well even in the non-dynamic baseline before turning to split--merge and refinement paths. Path 2 follows a split--merge pattern $2\to3\to3\to2$. Path 3 follows a nested refinement pattern $2\to2\to3\to4$. We consider $q\in\{18,36,60\}$ and $T\in\{200,500,1000,2000\}$. For each combination of dimension, path, and type, we generate one seasonal coefficient system and run $200$ replications with Gaussian innovations, diagonal innovation variance $0.5$, and burn-in length $500$. Performance is summarized by the mean spectral-norm error, mean accuracy, and mean ARI, where the four season-specific clustering results are averaged within each replication.
We compare OLS and lasso estimation results. For lasso estimation, we use the accelerated proximal gradient method beck2009fast with a common penalty multiplier selected by block cross-validation baek2021sparse. For each design, one pilot series is simulated. Let $N_{\mathrm{eff},m}$ be the usable sample size in season $m$. The baseline penalty is $\lambda_m^{\mathrm{base}}=\{\log(sq^2)/N_{\mathrm{eff},m}\}^{1/2}$ and we search over $c_\lambda\in\{0.10,0.15,0.20,\ldots,1.00\}$. For each candidate $c_\lambda$, we fit the seasonal lasso regressions with $\lambda_m=c_\lambda\lambda_m^{\mathrm{base}}$ and evaluate them by $10$-fold block cross-validation, summing the validation errors over the four seasons. The selected $c_\lambda$ is then fixed within the same setup. The PisCES smoothing parameter $\alpha_N$ is selected separately by a $5$-fold holdout criterion over the default grid in Section (ref).
Table (ref) shows the expected monotone effect of sample size: the spectral-norm error decreases and the accuracy increases in every setup. For instance, in $(q=18,\mathrm{path1},\mathrm{type1})$, the error drops and the accuracy rises as $T$ increases. Even in the most difficult setting $(q=60,\mathrm{path3},\mathrm{type2})$, the error still decreases and the accuracy increases. The relative difficulty ordering is stable throughout: larger $q$ is harder, Type 2 is harder than Type 1, and Path 3 is typically the most demanding, while Path 1 is the easiest. These patterns are consistent with the theory, because the monotone decline of the first-stage error with $T$ is accompanied by monotone improvement in community recovery.
The same ordering persists for the OLS estimation results shown in Appendix (ref). However, the rate of improvement is significantly slower. In particular, Figure (ref) shows that the ARI increases steadily for lasso results across all three paths, whereas the OLS curves remain substantially lower, especially for $q=36$ and $q=60$ and for Paths 2 and 3. Thus the gain from lasso estimation is not minor; it is what keeps spectral co-clustering informative once the problem becomes high-dimensional or the community path becomes more complex.
We keep the three horizon-specific coefficient matrices fixed within each setup. If node $i$ belongs to sending block $k$ and node $j$ belongs to receiving block $r$ at horizon $h\in\{S,M,L\}$, we generate
where $B_{(h)}$ is a horizon-specific block-probability matrix. We consider the following two specifications:
These correspond to Type 1 and Type 2, respectively, with all probabilities truncated to $[0.001,0.995]$. Type 1 serves as the benchmark design, with stronger own- and within-block signal and weaker between-block interactions. Type 2 is intentionally more difficult: it reduces the own- and within-block magnitudes while increasing the cross-block magnitudes, thereby weakening community separation. Specifically, Type 1 uses $(p_{\mathrm{diag},S},p_{\mathrm{off},S})=(0.95,0.02)$, $(p_{\mathrm{diag},M},p_{\mathrm{off},M})=(0.93,0.03)$, and $(p_{\mathrm{diag},L},p_{\mathrm{off},L})=(0.91,0.04)$; the corresponding Type 2 probabilities are $(0.92,0.90,0.88)$ for diagonal entries, $(0.05,0.06,0.07)$ for upper off-diagonal entries, and $(0.08,0.09,0.10)$ for lower off-diagonal entries across the short-, medium-, and long-horizon components.
Conditional on $A_{(h)}$, we normalize by sender-block size and define $[\widetilde{A}_{(h)}]_{ij}=[A_{(h)}]_{ij}/n_{h,k}^{y}$, where $n_{h,k}^{y}$ is the size of sending block $k$ at horizon $h$. This prevents larger sending blocks from generating stronger coefficients purely because they contain more nodes, and therefore keeps the block effect interpretable as an average per-sender effect. We then set $\widetilde{\Phi}_{(h)} = c_h\,\widetilde{A}_{(h)}$, $h\in\{S,M,L\}$ with $(c_S,c_M,c_L) = (0.34,\;0.28\sqrt{b_M},\;0.24\sqrt{b_L})$, and apply a common rescaling factor $a>0$ such that $\sigma_{\max}\bigl(a(\widetilde \Phi_{(S)}+\widetilde \Phi_{(M)}+\widetilde \Phi_{(L)})\bigr)=0.90$. A scale has been chosen so that the mean stationary marginal variance equals $0.5$.
Figure (ref) displays the three latent paths in the natural long-to-short order $\mathrm{LS}\rightarrow \mathrm{LR/MS}\rightarrow \mathrm{MR/SS}\rightarrow \mathrm{SR}$, where $\mathrm{LS}=Y_{(L)}$, $\mathrm{LR/MS}=Z_{(L)}=Y_{(M)}$, $\mathrm{MR/SS}=Z_{(M)}=Y_{(S)}$, and $\mathrm{SR}=Z_{(S)}$. Path 1 is static. Path 2 follows the split--merge pattern $2 \to 2 \to 3 \to 2$, whereas Path 3 follows the coarsening pattern $3 \to 3 \to 2 \to 2$.
We consider $q\in\{18,24,36\}$ and $T\in\{500,1000,2000,3000\}$. Each setup is run for 100 replications with burn-in length 300, $(b_M,b_L)=(3,10)$ and Gaussian innovations with diagonal variance of 0.5. We use both OLS and lasso estimation, but report only the PisCES-smoothed lasso results in the main text. The lasso penalty is chosen by block cross-validation over $c_\lambda\in\{0.10,0.15,0.20,\ldots,1.00\}$ while the PisCES smoothing parameter is selected separately based on Section (ref).
Table (ref) again shows monotone improvement as $T$ increases, although the VHAR design is clearly more challenging than the PVAR experiments. For example, in $(q=18,\mathrm{path1},\mathrm{type1})$, the spectral norm decreases and the accuracy increases relatively slowly while $T$ grows more rapidly. In the more difficult setting $(q=36,\mathrm{path2},\mathrm{type2})$, the spectral norm still decreases and the accuracy still increases, although more slowly. The same difficulty ordering is stable throughout the table: larger $q$ leads to poorer recovery, Type 2 is more difficult than Type 1, and the non-static paths are more difficult than Path 1. Path 1 is the easiest case because the block structure is unchanged across horizons, whereas Paths 2 and 3 require the method to track either a split--merge or a coarsening pattern across ordered horizons. This is also consistent with its role as a static benchmark, where the procedure effectively reduces to a directed analogue of the static spectral clustering setting and therefore should perform best.
The OLS results displayed in Appendix (ref) exhibit the same ordering, but with substantially larger spectral-norm errors and markedly lower ARIs. For instance, in $(q=36,\mathrm{path1}$, $\mathrm{type1},T=3000)$, the PisCES-smoothed lasso estimator achieves an ARI of $0.971$, whereas the corresponding OLS estimator achieves only $0.180$. Thus, sparse estimation is markedly more effective in the ScBM-VHAR models as well.
The nonfarm payroll employment series records paid U.S.\ workers in selected industries, excluding farmworkers, private household employees, and military personnel. We use the monthly U.S.\ Bureau of Labor Statistics series obtained from the Federal Reserve Economic Data (FRED).\footnote{Series descriptions and industry classifications are available at \url{https://fred.stlouisfed.org/} and \url{https://www.bls.gov/news.release/empsit.t17.htm}.} We aggregate the data to quarterly frequency from January 1990 to March 2020, which yields $T=120$ observations per series. To form a moderately high-dimensional and heterogeneous panel, we include major sectors together with selected third- and fourth-level subcategories, for a total of $q=22$ series; see Table (ref) in Appendix (ref). Throughout, we refer to sectors by the codes reported in that table.
Before estimation, we apply a log transformation and first differencing. Figure (ref) in Appendix (ref) shows time plots for the first eight series and sample ACFs and PACFs for Mining, Nondurable, Wholesale, and Retail. The ACFs display clear periodicity, while the higher-order PACFs are relatively sparse. This pattern supports a low-order periodic VAR specification. Cross-correlations, not reported here, also indicate substantial intersectoral dependence. At the same time, sectors within the same broad category do not necessarily move together. For example, the 2008 Subprime Mortgage crisis appears in most manufacturing- and trade-related series, whereas Utilities follows a distinct path and Manufacturing behaves differently from Mining and Construction. These features favor a community-based analysis over a grouping based only on broad sector labels.
We fit a lasso-based ScBM--PVAR model with $s=4$ seasons corresponding to Q1--Q4 and a common lag order $p_m=1$ for all $m=1,\ldots,4$. For each quarter, we estimate the transition matrix, construct the corresponding autoregressive matrix $\hat{\Phi}_m$, and inspect its singular values; see Figure (ref) in Appendix (ref). To enforce the cyclic community structure in Section (ref), we adopt the admissible configuration $$ \left(K_{y_1},K_{z_1},K_{y_2},K_{z_2},K_{y_3},K_{z_3},K_{y_4},K_{z_4}\right) = (2,3,3,3,3,2,2,2), $$ which yields a sending-side seasonal path $2 \to 3 \to 3 \to 2$ across Q1--Q4.
Figure (ref) reveals a clear cyclic seasonal topology. Q1 is organized around a broad two-community split between a business--trade--property block and a production--infrastructure--public-service block. In Q2, this coarse partition expands into three groups, separating a consumer-, property-, and public-demand block, a business-coordination block, and an infrastructure--information--health block. Q3 retains a three-way structure but with substantial recomposition of memberships, now distinguishing a local-demand and public block, a business--trade--information block, and a production--infrastructure--human-capital block. In Q4, the system coarsens back to a two-community structure. The annual dependence pattern is therefore characterized by mid-year differentiation followed by end-of-year recomposition.
A few sectors remain highly stable, most notably Accommodation and Arts, whereas Wholesale, Management, Transportation, and several local-demand and public-service sectors account for much of the seasonal reallocation. Economically, this suggests a recurrent business-centered core together with broader seasonal reshuffling in more mobile sectors. For comparison, Figure (ref) in Appendix reports the alternative $2 \to 2 \to 3 \to 2$ specification, which yields a coarser Q2 split but a broadly similar annual pattern.
Realized volatility (RV) is an ex-post measure of return variation, typically constructed from the sum of squared intraday returns. We compute RV from 5-minute returns and follow the cleaning and aggregation procedure in Section 4.3 of corsi2009simple; see also andersen2003modeling. Our panel contains $q=29$ stock-index RV series for major equity markets. Since the original Oxford-Man Institute data are no longer publicly available, we use the processed dataset of baek2021sparse. The sample spans January 3, 2010 to December 31, 2019, giving $T=2618$ observations and ending before the COVID-19 pandemic. The stock indices, MSCI classifications msci2025, and regional information are listed in Table (ref) in Appendix (ref).
Figure (ref) in Appendix (ref) reports time plots and sample ACF/PACF functions for four representative indices: FTSE 100, Nikkei 225, KOSPI Composite, and IPC Mexico. All four series exhibit persistent autocorrelation, consistent with long-memory-type volatility dynamics, but their volatility spikes differ markedly across markets. For instance, FTSE 100, KOSPI, and IPC Mexico show sharp increases around the third quarter of 2011, whereas Nikkei 225 peaks earlier and more gradually. KOSPI also appears closer to FTSE 100 than to Nikkei 225, which suggests that RV comovement is not driven solely by geographic proximity and motivates a network-based analysis of latent community structure.
We estimate the generalized ScBM--VHAR model under the standard specification $(b_M,b_L)=(5,22)$, so that the short-, medium-, and long-horizon components correspond to daily, weekly, and monthly RV aggregates. Scree plots of the estimated transition matrices, reported in Figure (ref) in Appendix (ref), suggest three dominant long-horizon components, a somewhat richer middle-horizon structure when each horizon is treated separately, and three short-horizon components. Under the imposed long-to-short restriction (ref), however, the ranks cannot be selected independently across horizons. We therefore adopt the admissible configuration
which yields three effective communities at each horizon and gives a parsimonious baseline representation of the cross-horizon structure. In what follows, the three VHAR components are interpreted as long-, medium-, and short-horizon volatility spillovers. We note that an alternative admissible specification producing a $(3,3,4,2)$ horizon pattern leads to a broadly similar picture, although it yields a somewhat more dynamic short-horizon reallocation.
Figure (ref) should be read from left to right, from the long horizon to the short horizon. It reveals a clear cross-horizon reorganization of the realized-volatility network. At an aggregate level, 17 of the 29 indices move across two distinct communities, 10 remain on a single aligned path, and only KS11 and SSEC pass through all three communities. The dependence structure is therefore not well captured by a single static partition. It is more naturally described as a sequence of horizon-specific communities with substantial but structured reallocation.
At the long horizon, one compact block consists of DJI, IXIC, SPX, and N225, forming a distinct U.S.-centered cluster with Nikkei 225 attached to it. A second and much larger block is Europe-heavy, though it also contains several emerging and peripheral markets. The remaining block is smaller and includes AORD, HSI, KS11, KSE, OMXHPI, OSEAX, and RUT. Thus the long-horizon structure is not a simple regional partition, but separates a compact U.S.-centered block, a broad Europe-heavy core, and a smaller Asia--Pacific and peripheral block.
The middle horizon yields the clearest segmentation. DJI, IXIC, SPX, and N225 remain together and are joined by KS11, KSE, and SSEC, forming a distinct U.S.--East Asia block. At the same time, AEX, BFX, FCHI, FTMIB, GDAXI, IBEX, OMXHPI, OMXSPI, SMSI, and STOXX50E form a more cohesive Europe-heavy developed core, with BVSP attached to that group. The remaining indices form a broader mixed peripheral block. The weekly scale is therefore where the global volatility network is most differentiated.
At the short horizon, the structure becomes more dynamic. The long- and middle-horizon U.S.-centered block no longer moves as a single unit: DJI joins a broad developed-market cluster, IXIC and SPX remain together in a separate short-run block, and N225 shifts into a broader peripheral group. A relatively stable developed-market trajectory nevertheless remains visible. In particular, the path $3 \to 2 \to 3$ is shared by AEX, BFX, FCHI, FTMIB, IBEX, SMSI, and BVSP, whereas KS11 and SSEC act more like bridge markets, following the paths $2 \to 1 \to 3$ and $3 \to 1 \to 3$, respectively. Figure (ref) reports the 3--4--3 specification for comparison. The broad long-horizon partition remains similar, but the middle horizon becomes more finely segmented and the short-horizon reallocation becomes more dynamic.
The lasso-based ScBM--VHAR results therefore point to a layered and highly dynamic dependence structure. A compact U.S.-centered block and a Europe-heavy developed core emerge at the long horizon. The middle horizon sharpens this separation, while the short horizon produces the strongest internal reorganization. This long-to-short progression suggests that near-term volatility spillovers are shaped less by simple geographic proximity than by heterogeneous market roles and transmission channels.
This paper establishes a framework for latent community paths in high-dimensional VAR-type models. By combining degree-corrected stochastic co-blockmodels, directed spectral co-clustering, and eigenvector smoothing, the proposed method tracks how directional groups persist, split, merge, and recompose across seasons or dependence horizons. The framework covers both ScBM--PVAR and generalized ScBM--VHAR models. Its theoretical contribution is to provide non-asymptotic guarantees for singular-vector perturbation and community misclassification through a modular link from first-stage estimation error to clustering error.
The empirical applications illustrate the value of this perspective. In U.S. nonfarm payrolls, the estimated paths distinguish a recurrent business-centered core from more mobile sectors with stronger seasonal reallocation. In global realized volatilities, the estimated paths reveal a compact U.S.-centered long-horizon block, a Europe-heavy developed core, and substantial short-horizon reorganization. These findings suggest that latent community paths provide an interpretable summary of dynamic dependence that is difficult to obtain from entrywise coefficient inspection alone.
CB was supported by the National Research Foundation of Korea grant funded by the Korean government (MSIT) (RS-2025-00519717). The authors thank Drs. Francis X. Diebold, Majid Al-Sadoon, and Vladas Pipiras for their comments at the 2025 NBER–NSF Time Series Conference, which substantially improved the quality of this paper. The authors are also grateful to Drs. Gu{\dh}mundur Stef{'a}n Gu{\dh}mundsson and Christian Brownlees for insightful discussions regarding the models. The R code for the simulation study and data applications is available at \url{https://github.com/crbaek/dynamic-network-clustering}.
{ {4pt} }