Testing Clustered Equal Predictive Ability with Unknown Clusters

Testing Clustered Equal Predictive Ability with Unknown Clusters

\def\spacingset#1{ {#1}} \spacingset{1}

\if11 \fi

\if01 {

} \fi

{\bf Keywords:} Forecast Evaluation; Hypothesis Testing; Panel Data; Sample Splitting; Selective Inference. \newline {\bf JEL classification}: C12, C23, C52, C53, C55.

\spacingset{1.8}

{5pt} {5pt}

Introduction

Despite a rich literature on testing equal predictive ability (EPA) in time series---see clark13 and rossi_forecasting_2021 for reviews---EPA testing in panel settings has only recently attracted attention. The main contributions are Akgun, Pirotte, Urga & Yang (akgun24, APUY) and Qu, Timmermann & Zhu (qu23b, QTZ), who study two null hypotheses: overall EPA (O-EPA), which states forecast equivalence on average across time and units, and clustered EPA (C-EPA), which states equivalence across $K \geq 2$ known clusters.

In many applied forecasting contexts, predictive performance varies across units such as countries or firms. For instance, dreher2008political show that IMF forecast quality differs significantly depending on whether countries received IMF assistance or were aligned with major donors in international platforms. More generally, forecasting accuracy may vary systematically across groups defined by income level, geography, political alignment, or development status. This implies heterogeneity across clusters, often unobserved by the researcher. Testing for EPA in such cases must account for clustered heterogeneity without prior knowledge of cluster structure.

The primary contribution of this paper is the development of conditional C-EPA tests for panel data with unknown cluster structure. Our framework extends APUY and QTZ in several directions. First, inspired by giacomini06, we allow for conditioning variables, offering a more flexible setup. Second, we estimate clusters using the Panel Kmeans Algorithm, which generalizes classical Kmeans by exploiting time variation. Third, to ensure valid post-clustering inference, we develop a selective conditional inference framework based on the polyhedral method lee_exact_2016. We propose a Wald-type test for pairwise homogeneity of cluster centers and derive its truncated $\chi$-variate asymptotic distribution conditional on the estimated clusters, along with an analytical characterization of the truncation region under Panel Kmeans. Fourth, we prove that information criterion (IC)-based selection of $K$ preserve validity without additional conditioning. Finally, rather than using a Wald test for joint C-EPA—which may be anti-conservative when many constraints are tested—we aggregate the evidence from all pairwise tests and the O-EPA test using a $p$-value combination approach that controls Type I error.

The main theoretical challenge lies in valid inference on cluster centers after estimating the unknown clusters. While methods such as hierarchical clustering and Kmeans are common, we focus on the Panel Kmeans Estimator, widely used in econometrics bonhomme15a,bonhomme2022discretizing,patton23. When predictive ability differences vary across but not within clusters, Panel Kmeans consistently recovers cluster structure under cluster separation. Under the C-EPA null, this assumption fails and all units belong to a single cluster, giving rise to the double dipping problem kriegeskorte09, where the same data are used for both clustering and inference. A common remedy is sample splitting: in cross-sections, gao24 show it does not yield valid inference, while in panels, patton23 propose a Split Sample test exploiting the time dimension. However, although sample splitting remains a natural way to deal with double dipping, the past literature highlighted some limitations of this approach: splits are often arbitrary hansen2012choice, structural breaks can invalidate the design, and dependence may compromise validity kuchibhotla22. patton23 offer an effective solution to the latter by discarding some periods between training and test sets. This reduces dependence but may also reduce power.

We propose an alternative selective inference framework that uses the full sample to estimate unknown clusters and conduct inference on their centers. This builds on the growing literature on the polyhedral method for post-selection inference lee_exact_2016,gao24,chen23. Our main motivation comes from gao24 and chen23, who compute selective $p$-values for testing equality of two cluster means in cross-sectional settings. Extending their methods to panels poses several nontrivial challenges. Unlike their pairwise focus, we test a joint null that all cluster means are zero. A recent generalization by yun2024selective considers joint equality across clusters, but applying it in our setting would require testing many constraints simultaneously, likely leading to poor small-sample performance.

Our methodology proceeds in three steps. First, we estimate cluster memberships and centers using a panel version of Lloyd’s Kmeans algorithm lloyd82, following bonhomme15a. The estimated centers capture average forecast performance differences within clusters. Second, we construct a test statistic based on the square root of a Wald statistic to measure forecast loss differences across clusters. As standard $\chi$ critical values are invalid, we condition on the estimated clusters, leading to a truncated $\chi$ distribution with analytically derived truncation sets. Third, we decompose the C-EPA null into $n_p = K(K-1)/2$ unique pairwise equality tests and an O-EPA test, then combine the resulting $p$-values using a combination method spreng23,vovk20,vovk2022admissible,gasparin2024combining.

Unlike much of the selective inference literature, which relies on strong assumptions such as normality, homoskedasticity, and independence gao24,chen23, our asymptotic theory accommodates heteroskedastic, dependent, and non-Gaussian panel data. We adopt a HAC variance estimator following sun13,sun14a, applied to cross-sectional averages of loss differentials. This yields test statistics robust to arbitrary forms and strengths of cross-sectional dependence (CD) driscoll98. We show that the tests are correctly sized and consistent under general alternatives, and that Panel Kmeans remains consistent under strong CD—extending beyond the weak dependence settings of bonhomme15a and patton23.

We assess the small sample properties of our tests through Monte Carlo simulations, comparing them to Split Sample statistics. The results show that our tests perform optimally even in very small samples, with negligible size distortions and substantial power under weak deviations from the C-EPA null.

We illustrate the empirical relevance of our method with an exchange rate forecasting application, comparing traditional time series models to modern machine learning approaches. Using a large panel of bilateral exchange rates against the U.S. dollar, we evaluate performance relative to an AR(1) benchmark. The results show substantial cluster heterogeneity and indicate that nonlinear models with macroeconomic fundamentals significantly outperform standard benchmarks. These findings are consistent with recent evidence in spreng23 and hillebrand2023exchange.

Section (ref) introduces the null and alternative hypotheses along with three motivating examples. Section (ref) develops the test statistics, while Section (ref) establishes their asymptotic properties. Section (ref) presents simulation results, while Section (ref) presents the empirical application. Section (ref) concludes. Additional material and proofs are reported in the Online Appendix (OA).

Setup and Motivating Examples

Testing Framework and Hypotheses

