Estimation of Large-Dimensional Covariance Matrices via Second-Order Stein-Type Regularization

Abstract

This paper tackles the problem of estimating the covariance matrix in large-dimension and small-sample-size scenarios. Inspired by the well-known linear shrinkage estimation, we propose a novel second-order Stein-type regularization strategy to generate well-conditioned covariance matrix estimators. We model the second-order Stein-type regularization as a quadratic polynomial concerning the sample covariance matrix and a given target matrix, representing the prior information of the actual covariance structure. To obtain available covariance matrix estimators, we choose the spherical and diagonal target matrices and develop unbiased estimates of the theoretical mean squared errors, which measure the distances between the actual covariance matrix and its estimators. We formulate the second-order Stein-type regularization as a convex optimization problem, resulting in the optimal second-order Stein-type estimators. Numerical simulations reveal that the proposed estimators can significantly lower the Frobenius losses compared with the existing Stein-type estimators. Moreover, a real data analysis in portfolio selection verifies the performance of the proposed estimators.

1. Introduction

As a fundamental problem in modern multivariate statistics and various practical applications, estimating the covariance matrix of a large-dimensional random vector has attracted significant interest in the last two decades [1,2]. The traditional sample covariance matrix (SCM) becomes unstable and ill-conditioned when the dimension increases proportionally with the sample size. The algorithm, which still employs SCM as the covariance matrix estimator, will result in a drop in performance or failure [3,4]. Although remaining unbiased, the SCM is not a satisfactory estimator of the actual covariance matrix anymore [5,6]. Therefore, it is of great concern to develop well-conditioned covariance matrix estimators in large-dimensional scenarios [7,8,9].

A comprehensive point of view to obtain well-conditioned estimators is by improving the SCM [10]. In the Stein-type regularization (linear shrinkage estimation), a target matrix is preset by rectifying the SCM according to prior information on the covariance structure. For example, the spherical target is a scalar matrix, with the coefficient being the average of the diagonal elements of the SCM [11]. The diagonal target is a diagonal matrix which retainsthe diagonal elements of the SCM [12,13]. Moreover, the Toeplitz-structured target is formed by averaging each diagonal of the SCM [14]. To some extent, the target matrices are covariance matrix estimators. Despite being biased, they usually enjoy low variances. The Stein-type estimator combines the SCM and the target matrix to balance the bias and variance. This method generates a well-conditioned estimator for the spherical target matrix by retaining the sample eigenvectors and shrinking the sample eigenvalues toward their grand mean. Moreover, the Stein-type regularization can be translated as a weighted average between the SCM and the target matrix. Many Stein-type estimators have been developed for various target matrices. Additionally, it is worth mentioning that the optimal Stein-type estimator can be expressed in closed form and significantly outperforms the traditional SCM under some appropriate criteria [15,16,17].

In addition, the nonlinear shrinkage estimation is proposed based on the random matrix theory [18]. By taking spectral decomposition on the SCM, this method retains the sample eigenvectors and estimates the actual eigenvalues by taking a nonlinear transformation on the sample eigenvalues. Then, the nonlinear shrinkage estimator is obtained by assembling the estimated eigenvalues and the sample eigenvectors. In the mean squared error (MSE) sense, the resulting nonlinear shrinkage estimator enjoys a significant advantage over the SCM. It also outperforms the Stein-type estimator for the spherical target. However, both the sample eigenvalues and sample eigenvectors have serious deficiencies in a high-dimensional case [19]. Hence, the existing nonlinear shrinkage strategy, which modifies the sample eigenvalues whileretaining the sample eigenvectors, has some limitations in improving the SCM performance. Moreover, the method can hardly handle the prior structure information employed in the Stein-type regularization. Hence, developing a new nonlinear shrinkage technique is essential to generate outperformed covariance matrix estimators.

This paper combines the SCM and the target matrix via a nonlinear shrinkage strategy to obtain well-conditioned estimators of a large-dimensional covariance matrix. The main contributions are the following:

• The second-order Stein-type estimator is modeled as a quadratic polynomial concerning the SCM and an almost surely (a.s.) positive definite target matrix. For the spherical and diagonal target matrices, the MSEs between the second-order Stein-type estimator and the actual covariance matrix are unbiasedly estimated under Gaussian distribution.
• We formulate the second-order Stein-type estimators for the two target matrices as convex quadratic programming problems. Then, the optimal second-order Stein-type estimators are immediately obtained.
• Some numerical simulations and application examples are provided for comparing the proposed second-order Stein-type estimators with the existing linear and nonlinear shrinkage estimators.

The outline of this paper is as follows. Section 2 proposes the second-order Stein-type estimator based on the Stein-type regularization. In Section 3, the spherical and diagonal matrices are employed as the target matrices. We obtain unbiased estimates of the MSEs between the second-order Stein-type estimators and the actual covariance matrix. The optimal second-order Stein-type estimators are obtained by solving the corresponding optimization problems. Section 4 provides some numerical simulations and two examples to discover the performance of the proposed estimators in large-dimensional scenarios. Section 5 concludes the major work of this paper.

2. Notation, Motivation, and Formulation

The symbol ${\mathbb{R}}^p$ denotes the set of entire p-dimensional real column vectors, ${\mathbb{R}}^{m\times n}$ denotes the set of entire $m\times n$ real matrices, and ${\mathbb{S}}^p$ denotes the set of entire $p\times p$ real symmetric matrices. The bold symbol $\mathbf{E}$ denotes the squared matrix having all entries 1 with appropriate dimensions. The symbol ${\mathbf{I}}_p$ denotes the $p\times p$ identity matrix. For a matrix $\mathbf{A}$, ${\mathbf{A}}^T$, $tr\left(\mathbf{A}\right)$, and $\parallel \mathbf{A}\parallel$, we denote its transpose, trace, and Frobenius matrix norm, respectively. For two matrices $\mathbf{A}$ and $\mathbf{B}$, $\mathbf{A}\circ \mathbf{B}$ means their Hadamard (element-wise) product.

Assume that $x_1,x_2,\dots ,x_n\in {\mathbb{R}}^p$ is an independent and identically distributed (i.i.d.) sample drawn from a certain distribution with mean 0 and covariance matrix $\mathsf{\Sigma }$. The SCM $\mathbf{S}$ is defined by

$$\mathbf{S}={\left(s_{ij}\right)}_{p\times p}=\frac{1}{n}\sum _{m=1}^n{x}_m{x}_m^T.$$

(1)

As is widely known, the SCM is ill conditioned in large-dimension scenarios and is even singular when $p>n$. The Stein-type regularization can produce a well-conditioned covariance matrix estimator based on the SCM [20,21,22,23].

For an a.s. positive definite target matrix $\mathbf{T}$ representing the prior information of the covariance structure, the Stein-type estimator combines the SCM and the target matrix with a linear function

$$f(\mathbf{S},\mathbf{T})=(1-w)\mathbf{S}+w\mathbf{T}.$$

(2)

Through an equivalent transformation, the expression given by (2) can be recombined as

$$f(\mathbf{S},\mathbf{T})=\mathbf{S}+w(\mathbf{T}-\mathbf{S}),$$

(3)

where $\mathbf{T}-\mathbf{S}$ is a regularization term and w is a tuning parameter. Moreover, the tuning parameter w lies in $(0,1]$ so as to keep the Stein-type estimator a.s. positive definite even when $n<p$. An interesting fact is that the matrix ${(\mathbf{T}-\mathbf{S})}^2$ is still symmetric and positive definite, which motivates us to consider further a quadratic function of the SCM and the target matrix.

For an a.s. positive definite target matrix $\mathbf{T}$, we model the second-order Stein-type estimator of the covariance matrix $\mathsf{\Sigma }$ as

$$\widehat{\mathsf{\Sigma }}=\mathbf{S}+w_1(\mathbf{T}-\mathbf{S})+w_2{(\mathbf{T}-\mathbf{S})}^2,$$

