Regularizing stock return covariance matrices via multiple testing of correlations

{{\sc{Regularizing stock return covariance matrices via multiple testing of correlations}} }

{Richard Luger}\footnotemark[1]

{Universit{\'e} Laval, Canada}

\footnotetext[1]{Correspondence to: Department of Finance, Insurance and Real Estate, Laval University, Quebec City, Quebec G1V 0A6, Canada.

E-mail address: {[email removed]}.

}

{\bf Abstract:} This paper develops a large-scale inference approach for the regularization of stock return covariance matrices. The framework allows for the presence of heavy tails and multivariate GARCH-type effects of unknown form among the stock returns. The approach involves simultaneous testing of all pairwise correlations, followed by setting non-statistically significant elements to zero. This adaptive thresholding is achieved through sign-based Monte Carlo resampling within multiple testing procedures, controlling either the traditional familywise error rate, a generalized familywise error rate, or the false discovery proportion. Subsequent shrinkage ensures that the final covariance matrix estimate is positive definite and well-conditioned while preserving the achieved sparsity. Compared to alternative estimators, this new regularization method demonstrates strong performance in simulation experiments and real portfolio optimization.

{\bf JEL classification:} C12; C15; C58; G11

{\bf Keywords:} Regularization; Multiple testing; Sign-based tests; Generalized familywise error rate; False discovery proportion

{

{This paper is forthcoming in the Journal of Econometrics:}

{Luger, R. (2024). Regularizing stock return covariance matrices via multiple testing of correlations. {\it Journal of Econometrics}, https://doi.org/10.1016/j.jeconom.2024.105753. } }

\thispagestyle{empty}

\pagenumbering{arabic}

\setcounter{footnote}{0}

Introduction

Estimating covariance matrices is a fundamental problem in multivariate statistical analysis, with wide-ranging applications in fields such as finance, economics, meteorology, climate research, spectroscopy, signal processing, pattern recognition, and genomics. In the realm of finance, accurate covariance matrix estimates are essential for capturing dependencies between asset returns -- a crucial input for portfolio optimization and risk management. Furthermore, many other statistical techniques rely on covariance matrix estimates, including regression analysis, discriminant analysis, principal component analysis, and canonical correlation analysis

Traditional estimation of the sample covariance matrix is known to perform poorly when the number of variables, $N$, is large compared to the number of observations, $T$. As the concentration ratio $N/T$ grows, there are simply too many parameters relative to the available data points and the eigenstructure of the sample covariance matrix gets distorted in the sense that the sample eigenvalues are more spread out than the population ones Johnstone:2001. The most egregious case occurs as $N/T > 1$, which causes the sample covariance matrix to become singular (non-invertible). By continuity, this matrix becomes ill-conditioned (i.e., its inverse incurs large estimation errors) as $N$ gets closer to $T$.

In such situations it is desirable to find alternative estimates that are more accurate and better conditioned than the sample covariance matrix. Regularization methods for large covariance matrices can be divided into two broad categories: (i) methods that aim to improve efficiency and obtain well-conditioned matrices, and (ii) methods that introduce sparsity (off-diagonal zeros) by imposing special structures on the covariance matrix or its inverse (the precision matrix). The first group includes linear shrinkage Ledoit-Wolf:2003, Ledoit-Wolf:2004, non-linear shrinkage Ledoit-Wolf:2012, condition-number regularization Won-Lim-Kim-Rajaratnam:2013, and split-sample regularization Abadir-Distaso-Zikes:2014. Methods that impose special structure include banding or tapering Bickel-Levina:2008b, Wu-Pourahmadi:2009 and thresholding Bickel-Levina:2008a, ElKaroui2008, Rothman-Levina-Zhu:2009, Cai-Liu:2011, which involves setting to zero the off-diagonal entries of the covariance matrix that are in absolute value below a certain data-dependent threshold.

\citet*{Bailey-Pesaran-Smith:2019}, hereafter BPS, develop an alternative thresholding approach using a multiple hypothesis testing procedure to assess the statistical significance of the elements of the sample correlation matrix; see also ElKaroui2008 who suggests a similar approach. The idea is to test all pairwise correlations simultaneously, and then to set to zero the elements that are not statistically significant. As with other thresholding methods, this multiple testing approach preserves the symmetry of the correlation matrix but it does not ensure its positive definiteness. BPS resolve this issue with an additional linear shrinkage step, whereby the correlation matrix estimator is shrunk towards the identity matrix to ensure positive definiteness. It must be emphasized that the BPS approach for reducing the number of spurious correlations is also of interest in the classical “{low $N,$ large $T$}" setting.

The simultaneous testing of all pairwise correlations gives rise to a multiple comparisons problem. Indeed if the multiplicity of inferences is not taken into account, then the probability that some of the true null hypotheses (of zero pairwise correlation) are rejected by chance alone may be unduly large. A traditional objective in multiple testing is to control the familywise error rate (FWER), defined as the probability of rejecting at least one true null hypothesis. BPS use ideas from the multiple testing literature, but from the get-go they state in their introduction (p. 508) that they “will not be particularly concerned with controlling the overall size of the joint $N(N-1)/2$ tests” of zero pairwise correlations. The simulation evidence presented in this paper reveals that the empirical FWER with BPS thresholding can be severely inflated, resulting in far too many erroneous rejections of the null hypothesis of zero correlation. This over-rejection problem is greatly exacerbated by the presence of heavy tails, which obviously defeats the purpose of achieving sparsity.

Resampling techniques can be exploited to control the flood of {Type I} errors that arise when many hypothesis tests are performed simultaneously. In particular, such techniques can be used to account for the joint distribution of the test statistics and obtain multiple testing procedures which achieve control of the FWER and other false positive error rate measures; see Westfall-Young:1993, Romano-Wolf:2005, and Romano-Shaikh-Wolf:2008. The primary goal of this paper is to extend the BPS multiple testing regularization approach so that it is applicable to financial stock returns. Maintaining proper control of the FWER in this context means that more spurious correlations are detected and greater sparsity is induced. More specifically, this paper makes two main contributions.

{\it First,} a sign-based Monte Carlo resampling technique Dwass:1957, Barnard:1963, Birnbaum:1974 is proposed to test pairwise correlations among stock returns. The theory in BPS rules out the possibility of time-varying conditional variances and covariances -- a well-known feature of financial returns Cont:2001. In turn, the presence of such effects gives rise to heavy tails and potential outliers in the distribution of returns. The procedures in this paper are developed in a general framework that allows for the presence of heavy tails and multivariate GARCH-type effects of unknown form. Indeed, the Monte Carlo resampling scheme proceeds conditional on the absolute values of the centered returns, since only their signs are randomized. The Lehmann-Stein:1949 impossibility theorem shows that such sign-based tests are the only ones that yield valid inference in the presence of non-normalities and heteroskedasticity of unknown form; see Dufour:2003 for more on this point.

{\it Second,} the resampling approach is used to develop both single-step and step-down test procedures that ensure control of two error rate measures advocated by Lehmann-Romano:2005. The first of these measures is the $k$-FWER, defined as the probability of committing $k$ or more {Type I} errors, which are commonly referred to as false positives or false discoveries. Setting $k=1$ yields procedures controlling the traditional FWER used in confirmatory research, where the goal is to test a set of hypotheses while rigorously controlling the probability of making at least one false discovery. Simulation evidence reveals that the power of the step-down $1$-FWER test procedure is on par with that of the (FWER-adjusted) BPS tests. The second error rate measure is the false discovery proportion (FDP), which represents the proportion of rejected null hypotheses that are erroneously rejected. Using the $k$-FWER (with $k > 1$) and the FDP as multiple testing criteria allows for a more lenient control over false rejections, thereby enhancing the ability to detect false null hypotheses. Unlike the traditional FWER, the $k$-FWER and the FDP offer a more nuanced approach to hypothesis testing in exploratory research. In these contexts, rejected hypotheses are generally not meant to be reported as end results, but are to be followed up with further validation experiments, such as a subsequent out-of-sample performance evaluation; see Goeman-Solari:2011.

Of course, the benefit of conducting multiple testing of correlations in terms of covariance matrix regularization is expected to increase with the true degree of sparsity in the population covariance matrix. In the context of stock portfolios, it is worth noting that this regularization approach provides a measure of diversification when returns tend to be positively correlated. Indeed the induced sparsity is expected to be proportional to the level of portfolio diversification, since in this case a well-diversified stock portfolio is precisely one in which the constituent assets demonstrate little or no correlation.

The rest of this paper is organized as follows. Section 2 presents the BPS multiple testing regularization approach. Section 3 establishes the financial context and Section 4 develops the multiple testing procedures. Section 5 presents the results of simulation experiments that compare the performance of the new regularization method to BPS and other covariance matrix estimators. Section 6 further illustrates the large-scale inference approach with an out-of-sample evaluation of portfolio optimization strategies. Section 7 offers some concluding remarks. All proofs, additional numerical results, and computing times for the multiple testing procedures are provided in the Supplementary material.

Multiple testing regularization

Consider a sample covariance matrix $\hat {\boldsymbol \Sigma} = [ \hat \sigma_{ij} ]_{N \times N}$ based a data sample of size $T$, and let $\hat {\boldsymbol \Gamma} = [ \hat \rho_{ij} ]_{N \times N}$ denote the corresponding correlation matrix with typical element $ \hat \rho_{ij}= \hat{\sigma}_{ij}/\sqrt{\hat{\sigma}_{ii} \hat{\sigma}_{jj} } $. As usual the sample covariance and correlation matrices are related via $\hat {\boldsymbol \Gamma} = \hat {\mathbf D}^{-1/2} \hat {\boldsymbol \Sigma} \hat {\mathbf D}^{-1/2}$, where $\hat {\mathbf D} = \text{diag}(\hat \sigma^2_{1},\ldots,\hat \sigma^2_{N} )$ with $\hat \sigma^2_{i}=\hat \sigma_{ii}$. The BPS regularization strategy aims to improve $\hat {\boldsymbol \Sigma} $ by testing the family of $M= N(N-1)/2$ individual hypotheses $ H_{ij}$ in the two-sided setting

for $j=1,\ldots,N-1$ and $i=j+1,\ldots,N$, while controlling the FWER. The elements that are found to be statistically insignificant are then set to zero. Instead of covariances, BPS prefer to base inference on the sample correlations since they are all on the same scale. This leads to multiple testing procedures that are balanced in the sense that all constituent tests have about the same power Westfall-Young:1993. Note that the entries of the sample correlation matrix are intrinsically dependent even if the original observations are independent.

There are two types of FWER control. To introduce these, define an index $\ell$ taking values in the set $\mathcal M =\{1,\ldots, M \}$ as $j=1,\ldots,N-1$ and $i=j+1,\ldots,N$ so that $ H_{1} = H_{2,1}, \ldots, H_{N-1} = H_{N,1},H_{N} = H_{3,2}, \ldots, H_{M} = H_{N,N-1} $. Of these $M$ hypotheses tested, let $R$ denote the number of hypotheses rejected. Furthermore, let $\mathcal M_0 = \{\ell : H_{\ell} \text{ is true} \}$ denote the index set of true hypotheses. The number of false positive decisions (i.e., the number of {Type I} errors) is denoted by $F$. Given the nominal significance level $\alpha \in (0,1)$, the FWER is said to be controlled in the weak sense when \sloppy $\Pr \left( F \ge 1 \, \vert \, \bigcap_{\ell \in \mathcal M} H_{\ell} \right) = \Pr \left( \text{Reject at least one } H_{\ell} \, \vert \, \bigcap_{\ell \in \mathcal M} H_{s} \right) \le \alpha$, where the conditioning is on the complete null hypothesis that $\mathcal M_0 =\mathcal M$ or $|\mathcal M_0| = M$.\footnote{The notation $\Pr(\text{event} \, \vert \, H)$ refers to the probability of the event occurring when $H$ holds true.} Strong control is achieved when $\Pr \left( F \ge 1 \, \vert \, \bigcap_{\ell \in \mathcal M_0} H_{\ell} \right) = \Pr \left( \text{Reject at least one } H_{\ell}, \, \ell \in \mathcal M_0 \, \vert \, \bigcap_{\ell \in \mathcal M_0} H_{\ell} \right) \le \alpha$, regardless of the partial null hypothesis (i.e., the particular intersection of hypotheses that happens to be true).

The BPS thresholding estimator, denoted here by $\hat {\boldsymbol \Gamma}_{\text{BPS}} $, has off-diagonal entries computed as

wherein $ \mathds{1}\{\cdot\}$ denotes the indicator function. The critical value $c_{\alpha}(N)$ appearing in ((ref)) is given by

where $\Phi^{-1}(\cdot)$ is the quantile function of a standard normal variate and $f(N)$ is a general function of $N$ chosen to ensure $\text{FWER} \le \alpha$ in the strong sense. Observe that the term $T^{-1/2} c_{\alpha}(N)$ in ((ref)) is a `universal' threshold value in the sense that a single value is used to threshold all the off-diagonal elements of the correlation matrix $\hat {\boldsymbol \Gamma} $. Among other asymptotic properties, BPS show that $\hat {\boldsymbol \Gamma}_{\text{BPS}} $ converges to the true $ {\boldsymbol \Gamma} $ as the sample size grows. As one expects, the payoff in terms of noise reduction with this approach increases with the actual number of zeros in $ {\boldsymbol \Gamma} $.

When the number $M$ of tested hypotheses is very large and/or when strong evidence is not required, control of the FWER at conventional levels can be too stringent. In such cases one may wish to control the $k$-FWER, with $k$ set to a value greater than $1$ Lehmann-Romano:2005. Control of the $k$-FWER can also be either weak or strong. Given $\alpha \in (0,1)$ and $k \in \{1,\ldots, M\}$, weak control occurs whenever $\Pr \left( F \ge k \, \vert \, \bigcap_{\ell \in \mathcal M} H_{\ell} \right) \le \alpha$ and strong control is achieved when $\Pr \left( F \ge k \, \vert \, \bigcap_{\ell \in \mathcal M_0} H_{\ell} \right) \le \alpha$.

The choice of $k$ should be guided by the research objectives and the balance between minimizing the risk of false positives and maximizing the discovery of true effects. If the primary objective is to maintain strict control over {Type I} errors, then a constant, relatively low value of $k$ might be chosen regardless of the number of tested hypotheses. This ensures a consistent level of control over false positives. Conversely, when there is room for a controlled number of {Type I} errors and a desire to uncover more findings, a higher value of $k$ might be selected, for instance by letting $k$ increase with the number of tested hypotheses.

Another possibility considered in Lehmann-Romano:2005 is to maintain control of the FDP, defined as

The FDP thus represents the number of false rejections relative to the total number of rejections (and equals zero if there are no rejections at all). FDP is a useful criterion when the focus is on controlling the rate at which null hypotheses are incorrectly rejected. Note that the FDP in ((ref)) is a random variable, and FDP control focuses on its tail probabilities. For given values of $\gamma \in [0,1)$ and $\alpha \in (0,1)$ specified by the user, strong control (of the tail probability) of the FDP means that $\Pr \left( \text{FDP} > \gamma \, \vert \, \bigcap_{\ell \in \mathcal M_0} H_{\ell} \right) \le \alpha$. The interpretation is that, among the hypotheses that are rejected, the proportion of false discoveries may exceed $\gamma$ with probability no larger than $\alpha.$

Control of the FDP is a reasonable practice especially in non-confirmatory settings, where a certain proportion of false discoveries is considered acceptable. The choice of $\gamma$, which determines the acceptable rate of false discoveries, involves a trade-off. A lower $\gamma$ reduces the risk of false discoveries but potentially missing some true effects. Indeed, setting $\gamma=0$ is equivalent to controlling the conservative $1$-FWER criterion. Conversely, increasing the value of $\gamma$ allows for more discoveries but at the cost of accepting a higher rate of false discoveries.

Positive definiteness

As with other thresholding methods, the matrix $\hat {\boldsymbol \Gamma}_{\text{BPS}}$ obtained via ((ref)) is not necessarily well-conditioned nor even positive definite. BPS solve this problem by shrinking $\hat {\boldsymbol \Gamma}_{\text{BPS}}$ towards $ \mathbf I_N$, the $N \times N$ identity matrix. Let $\lambda_{\min}(\hat {\boldsymbol \Gamma}_{\text{BPS}})$ denote the minimum eigenvalue of $\hat {\boldsymbol \Gamma}_{\text{BPS}}$ and set a limit $\epsilon > 0$ to avoid solutions that are too close to being singular.\footnote{In the simulation experiments and empirical application, the limit was set as $\epsilon = 0.01.$} The `shrinkage upon multiple testing' correlation matrix estimator is given by

with shrinkage parameter $\xi \in [\xi_0,1]$, where $\xi_0 = \big(\epsilon- \lambda_{\min}(\hat {\boldsymbol \Gamma}_{\text{BPS}})\big)/\big(1-\lambda_{\min}(\hat {\boldsymbol \Gamma}_{\text{BPS}})\big)$ if $\lambda_{\min}(\hat {\boldsymbol \Gamma}_{\text{BPS}}) < \epsilon$, and $\xi_0 = 0$ if $\lambda_{\min}(\hat {\boldsymbol \Gamma}_{\text{BPS}}) \ge \epsilon$. Note that by shrinking towards the identity matrix, the resulting correlation matrix estimate preserves the zeros (sparsity) achieved by $\hat {\boldsymbol \Gamma}_{\text{BPS}}$ and its diagonal elements do not deviate from unity. The computation of ((ref)) is made operational by replacing $\xi$ by $\xi^*$, which is found numerically as \[ \xi^* = \arg \min_{\xi_0 \le \xi \le 1} \left\lVert \hat{\boldsymbol \Gamma}^{-1}_0 - \hat {\boldsymbol \Gamma}_{\text{BPS}}^{-1}(\xi) \right\rVert^2_F, \] where $\hat{\boldsymbol \Gamma}_0$ is a reference matrix and $\lVert \mathbf A \rVert_F$ denotes the Frobenius norm of $\mathbf A $.\footnote{Recall that for a matrix $\mathbf A$ with elements $a_{ij}$, its Frobenius norm is defined as $\lVert \mathbf A \rVert_F = \sqrt{ \sum_i \sum_j a_{ij}^2 }$. In the implementation, the value of $\xi^*$ was found by grid search with a step size of $\epsilon/2$.}

Following Schafer-Strimmer:2005, the reference matrix $\hat{\boldsymbol \Gamma}_0$ is found by applying the linear shrinkage approach of Ledoit-Wolf:2004 to the sample correlation matrix, which yields \[ \hat{\boldsymbol \Gamma}_0 = \hat \theta^* \mathbf I_N + (1- \hat \theta^*) \hat {\boldsymbol \Gamma}, \] where \[ \hat \theta^* = 1- \frac{\underset{i \neq j}{\sum \sum} \hat \rho_{ij} \left[ \hat \rho_{ij} - \frac{\hat \rho_{ij} (1-\hat \rho_{ij}^2 )}{2T} \right] }{\frac{1}{T} \underset{i \neq j}{\sum \sum} (1- \hat \rho_{ij}^2 )^2 + \underset{i \neq j}{\sum \sum}\left[ \hat \rho_{ij} - \frac{\hat \rho_{ij} (1-\hat \rho_{ij}^2 )}{2T} \right]^2}, \] with the proviso that if $\hat \theta^* < 0$ then $\hat \theta^*$ is set to $0$, and if $\hat \theta^* > 1$ then it is set to $1$.\footnote{ The analytical expression for $\hat \theta^*$ is an estimate of the optimal value of the shrinkage parameter that minimizes $\mathbb E (\lVert \hat{\boldsymbol \Gamma}_0 - {\boldsymbol \Gamma} \rVert^2_F)$, assuming $\mathbb E(\hat \rho_{ij})$ can be approximated by $\hat \rho_{ij} - \hat \rho_{ij}(1-\hat \rho^2_{ij})/2T$ and $\text{Var}(\hat \rho_{ij})$ by $(1-\hat \rho_{ij}^2)^2/T$ Soper-Young-Cave-Lee-Pearson:1917. As Ledoit-Wolf:2003 explain, most shrinkage covariance matrix estimators are based on such first-cut assumptions. } The resulting covariance matrix estimate is given by $\hat {\boldsymbol \Sigma}_{\text{BPS}}(\xi^*) = \hat {\mathbf D}^{1/2} \hat {\boldsymbol \Gamma}_{\text{BPS}}(\xi^*) \hat {\mathbf D}^{1/2},$ wherein $ \hat {\boldsymbol \Gamma}_{\text{BPS}}(\xi^*)$ corresponds to ((ref)) evaluated with $\xi^*$.

Financial context

Consider $N$ financial assets with time-$t$ returns $\mathbf r_t =(r_{1,t},\ldots,r_{N,t})^{\prime}$, for $t=1,\ldots,T$, decomposed as

where $\boldsymbol \mu = \mathbb E(\mathbf r_t) = \big(\mu_{1}, \ldots, \mu_{N}\big)^{\prime}$ is a vector of location parameters. The error $ \boldsymbol \varepsilon_t = (\varepsilon_{1,t},\ldots,\varepsilon_{N,t})^{\prime}$ in ((ref)) consists of an innovation vector $\mathbf z_t$ satisfying Assumption (ref) below, and an unspecified $N \times N$ “square root" matrix $\boldsymbol \Sigma_t^{1/2} $ such that $\boldsymbol \Sigma_t^{1/2} (\boldsymbol \Sigma_t^{1/2})^{\prime} = \boldsymbol \Sigma_t = \mathbb E(\boldsymbol \varepsilon_t \boldsymbol \varepsilon_t^{\prime} \, \vert \, \mathcal F_{t-1})$, where $\mathcal F_{t-1} = (\mathbf r_{t-1}^{\prime},\mathbf r_{t-2}^{\prime}, \ldots)^{\prime}$. This framework is compatible with several popular models of time-varying covariances, such as multivariate GARCH models Silvennoinen-Terasvirta:2009, Boudt-Galanos-Payseur-Zivot:2019 and multivariate stochastic volatility models Chib-Omori-Asai:2009.

This assumption means that $\mathbf z_t$ admits the stochastic representation $\mathbf z_t \stackrel{d}{=} \mathbf H \mathbf z_t $, where the symbol $\stackrel{d}{=}$ stands for an equality in distribution and $\mathbf H $ is any $N \times N$ orthogonal matrix such that $\mathbf H^{\prime} \mathbf H = \mathbf H \mathbf H^{\prime} = \mathbf I_N$. This class includes the multivariate versions of the standardized normal, Student $t$, logistic, and {Laplace} distributions, among many others Fang-Kotz-Ng:1990.

When $\mathbf z_t$ has a well-defined density, then Assumption (ref) is equivalent to assuming that the conditional distribution of $\mathbf r_t $ is elliptically symmetric, meaning that its density has the form $|\boldsymbol \Sigma_t^{-1/2}| f\big( (\mathbf r_t - \boldsymbol \mu )^{\prime} \boldsymbol \Sigma_t^{-1} (\mathbf r_t - \boldsymbol \mu ) \big)$ for some non-negative scalar function $f(\cdot)$. Elliptically symmetric distributions play a very important role in mean-variance analysis (cf. Section 6) because they guarantee full compatibility with expected utility maximization regardless of investor preferences Berk:1997, Chamberlain:1983, Owen-Rabinovitch:1983.

In the context of ((ref)), the complete null hypothesis is formally stated as

for $t=1,\ldots,T$, where $\mathbf D_t^{1/2}$ is a diagonal matrix (i.e., with zeros outside the main diagonal). Observe that conditional heteroskedasticity is permitted under $H_0$; i.e., the diagonal elements of $\mathbf D_t^{1/2}$ may be time-varying. It is easy to see that when Assumption (ref) holds and $H_0$ is true, the error vector $\boldsymbol \varepsilon_t =(\varepsilon_{1,t},\ldots,\varepsilon_{N,t})^{\prime}$ becomes sign-symmetric Serfling:2006 in the sense that \[ \boldsymbol \varepsilon_t \stackrel{d}{=} \mathbb S_t \boldsymbol \varepsilon_t, \] for $t=1,\ldots,T$ and for all $N \times N$ diagonal matrices $ \mathbb S_t $ with $\pm 1$ on the diagonal.

With this assumption the sign-symmetry condition $(\varepsilon_{i,t}, \varepsilon_{j,t}) \stackrel{d}{=} (\pm \varepsilon_{i,t}, \pm \varepsilon_{j,t}) $ implies $\sigma_{ij}=0$, for $i \neq j$, where $\sigma_{ij}$ is the $(i,j)^\text{th}$ element of $\boldsymbol \Sigma$ Randles-Wolfe:1979.

Multiple testing procedures

Assume momentarily that the true value ${\boldsymbol \mu}^{\ast}$ of $\boldsymbol \mu$ in ((ref)) is known. For instance, with daily returns it is often reasonable to assume that ${\boldsymbol \mu}^{\ast} = \boldsymbol 0$. The case of unknown location parameters will be dealt with in Section 4.6.

Given the value of ${\boldsymbol \mu}^{\ast}$, centered returns can then be defined as $\mathbf y_t = \mathbf r_t - {\boldsymbol \mu}^{\ast} = (y_{1,t},\ldots,y_{N,t})^{\prime}$, for $t=1,\ldots,T$, and these have the same properties as $ \boldsymbol \varepsilon_t$. The time series of centered returns are collected into the $T \times N$ matrix $\mathbf Y = [\mathbf y_1,\ldots,\mathbf y_T ]^{\prime}$. Following BPS, inference is based on the pairwise correlations $ \hat \rho_{ij} = { \hat \sigma_{ij} }/{\sqrt{\hat \sigma_{ii} \hat \sigma_{jj}} } $ that constitute the matrix $\hat {\boldsymbol \Gamma} = [ \hat \rho_{ij} ]_{N \times N}$. This matrix can be obtained from the familiar relationship $\hat {\boldsymbol \Gamma} = \hat {\mathbf D}^{-1/2} \hat {\boldsymbol \Sigma} \hat {\mathbf D}^{-1/2}$ with $\hat {\mathbf D} = \text{diag}(\hat \sigma^2_{1},\ldots,\hat \sigma^2_{N} )$ and $\hat \sigma^2_{i}=\hat \sigma_{ii}$, where $\hat {\boldsymbol \Sigma}$ now comprises the variances and covariances about the origin, computed as $ \hat \sigma_{ij} = T^{-1} \sum_{t=1}^T y_{i,t} y_{j,t}$, for $i,j=1,\ldots,N$.

Let $\tilde{ \mathbb S}_t = \text{diag}(\tilde s_{1,t},\ldots,\tilde s_{N,t})$, for $t=1,\ldots,T$, where $\tilde s_{i,t} $ are independent Rademacher random draws such that $\Pr(\tilde s_{i,t} =1) = \Pr(\tilde s_{i,t} =-1 ) =1/2$, for each $i, t$. An artificial sample $\tilde{\mathbf Y} = \left[\tilde{\mathbf y}_1,\ldots,\tilde{\mathbf y}_T \right]^{\prime}$ with $\tilde { \mathbf y}_t = ( \tilde y_{1,t},\ldots, \tilde y_{N,t})'$ is then defined as

If Assumption (ref) holds and $H_0$ is true, then $\mathbf Y \stackrel{d}{=} \tilde{\mathbf Y} $, for each of the $2^{TN}$ possible matrix realizations of $ \tilde{\mathbf Y} $, given $|\mathbf Y |$. Here $|\mathbf Y |$ is the matrix of entrywise absolute values of $\mathbf Y$. For a given artificial sample $\tilde { \mathbf Y}$, let $\tilde {\boldsymbol \Gamma} = [ \tilde \rho_{ij} ]_{N \times N}$ denote the associated correlation matrix comprising the pairwise correlations about the origin $\tilde \rho_{ij} = { \tilde \sigma_{ij} }/{\sqrt{\tilde \sigma_{ii} \tilde \sigma_{jj}} }$, where $\tilde \sigma_{ij} = {T}^{-1} \sum_{t=1}^T \tilde y_{i,t} \tilde y_{j,t}$.

Proposition (ref) shows that $\hat {\boldsymbol \Gamma} $ is conditionally pivotal under $H_0$, meaning that its sign-randomization distribution does not depend on any nuisance parameters. In principle, critical values could be found from the conditional distribution of $\hat {\boldsymbol \Gamma} $ derived from the $2^{TN}$ equally likely values represented by $ \tilde {\boldsymbol \Gamma} $. Determination of this distribution from a complete enumeration of all possible realizations of $\tilde {\boldsymbol \Gamma} $ is obviously impractical. Following Zhu-Neuhaus:2000, the algorithms developed next use a non-parametric Monte Carlo test technique in order to circumvent this problem and still obtain exact $p$-values.

Unadjusted $p$-values

It is useful to first describe how to obtain the Monte Carlo $p$-values without multiplicity adjustments, even if they are not used directly for multiple testing regularization. Indeed, these raw $p$-values are the foundational blocks for the development of multiplicity-adjusted $p$-values. Some additional notation will facilitate the explanation of the Monte Carlo test technique. With a correlation matrix ${\boldsymbol \Gamma} = [\rho_{ij}]_{N \times N}$ as input, let $\boldsymbol \rho = \text{vechs}({\boldsymbol \Gamma} ) = (\rho_{2,1},\ldots,\rho_{N,1}, \rho_{3,2}, \ldots, \rho_{N,N-1})^{\prime}$ be the $M \times 1$ vector resulting from its strict half-vectorization and let $\text{fill}(\cdot)$ denote the inverse function such that ${\boldsymbol \Gamma} = \text{fill}(\text{vechs}({\boldsymbol \Gamma} ))$. Furthermore, let $ (\hat \rho_{1},\ldots, \hat \rho_{M})^{\prime} = \text{vechs}(\hat {\boldsymbol \Gamma} )$ so that $|\hat \rho_{\ell}|$ is the statistic for $H_{\ell}$, $\ell=1,\ldots,M$.

Note that sampling according to ((ref)) yields a discrete distribution of $\tilde {\boldsymbol \Gamma}$ values, which means that ties among the resampled values can occur, at least theoretically. Following Dufour:2006, these are dealt with by working with lexicographic (tie-breaking) ranks. Algorithm (ref) details the steps to obtain the unadjusted Monte Carlo $p$-values with $B-1$ resampling draws chosen so that $\alpha B$ is an integer, where $\alpha \in (0,1)$ is the desired significance level.

The Monte Carlo $p$-values obtained from Algorithm (ref) have the usual interpretation: $ \tilde p_\text{U}(| \hat \rho_{\ell} | )$ is the proportion of $| \tilde \rho_{\ell} |$ values as extreme or more extreme than the observed $| \hat \rho_{\ell} |$ value in its resampling distribution.

Single-step adjusted $p$-values for $k$-FWER control

Westfall-Young:1993 propose several resampling-based methods to adjust $p$-values so as to account for multiplicity. Adjusted $p$-values are defined as the smallest significance level for which one still rejects an individual hypothesis $H_{\ell}$, given a particular multiple test procedure. Let $\text{$k$-max}(|\boldsymbol \rho|)$ denote the $k^\text{th}$ largest value of $|\boldsymbol \rho|$. So if the elements of $|\boldsymbol \rho|$ are sorted in decreasing order as $|\rho_{(1)}| \ge |\rho_{(2)}| \ge \ldots \ge |\rho_{(M)}|$, then $\text{$k$-max}(|\boldsymbol \rho|)=|\rho_{(k)}|$.

A straightforward extension of Westfall and Young's single-step (SS) maxT adjusted p-values to the present context for $k$-FWER control yields the definition

for $\ell=1,\ldots,M$, where $\tilde{\boldsymbol \rho} = \text{vechs}(\tilde {\boldsymbol \Gamma} )$ and $H_0$ is the complete null hypothesis in ((ref)). When $k=1$, ((ref)) reduces to Definition (2.8) in Westfall-Young:1993. In words, this says that the SS adjusted $p$-value is the probability that the $k^\text{th}$ largest absolute correlation in the artificial data is greater than the observed absolute correlation in the actual data.

Lemma (ref) paves the way for the computation of the Monte Carlo version of ((ref)) as described next. The number of resampling draws $B-1$ is assumed to be chosen so that $\alpha B$ is an integer, where $\alpha \in (0,1)$ is the desired $k$-FWER.

The proof in Westfall-Young:1993 that their SS adjusted $p$-values control the FWER in the strong sense relies heavily on the assumption of subset pivotality. That is, they assume that the joint distribution of unadjusted $p$-values under any partial null hypothesis is identical to that under the complete null hypothesis. As noted by Westfall-Young:1993 and Romano-Wolf:2005, this assumption fails in the context of testing pairwise correlations. However, as the next proposition shows, subset pivotality is not a necessary condition for strong control; see also Romano-Wolf:2005, Westfall-Troendle:2008, and Goeman-Solari:2010.

The proof of this result (given in the Supplementary material) makes clear that the SS procedure becomes more conservative as the cardinality of $\mathcal M_0$ decreases.

Step-down adjusted $p$-values for $k$-FWER control

A disconcerting feature of ((ref)) is that all the $p$-values are adjusted according to the distribution of the $k^\text{th}$ largest absolute correlation. Potentially less conservative $p$-values may be obtained from step-down adjustments that result in uniformly smaller $p$-values, while retaining the same protection against {Type I} errors. The underlying idea is analogous to Holm:1979's (Holm:1979) sequential refinement of the Bonferroni adjustment, which eliminates from consideration any null hypotheses that are rejected at a previous step.\footnote{Romano-Wolf:2005 discuss the idealized step-down method; see also Romano-Wolf:2016. That method is not feasible with the resampling scheme developed here, because it is not possible to generate artificial data that obey the null hypothesis for each possible intersection of true null hypotheses.}

With the absolute correlations $|\hat \rho_{1}|,\ldots,|\hat \rho_{M}|$, let the ordered test statistics have index values $\pi_1,\ldots,\pi_M$ so that $|\hat \rho_{\pi_1}| \ge |\hat \rho_{\pi_2}| \ge \ldots \ge |\hat \rho_{\pi_M}|$. To control the $k$-FWER, the definition of Westfall and Young's step-down (SD) maxT adjusted p-values can be extended as follows:

wherein the sequence of index values $\pi_1,\ldots,\pi_M$ is held fixed.\footnote{That is, the adjustments are made by “stepping down” from the largest test statistic to the smallest.}

Instead of adjusting all $p$-values according to the distribution of the $k^\text{th}$ largest absolute correlation, this approach only adjusts the $p$-values for $|\hat \rho_{\pi_1}|, |\hat \rho_{\pi_2}|, \ldots, |\hat \rho_{\pi_k}|$ using this distribution. The remaining $p$-values in steps $k+1,\ldots, M$ are then adjusted according to the distributions of smaller and smaller sets of maximum absolute correlations, where the $\min$ operator is used to ensure that $ \tilde m^k_{\ell} \le \tilde m^k_{\ell-1} \le \text{$k$-max}\bigl(|\tilde{ \boldsymbol \rho}|\bigr),$ for $\ell=k+1,\ldots, M.$

Note that the SD adjusted $p$-values have the same step-down monotonicity as the original test statistics; i.e., smaller $p$-values are associated with larger values of the $|\hat \rho_{\ell}|$ test statistics. This is obvious for $p^k_{\text{SD}}(|\hat \rho_{\pi_1}|), \ldots, p^k_{\text{SD}}(|\hat \rho_{\pi_k}|)$, since these $p$-values follow a single-step adjustment. For $\ell=k+1,\ldots, M$, it is the application of the $\max$ operator (outside the square brackets) at each of those steps that guarantees the remaining step-down monotonicity. This approach can yield power improvements since the SD adjusted $p$-values are uniformly no larger than their SS counterparts. When $k=1$, Definition ((ref)) simplifies to Westfall and Young's definition of SD adjusted $p$-values Ge-Dudoit-Speed:2003.

Extending Westfall-Young:1993 and Ge-Dudoit-Speed:2003, Algorithm (ref) below shows how to compute the Monte Carlo version of the SD $p$-values defined in ((ref)). Here again the number of resampling draws $B-1$ is assumed to be chosen so that $\alpha B$ is an integer, where $\alpha \in (0,1)$ is the user's desired $k$-FWER.

Adjusted $p$-values for FDP control

Following Romano-Wolf:2007, control of the FDP in ((ref)) can be achieved by sequentially applying a $k$-FWER-controlling procedure, for $k=1,2,\ldots$, until a stopping rule indicates termination. Whether Algorithm (ref) or Algorithm (ref) is employed, it is important that the same underlying $\tilde{\mathbf Y}_b$, $b=1,\ldots,B-1$, and $u_b$, $b=1,\ldots,B$, be used for each $k$ to ensure coherence of the resulting $p$-values. This can be done simply by resetting the seed of the random number generator to the same value each time $k$ is changed, as illustrated in Algorithm (ref) below.

The user first specifies the FDP exceedance threshold $\gamma \in [0,1)$ and chooses $B$ so that $\alpha B$ is an integer, where $\alpha \in (0,1)$ is the desired probability level. In the following, $ \tilde p^{k}_{\diamond}(| \hat \rho_{\ell} |) $ represents either $ \tilde p^k_\text{SS}(| \hat \rho_{\ell} |)$ or $ \tilde p^k_\text{SD}(| \hat \rho_{\ell} |) $ computed according to Algorithms (ref) and (ref), respectively. Let {\tt seed} be a value that will be used to set the seed of the random number generator.\footnote{In the R programming language this can be done with the command {\tt set.seed({\tt seed})} where {\tt seed}=8032, for example.}

For large $M$, Algorithm (ref) can be exceedingly slow, since it finds $k^*$ by progressing sequentially starting with $k=1$, then $2$, and so on. The monotonicity of the underlying $k$-FWER-controlling procedure can be further exploited to achieve important speed gains via a bisection method, as follows.

The bisection method, known for its numerical robustness, is particularly well-suited to handle the discontinuities within the underlying step function $R_k$. Steps 1 and 2 isolate the largest value $k_l$ such that $1 \le k_l \le k^* \le M$ and $\gamma(R_{k_l} + 1) = \gamma(R_{k^{\ast}} + 1)$. This is achieved by successively narrowing down the interval in which $k_l$ lies through a process of halving at each iteration of Step 2. The call in Step 3 to Algorithm (ref), initialized with $k_l$, then returns $\tilde p^{\gamma}_{\diamond}(| \hat \rho_{\ell} |)$, $\ell=1,\ldots,M,$ with at most two iterations of the while loop in that algorithm.

Covariance matrix estimators

The next step in the construction of the proposed multiple-testing regularized estimators is to set to zero the statistically insignificant entries of the sample correlation matrix $\hat {\boldsymbol \Gamma} = [ \hat \rho_{ij} ]_{N \times N}$, defined previously. Let $ \tilde p^{\diamond}_{\diamond}(| \hat \rho_{ij} |) $ represent either $ \tilde p^k_\text{SS}(| \hat \rho_{ij} |)$, $ \tilde p^k_\text{SD}(| \hat \rho_{ij} |) $, $ \tilde p^\gamma_\text{SS}(| \hat \rho_{ij} |)$, or $ \tilde p^\gamma_\text{SD}(| \hat \rho_{ij} |) $. A correlation $p$-value matrix corresponding directly to $\hat {\boldsymbol \Gamma} $ is given by $\text{fill}\big(\tilde p^{\diamond}_{\diamond}(| \hat \rho_{1} |), \ldots, \tilde p^{\diamond}_{\diamond}(| \hat \rho_{M} |) \big)$ with the diagonal elements set to zero.\footnote{These zero diagonal values are irrelevant since the diagonal elements of ${\boldsymbol \Sigma}$ are not tested.}

The associated correlation matrix estimator $\hat {\boldsymbol \Gamma}^{\diamond}_{\diamond} = [ \hat \rho^{\diamond}_{\diamond, ij} ]_{N \times N}$ has entries given by

These adjustments to the sample correlation matrix are made for $j=1,\ldots,N-1$ and $i=j+1,\ldots,N$; the diagonal elements of $\hat {\boldsymbol \Gamma}^{\diamond}_{\diamond} $ are obviously set as $\hat \rho^{\diamond}_{\diamond,ii} =1$; and symmetry is imposed by setting $\hat \rho^{\diamond}_{\diamond,ij} = \hat \rho^{\diamond}_{\diamond,ji}.$ In contrast to the universal threshold critical value used in ((ref)) for each of the pairwise correlations, note that the $p$-values used in ((ref)) are fully data-driven and adapt to the variability of individual entries of the correlation matrix.

Proceeding to the shrinkage step in ((ref)) of the BPS approach with $\hat {\boldsymbol \Gamma}^{\diamond}_{\diamond} $ instead of $\hat {\boldsymbol \Gamma}_{\text{BPS}} $ yields the positive definite correlation matrix estimator $\hat {\boldsymbol \Gamma}^{\diamond}_{\diamond}(\xi^*) = \xi^* \mathbf I_N + (1-\xi^*) \hat {\boldsymbol \Gamma}^{\diamond}_{\diamond},$ where $\xi^* = \arg \min_{\xi_0 \le \xi \le 1} \left\lVert \hat{\boldsymbol \Gamma}^{-1}_0 - \hat {\boldsymbol \Gamma}_{\diamond}^{\diamond^{ -1}}(\xi) \right\rVert^2_F, $ and the associated covariance matrix estimator $\hat {\boldsymbol \Sigma}^{\diamond}_{\diamond} (\xi^*) = \hat {\mathbf D}^{1/2} \hat {\boldsymbol \Gamma}^{\diamond}_{\diamond}(\xi^*) \hat {\mathbf D}^{1/2}.$

Unknown location parameters

Observe that the null distribution generated according to ((ref)) depends on ${\boldsymbol \mu}^{\ast}$ in ((ref)) only through the subtractive transformations $\mathbf y_t = \mathbf r_t - \boldsymbol \mu^{\ast}$, $t=1,\ldots,T$. When the values comprising ${\boldsymbol \mu}^{\ast}$ are unknown, it will be assumed that they can be estimated consistently. Let $ \hat {\boldsymbol \mu}_T $ denote the equation-by-equation estimator of ${\boldsymbol \mu}^{\ast}$ obtained from a sample of size $T$. For instance, this could be the sample mean. With $ \hat {\boldsymbol \mu}_T $ in hand, the Monte Carlo procedures can proceed as before except that $\hat{\mathbf y}_t = \mathbf r_t - \hat{\boldsymbol \mu}_T $ replaces $\mathbf y_t = \mathbf r_t - {\boldsymbol \mu}^{\ast}$.

Let $\hat \rho_{ij}( \hat {\boldsymbol \mu}_T )$, $1 \le j < i \le N$, denote the correlations estimated on the basis of $\hat{\mathbf Y} = \left[\hat{\mathbf y}_1,\ldots,\hat{\mathbf y}_T \right]^{\prime}$. Since $ \hat {\boldsymbol \mu}_T \stackrel{p}{\rightarrow} {\boldsymbol \mu}^{\ast}$, Slutsky's theorem ensures that $\hat{ \mathbf Y} \stackrel{d}{\rightarrow} \mathbf Y \stackrel{d}{=} \tilde{\mathbf Y}$, under $H_0$, as $T \rightarrow \infty$. From Zhu-Neuhaus:2000 and Toulis-Bean:2021, it then follows that the Monte Carlo $p$-values computed according to Algorithm (ref) are asymptotically valid in the sense that $\Pr( \tilde p_\text{U}(| \hat \rho_{ij}( \hat {\boldsymbol \mu}_T ) |) \le \alpha \, \vert \, H_0 ) = \alpha + o_p(1)$, for $1 \le j < i \le N$, as $T \rightarrow \infty$. An immediate consequence is that the adjusted $p$-values computed according to Algorithms (ref)--(ref) are also asymptotically valid, since they all rest on Algorithm (ref).

Simulation experiments

This section examines the performance of the proposed Monte Carlo regularized covariance estimators, with the BPS approach serving as the natural benchmark for comparisons. The simulation experiments are designed to resemble the empirical application presented in the next section. For this purpose the data-generating process for daily returns $\mathbf r_t=(r_{1,t},\ldots,r_{N,t})^{\prime}$ is specified as a CCC model Bollerslev:1990 of the form

where $\mathbf D_t = \text{diag}(\sigma_{1,t}^2,\ldots, \sigma_{N,t}^2)$ is an $N \times N$ diagonal matrix comprising the time-$t$ conditional variances, and $ \boldsymbol \Gamma$ is a constant conditional correlation (CCC) matrix. The vector $ \boldsymbol \mu $ is set to zero, but this is not assumed known and the multiple testing procedures are applied with $\mathbf y_t = \mathbf r_t - \bar{\mathbf r}$, where $\bar{\mathbf r}$ is the vector of sample means (cf. Section 4.6). The conditional variances appearing in $\mathbf D_t$ evolve according to the standard $\text{GARCH}(1,1)$ model $ \sigma^2_{i,t} = \theta_0 + \theta_1 r_{i,t-1}^2 + \theta_2 \sigma^2_{i, t-1} $ with common parameters across assets set as $\theta_0 = 0.01$, $\theta_1=0.1$, and $\theta_2=0.85$. These are typical values found with daily returns data. The innovation terms $ z_{i,t}$ are i.i.d. according either to the standard normal distribution, or the standardized $t$-distribution with 12 or 6 degrees of freedom.

The correlation structure is defined as follows. Given a value $\delta$, $0 \le \delta \le 1$, the vector $\mathbf c = (c_1,\ldots,c_N)^{\prime}$ is filled in with $N_c= \lfloor \delta N \rfloor$ non-zero elements drawn from the triangular distribution on $[0,1]$ with mode at $1$, and the remaining $N-N_c$ elements are set to zero. The positions of the zero and non-zero elements within $\mathbf c $ are random. Following BPS, this vector is then used to obtain a well-defined correlation matrix as $\boldsymbol \Gamma = \mathbf I_N + \mathbf c \mathbf c^{\prime} - \text{diag}( \mathbf c \mathbf c^{\prime} ).$ The complete null hypothesis in ((ref)) is thus represented by $\delta=0$ and increasing the value of $\delta$ towards 1 results in a smaller set $\mathcal M_0$ of true hypotheses. With $\boldsymbol \Gamma$ in hand, the unconditional covariance matrix is then found as ${\boldsymbol \Sigma} = \mathbf D^{1/2} \boldsymbol \Gamma \mathbf D^{1/2},$ where $\mathbf D = \text{diag}(\sigma_{1}^2,\ldots, \sigma_{N}^{2})$ comprises the GARCH-implied unconditional variances given by $\sigma_i^2=\theta_0/(1-\theta_1-\theta_2)$, for $i=1,\ldots,N.$ Vector sizes $N=25$, $100$, and $500$, and sample sizes $T=63$, $126$, and $252$ are considered.

Following BPS, two choices are used to complete the definition in ((ref)) of the universal critical value: (i) $f(N) = N^2$, which results in a more conservative adjustment; and (ii) $f(N) = N(N-1)/2$, representing the well-known Bonferroni rule. The results based on these choices are labelled $\text{BPS}_a$ and $\text{BPS}_b$, respectively.

In the application of the $k$-FWER procedures, reported as $\text{SS}_{k}$ and $\text{SD}_{k}$, two heuristic rules are considered: (i) $k=\lfloor \log{M} \rfloor$ and (ii) $k=\lfloor \sqrt{M} \rfloor$. These choices scale the allowable number of {Type I} errors with the number of hypotheses being tested, allowing for a more permissive approach to error control. The exceedance threshold in the FDP procedures, whose results are reported as $\text{SS}_{\gamma}$ and $\text{SD}_{\gamma}$, is set at $\gamma=0.1.$ The number of replications for the Monte Carlo procedures is established with $B=100$ and all the tests are conducted with $\alpha = 0.05.$

Numerical results

Table 1 reports the empirical rejection rates of the FWER procedures. Panels A and B report the empirical FWER and Panel C reports the empirical average power, defined as the average expected number of correct rejections among all false null hypotheses Bretz-Hothorn-Westfall:2010. Under the complete null hypothesis ($\delta=0$) the SS and SD test procedures have identical rejection rates, indicated on the lines labelled SS/SD in Panel A. It is seen that the SS (SD) procedure does a good job at keeping the FWER close to the desired $5\%$ value when $\delta=0$. This case provides a comparison benchmark for the FWER results in Panel B when $\delta$ is increased to $0.9$. These tests are seen in Panel B to also maintain control of the FWER in the strong sense. Comparing Panels A and B shows that the SS and SD test procedures become more conservative when the number of true null hypotheses diminishes ($\delta$ increasing), as expected from the developed theory.

The BPS method also achieves good control of the FWER under normality. And just like the SS and SD procedures, the BPS tests are seen in Panel B to also become conservative as $\delta$ reaches 0.9. However when the error terms are non-normal, the BPS approach tends to spuriously over-reject. Table 1 reveals that the BPS over-rejection problem worsens as: (i) the degree of tail heaviness increases (from $t_{12}$ to $t_6$ errors), (ii) $N$ increases, and (iii) $T$ increases. The most egregious instance occurs with $T=252$ where the FWER with $\text{BPS}_a$ and $\text{BPS}_b$ attains nearly $100\%$ when $N=500$ under $t_6$ errors. This makes clear that finite-sample FWER control with the BPS thresholding estimator based on $c_{\alpha}(N)$ in ((ref)) is heavily dependent on normality, even as the sample size $T$ increases. Therefore in order to ensure a fair comparison, all the power results for the BPS approach are based on strong FWER-adjusted critical values.\footnote{ Of course, the FWER-adjusted BPS method is not feasible in practice. It is merely used here as a benchmark for the SS and SD procedures.}

Panel C of Table 1 shows the power of the three FWER procedures when $\delta=0.9$. Note that $\text{BPS}_a$ and $\text{BPS}_b$ have identical size-adjusted power, reported on the lines labelled $\text{BPS}_a/\text{BPS}_b$. The results in Table 1 show that SS tends to be less powerful than SD, whose average power is quite close to that of BPS, especially as the sample size grows and tail heaviness increases. As expected, power is seen to increase with $T$, and to decrease as $N$ grows large and tail heaviness increases.

Table 2 reports the empirical rejection rates of the $k$-FWER procedures when $\delta=0.9.$ In Panel A, the empirical probability of making $k$ or more false rejections is effectively zero, corroborating Propositions (ref) and (ref). Table 2 reveals that the average power of the $k$-FWER procedures is higher when compared to the SS and SD results presented in Table 1. In fact, in Panel B, both $\text{SS}_{k}$ and $\text{SD}_{k}$ exhibit increasing power, with the gap between them narrowing as $k$ increases (cf. Remark (ref)).

The empirical exceedance probabilities and average power of the FDP procedures with $\gamma=0.1$ are reported in Table 3, again for the case where $\delta=0.9.$ Panel A shows that the procedures effectively control the FDP criterion by ensuring that the empirical probability of having a proportion of false discoveries exceed $\gamma$ is zero. From Panel B, it is clear that $\text{SS}_{\gamma}$ and $\text{SD}_{\gamma}$ tend to have the same average power as $M$ grows large. Compared to the $k$-FWER results in Table 2, the FDP procedures appear to reject more false hypotheses, particularly when there is a large number $M$ of hypotheses being tested.

Application to portfolio optimization

The proposed multiple testing procedures are further illustrated in this section with an application to Markowitz portfolio optimization using daily stock return data downloaded from the Center for Research in Security Prices (CRSP) for the period starting on January 1, 2004 and ending on December 31, 2021 (4531 observations).\footnote{Only U.S.-based common stocks as identified via CRSP share codes 10 and 11 from the NYSE, NYSE MKT, and NASDAQ exchanges are considered. Scheuch-Voigt-Weiss:2023 show how to access, filter, and download daily CRSP data in small batches through a Wharton Research Data Services (WRDS) connection using {R}.} The considered problem is that of an investor whose objective is to hold a global minimum variance (GMV) portfolio based on $N \in \{25, 100, 500 \}$ stocks over the next $21$ days, which is considered a “month." As Engle-Ledoit-Wolf:2019 and DeNard-Engle-Ledoit-Wolf:2022 argue in a similar context, the accuracy of the covariance matrix estimators should be evaluated primarily by the out-of-sample standard deviation achieved by the GMV portfolios. Nonetheless, there is also a significant interest in exploring additional performance measures. This interest persists even when considering that more realistic portfolio management scenarios would involve limitations on maximum position, industry position, factor exposure, etc. DeNard-Engle-Ledoit-Wolf:2022.

Let $t_{\text{b}}$ refer to the day when the portfolio is initially formed and subsequent rebalancing days. On those days the investor uses the returns from the past $L \in \{ 63, 126, 252 \}$ days to obtain $\hat{ \boldsymbol \Sigma }_{t_{\text{b}}}$, an estimate of the covariance matrix. The investor then uses $\hat{ \boldsymbol \Sigma }_{t_{\text{b}}}$ to find the GMV portfolio weights $ \hat {\boldsymbol \omega}_{t_{\text{b}} } = (\omega_{1,t_{\text{b}}},\ldots,\omega_{N,t_{\text{b}}})^{\prime}$ by solving the problem

with the assignment $\boldsymbol \Sigma \leftarrow \hat{ \boldsymbol \Sigma }_{t_{\text{b}}}$, and where $\boldsymbol \iota$ denotes an $N \times 1$ vector of ones. The second constraint in ((ref)) restricts each $\omega_{i,t_{\text{b}}}$ to be non-negative, meaning that short sales are prohibited. When short sales are allowed, the analytical solution to the GMV problem is $ \hat{ \boldsymbol \Sigma }_{t_{\text{b}} } ^{-1} \boldsymbol \iota / ( \boldsymbol \iota^{\prime} \hat{ \boldsymbol \Sigma }_{t_{\text{b}} } ^{-1} \boldsymbol \iota) $. With the short-selling restriction, the GMV problem is a quadratic optimization problem that can be solved numerically.\footnote{See Scheuch-Voigt-Weiss:2023 for a more detailed description of constrained portfolio optimization along with R code.} Restricting short sales is equivalent to shrinking towards zero the larger elements of the covariance matrix that would otherwise imply negative weights Jagannathan-Ma:2003, offering an interesting point of comparison.

With each choice of estimation window length $L$, the initialization $t_{\text{b}} \leftarrow 252$ corresponds to January 3, 2005, when the portfolio is first formed. This portfolio is then held for $21$ days and the resulting realized out-of-sample portfolio returns are $\hat {\boldsymbol \omega}_{t_{\text{b}} }^{\prime} \mathbf r_{\tau} $, for $\tau = t_{\text{b}}+1,\ldots,t_{\text{b}} + 21.$ After its initial formation, the portfolio is rebalanced on days $t_{\text{b}} \leftarrow t_{\text{b}} + 21$. This process yields 203 $t_{\text{b}}$'s and 4263 out-of-sample returns. Rebalancing consists of finding new GMV weights by solving ((ref)) with the updated covariance matrix estimate based on the returns from the last $L$ days. Following Engle-Ledoit-Wolf:2019, the investment universe on days $t_{\text{b}}$ is obtained by finding the $N$ largest capitalization stocks with a complete return history over the most recent $L$ days and a complete return future over the next 21 days. This way the composition of the investment universe evolves with each rebalancing day.

To obtain $\hat{ \boldsymbol \Sigma }_{t_{\text{b}}}$, it is natural to first consider the sample covariance matrix. Of course, the other choices considered include the approaches based on FWER control ($\text{BPS}_a$, $\text{BPS}_b$, $\text{SS}$, $\text{SD}$), the approaches based on $k$-FWER control ($\text{SD}_{\lfloor \log{M} \rfloor}$ and $\text{SD}_{\lfloor \sqrt{M} \rfloor}$), and the approaches based on FDP control ($\text{SS}_{\gamma}$ and $\text{SD}_{\gamma}$, with $\gamma=0.1$), which proceed by testing the significance of the $M \in \{300, 4950, 124750 \}$ distinct covariances in the rolling-window scheme.\footnote{The nominal levels are set to $\alpha = 5\%$ and $B-1=99$ resampling draws are used in the computation of the Monte Carlo $p$-values.} The effectiveness of those approaches is further assessed by comparing them to portfolios based on two shrinkage covariance matrix estimators: (i) the linear shrinkage (LS) estimator of Ledoit-Wolf:2004, which shrinks the sample covariance matrix towards the identity matrix; and (ii) the non-linear shrinkage (NLS) estimator proposed by Ledoit-Wolf:2015.\footnote{Specifically, the Ledoit-Wolf shrinkage covariance matrix estimates are computed with the {\tt linshrink_cov} and {\tt nlshrink_cov} commands available with the R package `nlshrink' Ramprasad:2016.}

The performance evaluation also includes the equally weighted (EW) portfolio, which bypasses ((ref)) altogether and simply sets $\hat \omega_{i,t_{\text{b}}}=1/N,$ $i=1,\ldots, N.$ This naive portfolio is a standard benchmark for comparisons; see DeMiguel-Garlappi-Uppal:1999 and Kirby-Ostdiek:2012, among others. The final strategy considered is the volatility timing (VT) portfolio with weights $\hat \omega_{i,t_{\text{b}}} = \hat\sigma^{-2}_{i,t_{\text{b}}} / \sum_{i=1}^N \hat\sigma^{-2}_{i,t_{\text{b}}},$ $i=1,\ldots, N,$ which was suggested by Kirby-Ostdiek:2012 as a competitor to EW. The VT portfolio can be seen as an aggressive form of shrinkage that sets to zero all the off-diagonal elements of the sample covariance matrix.

Note that there is no estimation risk associated with the EW strategy, which helps reduce the portfolio turnover. On the contrary, active strategies that generate high turnover will suffer more in the presence of transaction costs. To see this, note that for every dollar invested in the portfolio at time $t_{\text{b}}$, there are $ \hat \omega_{i,t_{\text{b}}} \prod_{\tau=t_{\text{b}}+1}^{t} (1+ r_{i,\tau}) $ dollars invested in asset $i$ at time $t$, for $t > t_{\text{b}}$ and as long as the portfolio is not rebalanced. Hence, at any time $t$ until the the next rebalancing occurs, the actual weight of asset $i$ in the portfolio is \[ \omega^{\ast}_{i,t} = \frac{ \hat \omega_{i,t_{\text{b}}}\prod_{\tau=t_{\text{b}}+1}^{t} (1+ r_{i,\tau}) } { \sum_{i=1}^N \hat \omega_{i,t_{\text{b}}} \prod_{\tau=t_{\text{b}}+1}^{t} (1+ r_{i,\tau}) }. \] When rebalancing occurs, the portfolio turnover can be defined as $\text{TO}_{t} = \sum_{i=1}^N \vert \hat \omega_{i, t } - \omega^{\ast}_{i,t} \vert, $ where $\hat \omega_{i, t }$ is the updated weight for asset $i$ at rebalancing time $t= t_{\text{b}}$. Denoting by $\kappa$ the transaction cost in proportion to the amount of wealth invested, the proportional cost of rebalancing all the portfolio positions is $\text{TC}_{t} = \kappa \text{TO}_{t}$ when $t= t_{\text{b}}$. Therefore, starting with $t_{\text{b}} \leftarrow 252$ and a normalized initial wealth of $W_{t_{\text{b}}} = \$1$ on that first day, wealth subsequently evolves according to \[ W_{t+1} = \left\{

\right. \] for $t = 252,\ldots,T-1$, with the updating rule $t_{\text{b}} \leftarrow t_{\text{b}} + 21$ to determine the rebalancing days. The out-of-sample portfolio return net of transaction costs is then given by $(W_{t+1}- W_{t})/W_{t}$. Following DeNard-Engle-Ledoit-Wolf:2022, $\kappa$ is set to 5 basis points to account for transaction costs.\footnote{French:2008 estimates the cost of trading stocks listed on the NYSE, AMEX, and NASDAQ, including total commissions, bid-ask spreads, and other costs investors pay for trading services. He finds that these costs have dropped significantly over time “from 146 basis points in 1980 to a tiny 11 basis points in 2006.”}

For each portfolio strategy, the following out-of-sample performance metrics are computed: (i) AV, the average return (annualized by multiplying by 252) in percent; (ii) SD, the standard deviation of returns (annualized by multiplying by $\sqrt{252}$) in percent; (iii) IR, the information ratio given as the annualized mean divided by the annualized standard deviation; (iv) TO, the average turnover; (v) MDD, the maximum drawdown in percent over the trading period;\footnote{A drawdown refers to the loss in the investment's value from a peak to a trough, before a new peak is attained. MDD is thus an indicator of downside risk.} and (vi) TW, the terminal wealth in dollars at the end of the trading period. Tables 4 and 5 present the results for $N=100$ and $500$, respectively. (The results with $N=25$ are in the Supplementary material.) When $N > L$, the sample covariance matrix is singular. Those cases are indicated with n/a in the tables.

Since a major objective of this paper is to show that the BPS approach can be improved upon, the statistical significance of the standard deviation differential between a given portfolio strategy and $\text{BPS}_b$ is assessed using the two-sided $p$-value of the prewhitened $\text{HAC}_{\text{PW}}$ method described in Ledoit-Wolf:2011 for testing the null hypothesis of equal standard deviations.\footnote{Note that these $p$-values assess the significance of pairwise differences; they no dot account for the multiplicity of tests.} This is done for each $(N, L)$ combination and whether short sales are allowed or not. Observe that the performance of the EW portfolio is influenced by the choice of $L$, since the composition of the investment universe on rebalancing days depends on the last $L$ and next 21 days. The main findings with respect to the out-of-sample standard deviation can be summarized as follows:

While minimizing volatility is a primary goal for a GMV investor, examining the other metrics in Tables 4 and 5 offers a more comprehensive view of risk, return, and overall performance. For instance, the TW columns reveal that wealth accumulated with LS and NLS strategies tends to be lower, especially when short selling is allowed. This can also be appraised from the AV columns in Tables 4 and 5, which show the tendency of LS and NLS to yield lower mean returns in comparison to the other portfolio strategies.

The EW portfolio generally results in the greatest accumulated wealth, owing to its near zero turnover. As expected, however, this payoff involves greater risk as can be seen from the SD and MDD columns. Indeed, the standard deviation under the EW strategy is seen everywhere in Tables 4 and 5 to be statistically larger at the 1% level than under $\text{BPS}_b$. And the MDD, which indicates the most significant loss sustained during the trading period, is exacerbated under the EW strategy.

The multiple testing strategies appear as “Goldilocks" solutions, striking a balance between risk and reward. In particular, when $N=100$ and $N=500$ and short selling is prohibited (Panel B in Tables 4 and 5), the stricter SS strategy tends to result in greater end-of-period wealth when compared to the other multiple testing strategies, as seen in the TW columns. This is further corroborated by comparing SS with VT in Panel B of Table 5, where terminal wealth with the even more stringent VT portfolio is higher than with the SS portfolio. On the other hand, the MDD columns reveal that the more lenient FDP procedures ($\text{SS}_{0.1}$ and $\text{SD}_{0.1}$) generally offer better protection against downside risk.

Further insight into these results can be gleaned from Figures 1--3, which show the proportion of correlations declared statistically significant by the multiple testing procedures each time the portfolio is rebalanced, given an estimation window of length $L=252$ days and $N=500$ assets. Figure 1 shows the results for the FWER procedures, Figure 2 for the $k$-FWER procedures, and Figure 3 for the FDP procedures. The solid line in Figure 1 shows that among the FWER procedures, $\text{BPS}_b$ declares the greatest number of non-zero correlations, followed by SD (dotted line), and then SS (dashed line). A comparison with Figure 2 reveals that the number of significant correlations found with $\text{BPS}_b$ generally falls in between $\text{SS}_{11}$ and $\text{SS}_{353},$ except for the period from March 2020 to March 2021 when $\text{BPS}_b$ is above $\text{SS}_{353}$. And from Figure 3 it is clear that $\text{SS}_{0.1}$ is most lenient in declaring non-zero correlations.

Finally note that these figures provide a gauge of the portfolio's overall diversification profile under a given multiple testing criterion. An improved profile is indicated when the lines dip; i.e., in periods when there are relatively fewer significant correlations among these large-cap stocks whose prices typically move in tandem. It is interesting to observe that $\text{BPS}_b$ surges above $\text{SS}_{353}$ in March 2020, precisely when the World Health Organization declared the COVID-19 outbreak a pandemic and stock market volatility soared Baker-Bloom-Davis-Kost-Sammon-Viratyosin:2020. This surge is likely due to the BPS approach's tendency to produce spurious rejections when confronted with non-Gaussian data conditions.

Concluding remarks

This paper has developed a sign-based Monte Carlo resampling method to regularize stock return covariance matrices. Following BPS, the method begins by testing the significance of pairwise correlations and then sets to zero the sample correlations whose multiplicity-adjusted $p$-values fall above the specified threshold level. A subsequent shrinkage step ensures that the final covariance matrix estimate is positive definite and well-conditioned, while preserving the zero entries achieved by thresholding.

The multiple testing procedures developed in this paper extend the BPS approach by offering strong control of the traditional FWER, the $k$-FWER, or the FDP, even in the presence of unknown non-normalities and heteroskedasticity. While the conservative FWER is arguably the most appropriate error rate measure for confirmatory purposes, the more lenient $k$-FWER (with $k>1$) and the FDP may be preferred in exploratory settings where less control is desired in exchange for more power. This was illustrated in an application to portfolio optimization where the goal of multiple testing was not the rigorous validation of candidate hypotheses. Instead, it aimed to regularize stock return covariance matrices, which were then evaluated for their influence on the out-of-sample performance of GMV portfolios.

{\bf Code availability.} An open-access implementation of the procedures developed in this paper is available on GitHub at {https://github.com/richardluger/CovRegMT}. The R code produces multiplicity-adjusted $p$-values and the associated covariance matrix estimates.

Acknowledgements

I wish to thank the editor Serena Ng, an associate editor, and two anonymous reviewers whose constructive feedback greatly improved the final paper. This work draws on research supported by {l’Autorit{\'e} des march{\'e}s financiers (AMF Qu{\'e}bec)} and the Social Sciences and Humanities Research Council of Canada.

\thispagestyle{empty}

\thispagestyle{empty}

\thispagestyle{empty}

\footnotetext[1]{Correspondence to: Department of Finance, Insurance and Real Estate, Laval University, Quebec City, Quebec G1V 0A6, Canada.

E-mail address: {[email removed]}. }

\pagenumbering{arabic}

\setcounter{footnote}{0}