Let $\widehat{Y}_{a,it}$ denote the $\tau$-steps-ahead forecast, $\tau \geq 1$, of agent $a=1,2$ for the target variable $Y_{it}$, made at time $t-\tau$, for $t=1,\dots,T$ and $i=1,\dots,N$. The index $a$ represents a forecasting agent, such as the IMF or OECD (as in APUY and QTZ), or a forecasting model. To the best of our knowledge, there is no study deriving the asymptotic properties of the tests for comparing the out-of-sample forecasts made by panel data models in a theoretical level, though the corresponding time series literature is extensive west1996asymptotic,clark2001tests,clark13,clark2014tests,clark2015nested,giacomini06. Let $L(\cdot,\cdot)$ denote a generic loss function, which may be quadratic, absolute, or not necessarily in forecast error form gneiting11. Define the loss differentials as $\Delta L_{it} = L(\widehat{Y}_{1,it},Y_{it}) - L(\widehat{Y}_{2,it},Y_{it})$, where all variables are defined on a complete probability space $(\mathit{\Omega},\mathcal{E},\mathbb{P})$.

The null hypothesis of interest is the generalized C-EPA hypothesis, where “generalized” refers to the inclusion of conditioning variables—unlike the unconditional nulls in APUY and QTZ. It is stated as

where $\mathcal{F}_{t} \subseteq \mathcal{E}$ is a conditioning set, and $\mathcal{C}_k = \{i : k_i = k\}$, with $k_i \in \{1,\dots,K\}$ indicating cluster membership. The clusters are mutually exclusive and exhaustive: $\mathcal{C}_k \cap \mathcal{C}_g = \emptyset$ for $k \neq g$ and $\bigcup_{k=1}^K \mathcal{C}_k = \{1,\dots,N\}$. The alternative is

We implicitly assume the conditional expectations are time invariant almost surely. With more complex notation, one could instead consider time-averaged expectations, but this may require alternative variance estimation harvey2024testing or clustering methods.