(4)

where $w_1,w_2$ are the tuning parameters. It is easy to find out that $\widehat{\mathsf{\Sigma }}\in {\mathbb{S}}^p$. In the same manner, we further assume $w_1\in (0,1]$ and $w_2\ge 0$ to keep the covariance estimator given by (4) to be a.s. positive definite. We note that the constraint on the tuning parameters is an easy-to-implement condition but is not necessary. One can also consider alternative assumptions, such as the condition number constraint, to obtain positive definite estimators of the large-dimensional covariance matrix [24].

Next, we choose the optimal tuning parameters in (4). In the Stein-type regularization, the MSE between the actual covariance matrix and its estimator is the most commonly used loss function. It includes unknown scalars concerning the expectation operator and the actual covariance matrix. One practical way is to make the MSE available by estimating the unknown scalars and obtaining the optimal shrinkage intensity by minimizing the available MSE [15]. Therefore, we still follow the above steps to find optimal tuning parameters in the second-order Stein-type estimator. To be specific, the loss function of the second-order Stein-type estimator $\widehat{\mathsf{\Sigma }}$ is defined by

$${\mathcal{M}}_{\mathbf{T}}\left(\mathbf{w}\right|\mathsf{\Sigma })=\mathbb{E}[\parallel \widehat{\mathsf{\Sigma }}-\mathsf{\Sigma }{\parallel }^2],$$

(5)

where $\mathbf{w}={(w_1,w_2)}^T$. Through substituting the expression of $\widehat{\mathsf{\Sigma }}$ into (5), we can obtain

$${\mathcal{M}}_{\mathbf{T}}\left(\mathbf{w}\right|\mathsf{\Sigma })=\frac{1}{2}{\mathbf{w}}^T\mathbf{H}\left(\mathbf{T}\right)\mathbf{w}-2{\mathbf{w}}^T\mathbf{b}\left(\mathbf{T}\right)+c,$$

(6)

where

$$\begin{array}{cc}\hfill \mathbf{H}\left(\mathbf{T}\right)& =2\left[\begin{array}{cc}\hfill \mathbb{E}[tr{(\mathbf{S}-\mathbf{T})}^2]& \hfill -\mathbb{E}[tr{(\mathbf{S}-\mathbf{T})}^3]\\ \hfill -\mathbb{E}[tr{(\mathbf{S}-\mathbf{T})}^3]& \hfill \mathbb{E}[tr{(\mathbf{S}-\mathbf{T})}^4]\end{array}\right],\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \mathbf{b}\left(\mathbf{T}\right)& =\left[\begin{array}{c}\mathbb{E}[tr(\mathbf{S}-\mathsf{\Sigma }\left)\right(\mathbf{S}-\mathbf{T}\left)\right]\\ -\mathbb{E}[tr(\mathbf{S}-\mathsf{\Sigma }){(\mathbf{S}-\mathbf{T})}^2]\end{array}\right],\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill c& =\mathbb{E}[tr{(\mathbf{S}-\mathsf{\Sigma })}^2].\hfill \end{array}$$

(9)

Therefore, the second-order Stein-type regularization can be modeled as the following optimization problem:

$$\begin{array}{cc}\hfill min& \frac{1}{2}{\mathbf{w}}^T\mathbf{H}\left(\mathbf{T}\right)\mathbf{w}-2{\mathbf{w}}^T\mathbf{b}\left(\mathbf{T}\right)+c\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& 0\le {\mathbf{e}}_1\mathbf{w}\le 1,\hfill \\ \hfill & {\mathbf{e}}_2\mathbf{w}\ge 0,\hfill \end{array}$$

(10)

where ${\mathbf{e}}_1=(1,0)$ and ${\mathbf{e}}_2=(0,1)$. It is easy to see that the loss function in problem (10) is a binary quadratic polynomial function about $\mathbf{w}$, and $\mathbf{H}\left(\mathbf{T}\right)$ is the Hessian matrix. By the Cauchy–Schwarz inequality, we have

$${\left(\mathbb{E}[tr{(\mathbf{S}-\mathbf{T})}^3]\right)}^2\le \mathbb{E}[tr{(\mathbf{S}-\mathbf{T})}^2]\mathbb{E}[tr{(\mathbf{S}-\mathbf{T})}^4].$$

(11)

Therefore, the Hessian matrix $\mathbf{H}\left(\mathbf{T}\right)$ is positive definite. Then, the optimization problem (10) is a convex quadratic program. However, it cannot be solved because the quantities $\mathbf{H}\left(\mathbf{T}\right),\mathbf{b}\left(\mathbf{T}\right)$, and c in the objective function are unknown. When the underlying distribution is given and the target matrix is prespecified, we can estimate the unknown quantities $\mathbf{H}\left(\mathbf{T}\right),\mathbf{b}\left(\mathbf{T}\right)$, and c. Then, the optimization problem (10) turns out to be available based on plug-in strategy and can be effectively solved.

Remark 1.

The unknown quantity c does not affect the choice of optimal tuning parameter $\mathbf{w}$ in the optimization problem (10). Moreover, it is the theoretical MSE between the actual covariance matrix and the classic SCM and plays an important role in evaluating the performance of improved covariance estimators based on the SCM.

3. Optimal Second-Order Stein-Type Estimators

In this section, as the target matrix is specified, we estimate the corresponding unknown quantities under Gaussian distribution, then establish the available version of the optimization problem to obtain the optimal second-order Stein-type estimator.

3.1. Target Matrices

As mentioned before, the target matrix represents the prior information of the actual covariance structure. In the Stein-type regularization, the commonly used target matrices include the spherical target, the diagonal target, the Toeplitz-structured target, and the tapered SCM. Among these, the spherical target and the diagonal target are a.s. positive definite, whereas both the Toeplitz-structured target and the tapered SCM are unnecessary. Thereby, in the second-order Stein-type regularization, we employ the spherical target and the diagonal target, given by

$${\mathbf{T}}_1=\frac{tr\left(\mathbf{S}\right)}{p}{\mathbf{I}}_p,{\mathbf{T}}_2=diag(s_{11},\dots ,s_{pp}).$$

(12)

The diagonal target ${\mathbf{T}}_2$ is also denoted as ${\mathbf{D}}_{\mathbf{S}}$ because it consists of the diagonal elements of the SCM.

3.2. Available Loss Functions

For the target matrices ${\mathbf{T}}_1$ and ${\mathbf{T}}_2$, we unbiasedly estimate the loss functions given by (6) through plugging in the estimates of unknown quantities $\mathbf{H}\left({\mathbf{T}}_i\right),\mathbf{b}\left({\mathbf{T}}_i\right),i=1,2$ and c under Gaussian distribution.

First of all, by directly removing the expectation operator, the Hessian matrices $\mathbf{H}\left({\mathbf{T}}_i\right)$ can be estimated by

$$\widehat{\mathbf{H}}\left({\mathbf{T}}_i\right)=2\left[\begin{array}{cc}\hfill tr{(\mathbf{S}-{\mathbf{T}}_i)}^2& \hfill -tr{(\mathbf{S}-{\mathbf{T}}_i)}^3\\ \hfill -tr{(\mathbf{S}-{\mathbf{T}}_i)}^3& \hfill tr{(\mathbf{S}-{\mathbf{T}}_i)}^4\end{array}\right].$$

(13)

Furthermore, $\widehat{\mathbf{H}}\left({\mathbf{T}}_i\right),i=1,2$ are, respectively, unbiased estimates of $\mathbf{H}\left({\mathbf{T}}_i\right),i=1,2$.

Next, for $i=1,2$, we decompose the unknown vectors $\mathbf{b}\left({\mathbf{T}}_i\right)$ into two terms,