Two special cases of the null hypothesis (ref) and its alternative are of particular interest. The first is the unconditional C-EPA hypothesis, obtained when $\mathcal{F}_{t} = \{\emptyset, \mathit{\Omega}\}$. For predetermined clusters, tests for this null have been developed by APUY and QTZ under various assumptions on autocorrelation and CD in loss differentials. The second is the conditional C-EPA hypothesis, which includes two useful sub-cases. First, let $\mathcal{F}_{t} = \sigma(\{W_{is}\}_{i=1}^N, s \leq t)$, where $W_{it} = (Y_{it},X'_{it})'$ includes external predictors $X_{it}$ used for $\widehat{Y}_{a,it}$. This yields a meaningful conditional null of the form (ref). Second, set $\mathcal{F}_{t} = \sigma(F_s, s \leq t)$, where $F_t$ denotes measurable-$\mathcal{E}$ common factors, such as dummies for the global financial crisis or COVID-19. Properly chosen, these factors allow detection of local differences in predictive ability.

The two conditioning schemes—on observed covariates and on common factors—are not mutually exclusive. In practice, forecast errors may arise from panel models that include both external predictors and common factors, with residuals exhibiting spatial or network dependence. Such models capture strong CD via factors and weak CD via spatial interactions chudik2011weak. As a result, loss differentials may reflect multiple CD sources due to model differences. While our framework accommodates general CD, explicitly modeling the CD structure could improve inference power (see APUY).

The null hypothesis $\mathcal{H}_{0}$ implies $|\mathcal{C}_k|^{-1} \sum_{i \in \mathcal{C}_k} \mathbb{E}(\widetilde{H}_{i,t-\tau} \Delta L_{it}) = 0$ for any measurable-$\mathcal{E}$ vector $\widetilde{H}_{it}$ giacomini06. Taking expectations with respect to $\widetilde{H}_{i,t-\tau}$ yields an unconditional moment condition. Let $H_{it}$ denote such a $P \times 1$ vector (a “testing function” in giacomini06), and $Z_{it} = H_{i,t-\tau}\Delta L_{it}$ with $\mu_{i}^0 = \mathbb{E}(Z_{it})$. Define $\theta^0_{k}(\mathcal{C}) = |\mathcal{C}_k|^{-1} \sum_{i \in \mathcal{C}_k} \mu_i^0$ where $\mathcal{C} = \{ \mathcal{C}_1,\dots,\mathcal{C}_K \}$. The null then implies

This transformation, standard in forecast evaluation and GMM settings, enables inference without explicitly modeling the $\sigma$-field $\mathcal{F}_{t-\tau}$. Although it does not preserve the full conditional distribution of $\Delta L_{it}$, it retains enough structure for testing, provided the test function is informative. In practice, the choice of $H_{i,t-\tau}$—e.g., lagged loss differentials, regressors, or common factors—affects both power and interpretation.

Examples

We present three examples illustrating the importance of accounting for unknown clusters in C-EPA testing.

Example 1: Time series forecasting. In time series forecasting, benchmark models, e.g. AR(1), are often compared to more flexible alternatives. For example, marcellino2006comparison compare direct and iterated AR forecasts across various macro series.

Consider $N$ bivariate time series $\{Y_{it},X_{it}\}_{t=0}^T$ generated from one of two latent clusters: \[ Y_{it} = \left\{

\right. \] with $U_{it} \sim iid(0, \sigma^2)$, and predictors being fixed quantities. Two forecasters have imperfect knowledge of the DGP and make the following forecasts:

The least squares estimators $\hat{\alpha}_i$, $\hat{\beta}_i$, and $\tilde{\beta}_i$ are computed from a fixed estimation window and are therefore subject to sampling variability. Each forecaster performs well on one cluster and poorly on the other. For $\mathcal{C}_1$, Forecaster 1 includes the correct intercept, while Forecaster 2 omits it and is biased. For $\mathcal{C}_2$, the true DGP has no intercept, so Forecaster 2 is correct, and Forecaster 1 overfits with an unnecessary constant.

This setup yields systematic differences in forecast accuracy across clusters. In Section (ref) of the OA, we derive the expected quadratic loss differential between the two forecasters, $\Delta L_{it} = \mathbb{E}[(\widehat{Y}^{(1)}_{i,T+1} - Y_{i,T+1})^2] - \mathbb{E}[(\widehat{Y}^{(2)}_{i,T+1} - Y_{i,T+1})^2]$, which is given by \begingroup

\endgroup where $\Delta_i = [ \mathbb{V}(\hat{\beta}_i) - \mathbb{V}(\tilde{\beta}_i) + \mathbb{B}(\hat{\beta}_i)^2 - \mathbb{B}(\tilde{\beta}_i)^2 ] X_{i,T}^2 + 2 X_{i,T} \text{Cov}(\hat{\alpha}_i, \hat{\beta}_i)$ with $\mathbb{B}(\cdot)$ denoting the bias of an estimator.

This decomposition highlights how heterogeneity in specification and precision drives cross-cluster performance gaps, motivating the C-EPA hypothesis as a testable implication of latent structure in forecast accuracy.

Example 2: Panel data forecasting. Latent group structures became popular in panel data analysis in the last decade bonhomme15a,su16,ando2017clustering,lumsdaine2023estimation. Suppose that two forecasters are interested in a variable $Y_{it}$ whose DGP is given by \[ Y_{it} = \beta'_{k_i} X_{i,t-1} + U_{it}, \quad U_{it} \sim iid(0, \sigma^2), \quad k_i \in \{1, \dots, K\}. \] We assume that the vector of predictors $X_{i,t-1}$ is known and fixed, and that the forecast errors $U_{it}$ are independent of all regressors. Two forecasters make the following two forecasts:

While the pooled estimator $\hat{\beta}$ suffers from misspecification bias if $\beta_{k_i} \neq \beta$, the individual estimator $\hat{\beta}_i$ is unbiased but suffers from increased variance due to limited time series observations. Let $\Delta L_{it} = \mathbb{E}[(\widehat{Y}^{\text{pooled}}_{i,T+1} - Y_{i,T+1})^2] - \mathbb{E}[(\widehat{Y}^{\text{het}}_{i,T+1} - Y_{i,T+1})^2]$. Under standard regularity conditions, we have

where $\Sigma_X = |\mathcal{C}_k|^{-1} \sum_{i \in \mathcal{C}_k} X_{i,T} X_{i,T}'$ is the empirical second moment matrix of regressors in cluster $\mathcal{C}_k$, and $\overline{\mathbb{V}(\hat{\beta}_i)} = |\mathcal{C}_k|^{-1} \sum_{i \in \mathcal{C}_k} \mathbb{V}(\hat{\beta}_i)$ is the average variance of unit-specific estimators. The proof is given in Section (ref) of the OA. This shows how strong group-level heterogeneity leads to systematic differences in forecast performance across units.

Example 3: Forecasting with machine learning methods. Machine learning methods are increasingly popular in economics athey2018impact, haghighi2025machine. In high-dimensional forecasting, researchers often compare linear approaches like LASSO to nonlinear ones such as random forests (RF). For example, goulet2022machine examine various data-rich and data-poor models, finding that ML methods have the advantage of capturing nonlinearities linked to uncertainty, financial stress, and housing bubbles. Suppose that two methods are trained and evaluated using validation MSE:

When only some units exhibit nonlinear patterns, averaging MSE across units can obscure performance differences. To address this, one might apply a second ML tool—clustering—on forecast loss differentials. Testing the C-EPA null then helps reveal cluster-specific model dominance. If it were not already in use, we would label this usage of our proposed method “double machine learning.”

Test Statistics

We begin by decomposing the C-EPA hypothesis into two components: homogeneity and O-EPA. The null hypothesis (ref) can be written as $\mathcal{H}'_{0}: \mathcal{H}^{homo}_{0} \cap \mathcal{H}^{oepa}_{0}$, where

is the homogeneity hypothesis, and

is the O-EPA hypothesis, where overall predictive performance difference is a weighted average of cluster means. The O-EPA parameter is invariant to the specific clustering used.

Both $\mathcal{H}^{homo}_{0}$ and $\mathcal{H}^{oepa}_{0}$ are empirically relevant. Tests of the unconditional O-EPA hypothesis with known clusters have been analyzed by APUY under various CD assumptions. Testing $\mathcal{H}^{homo}_{0}$ is important beyond EPA contexts; see patton23. In Section (ref) of the OA, we develop a test for $\mathcal{H}^{homo}_{0}$.

Testing Pairwise Equality with Unknown Clusters

We begin by introducing the Panel Kmeans Estimator of the clusters. When no prior information is available on the clusters $\mathcal{C}_k$, $k = 1,\dots,K$, one may estimate them using Panel Kmeans applied to the stacked panel $Z = (Z'_{11},Z'_{12},\dots,Z'_{NT})'$, denoted $\mathcal{C}(Z)$. For a given $K$, the cluster memberships and centers are defined by:

where $\widehat{\mathcal{C}}_k = \{i : \hat{k}_i(Z) = k\}$, with $\hat{k}_i(Z) \in \{1,\dots,K\}$ indicating the estimated cluster membership. This optimization is typically solved by an iterative algorithm lloyd82,hartigan75. Algorithm (ref) in Section (ref) of the OA implements a generalized version of Lloyd’s method for computing these estimates where we also discuss practical aspects of the algorithm. The theoretical properties of the Panel Kmeans Estimator are presented in Section (ref).

We now develop a test for each pairwise sub-hypothesis in (ref). The homogeneity null $\mathcal{H}^{homo}_{0}$ is the intersection of $n_p = K(K-1)/2$ distinct pairwise equalities. For each $k,g \in \{1,\dots,K\}$, $k\neq g$, we define the test statistic $D_{k,g}(\widehat{\mathcal{C}})$ as the square root of the corresponding Wald statistic:

where $ \widehat{\Sigma}_{k,g}(\widehat{\mathcal{C}}) = \widehat{\omega}_{k,k}(\widehat{\mathcal{C}}) + \widehat{\omega}_{g,g}(\widehat{\mathcal{C}}) - 2\widehat{\omega}_{k,g}(\widehat{\mathcal{C}}), $ and $\widehat{\omega}_{k,g}(\widehat{\mathcal{C}})$ is the $\{k,g\}$th $P \times P$ block of $\widehat{\Omega}(\widehat{\mathcal{C}})$, an orthonormal series (OS) variance estimator given by

with $\bar{Z}_t(\widehat{\mathcal{C}}) = [\bar{Z}'_{1,t}(\widehat{\mathcal{C}}),\dots,\bar{Z}'_{K,t}(\widehat{\mathcal{C}})]'$, $\bar{Z}_{k,t}(\widehat{\mathcal{C}}) = |\widehat{\mathcal{C}}_k|^{-1} \sum_{i \in \widehat{\mathcal{C}}_k} Z_{it}$ and $\widehat{\theta}(\widehat{\mathcal{C}}) = [\widehat{\theta}^{\prime}_{1}(\widehat{\mathcal{C}}),\dots,\widehat{\theta}^{\prime}_{K}(\widehat{\mathcal{C}})]'$. We use the square root of the Wald statistic because its decomposition into its norm and direction is a linear function of the data, which is essential for deriving the truncation region in its conditional distribution (see Equation (ref) below). We refer to the discussion in Section (ref) for desired properties of the OS estimator. Under regularity conditions, $D_{k,g}(\mathcal{C}) \overset{d}{\longrightarrow} \chi_{\scriptscriptstyle P}$ as $(T,N) \to \infty$ for fixed $\mathcal{C}$, where $\overset{d}{\longrightarrow}$ denotes convergence in distribution. However, as discussed in the introduction, critical values from this limiting distribution are invalid when clusters are estimated. We therefore define the asymptotic selective Type I error rate as the basis for valid testing under unknown clusters.

A valid test of $\mathcal{H}^{k,g}_0$ controls the selective Type I error at level $q$, conditional on the clustering produced by the Panel Kmeans Algorithm. Specifically, the conditioning event in (ref) implies that $\mathcal{H}^{k,g}_0$ is rejected if the probability of observing a test statistic at least as large as the realized one does not exceed $q$ over all $Z$ yielding the same clustering as $z$.

As noted by chen23, directly characterizing the conditioning set is nontrivial. Instead, we condition on the cluster assignments obtained at each iteration $m = 1,\dots, M$ of the algorithm. Two additional conditioning terms emerge from a decomposition of $Z$ into components aligned with and orthogonal to the test statistic $D_{k,g}(\widehat{\mathcal{C}})$:

where $\hat{J}_Z = \mathrm{dir}[\widehat{\Sigma}^{-1/2}_{k,g}(\widehat{\mathcal{C}})Z'\hat{\nu}_{k,g}]$ and \[ \widehat{\Pi}_{k,g} = I - \frac{\hat{\nu}_{k,g} \hat{\nu}_{k,g}'}{\lVert \hat{\nu}_{k,g} \rVert^2}, \quad \hat{\nu}_{k,g,i} = \iota_T \hat{\delta}_{k,g,i}, \quad \hat{\delta}_{k,g,i} = \frac{\mathbf{1} \{ \hat{k}_i(Z) = k \}}{|\widehat{\mathcal{C}}_k|} - \frac{\mathbf{1} \{ \hat{k}_i(Z) = g \}}{|\widehat{\mathcal{C}}_g|}, \] and $\iota_T$ is a $T \times 1$ vector of ones. This decomposition is derived in Section (ref) of the OA and forms the basis for characterizing the conditional distribution of $D_{k,g}(\widehat{\mathcal{C}})$ given $\widehat{\mathcal{C}}$. Namely, it decomposes the observed data $Z$ into two orthogonal components. First one determines the value of the test statistic $D_{k,g}(\widehat{\mathcal{C}})$, and the second one remains invariant under perturbations of the test statistic in its direction. By conditioning on both the orthogonal projection $\widehat{\Pi}_{k,g} Z$ and the direction $\hat{J}_Z$, we are able to hold fixed the information that does not affect clustering. This in turn enables us to characterize the truncated distribution of $D_{k,g}(\widehat{\mathcal{C}})$ conditional on the clustering outcome $\widehat{\mathcal{C}}$, as detailed in Section (ref) of the OA.

Following this discussion, we define the asymptotic $p$-value for testing $\mathcal{H}^{k,g}_0$ as

for $k,g \in \{1,\dots,K\}$, where the conditioning set is defined as

with $\hat{J}_z = \mathrm{dir}[\widehat{S}^{-1/2}_{k,g}(\widehat{\mathcal{C}})z'\hat{\nu}_{k,g}]$, $\widehat{S}_{k,g}(\widehat{\mathcal{C}})$ denoting the realization of $\widehat{\Sigma}_{k,g}(\widehat{\mathcal{C}})$ associated with $z$. The first condition in $\mathcal{A}$ is central to the selective conditional inference framework: it requires that each unit’s cluster assignment at every iteration $m$ of the Panel Kmeans Algorithm using $Z$ matches that from the observed realization $z$, i.e., $k^{(m)}_i(Z) = k^{(m)}_i(z)$. This ensures we condition on the event that $Z$ yields the same clustering as $z$, as required by Definition (ref). The remaining two conditions remove the nuisance terms $\widehat{\Pi}_{k,g} Z$ and $\hat{J}_Z$ in (ref), which would otherwise make the conditional distribution of $D_{k,g}(\widehat{\mathcal{C}})$ intractable. These are standard in the selective inference literature gao24,chen23.

The asymptotic $p$-value $p_{\infty} [ d_{k,g}(\widehat{\mathcal{C}}) ]$ is based on the selective inference methodology of chen23 but it generalizes it in several ways. First of all, here, we have double indexed random variables $Z_{it}$, $i=1,\dots,N$, $t=1,\dots,T$. Second, their study does not allow for dependencies between $Z_{it}$ and $Z_{js}$, for either $i \neq j$ or $t \neq s$, but only across different variables of the same observation, i.e. between $Z_{p,it}$ and $Z_{c,it}$, the $p$-th and the $c$-th elements of $Z_{it}$. Whereas, we allow for arbitrary autocorrelation and CD as well as dependencies between different elements of $Z_{it}$. Third, their method depends crucially on the normality of the data generating process, whereas we make use of a CLT (see Lemma (ref) below) by exploiting the time series dimension of the data.

Next proposition shows how to calculate a $p$-value in observed samples following this definition under standard assumptions, which we present in Section (ref).

The vector $z(\phi)$ defines a perturbation of the original data $z$. Varying $\phi$ moves clusters $k$ and $g$ closer or farther apart along the direction $\widehat{S}^{-1/2}_{k,g}(\widehat{\mathcal{C}})z'\hat{\nu}_{k,g}$. When $\phi = d_{k,g}(\widehat{\mathcal{C}})$, $z(\phi) = z$; for $\phi > d_{k,g}(\widehat{\mathcal{C}})$, the clusters are pulled apart; and for $\phi < d_{k,g}(\widehat{\mathcal{C}})$, they are pushed together—with $\phi = 0$ implying identical centers. Thus, $\phi$ measures the degree of perturbation chen23. Switching from Kmeans to Panel Kmeans alters the geometry of the selection region, requiring new derivations for truncation sets. We outline the steps for computing the selective $p$-value in Section (ref) of the OA via a characterization of the truncation set $\mathcal{T}$ for Panel Kmeans.

The O-EPA Test

The second sub-hypothesis of the C-EPA hypothesis (ref), namely $\mathcal{H}^{oepa}_{0}$, states that the two forecasts are equally good on average given past information. To test this sub-hypothesis, consider the test statistic $$ W_{oepa} = a_{B} T\bar{Z}'_{o} \widehat{{\Omega}}^{-1}_{o} \bar{Z}_{o},$$ where $a_{B} = (B-P+1)/(PB)$, $\bar{Z}_{o} = T^{-1} \sum_{t=1}^T \bar{{Z}}_{t}$, $\bar{{Z}}_{t} = N^{-1} \sum_{i=1}^N Z_{it}$, and $\widehat{{\Omega}}_{o}$ is given by $$\widehat{\Omega}_{o} = B^{-1} \sum_{j=1}^B \widehat{\Lambda}_{o,j} \widehat{\Lambda}_{o,j}', \quad \widehat{\Lambda}_{o,j} = \sqrt{2/T} \sum_{t=1}^T [\bar{{Z}}_{t} - \bar{Z}_{o}] \cos \left[ \pi j ( t-1/2 )/T \right].$$ The test rejects the O-EPA null if $p(w_{oepa}) = \mathbb{P}_{\mathcal{H}_0} \left[ \mathbb{F}_{P,B-P+1} \geq w_{oepa} \right] \leq q$, where $q \in (0,1)$ is the nominal Type I error rate. When $B = T$ and $P = 1$, the statistic reduces to a Wald-type test that is robust to cross-sectional dependence but ignores autocorrelation. This corresponds to the $S^{(3)}$ test of APUY with a bandwidth set to zero.

The C-EPA Test with Unknown Clusters

We now introduce the main test statistic for the C-EPA null $\mathcal{H}_0$. Based on the results of the previous sections, we define a $p$-value combination statistic that aggregates the $n_p$ pairwise tests and the O-EPA test which is given by

where $r \in [-\infty, -1)$.

This test statistic belongs to the class of precise merging functions, satisfying both monotonicity and sharpness properties under arbitrary dependence of the input $p$-values. The normalization factor $[r/(r+1)] (n_p+1)^{1+1/r}$ guarantees that the statistic in (ref) defines a valid $p$-value under the global null hypothesis. This is shown in Theorem 2 of vovk20 and generalized in Theorem 3 of vovk2022admissible, where the authors establish the admissibility and optimality of such M-family-based merging functions. In particular, the proposed $F_{SI,r}$ controls the family-wise Type I error under any form of dependence between the constituent $p$-values.

Unlike Fisher's method fisher32, which assumes independence, or Bonferroni's $p$-merging function, which is conservative, this choice of merging function maintains optimal Type I control under general dependence structures.

A similar $p$-merging function was recently used by spreng23 in a multiple forecast comparison setting. The difference between our proposal and that of the authors lies on the choice of the calibration constant $b_{r,n_p}$, using the notation of vovk20. While spreng23 sets $b_{r,n_p} = r/(r+1)$, we follow exactly the constant suggested by Proposition 5 of vovk20 and set $b_{r,n_p} = [r/(r+1)] (n_p+1)^{1+1/r}$. We found that this choice results in smaller size distortions in our particular framework with a small number of $p$-values combined.

Asymptotic Theory

Assumptions and Two Useful Lemmata

We state the assumptions and two preliminary results underlying the asymptotic theory of the proposed tests. We introduce some new notation: $C$ denotes a generic positive constant, and $(T, N) \to \infty$ refers to joint divergence with $N = N(T)$ growing as $T \to \infty$. Let $V_{it} = Z_{it} - \mu_i^0$, and denote its $p$th element by $V_{p,it}$ for $p = 1,\dots,P$. The first three assumptions below (G$\#$) are generic, required for both size and power; the last three (S$\#$) are specific to power under the alternative hypothesis $\mathcal{H}_1$.

Assumptions (ref)(ref) and (ref)(ref) ensure well-defined cluster centers and finite moments up to the fourth, so that means and variances are consistently estimable under regularity. Assumption (ref)(ref) restricts time dependence. No restriction is imposed on CD, which may be weak or strong (see discussion after Lemma (ref)).

Assumption (ref) controls cluster sizes asymptotically. It is standard in the clustering literature (e.g., Assumption 2(a) of bonhomme15a, A1(vii) of su16) and requires each cluster to have non-negligible mass. This could be relaxed at the cost of more complex notation.

Assumption (ref) imposes mixing conditions. The matrix $\Omega_i$ is assumed positive definite—a requirement for Diebold-Mariano-type EPA tests west1996asymptotic. It holds when forecasts come from non-nested models or nested models under conditions in giacomini06, such as fixed or rolling estimation windows. Expanding windows are excluded for nested comparisons clark2015nested,mccracken2020diverging,zhu2022can.

Assumption (ref) requires identical means within clusters but different across them. Assumption (ref) imposes a lower bound on inter-cluster distances, ensuring well-separated centers and thus violation of $\mathcal{H}_0$. While not required, this guarantees test power. Notably, even when $K^0 = 1$, the tests may reject $\mathcal{H}_0$ if the overall mean differs from zero, as shown below.

Assumption (ref) strengthens dependence and tail conditions on $V_{it}$ beyond Assumptions (ref) and (ref), ensuring consistent estimation of cluster memberships and asymptotic equivalence between Panel Kmeans and oracle estimators.

We now state two lemmata essential for the theoretical analysis of the test statistics. Define the $KP \times 1$ vectors $\widehat{\theta}(\mathcal{C}) = [\widehat{\theta}^{\prime}_{1}(\mathcal{C}),\dots,\widehat{\theta}^{\prime}_{K}(\mathcal{C})]'$, $\theta^0(\mathcal{C}) = [\theta^{0\prime}_{1}(\mathcal{C}),\dots,\theta^{0\prime}_{K}(\mathcal{C})]'$ and let $ \Omega(\mathcal{C}) = \mathbb{V} \{ \sqrt{T} [\hat{\theta}(\mathcal{C}) - \theta^0(\mathcal{C})] \}$, $\mathcal{N}(\mathcal{C}) = \mathrm{diag}(|\mathcal{C}_1|,\dots,|\mathcal{C}_K|) \otimes I_P. $ The following result gives the standard properties of sample means for a fixed clustering $\mathcal{C}$. This remains useful even when clusters are estimated, but inference is conditional on them, as will be in our case.

Part (ref) of Lemma (ref) establishes consistency of sample means for fixed cluster assignments under Assumptions (ref)--(ref). Part (ref) is a CLT. The scalar $\epsilon \in [1/2,1]$ captures the degree of CD: $\epsilon = 1$ corresponds to strong CD (e.g., factor models), while $\epsilon \in [1/2,1)$ covers weak CD (e.g., spatial models or independence). See chudik2011weak for a thorough discussion, and bailey2016exponent for methods to estimate $\epsilon$. The parameter $\epsilon$ allows for a unified treatment of strong and weak CD. While Lemma (ref) applies to fixed clusters, we use it in a conditional framework to analyze tests with estimated clusters.

Define $\theta^0 := \theta^0(\mathcal{C}^0)$, that is, the true centers of the true clusters of the population. The following lemma establishes the properties of the Panel Kmeans Estimators when clusters are well separated, Assumption (ref) in particular.

Lemma (ref) establishes the properties of Panel Kmeans when the clusters are well-separated. Based on this result, a naive test of C-EPA would estimate clusters using Panel Kmeans and plug them into a Wald statistic resulting in $W(\widehat{\mathcal{C}})$. The test rejects the null if $p[w(\widehat{\mathcal{C}})] \leq q$ for some $q \in (0,1)$. However, this approach is invalid under $\mathcal{H}_0$: the clusters are homogeneous and they are estimated from the same data used for testing.

Recent work patton23,chen23,gao24 shows that testing for homogeneity after clustering yields anti-conservative tests. Clustering under the null typically produces artificially separated group means, inflating Type I error rates unless the selection step is accounted for. The null hypotheses in these studies are nested within ours, so their critique applies here. We demonstrate the failure of this naive approach through simulations in Section (ref).

Main Results

This section establishes the asymptotic properties of the proposed test statistics. The first result is on the asymptotic validity of $p[D_{k,g}(\widehat{\mathcal{C}})]$ for testing the pairwise homogeneity null $\mathcal{H}^{k,g}_0$ defined in Definition (ref).

Part (ref) shows that $p[D_{k,g}(\widehat{\mathcal{C}})]$ is asymptotically a $p$-variable in the sense of vovk20 under the null of pairwise cluster equality. Following the convention, we refer to both $p[D_{k,g}(\widehat{\mathcal{C}})]$ and its realization $p[d_{k,g}(\widehat{\mathcal{C}})]$ as $p$-values. Part (ref) establishes the consistency of $D_{k,g}(\widehat{\mathcal{C}})$ when $\mathcal{H}^{k,g}_0$ fails, assuming $K^0$ is known, i.e. $K = K^0$. This assumption is relaxed in Section (ref) of the OA, where we introduce an IC as well as a crosss-validation (CV) method to estimate $K^0$.

Next, we establish the asymptotic properties of the O-EPA test statistic which is the second main component of our proposed test of C-EPA.

Part (ref) of the theorem shows that the limiting distribution of the test statistic is an $\mathbb{F}_{P,B-P+1}$ variate for fixed $B$. When $B \longrightarrow \infty$, we have $W_{oepa}/a_B \overset{d}{\longrightarrow} \chi^2_{\scriptscriptstyle KP}$ which follows as a corrollary to the theorem. The results of sun13 show that when $B$ is not large, using the $\mathbb{F}_{KP,B-KP+1}$ critical values instead of (scaled) $\chi^2_{\scriptscriptstyle KP}$ critical values results in better size properties. Part (ref) of the theorem shows that the test statistic is consistent.

Having established the properties of the pairwise Homogeneity and O-EPA tests, we now turn to those of the proposed C-EPA test. The following result summarizes the desired asymptotic properties of (ref).

The asymptotic result shows that the proposed selective inference test successfully controls the Type I error rate and it is consistent as its power approaches one when either $\mathcal{H}^{homo}_{0}$ or $\mathcal{H}^{oepa}_{0}$ fails. The finite sample properties of the test statistic are investigated in Section (ref) where the simulation results confirm these theoretical expectations.

Monte Carlo Study

We study the finite sample size and power properties of the test statistics. In Section (ref) we describe the Monte Carlo design and in Section (ref) we report and comment on the results.

Design

To investigate the finite sample properties of the testing procedures, we generate observations from a panel AR(1) process given by:

This DGP, as well as our setup that we describe below, is similar to that of hoga2023testing except that their focus is on measurement errors in the target variable whereas ours is on clustered heterogeneity.

Two forecasters, indexed by \( a = 1,2 \), aim to construct one-step-ahead forecasts of \( Y_{it} \) without observing the true data-generating process. Forecaster 1 includes an intercept but adds noise, while Forecaster 2 omits the intercept. Their models are:

for \( t = 1,\dots,T \) and \( i = 1,\dots,N \), where \( k_i \in \{1,2,3\} \) indicates latent cluster membership. Following hoga2023testing, we assume both forecasters use the true slope and, if applicable, the true intercept. This is justified by noting that the noise in Forecaster 1 may reflect overfitting to heterogeneity, while Forecaster 2’s misspecification omits the intercept.

The noise term \( \varepsilon_{it} \) is constructed to have zero mean and cluster-specific forecast variance, and evolves as a stationary process: \[ \varepsilon_{it} = \phi \varepsilon_{i,t-1} + \lambda F_t + \sqrt{ \sigma^2_{\varepsilon,k_i} (1 - \phi^2) - \lambda^2 } \cdot \xi_{it}, \qquad \xi_{it} \sim iid \, \mathbb{N}(0,1), \] where \( F_t \sim iid \, \mathbb{N}(0,1) \) is a common factor independent of \( \xi_{it} \). The parameter \( \phi \in (-1,1) \) governs AR(1) persistence, and \( \lambda \) controls the strength of cross-sectional dependence (CD) via \( F_t \). The forecast variance for Forecaster 1 in cluster \( k_i \) is \( \sigma^2_{\varepsilon,k_i} = \alpha^2 (1 - \rho_{k_i})^2 + \psi_{k_i} \).

We implement both unconditional and conditional EPA tests, corresponding to $H_{i,t-1} = 1$ and $H_{i,t-1} = (1, Y_{i,t-1})'$, respectively. Let $\Delta L_{it} = (Y_{it} - \widehat{Y}_{it}^{(1)})^2 - (Y_{it} - \widehat{Y}_{it}^{(2)})^2$. By straightforward calculations hoga2023testing, we have: \begingroup \[ \mathbb{E} (H_{i,t-1} \Delta L_{it}) = \left\{

\right. \] \endgroup Thus, the expected loss differential depends solely on the noise variance $\psi_{k_i}$ in the unconditional case, and on both $\psi_{k_i}$ and the unconditional mean $\mu$ in the conditional case.

In all experiments, we set $\mu = 1$, $\phi = 0.2$, and $\lambda = 0.2$. For the AR(1) process of $Y_{it}$, panel units are divided into three latent clusters of unequal sizes, aligned with the structure of the loss differentials: \begingroup

\endgroup with $(\rho_1, \rho_2, \rho_3) = (0.1, 0.2, 0.3)$, so that Cluster 3 is twice as large as Cluster 1 and Cluster 2. To assess size, we set $(\psi_1, \psi_2, \psi_3) = (0, 0, 0)$. Power is examined under two alternatives with $K^0 = 3$: Case 1 — O-EPA fails: $(\psi_1, \psi_2, \psi_3) = \psi/2 + \psi \cdot (-1.2, -0.8, 1)$, Case 2 — O-EPA holds: $(\psi_1, \psi_2, \psi_3) = \psi \cdot (-1.2, -0.8, 1)$. The parameter $\psi$ governs deviation from the null, with values $\psi \in \{0.125, 0.25, 0.375, 0.5\}$. We assess size across all $(T, N)$ combinations with $N \in \{80,120,160\}$ and $T \in \{20,50,100,200\}$. Due to the computational cost of the proposed procedures, power analysis is restricted to $N = 80$ and $T \in \{50,200\}$. As the loss differentials exhibit strong CD, increasing $N$ has little to no effect on power. All results are based on 1000 replications.

We implement four types of tests: Predetermined, Naive, Split Sample, and Selective Inference. Each is conducted under both unconditional and conditional specifications. Implementation details are as follows: Predetermined: As in Section (ref) of the OA, with $k_i = k_i^0$ for all $i = 1,\dots,N$. Naive: As in Section (ref) of the OA, with $k_i = \hat{k}_i(Z)$ from Algorithm (ref). Split Sample: As in Section (ref) of the OA, using $\mathcal{S}_1 = \{1,\dots,0.2 \cdot T\}$ for training and $\mathcal{S}_2 = \{0.2 \cdot T + 1 + l,\dots,T\}$ for testing, with $l = \lfloor \sqrt{0.2 \cdot T} \rfloor$; clusters are $k_i = \hat{k}_i(Z_{\mathcal{S}_1})$, i.e. output of Algorithm (ref) with input $Z_{\mathcal{S}_1}$, data corresponding to the training sample $\mathcal{S}_1$. Selective Inference: As in Section (ref), with $k_i = \hat{k}_i(Z)$ from Algorithm (ref). All tests are robust to arbitrary autocorrelation and CD. The number of cosines in the LRV estimator is set as $B = \min(\lfloor P T^{2/3} \rfloor, T)$ for full-sample tests, and $B = \min(\lfloor P |\mathcal{S}_2|^{2/3} \rfloor, |\mathcal{S}_2|)$ for Split Sample tests. Since latent clustering is central to our framework, all tests—except Predetermined—are implemented using $\widehat{K}_{IC}$ given in Section (ref) of the OA. When applicable, Algorithm (ref) is run with 10 random initializations and a maximum of 100 iterations.

Results

We report the results in two parts, size and power properties, respectively. A robustness check for structural breaks in the process is reported in Section (ref) of the OA.

Table (ref) reports rejection rates of the four C-EPA tests under the null, evaluated at the 5% level, separately for unconditional and conditional versions. The Naive test, which treats estimated clusters as known, rejects 100% of the time in all configurations. This highlights the risk of ignoring model selection when clusters are data-driven.

In contrast, the Predetermined test—using fixed, exogenous clusters—yields rejection rates near the nominal level, ranging from 0.04 to 0.07. For instance, with $N = 120$ and $T = 100$, rejection rates are 0.05 (unconditional) and 0.06 (conditional). While a useful benchmark, its reliance on known cluster structure limits practical use.

The Split Sample test also shows reasonable size control, with rejection rates between 0.03 and 0.11. For example, with $N = 160$ and $T = 20$, the rates are 0.07 (unconditional) and 0.09 (conditional)—slightly above nominal but acceptable in small samples. By using disjoint subsamples, it reduces selection bias but sacrifices power due to smaller samples.

The Selective Inference test, which adjusts for cluster estimation via truncation-based conditioning, consistently achieves accurate size control. Rejection rates stay close to the 5% level—for instance, 0.05 (unconditional) and 0.06 (conditional) at $N = 80$, $T = 50$. This confirms that the method effectively corrects for data-driven clustering without requiring sample splitting or external information.

To sum up, Naive test leads to severe over-rejection, while Split Sample and Selective Inference maintain valid size. Among feasible methods, Selective Inference offers the most reliable size performance across a wide range of settings.

First part of Table (ref) reports rejection rates for the four C-EPA tests under the alternative where the O-EPA hypothesis fails. As expected, all tests gain power as $\psi$ and $T$ increase, though at different rates. The Naive test rejects nearly 100% of the time, regardless of sample size or effect strength. The Predetermined test, which uses true cluster assignments, performs well—e.g., at $T = 50$, $\psi = 0.125$, power is 87% (unconditional) and 74% (conditional); with $T = 200$, it reaches 100% in all cases.

The Split Sample test shows lower power for small $T$ and weak signals—only 20% (unconditional) and 15% (conditional) at $T = 50$, $\psi = 0.125$—but improves with larger $T$, reaching 100% at $T = 200$, $\psi = 0.25$. This reflects the efficiency-size trade-off of data splitting.

The Selective Inference test behaves similarly but often outperforms Split Sample in conditional settings. At $T = 50$, $\psi = 0.125$, power is 19% (unconditional) and 16% (conditional); at $T = 200$, $\psi = 0.25$, power reaches 100%.

In summary, under O-EPA violations, both Split Sample and Selective Inference control size and achieve high power as signal strength grows. The Predetermined test sets an upper bound, while Selective Inference offers a robust alternative that avoids over-rejection.

Second part of Table (ref) reports rejection rates under the alternative where O-EPA holds but heterogeneity exists within clusters. This setting evaluates whether tests can detect within-cluster predictive differences despite similar overall performance.

The Predetermined test sets a power benchmark. Rejection rates are high in all cases—even in small samples and weak deviations (e.g., $T = 50$, $\psi = 0.125$, power is 0.78 unconditional and 0.65 conditional)—confirming that signals are detectable under ideal clustering.

As before, the Naive test always rejects (power = 1.00) regardless of signal strength. The Split Sample test performs well: power is low for weak signals and short panels (e.g., 0.07 at $T = 50$, $\psi = 0.125$), but increases rapidly. At $T = 50$, $\psi = 0.375$, power reaches 80% (unconditional) and 57% (conditional); for $T = 200$ and $\psi \geq 0.25$, rejection exceeds 95%.

The Selective Inference test shows lower power in this setting. For $T = 50$, $\psi = 0.125$, rejection is near nominal (0.06). Power rises gradually: at $T = 200$, $\psi = 0.375$, power is 31% (unconditional) and 53% (conditional); for $\psi = 0.5$, it improves to 64% and 67%. This reflects two factors: (i) additional conditions due to the nuisances in the conditional distribution limit power; and (ii) inclusion of the O-EPA test in the $p$-value combination reduces sensitivity when O-EPA holds.

Despite lower power when O-EPA holds, the selective test is the only viable C-EPA method in most empirical settings. It controls false positives but may under-reject when deviations are subtle. However, in many realistic empirical settings Split Sample statistics may fail while Selective Inference keeps its validity (see Section (ref) of the OA).

Empirical Illustration

This section implements alternative forecasting techniques for monthly exchange rate returns to compare the performance of machine learning techniques with the AR(1) benchmark.

Data

The empirical analysis uses monthly bilateral exchange rates from the IMF (1999–2023) and macroeconomic predictors from FRED-MD, resulting in a balanced panel of 131 series after standard filtering. Forecasts are constructed recursively using a fixed 60-month window, yielding $T = 238$ one-step-ahead forecast errors. To assess model performance, we compute quadratic loss differentials relative to an AR(1) benchmark across a wide range of models. Descriptive results reveal that while the AR(1) model is difficult to beat uniformly, more flexible methods—particularly XGBoost—deliver substantial gains in specific environments, especially where the benchmark model performs poorly. Regularized linear models like Elastic Net (EN) offer smaller but more stable improvements. Full details on data handling, forecast design, and summary statistics are provided in Section (ref) of the OA.

Results

Table (ref) reports the $p$-values from a series of C-EPA tests applied to loss differentials between five forecasting models and the AR(1) benchmark. The aim is to detect whether the models improve predictive accuracy overall or within specific clusters of currency pairs.

We first look at the O-EPA test results. We see that for all models but SVM, the O-EPA hypothesis is rejected at least at the 10% level in all settings. It is seen in the summary statistics reported in Table (ref) of Section (ref) of the OA that AR($p$), XGBoost and RF perform better overall with respect to AR(1), whereas EN is worse. Hence, in an unconditional setting, the superiority of the first three methods and the inferiority of the last, against AR(1), are confirmed by the O-EPA test results.

Across all settings, SVM stands out as the only method consistently associated with very high $p$-values in the O-EPA test (e.g., 0.90, 0.95, 0.97), indicating no statistically significant improvement over AR(1) on average over all units and time periods. However, these high $p$-values do not imply poor performance; rather, they reflect that gains are not homogeneous across all cross-sectional units. This interpretation is supported by the rejection of the Homogeneity test at the 10% level in the conditional test with the lagged target and when the number of clusters is chosen by CV ($p$-value = 0.07). This suggests that SVM's performance is heterogeneous conditional on the past realization of the target variable. Moreover, the selective inference C-EPA test is significant at the 10% level ($p$-value = 0.09).

More generally, the rejection of the homogeneity null in several cases justifies the use of our Selective Inference C-EPA testing procedure. For example, when conditioning on the lagged target variable, the Homogeneity test rejects for SVM and XGBoost depending on the clustering method, and in many cases selective C-EPA $p$-values very low (e.g., RF yields a $p$-value of 0.00 in all settings.). These results confirm that forecast gains may vary across clusters, making clustered tests essential to discover such patterns.

Overall, these results highlight that clustered inference can detect model improvements that are missed by aggregate tests, and that conditioning and clustering are both essential tools in evaluating forecast performance in panel settings with heterogeneous effects.

Conclusion

This paper developed a statistical framework for testing hypotheses on the cluster centers of a panel process after estimating the clusters via Panel Kmeans. We applied this framework to conditional C-EPA testing to compare forecast performance across agents or models. To address the “double dipping” problem, we proposed a conditional testing procedure based on advances in selective inference. The method computes a $p$-value for the C-EPA hypothesis interpreted as the rejection frequency under the null across realizations yielding the same clustering. We compared its performance to that of simpler Split Sample tests, both theoretically and through Monte Carlo simulations.

Simulations show both methods perform well in small samples: they are correctly sized and have power against relevant alternatives. Selective inference tests, in particular, perform strongly and emerge as the preferred method given their theoretical and practical advantages.

Finally, using a large panel of exchange rates, we compared alternative time series and machine learning models to an AR(1) benchmark. The results show that accounting for latent clusters in forecast loss differentials can substantially improve predictive performance.

\if11 {

Acknowledgements

We thank Lucy L. Gao, Antonio Montañés, Ryo Okui, Hashem Pesaran, Esther Ruiz-Ortega, and the participants of the Centre for Econometric Analysis Occasional Econometrics Seminar at Bayes Business School (London, 27 Oct. 2023), the Annual Spatial Econometrics Association Conference (San Diego, 16–17 Nov. 2023), the 29th International Panel Data Conference (Orléans, 3–5 July 2024), and the GIAM Seminar at Galatararay University (Istanbul, 29 April 2025), for their helpful comments and suggestions.

Disclosure and Data Availability

The authors report there are no competing interests to declare. All methods are implemented in the clusteredEPA and PanelKmeansInference R packages, which, together with replication materials and data, are in \url{https://github.com/akoguzhan/}. } \fi

\if01 {

} \fi

\setcounter{section}{0} \setcounter{subsection}{0} \setcounter{subsubsection}{0} \setcounter{equation}{0} \setcounter{figure}{0} \setcounter{table}{0} \setcounter{footnote}{0}