$$\begin{array}{cc}\hfill \mathbf{b}\left({\mathbf{T}}_i\right)& =\left[\begin{array}{c}\mathbb{E}[tr(\mathbf{S}-\mathsf{\Sigma })(\mathbf{S}-{\mathbf{T}}_i)]\\ -\mathbb{E}[tr(\mathbf{S}-\mathsf{\Sigma }){(\mathbf{S}-{\mathbf{T}}_i)}^2]\end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{c}\mathbb{E}[tr\left(\mathbf{S}(\mathbf{S}-{\mathbf{T}}_i)\right)]\\ -\mathbb{E}[tr\left(\mathbf{S}{(\mathbf{S}-{\mathbf{T}}_i)}^2\right)]\end{array}\right]+\left[\begin{array}{c}-\mathbb{E}[tr\left(\mathsf{\Sigma }(\mathbf{S}-{\mathbf{T}}_i)\right)]\\ \mathbb{E}[tr\left(\mathsf{\Sigma }{(\mathbf{S}-{\mathbf{T}}_i)}^2\right)]\end{array}\right]\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & \triangleq \mathbf{u}\left({\mathbf{T}}_i\right)+\mathbf{v}\left({\mathbf{T}}_i\right).\hfill \end{array}$$

(15)

Similar to the Hessian matrices $\mathbf{H}\left({\mathbf{T}}_1\right)$ and $\mathbf{H}\left({\mathbf{T}}_2\right)$, the first term $\mathbf{u}\left({\mathbf{T}}_i\right)$ can be unbiasedly estimated by

$$\widehat{\mathbf{u}}\left({\mathbf{T}}_i\right)=\left[\begin{array}{c}tr\left(\mathbf{S}(\mathbf{S}-{\mathbf{T}}_i)\right)\\ -tr\left(\mathbf{S}{(\mathbf{S}-{\mathbf{T}}_i)}^2\right)\end{array}\right].$$

(16)

Therefore, we only need to estimate the second term $\mathbf{v}\left({\mathbf{T}}_i\right)$. It is challenging to estimate $\mathbf{v}\left({\mathbf{T}}_i\right)$ unbiased because it includes both the expectation operator and the actual covariance matrix $\mathsf{\Sigma }$. We need the following moment properties about the Wishart distribution [25].

Lemma 1.

Denote $\mathbf{A}$ and $\mathbf{B}$ as arbitrary symmetric nonrandom matrices, $\mathbf{S}$ is the sample covariance matrix given by (1), and then the following equalities hold under Gaussian distribution:

$$\begin{array}{cc}\hfill \mathbb{E}\left[\mathbf{S}\mathbf{A}\mathbf{S}\right]& =\frac{n+1}{n}\mathsf{\Sigma }\mathbf{A}\mathsf{\Sigma }+\frac{1}{n}tr\left(\mathsf{\Sigma }\mathbf{A}\right)\mathsf{\Sigma },\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill \mathbb{E}[tr(\mathbf{A}\mathbf{S}\left)\mathbf{S}\right]& =tr\left(\mathbf{A}\mathsf{\Sigma }\right)\mathsf{\Sigma }+\frac{2}{n}\mathsf{\Sigma }\mathbf{A}\mathsf{\Sigma },\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill \mathbb{E}[tr(\mathbf{A}\mathbf{S})tr(\mathbf{B}\mathbf{S}\left)\right]& =tr\left(\mathbf{A}\mathsf{\Sigma }\right)tr\left(\mathbf{B}\mathsf{\Sigma }\right)+\frac{2}{n}tr\left(\mathbf{A}\mathsf{\Sigma }\mathbf{B}\mathsf{\Sigma }\right).\hfill \end{array}$$

(19)

By Lemma 1, letting $\mathbf{A}=\mathbf{B}={\mathbf{I}}_p$, we can obtain

$$\mathbb{E}[tr\left({\mathbf{S}}^2\right)]=\frac{n+1}{n}tr\left({\mathsf{\Sigma }}^2\right)+\frac{1}{n}{tr}^2\left(\mathsf{\Sigma }\right),\mathbb{E}\left[{tr}^2\left(\mathbf{S}\right)\right]={tr}^2\left(\mathsf{\Sigma }\right)+\frac{2}{n}tr\left({\mathsf{\Sigma }}^2\right).$$

(20)

Moreover, letting $\mathbf{A}=\mathsf{\Sigma }$ and $\mathbf{B}={\mathbf{I}}_p$, we have

$$\begin{array}{cc}\hfill \mathbb{E}[tr\left({\mathbf{S}}^2\mathsf{\Sigma }\right)]& =\frac{n+1}{n}tr\left({\mathsf{\Sigma }}^3\right)+\frac{1}{n}tr\left(\mathsf{\Sigma }\right)tr\left({\mathsf{\Sigma }}^2\right),\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill \mathbb{E}[tr(\mathbf{S}\mathsf{\Sigma })tr(\mathbf{S}\left)\right]& =tr\left(\mathsf{\Sigma }\right)tr\left({\mathsf{\Sigma }}^2\right)+\frac{2}{n}tr\left({\mathsf{\Sigma }}^3\right).\hfill \end{array}$$

(22)

Lemma 1 is very helpful to compute the second term $\mathbf{v}\left({\mathbf{T}}_i\right)$ in (14). For the spherical target matrix ${\mathbf{T}}_1$, we have

$$\mathbb{E}[tr\left(\mathsf{\Sigma }(\mathbf{S}-{\mathbf{T}}_1)\right)]=tr\left({\mathsf{\Sigma }}^2\right)-\frac{1}{p}{tr}^2\left(\mathsf{\Sigma }\right),$$

(23)

and

$$\begin{array}{cc}\hfill \mathbb{E}[tr\left(\mathsf{\Sigma }{(\mathbf{S}-{\mathbf{T}}_1)}^2\right)]& =\mathbb{E}[tr\left({\mathbf{S}}^2\mathsf{\Sigma }\right)]-\frac{2}{p}\mathbb{E}[tr\left(\mathbf{S}\right)tr\left(\mathbf{S}\mathsf{\Sigma }\right)]+\frac{1}{p^2}\mathbb{E}\left[{tr}^2\left(\mathbf{S}\right)\right]tr\left(\mathsf{\Sigma }\right)\hfill \\ \hfill & =\frac{np+p-4}{np}tr\left({\mathsf{\Sigma }}^3\right)+\frac{p^2-2np+2}{n{p}^2}tr\left(\mathsf{\Sigma }\right)tr\left({\mathsf{\Sigma }}^2\right)+\frac{1}{p^2}{tr}^3\left(\mathsf{\Sigma }\right).\hfill \end{array}$$

(24)

For the diagonal target matrix ${\mathbf{T}}_2$, we have

$$\mathbb{E}[tr\left(\mathsf{\Sigma }(\mathbf{S}-{\mathbf{T}}_2)\right)]=tr\left({\mathsf{\Sigma }}^2\right)-tr\left({\mathbf{D}}_{\mathsf{\Sigma }}^2\right),$$

(25)

where ${\mathbf{D}}_{\mathsf{\Sigma }}=diag({\mathsf{\Sigma }}_{11},\dots ,{\mathsf{\Sigma }}_{pp})$. Moreover, we can obtain

$$\begin{array}{cc}\hfill & \mathbb{E}[tr\left(\mathsf{\Sigma }{(\mathbf{S}-{\mathbf{T}}_2)}^2\right)]\hfill \\ \hfill & =\mathbb{E}[tr\left({\mathbf{S}}^2\mathsf{\Sigma }\right)]-2\mathbb{E}[tr\left(\mathbf{S}{\mathbf{T}}_2\mathsf{\Sigma }\right)]+\mathbb{E}[tr\left({\mathbf{T}}_2^2\mathsf{\Sigma }\right)]\hfill \\ \hfill & =\frac{n+1}{n}tr\left({\mathsf{\Sigma }}^3\right)+\frac{1}{n}tr\left(\mathsf{\Sigma }\right)tr\left({\mathsf{\Sigma }}^2\right)-\frac{2n+4}{n}tr\left({\mathbf{D}}_{\mathsf{\Sigma }}{\mathsf{\Sigma }}^2\right)+\frac{n+2}{n}tr\left({\mathbf{D}}_{\mathsf{\Sigma }}^3\right).\hfill \end{array}$$

(26)

Denote

$$\begin{array}{ccc}{a}_1=tr\left({\mathsf{\Sigma }}^2\right),\hfill & a_2={tr}^2\left(\mathsf{\Sigma }\right),\hfill & a_3=tr\left({\mathbf{D}}_{\mathsf{\Sigma }}^2\right),\hfill \\ b_1=tr\left({\mathsf{\Sigma }}^3\right),\hfill & b_2={tr}^3\left(\mathsf{\Sigma }\right),\hfill & b_3=tr\left(\mathsf{\Sigma }\right)tr\left({\mathsf{\Sigma }}^2\right),\hfill \\ c_1=tr\left({\mathbf{D}}_{\mathsf{\Sigma }}^3\right),\hfill & c_2=tr\left({\mathbf{D}}_{\mathsf{\Sigma }}{\mathsf{\Sigma }}^2\right),\hfill \end{array}$$

(27)

then the vectors $\mathbf{v}\left({\mathbf{T}}_i\right),i=1,2$ can be rewritten as

$$\begin{array}{cc}\hfill \mathbf{v}\left({\mathbf{T}}_1\right)& =\left[\begin{array}{c}\frac{1}{p}{a}_2-a_1\\ \frac{np+p-4}{np}{b}_1+\frac{1}{p^2}{b}_2+\frac{p^2-2np+2}{n{p}^2}{b}_3\end{array}\right],\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill \mathbf{v}\left({\mathbf{T}}_2\right)& =\left[\begin{array}{c}{a}_3-a_1\\ \frac{n+1}{n}{b}_1+\frac{1}{n}{b}_3+\frac{n+2}{n}{c}_1-\frac{2n+4}{n}{c}_2\end{array}\right].\hfill \end{array}$$

(29)

It is worth noting that each element in $\mathbf{v}\left({\mathbf{T}}_i\right),i=1,2$ is a linear combination of the quantities in (27). Therefore, we only need to find out the estimates of the quantities in (27).

Firstly, the unbiased estimates of the quantities $a_i,i=1,2,3$ were proposed in [12,26,27],

$$\begin{array}{cc}\hfill {\alpha }_1& ={\tau }_a\left(ntr\left({\mathbf{S}}^2\right)-{tr}^2\left(\mathbf{S}\right)\right),{\alpha }_2={\tau }_a\left((n+1){tr}^2\left(\mathbf{S}\right)-2tr\left({\mathbf{S}}^2\right)\right),\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\alpha }_3& ={\tau }_a(n-1)tr\left({\mathbf{D}}_{\mathbf{S}}^2\right),\hfill \end{array}$$

(31)

where ${\tau }_a=\frac{n}{(n-1)(n+2)}$.

Secondly, denote the matrix $\mathbf{W}$ as

$$\mathbf{W}=\left[\begin{array}{ccc}{n}^2& 2& -3n\\ 16& n^2+3n-2& -6(n+2)\\ -4n& -(n+2)& n^2+2n+4\end{array}\right],$$

(32)

Then, the unbiased estimates of $b_i,i=1,2,3$ can be obtained by the following theorem.

Theorem 1.

Under Gaussian distribution, the following equation holds when $n\ge 3$:

$$\mathbb{E}\left[\begin{array}{c}tr\left({\mathbf{S}}^3\right)\\ {tr}^3\left(\mathbf{S}\right)\\ tr\left(\mathbf{S}\right)tr\left({\mathbf{S}}^2\right)\end{array}\right]={\left({\tau }_b\mathbf{W}\right)}^{-1}\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right],$$

(33)

where ${\tau }_b=\frac{n^2}{(n-1)(n-2)(n+2)(n+4)}$.

Proof. 

The actual covariance matrix $\mathsf{\Sigma }$ has the spectral decomposition which is described as $\mathsf{\Sigma }={\mathsf{\Gamma }}^T\mathsf{\Lambda }\mathsf{\Gamma }$, where $\mathsf{\Lambda }=diag({\lambda }_1,\dots ,{\lambda }_p)$ is a diagonal matrix consisting of the eigenvalues and $\mathsf{\Gamma }$ is the corresponding unitary matrix. Define $\mathbf{F}={\mathsf{\Sigma }}^{\frac{1}{2}}$, and $\mathbf{F}$ is a symmetric matrix; then we have ${\mathbf{F}}^2=\mathsf{\Sigma }$. For $m=1,\dots ,n$, denote ${\mathbf{z}}_m={\mathbf{F}}^{-1}{x}_m$, then ${\mathbf{z}}_m$ is an i.i.d. sample and ${\mathbf{z}}_m\sim N(\mathbf{0},{\mathbf{I}}_p)$. Let $\mathbf{X}=(x_1,\dots ,x_n)$ and $\mathbf{Z}=({\mathbf{z}}_1,\dots ,{\mathbf{z}}_n)$; then, we have $\mathbf{X}=\mathbf{F}\mathbf{Z}$. Notice that the SCM is $\mathbf{S}=\frac{1}{n}{\sum }_{m=1}^n{x}_m{x}_m^T$; therefore, we have

$$n\mathbf{S}=\mathbf{X}{\mathbf{X}}^T=\mathbf{F}\mathbf{Z}{\mathbf{Z}}^T\mathbf{F}.$$

(34)

Moreover, we can obtain

$$tr\left(n\mathbf{S}\right)=tr\left(\mathbf{F}\mathbf{Z}{\mathbf{Z}}^T\mathbf{F}\right)=tr\left({\mathbf{Z}}^T\mathsf{\Sigma }\mathbf{Z}\right)=tr\left({\mathbf{Z}}^T{\mathsf{\Gamma }}^T\mathsf{\Lambda }\mathsf{\Gamma }\mathbf{Z}\right).$$

(35)

Define a matrix ${\mathbf{Q}}^T={\mathbf{Z}}^T{\mathsf{\Gamma }}^T=({\mathbf{q}}_1,\dots ,{\mathbf{q}}_p)$ and denote $v_{ii}={\mathbf{q}}_i^T{\mathbf{q}}_i$ and $v_{ij}={\mathbf{q}}_i^T{\mathbf{q}}_j={\mathbf{q}}_j^T{\mathbf{q}}_i$ for $i,j\in \{1,\dots ,p\}$, then the above equation can be rewritten as

$$tr\left(n\mathbf{S}\right)=tr\left({\mathbf{Q}}^T\mathsf{\Lambda }\mathbf{Q}\right)=\sum _{i=1}^p{\lambda }_i{v}_{ii}.$$

(36)

In a same manner, we can obtain

$$\begin{array}{c}\hfill tr{\left(n\mathbf{S}\right)}^2=\sum _{i,j=1}^p{\lambda }_i{\lambda }_j{v}_{ij}^2,{tr}^2\left(n\mathbf{S}\right)=\sum _{i,j=1}^p{\lambda }_i{\lambda }_j{v}_{ii}{v}_{jj},\end{array}$$

(37)

and

$$\begin{array}{cc}\hfill tr{\left(n\mathbf{S}\right)}^3& =\sum _{i,j,k=1}^p{\lambda }_i{\lambda }_j{\lambda }_k{v}_{ij}{v}_{ik}{v}_{jk},\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill {tr}^3\left(n\mathbf{S}\right)& =\sum _{i,j,k=1}^p{\lambda }_i{\lambda }_j{\lambda }_k{v}_{ii}{v}_{jj}{v}_{kk},\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill tr\left(n\mathbf{S}\right)tr{\left(n\mathbf{S}\right)}^2& =\sum _{i,j,k=1}^p{\lambda }_i{\lambda }_j{\lambda }_k{v}_{ij}^2{v}_{kk}.\hfill \end{array}$$

(40)

By the moment properties of random variables $v_{ii}$ and $v_{ij}$ in [26,28], we have

$$\begin{array}{c}\hfill \begin{array}{ccc}\mathbb{E}\left[v_{ii}\right]=n,\hfill & \mathbb{E}\left[v_{ii}^2\right]=n(n+2),\hfill & \mathbb{E}\left[v_{ij}^2\right]=n,\hfill \\ \mathbb{E}\left[v_{ii}{ij}^2\right]=n(n+2),\hfill & \mathbb{E}\left[v_{ij}{v}_{ik}{v}_{jk}\right]=n\hfill & \mathbb{E}\left[v_{ii}^3\right]=n(n+2)(n+4).\hfill \end{array}\end{array}$$

(41)

where $i,j,k$ are arbitrary mutually unequal numbers in {1, ..., p}. Next, we compute the mathematical expectation of Equations (38)–(40) based on the moment properties in (41). Denote ${\mu }_1={\sum }_{i=1}^p{\lambda }_i^3$, ${\mu }_2={\sum }_{i\ne j}^p{\lambda }_i^2{\lambda }_j$ and ${\mu }_3={\sum }_{i\ne j\ne k}^p{\lambda }_i{\lambda }_j{\lambda }_k$, where $i\ne j\ne k$ means that $i,j,k$ are mutually unequal, and we have

$$\begin{array}{cc}\hfill \mathbb{E}[tr{\left(n\mathbf{S}\right)}^3]& =n(n+2)(n+4){\mu }_1+3n(n+2){\mu }_2+n{\mu }_3,\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill \mathbb{E}\left[{tr}^3\left(n\mathbf{S}\right)\right]& =n(n+2)(n+4){\mu }_1+3n^2(n+2){\mu }_2+n^3{\mu }_3,\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill \mathbb{E}[tr\left(n\mathbf{S}\right)tr{\left(n\mathbf{S}\right)}^2]& =n(n+2)(n+4){\mu }_1+n{(n+2)}^2{\mu }_2+n^2{\mu }_3.\hfill \end{array}$$

(44)

Denote a matrix $\mathbf{D}$ as

$$\mathbf{D}=\left[\begin{array}{ccc}n(n+2)(n+4)& 3n(n+2)& n\\ n(n+2)(n+4)& 3n^2(n+2)& n^3\\ n(n+2)(n+4)& n{(n+2)}^2& n^2\end{array}\right],$$

(45)

then the above Equations (42)–(44) can be rewritten in the form of a matrix equation:

$$\left[\begin{array}{c}\mathbb{E}[tr{\left(n\mathbf{S}\right)}^3]\\ \mathbb{E}\left[{tr}^3\left(n\mathbf{S}\right)\right]\\ \mathbb{E}[tr\left(n\mathbf{S}\right)tr{\left(n\mathbf{S}\right)}^2]\end{array}\right]=\mathbf{D}\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ {\mu }_3\end{array}\right].$$

(46)

Furthermore, the unknown quantities $b_i,i=1,2,3$ can be decomposed as

$$b_1={\mu }_1,b_2={\mu }_1+3{\mu }_2+{\mu }_3,b_3={\mu }_1+{\mu }_2.$$

(47)

Therefore, we have

$$\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ {\mu }_3\end{array}\right]={\left[\begin{array}{ccc}1& 0& 0\\ 1& 3& 1\\ 1& 1& 0\end{array}\right]}^{-1}\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right].$$

(48)

By the Equations (46) and (48), we have

$$\begin{array}{cc}\hfill \left[\begin{array}{c}\mathbb{E}[tr{\left(n\mathbf{S}\right)}^3]\\ \mathbb{E}\left[{tr}^3\left(n\mathbf{S}\right)\right]\\ \mathbb{E}[tr\left(n\mathbf{S}\right)tr{\left(n\mathbf{S}\right)}^2]\end{array}\right]& =\mathbf{D}{\left[\begin{array}{ccc}1& 0& 0\\ 1& 3& 1\\ 1& 1& 0\end{array}\right]}^{-1}\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right].\hfill \end{array}$$

(49)

For $n\ge 3$, the following equality holds:

$$\mathbf{D}{\left[\begin{array}{ccc}1& 0& 0\\ 1& 3& 1\\ 1& 1& 0\end{array}\right]}^{-1}={\left({\tau }_b\mathbf{W}\right)}^{-1}.$$

(50)

Therefore, we can obtain

$$\begin{array}{c}\hfill \mathbb{E}\left[\begin{array}{c}tr{\left(n\mathbf{S}\right)}^3\\ {tr}^3\left(n\mathbf{S}\right)\\ tr\left(n\mathbf{S}\right)tr{\left(n\mathbf{S}\right)}^2\end{array}\right]={\left({\tau }_b\mathbf{W}\right)}^{-1}\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right].\end{array}$$

(51)

This completes the proof. □

Then, the unknown scalars $b_i,i=1,2,3$ can be unbiasedly estimated by

$$\left[\begin{array}{c}{\beta }_1\\ {\beta }_2\\ {\beta }_3\end{array}\right]={\tau }_b\mathbf{W}\left[\begin{array}{c}tr\left({\mathbf{S}}^3\right)\\ {tr}^3\left(\mathbf{S}\right)\\ tr\left(\mathbf{S}\right)tr\left({\mathbf{S}}^2\right)\end{array}\right].$$

(52)

Remark 2.

In a large sample scenario, the unknown scalars $b_i,i=1,2,3$ can be consistently estimated by $tr\left({\mathbf{S}}^3\right)$, ${tr}^3\left(\mathbf{S}\right)$ and $tr\left(\mathbf{S}\right)tr\left({\mathbf{S}}^2\right)$ [29]. Theorem 1 shows that these estimates are biased. Moreover, the biases become non-ignorable in high-dimensional situations [30]. Furthermore, by Theorem 1, the biases can be eliminated by the linear combinations of $tr\left({\mathbf{S}}^3\right)$, ${tr}^3\left(\mathbf{S}\right)$, and $tr\left(\mathbf{S}\right)tr\left({\mathbf{S}}^2\right)$.

Thirdly, denote the matrices $\mathbf{G}$, $\mathbf{R}$, and $\mathbf{K}$ as follows:

$$\begin{array}{cc}\hfill \mathbf{G}& =\left(g_{ij}\right)\mathrm{with}{g}_{ij}={\parallel x^i\circ x^j\parallel }^2,\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill \mathbf{R}& =\left(r_{ij}\right)\mathrm{with}{r}_{ij}={\parallel x^i\circ x^i\circ x^j\parallel }^2,\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill \mathbf{K}& =\left(k_{ij}\right)\mathrm{with}{k}_{ij}=(x^i\circ x^i){(x^i\circ x^j)}^T,\hfill \end{array}$$

(55)

where $x^i$ is the observations of i-th variable for $i=1,\dots ,p$, then the following theorem holds.

Theorem 2.

Under Gaussian distribution, the following equations hold when $n\ge 3$:

$$\begin{array}{cc}\hfill \mathbb{E}[n^{3}tr\left({\mathbf{D}}_{\mathbf{S}}^3\right)-3ntr\left(\mathbf{G}{\mathbf{D}}_{\mathbf{S}}\right)+2tr\left(\mathbf{R}\right)]& ={\tau }_c^{-1}{c}_1,\hfill \end{array}$$

(56)

$$\begin{array}{cc}\hfill \mathbb{E}[n^{3}tr\left({\mathbf{D}}_{\mathbf{S}}{\mathbf{S}}^2\right)-2ntr\left(\mathbf{K}\mathbf{S}\right)-ntr\left(\mathbf{G}{\mathbf{D}}_{\mathbf{S}}\mathbf{E}\right)+2tr\left(\mathbf{R}\mathbf{E}\right)]& ={\tau }_c^{-1}{c}_2,\hfill \end{array}$$

(57)

where ${\tau }_c=\frac{1}{n(n-1)(n-2)}$.

Proof. 

Let $\mathbf{F}=\left(f_{ij}\right)$, and $\mathbf{F}$ is a symmetric matrix; then, we have ${\mathbf{F}}^2=\mathsf{\Sigma }=\left({\sigma }_{ij}\right)$. Therefore, for arbitrary $i,j\in \{1,\dots ,p\}$, the equalities ${\sigma }_{ij}={\sum }_{k=1}^p{f}_{ik}{f}_{jk}$ and ${\sigma }_{ii}={\sum }_{k=1}^p{f}_{ik}^2$ hold. For $m=1,\dots ,n$, denote $x_m={(x_{m1},x_{m2},\dots ,x_{mp})}^T$ and ${\mathbf{z}}_m={\mathbf{F}}^{-1}{x}_m$, then ${\mathbf{z}}_m$ is an i.i.d. sample and ${\mathbf{z}}_m\sim N(\mathbf{0},{\mathbf{I}}_p)$. Let ${\mathbf{z}}_m={(z_{m1},z_{m2},\dots ,z_{mp})}^T$, then $z_{mk},m=1,\dots ,n,k=1,\dots ,p$ are mutually independent standard Gaussian random variables. For arbitrary $m\in \{1,\dots ,n\}$ and $i,j\in \{1,\dots ,p\}$, we have $x_{mi}={\sum }_{k=1}^p{f}_{ik}{z}_{mk}$ and $x_{mj}={\sum }_{k=1}^p{f}_{jk}{z}_{mk}$. Denote that SCM $\mathbf{S}=\left(s_{ij}\right)$, then $s_{ij}$ can be decomposed as follows:

$$s_{ij}=\frac{1}{n}\sum _{m=1}^n{x}_{mi}{x}_{mj}=\frac{1}{n}\sum _{m=1}^n\sum _{k_1,k_2=1}^p{f}_{i{k}_1}{f}_{j{k}_2}{z}_{m{k}_1}{z}_{m{k}_2}.$$

(58)

Then, for arbitrary $m\in \{1,\dots ,n\}$ and $i,j\in \{1,\dots ,p\}$, we have

$$\mathbb{E}\left[x_{mi}{x}_{mj}\right]=\sum _{k_1,k_2=1}^p{f}_{i{k}_1}{f}_{j{k}_2}\mathbb{E}\left[z_{m{k}_1}{z}_{m{k}_2}\right]=\sum _{k=1}^p{f}_{ik}{f}_{jk}={\sigma }_{ij}.$$

(59)

Especially when $i=j$, we have $\mathbb{E}\left[x_{mi}^2\right]={\sigma }_{ii}$. Then, we can obtain

$$\begin{array}{cc}\hfill & \mathbb{E}[n^{3}tr\left({\mathbf{D}}_{\mathbf{S}}^3\right)-3ntr\left(\mathbf{G}{\mathbf{D}}_{\mathbf{S}}\right)+2tr\left(\mathbf{R}\right)]=\sum _{m_1\ne m_2\ne m_3}^n\sum _{i=1}^p\mathbb{E}\left[x_{m_{1}i}^2{x}_{m_{2}i}^2{x}_{m_{3}i}^2\right]\hfill \\ \hfill & =\sum _{m_1\ne m_2\ne m_3}^n\sum _{i=1}^p\mathbb{E}\left[x_{m_{1}i}^2\right]\mathbb{E}\left[x_{m_{2}i}^2\right]\mathbb{E}\left[x_{m_{3}i}^2\right]=\sum _{m_1\ne m_2\ne m_3}^n\sum _{i=1}^p{\sigma }_{ii}^3\hfill \\ \hfill & =n(n-1)(n-2)\sum _{i=1}^p{\sigma }_{ii}^3={\tau }_c^{-1}tr\left({\mathbf{D}}_{\mathsf{\Sigma }}^3\right).\hfill \end{array}$$

(60)

Furthermore, we have

$$\begin{array}{cc}\hfill & \mathbb{E}[n^{3}tr\left({\mathbf{D}}_{\mathbf{S}}{\mathbf{S}}^2\right)-2ntr\left(\mathbf{K}\mathbf{S}\right)-ntr\left(\mathbf{G}{\mathbf{D}}_{\mathbf{S}}\mathbf{E}\right)+2tr\left(\mathbf{R}\mathbf{E}\right)]\hfill \\ \hfill & =\sum _{m_1\ne m_2\ne m_3}^n\sum _{i,j=1}^p\mathbb{E}\left[x_{m_{1}i}^2{x}_{m_{2}i}{x}_{m_{2}j}{x}_{m_{3}i}{x}_{m_{3}j}\right]=\sum _{m_1\ne m_2\ne m_3}^n\sum _{i,j=1}^p\mathbb{E}\left[x_{m_{1}i}^2\right]\mathbb{E}\left[x_{m_{2}i}{x}_{m_{2}j}\right]\mathbb{E}\left[x_{m_{3}i}{x}_{m_{3}j}\right]\hfill \\ \hfill & =\sum _{m_1\ne m_2\ne m_3}^n\sum _{i,j=1}^p{\sigma }_{ii}{\sigma }_{ij}^2=n(n-1)(n-2)\sum _{i,j=1}^p{\sigma }_{ii}{\sigma }_{ij}^2={\tau }_c^{-1}tr\left({\mathbf{D}}_{\mathsf{\Sigma }}{\mathsf{\Sigma }}^2\right).\hfill \end{array}$$

(61)

□

Then, the unknown scalars $c_1$ and $c_2$ can be unbiasedly estimated by

$$\begin{array}{cc}\hfill {\gamma }_1& ={\tau }_c[n^{3}tr\left({\mathbf{D}}_{\mathbf{S}}^3\right)-3ntr\left(\mathbf{G}{\mathbf{D}}_{\mathbf{S}}\right)+2tr\left(\mathbf{R}\right)],\hfill \end{array}$$

(62)

$$\begin{array}{cc}\hfill {\gamma }_2& ={\tau }_c[n^{3}tr\left({\mathbf{D}}_{\mathbf{S}}{\mathbf{S}}^2\right)-2ntr\left(\mathbf{K}\mathbf{S}\right)-ntr\left(\mathbf{G}{\mathbf{D}}_{\mathbf{S}}\mathbf{E}\right)+2tr\left(\mathbf{R}\mathbf{E}\right)].\hfill \end{array}$$

(63)

By plugging the estimates of quantities in (27) into (28) and (29), the unbiased estimates of $\mathbf{v}\left({\mathbf{T}}_1\right)$ and $\mathbf{v}\left({\mathbf{T}}_2\right)$ are given by

$$\begin{array}{cc}\hfill \widehat{\mathbf{v}}\left({\mathbf{T}}_1\right)& =\left[\begin{array}{c}\frac{1}{p}{\alpha }_2-{\alpha }_1\\ \frac{np+p-4}{np}{\beta }_1+\frac{1}{p^2}{\beta }_2+\frac{p^2-2np+2}{n{p}^2}{\beta }_3\end{array}\right],\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill \widehat{\mathbf{v}}\left({\mathbf{T}}_2\right)& =\left[\begin{array}{c}{\alpha }_3-{\alpha }_1\\ \frac{n+1}{n}{\beta }_1+\frac{1}{n}{\beta }_3+\frac{n+2}{n}{\gamma }_1-\frac{2n+4}{n}{\gamma }_2\end{array}\right].\hfill \end{array}$$

(65)

By the Equations (16), (64), and (65), the unbiased estimates of the vectors $\mathbf{b}\left({\mathbf{T}}_i\right),i=1,2$ are given by

$$\widehat{\mathbf{b}}\left({\mathbf{T}}_i\right)=\widehat{\mathbf{u}}\left({\mathbf{T}}_i\right)+\widehat{\mathbf{v}}\left({\mathbf{T}}_i\right),i=1,2.$$

(66)

In addition, the constant c in (10) can be further calculated under Gaussian distribution, which is

$$c=\frac{1}{n}{a}_1+\frac{1}{n}{a}_2.$$

(67)

Therefore, we can obtain that the unbiased estimate of c is given by

$$\widehat{c}=\frac{1}{n}{\alpha }_1+\frac{1}{n}{\alpha }_2={\tau }_a\left(\frac{n-2}{n}tr\left({\mathbf{S}}^2\right)+{tr}^2\left(\mathbf{S}\right)\right).$$

(68)

To sum up, by Equations (13), (66) and (68), we can obtain that the unbiased estimates of the loss functions ${\mathcal{M}}_{{\mathbf{T}}_i}\left(\mathbf{w}\right|\mathsf{\Sigma })$ are given by

$${\widehat{\mathcal{M}}}_{{\mathbf{T}}_i}\left(\mathbf{w}\right)=\frac{1}{2}{\mathbf{w}}^T\widehat{\mathbf{H}}\left({\mathbf{T}}_i\right)\mathbf{w}-2{\mathbf{w}}^T\widehat{\mathbf{b}}\left({\mathbf{T}}_i\right)+\widehat{c}.$$

(69)

3.3. Optimal Second-Order Stein-Type Estimators

For the target matrices ${\mathbf{T}}_i,i=1,2$, through replacing the objective function in (10) with its unbiased estimate given by (69), we further formulate the second-order Stein-type estimators as the following optimization problems:

$$\begin{array}{cc}\hfill min& \frac{1}{2}{\mathbf{w}}^T\widehat{\mathbf{H}}\left({\mathbf{T}}_i\right)\mathbf{w}-2{\mathbf{w}}^T\widehat{\mathbf{b}}\left({\mathbf{T}}_i\right)+\widehat{c}\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& 0\le {\mathbf{e}}_1\mathbf{w}\le 1,\hfill \\ \hfill & {\mathbf{e}}_2\mathbf{w}\ge 0.\hfill \end{array}$$

(70)

For $i=1,2$, the Hessian matrix of the objective function in (70) is $\widehat{\mathbf{H}}\left({\mathbf{T}}_i\right)$. By the following inequality

$${\left(tr{(\mathbf{S}-{\mathbf{T}}_i)}^3\right)}^2\le tr{(\mathbf{S}-{\mathbf{T}}_i)}^{2}tr{(\mathbf{S}-{\mathbf{T}}_i)}^4,$$

(71)

we can obtain that $\widehat{\mathbf{H}}\left({\mathbf{T}}_i\right)$ is positive definite. Therefore, the optimization problem (70) is a convex quadratic program. Furthermore, we can obtain the globally optimal solution by an efficient algorithm.

For the target matrices ${\mathbf{T}}_i,i=1,2$, by denoting the corresponding optimal tuning parameters as ${\mathbf{w}}^i={(w_1^i,w_2^i)}^T$, the optimal second-order Stein-type estimators can be expressed as

$${\widehat{\mathsf{\Sigma }}}_i=\mathbf{S}+w_1^i(\mathbf{T}-\mathbf{S})+w_2^i{(\mathbf{T}-\mathbf{S})}^2.$$

(72)

Remark 3.

The proposed second-order Stein-type estimators are both well conditioned. Moreover, by taking the spectral decomposition, the SCM can be expressed as

$$\mathbf{S}={\mathbf{U}}^T\Delta \mathbf{U},$$

(73)

where $\Delta =diag({\delta }_1,\dots ,{\delta }_p)$ is a diagonal matrix consisting of the eigenvalues and $\mathbf{U}$ is the corresponding unitary matrix. Then, the second-order Stein-type estimator ${\widehat{\mathsf{\Sigma }}}_1$ has the following spectral decomposition:

$${\widehat{\mathsf{\Sigma }}}_1={\mathbf{U}}^T\tilde{\Delta }\mathbf{U},$$

(74)

where $\tilde{\Delta }=diag({\tilde{\delta }}_1,\dots ,{\tilde{\delta }}_p)$ with

$${\tilde{\delta }}_i={\delta }_i+w_1^1(\theta -{\delta }_i)+w_2^1{(\theta -{\delta }_i)}^2.$$

(75)

where $\theta =\frac{tr\left(\mathbf{S}\right)}{p}$ is the mean of sample eigenvalues. Therefore, the proposed estimator ${\widehat{\mathsf{\Sigma }}}_1$ shrinks the sample eigenvalues by a nonlinear transformation whileretaining the sample eigenvectors.

4. Numerical Simulations and Real Data Analysis

This section presents numerical simulations and two application examples to discover the performance of the proposed second-order Stein-type estimators. The proposed covariance matrix estimators for the target matrices ${\mathbf{T}}_1$ and ${\mathbf{T}}_2$ are denoted as QS-T1 and QS-T2, respectively. The control estimators include the Stein-type estimator LS-T1 for ${\mathbf{T}}_1$ in [11], and the Stein-type estimator LS-T2 for ${\mathbf{T}}_2$ in [12] and the nonlinear shrinkage estimator NS developed in [18].

4.1. MSE Performance

We assume that the actual distribution is $\mathcal{N}(\mathbf{0},\mathsf{\Sigma })$, where the following models generate the covariance matrix:

(1)

Model 1: $\mathsf{\Sigma }={\left({\sigma }_{ij}\right)}_{p\times p}$ with ${\sigma }_{ii}=1$ and ${\sigma }_{ij}=0.1^{|i-j|}$ for $i\ne j$,

(2)

Model 2: $\mathsf{\Sigma }={\left({\sigma }_{ij}\right)}_{p\times p}$ with ${\sigma }_{ii}=max\{50-i,0\}+2$ and ${\sigma }_{ij}=0.6$ for $i\ne j$.

Under Model 1, the diagonal elements are equal to 1, and the off-diagonal elements are tiny. Therefore, the covariance matrix is close to a spherical matrix. Under Model 2, the diagonal elements are dispersive, and the off-diagonal elements correspond to weak correlations. Therefore, the covariance matrix is close to a diagonal matrix. We carry out random sampling in each Monte Carlo run and compute the Frobenius loss of each covariance matrix estimator. The MSE performance of each covariance matrix estimator is evaluated by averaging the Frobenius losses of $5\times {10}^3$ runs.

Figure 1 and Figure 2 report the logarithmic Frobenius loss of each estimator in large-dimensional scenarios where the dimension is 180 and the sample size varies from 10 to 100. Under Model 1, the Stein-type estimator LS-T1 and the second-order Stein-type estimator QS-T1 outperform the nonlinear shrinkage estimator NS and the shrinkage estimators LS-T2 and QS-T2, which employ the diagonal target matrix. Furthermore, the proposed second-order Stein-type estimator QS-T1 shows a significant advantage over the Stein-type estimator LS-T1, especially when the sample size is tiny. Similarly, under Model 2, the Stein-type estimator LS-T2 and the second-order Stein-type estimator QS-T2 perform better than the other three estimators. Moreover, the proposed estimator QS-T2 outperforms the corresponding Stein-type estimator LS-T2. Therefore, when the correct target matrix is employed, the second-order Stein-type estimators enjoy lower Frobenius losses than the linear and nonlinear shrinkage estimators in large-dimensional scenarios.

Figure 3 and Figure 4 report the logarithmic Frobenius loss of each estimator versus the dimension. The sample size is 30. We can find that the logarithmic Frobenius losses become more prominent as the dimension increases. In Figure 3, the actual covariance matrix is close to being spherical. The proposed second-order Stein-type estimator QS-T1 significantly outperforms the other estimators. In Figure 4, the actual covariance matrix is close to being diagonal. The proposed second-order Stein-type estimator QS-T2 enjoys lower Frobenius loss when the dimension exceeds 120.

The proposed second-order Stein-type estimators take significant advantage of the MSE performance over the linear and nonlinear shrinkage estimators, especially when the dimension is large compared to the sample size.

4.2. Portfolio Selection

In finance, assets with higher expected returns generally involve higher risks. Therefore, investors must constantly balance the expected return and the risk tolerance. The portfolio strategy is a popular way to reduce risk and enhance return. Therefore, portfolio selection plays a vital role in asset investment.

In 1952, Markowitz introduced the famous mean-variance optimization to determine the optimal portfolio weights [31]. Let $\mathbf{m}$ and $\mathsf{\Sigma }$ be the expectation and covariance matrix of the daily returns. For portfolio weight $\mathbf{k}$, the variance of the portfolio is defined as ${\sigma }^2={\mathbf{k}}^T\mathsf{\Sigma }\mathbf{k}$ in the Markowitz framework. As short selling is forbidden, the Markowitz portfolio optimization is formulated as the following mean-variance problem:

$$\begin{array}{cc}\hfill min& {\mathbf{k}}^T\mathsf{\Sigma }\mathbf{k}\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& {\mathbf{k}}^T\mathbf{m}=r,{\mathbf{k}}^T\mathbf{1}=1,\mathbf{k}\ge 0,\hfill \end{array}$$

(76)

where r is a given expected return. By, respectively, replacing $\mathbf{m}$ and $\mathsf{\Sigma }$ with their estimates $\widehat{\mathbf{m}}$ and $\widehat{\mathsf{\Sigma }}$, the optimal weight ${\mathbf{k}}_r$ can be solved by efficient quadratic programming algorithm. It is obvious that the Markowitz optimization only depends on the estimates of the first and second moments of the daily returns. The sample mean and the SCM in the classic portfolio perform well in the portfolio risk measurement; however, the SCM becomes unstable as the number of stocks is large, resulting in significant property loss [32]. Therefore, a well-performed covariance matrix estimator is important in current portfolio selection [2,33].

In practice, we consider a portfolio consisting of $p=95$ highly capitalized stocks from the New York stock exchange with ticker symbols AA, ABT, AIG, AIR, ALL, AMD, AP, APA, AXP, BA, BAC, BAX, BEN, BK, BMY, C, CAT, CCL, CHK, CL, COP, CPE, CVS, CVX, D, DB, DD, DE, DNR, DVN, EAT, EME, EMR, EXC, FCX, FDX, FNMA, GD, GE, GILD, GLW, HAL, HD, HIG, HON, HPQ, IBM, INTC, ITW, JNJ, JPM, KMB, KO, L, LLY, LMT, LOW, M, MCD, MDT, MMM, MO, MRK, MRO, MS, NBR, NC, NE, NEE, NL, NNN, ODC, OXM, OXY, PCG, PEP, PFE, PG, RIG, SLB, SO, T, TGT, TRV, TXN, UNH, UNP, USB, VLO, VZ, WFC, WMT, WWW, X, XOM. The dataset $\mathbf{X}$ contains $n=536$ daily close prices from 11 November 2016 to 31 December 2018 collected via Yahoo! Finance at https://au.finance.yahoo.com/lookup?s=DATA, accessed on 10 November 2022. For each stock i, the daily close price is preprocessed as the daily return by

$$\tilde{\mathbf{X}}(i,j)=\frac{\mathbf{X}(i,j+1)}{\mathbf{X}(i,j)}-1,j=1,\dots ,n-1.$$

(77)

The covariance matrix estimators LS-T1, LS-T2, NS, QS-T1, and QS-T2 are generated from the daily return $\tilde{\mathbf{X}}$. For an expected return r, the realized risk is defined as ${\sigma }_r=\sqrt{{\mathbf{k}}_r^T\widehat{\mathsf{\Sigma }}{\mathbf{k}}_r}$. Next, we employ the realized risks, one key index to evaluate the portfolio, to verify the performance of the covariance matrix estimators.

Figure 5 and Figure 6 plot the realized risks of three kinds of shrinkage estimators for different investment horizons. For a short-term investment of 44 trading days (2 months), the proposed QS-T2 has the lowest realized risk when the expected return is less than 0.5% and has the highest realized risk when the expected return exceeds 0.6%. The proposed QS-T1 and the Stein-type estimator LS-T1 have the lowest realized risk when the expected return exceeds 0.5%. For a long-term investment of 280 trading days (13 months), the expected return of the five estimators becomes similar for the short-term investment.

Figure 7 and Figure 8 plot the realized risks of shrinkage estimators for different expected return levels. The proposed estimator QS-T2 enjoys the lowest realized risk for a low expected return. The nonlinear shrinkage estimator NS has the highest realized risk. However, the proposed estimator QS-T2 performs worst when the expected return becomes high. The proposed estimator QS-T1, together with NS and LS-T1, performs best in this scenario.

The proposed second-order Stein-type estimator QS-T2 enjoys good portfolio selection for short-term investment and prudent return cases. The proposed second-order Stein-type estimator QS-T1 is recommended in long-term investment and high-return scenarios.

4.3. Discriminant Analysis

We further discover the performance of the second-order regularized estimators in small-sample-size situations. The Parkinson’s data are collected on the website https://archive-beta.ics.uci.edu/, accessed on 10 November 2022. $p=160$ biomedical voice attitudes are measured from $n_1$ patients and $n_2$ healthy individuals. Let $\widehat{\mathsf{\Sigma }}$ be the pooling estimator based on a certain estimation strategy. We use the following quadratic discriminant rule $M={(x_i-\overline{x})}^T{\widehat{\mathsf{\Sigma }}}^{-1}(x_i-\overline{x})$ to make a diagnosis for each $x_i$. M denotes the Mahalanobis distance between individual $x_i$ and the sample center $\overline{x}$. The individual $x_i$ is classified as Parkinson’s patient if $x_i$ is closer to the sample center of patients in the sense of Mahalanobis distance.

Table 1. The classification accuracy rate of Parkinson’s data.

${\mathit{n}}_2$	15	20	25	30	35	40	45	50
LS-T1	0.5165	0.5530	0.5646	0.5826	0.5870	0.5956	0.5981	0.6054
LS-T2	0.6262	0.6459	0.6607	0.6714	0.6796	0.6892	0.6918	0.6971
NS	0.4931	0.4842	0.4701	0.4688	0.4544	0.4458	0.4285	0.4304
QS-T1	0.5124	0.5465	0.5600	0.5761	0.5833	0.5938	0.5953	0.6025
QS-T2	0.6336	0.6539	0.6678	0.6747	0.6840	0.6902	0.6916	0.6957

reports the classification accuracy rates of different estimators. We can see that the Stein-type estimators LS-T1, LS-T2, QS-T1, and LS-T2 perform better than NS. Moreover, the accuracy rate of the Stein-type estimators becomes larger as the sample size increases. Moreover, the proposed second-order Stein-type estimator QS-T2 enjoys the largest accuracy rate when $n\le 40$, and the Stein-type estimator LS-T2 has the best performance when $n\ge 45$. Therefore, we can see that the proposed second-order Stein-type estimation performs better than the classic Stein-type estimation in this application.

5. Conclusions and Discussion

This paper investigated the problem of estimating a large-dimensional covariance matrix. Motivated by Stein’s strategy, we developed a novel strategy named the second-order Stein-type estimation. The proposed estimator is expected to be positive definite in the form of a quadratic binomial of the SCM and the target matrix. Firstly, we specified the spherical and diagonal targets in the second-order Stein-type regularization. The mean squared errors were, respectively, obtained for the two targets. Secondly, we unbiasedly estimated the two mean squared errors under the Gaussian distribution. Thirdly, the optimal parameters were obtained by solving the convex quadratic programming. The optimal second-order Stein-type estimators were obtained for the two target matrices. Finally, we verified the performance of the proposed estimators in numerical simulations and real data applications.

It is worth mentioning that the second-order Stein-type estimators were proposed under Gaussian distribution. In practical applications, the data may often deviate from the Gaussian distribution. Therefore, the problem of investigating the second-order Stein-type regularization under non-Gaussian distributions remains open and important.