Estimation of the Kronecker Covariance Model by Quadratic Form

TL;DR

None

Contribution

None

Abstract

We propose a new estimator, the quadratic form estimator, of the Kronecker product model for covariance matrices. We show that this estimator has good properties in the large dimensional case (i.e., the cross-sectional dimension $n$ is large relative to the sample size $T$). In particular, the quadratic form estimator is consistent in a relative Frobenius norm sense provided $\log^3n/T\to 0$. We obtain the limiting distributions of Lagrange multiplier (LM) and Wald tests under both the null and local alternatives concerning the mean vector $\mu$. Testing linear restrictions of $\mu$ is also investigated. Finally, our methodology performs well in the finite-sample situations both when the Kronecker product model is true, and when it is not true.

Tables

Table 1. Table 1: M T subscript 𝑀 𝑇 M_{T} , Σ ~ ~ Σ \tilde{\Sigma} , LW04 and LW17 stand for the sample covariance matrix, quadratic form estimator, Ledoit and Wolf ( 2004 ) ’s linear shrinkage estimator, and Ledoit and Wolf ( 2017 ) ’s direct nonlinear shrinkage estimator, respectively. MSE-1 and MSE-2 are the MSE in terms of Σ Σ \Sigma and Σ − 1 superscript Σ 1 \Sigma^{-1} , respectively. PRIAL-1 is the PRIAL in terms of Σ Σ \Sigma . T = 252 𝑇 252 T=252 and n = 2 10 = 1024 𝑛 superscript 2 10 1024 n=2^{10}=1024 . 0.000 means less than 0.001.

	$M_T$	$\tilde{\Sigma }$	LW04	LW17
	$\rho =0.5$
MSE-1	2.989	0.000	0.242	0.243
MSE-2	NA	0.000	0.311	0.308
PRIAL-1	0	1.000	0.919	0.919
size of LM	NA	0.051	1.000	1.000
size of Wald	NA	0.050	0.085	0.093
	$\rho =0.7$
MSE-1	1.760	0.000	0.429	0.430
MSE-2	NA	0.000	0.722	0.715
PRIAL-1	0	1.000	0.756	0.756
size of LM	NA	0.050	1.000	1.000
size of Wald	NA	0.051	0.158	0.164
	$\rho =0.85$
MSE-1	0.501	0.001	0.320	0.316
MSE-2	NA	0.002	0.980	0.980
PRIAL-1	0	0.998	0.360	0.370
size of LM	NA	0.051	1.000	1.000
size of Wald	NA	0.060	0.334	0.329

Table 2. Table 2: M T subscript 𝑀 𝑇 M_{T} , Σ ~ ~ Σ \tilde{\Sigma} , LW04 and LW17 stand for the sample covariance matrix, quadratic form estimator, Ledoit and Wolf ( 2004 ) ’s linear shrinkage estimator, and Ledoit and Wolf ( 2017 ) ’s direct nonlinear shrinkage estimator, respectively. MSE-1 and MSE-2 are the MSE in terms of Σ Σ \Sigma and Σ − 1 superscript Σ 1 \Sigma^{-1} , respectively. PRIAL-1 is the PRIAL in terms of Σ Σ \Sigma . ρ = 0.7 𝜌 0.7 \rho=0.7 . 0.000 means less than 0.001.

	$M_T$	$\tilde{\Sigma }$	LW04	LW17		$M_T$	$\tilde{\Sigma }$	LW04	LW17
	$n=2^9$, $T=252$					$n=2^9$, $T=504$
MSE-1	0.882	0.001	0.346	0.345		0.442	0.000	0.249	0.246
MSE-2	NA	0.001	0.676	0.656		NA	0.000	0.601	0.531
PRIAL-1	0	0.999	0.608	0.609		0	0.999	0.438	0.443
size of LM	NA	0.039	1.000	1.000		NA	0.053	1.000	1.000
size of Wald	NA	0.041	0.153	0.149		NA	0.058	0.148	0.151
	$n=2^{10}$, $T=252$					$n=2^{10}$, $T=504$
MSE-1	1.760	0.000	0.429	0.430		0.882	0.000	0.345	0.344
MSE-2	NA	0.000	0.722	0.715		NA	0.000	0.677	0.659
PRIAL-1	0	1.000	0.756	0.756		0	1.000	0.608	0.610
size of LM	NA	0.050	1.000	1.000		NA	0.059	1.000	1.000
size of Wald	NA	0.051	0.158	0.164		NA	0.062	0.168	0.168
	$n=2^{11}$, $T=252$					$n=2^{11}$, $T=504$
MSE-1	3.514	0.000	0.489	0.490		1.760	0.000	0.429	0.429
MSE-2	NA	0.000	0.747	0.744		NA	0.000	0.723	0.717
PRIAL-1	0	1.000	0.861	0.861		0	1.000	0.756	0.757
size of LM	NA	0.057	1.000	1.000		NA	0.057	1.000	1.000
size of Wald	NA	0.067	0.202	0.221		NA	0.060	0.169	0.181

Table 3. Table 3: M T subscript 𝑀 𝑇 M_{T} , Σ ~ ~ Σ \tilde{\Sigma} , LW04 and LW17 stand for the sample covariance matrix, the quadratic form estimator, the Ledoit and Wolf ( 2004 ) ’s linear shrinkage estimator, and the Ledoit and Wolf ( 2017 ) ’s direct nonlinear shrinkage estimator, respectively. T = 252 𝑇 252 T=252 . Powers are not size-adjusted.

	$M_T$	$\tilde{\Sigma }$	LW04	LW17		$M_T$	$\tilde{\Sigma }$	LW04	LW17
	$\rho =0.5$, $n=2^9$					$\rho =0.5$, $n=2^{10}$
power of LM	NA	0.890	0.730	0.866		NA	0.925	0.999	1.000
power of Wald	NA	0.905	0.689	0.734		NA	0.942	0.750	0.780
	$\rho =0.7$, $n=2^9$					$\rho =0.7$, $n=2^{10}$
power of LM	NA	1.000	1.000	1.000		NA	1.000	1.000	1.000
power of Wald	NA	1.000	0.742	0.833		NA	1.000	0.746	0.806
	$\rho =0.85$, $n=2^9$					$\rho =0.85$, $n=2^{10}$
power of LM	NA	1.000	1.000	1.000		NA	1.000	1.000	1.000
power of Wald	NA	1.000	0.990	1.000		NA	1.000	0.981	0.977

Table 4. Table 4: M T subscript 𝑀 𝑇 M_{T} , Σ ~ ~ Σ \tilde{\Sigma} , LW04 and LW17 stand for the sample covariance matrix, quadratic form estimator (factorisations given in parentheses), Ledoit and Wolf ( 2004 ) ’s linear shrinkage estimator, and Ledoit and Wolf ( 2017 ) ’s direct nonlinear shrinkage estimator, respectively. MSE-1 and MSE-2 are the MSE in terms of Σ Σ \Sigma and Σ − 1 superscript Σ 1 \Sigma^{-1} , respectively. PRIAL-1 and PRIAL-2 are the PRIAL in terms of Σ Σ \Sigma and Σ − 1 superscript Σ 1 \Sigma^{-1} , respectively. n = 20 , T = 40 , α 2 = 0.5 formulae-sequence 𝑛 20 formulae-sequence 𝑇 40 superscript 𝛼 2 0.5 n=20,T=40,\alpha^{2}=0.5 .

	$M_T$	$\tilde{\Sigma }$	$\tilde{\Sigma }$	$\tilde{\Sigma }$	$\tilde{\Sigma }$
		$\left(5\times 2\times 2\right)$	$\left(2\times 5\times 2\right)$	$\left(2\times 2\times 5\right)$	$\left(4\times 5\right)$
MSE-1	0.446	0.137	0.136	0.137	0.140
MSE-2	6.876	0.154	0.153	0.154	0.163
PRIAL-1	0	0.684	0.685	0.682	0.675
PRIAL-2	0	0.977	0.977	0.977	0.976
size of LM	0.004	0.043	0.050	0.038	0.038
size of Wald	0.690	0.092	0.087	0.081	0.094
	$\tilde{\Sigma }$	$\tilde{\Sigma }$	$\tilde{\Sigma }$	LW04	LW17
	$\left(5\times 4\right)$	$\left(2\times 10\right)$	$\left(10\times 2\right)$
MSE-1	0.139	0.189	0.188	0.113	0.129
MSE-2	0.163	0.293	0.288	0.122	0.148
PRIAL-1	0.679	0.570	0.571	0.738	0.702
PRIAL-2	0.976	0.957	0.958	0.982	0.978
size of LM	0.041	0.035	0.028	0.022	0.015
size of Wald	0.100	0.163	0.167	0.074	0.083

Table 5. Table 5: M T subscript 𝑀 𝑇 M_{T} , Σ ~ ~ Σ \tilde{\Sigma} , LW04 and LW17 stand for the sample covariance matrix, quadratic form estimator (factorisations given in parentheses), Ledoit and Wolf ( 2004 ) ’s linear shrinkage estimator, and Ledoit and Wolf ( 2017 ) ’s direct nonlinear shrinkage estimator, respectively. MSE-1 and MSE-2 are the MSE in terms of Σ Σ \Sigma and Σ − 1 superscript Σ 1 \Sigma^{-1} , respectively. PRIAL-1 and PRIAL-2 are the PRIAL in terms of Σ Σ \Sigma and Σ − 1 superscript Σ 1 \Sigma^{-1} , respectively. n = 20 , T = 40 formulae-sequence 𝑛 20 𝑇 40 n=20,T=40 .

	$M_T$	$\tilde{\Sigma }$	LW04	LW17		$M_T$	$\tilde{\Sigma }$	LW04	LW17
		$\left(2\times 5\times 2\right)$					$\left(2\times 5\times 2\right)$
	${\alpha }^2=0.25$					${\alpha }^2=0.50$
MSE-1	0.492	0.077	0.050	0.070		0.446	0.136	0.113	0.129
MSE-2	7.405	0.089	0.048	0.086		6.876	0.153	0.122	0.148
PRIAL-1	0	0.843	0.898	0.856		0	0.685	0.738	0.702
PRIAL-2	0	0.988	0.993	0.988		0	0.977	0.982	0.978
size of LM	0.004	0.042	0.035	0.020		0.004	0.050	0.022	0.015
size of Wald	0.690	0.083	0.064	0.066		0.690	0.087	0.074	0.083
	$M_T$	$\tilde{\Sigma }$	LW04	LW17		$M_T$	$\tilde{\Sigma }$	LW04	LW17
		$\left(2\times 5\times 2\right)$					$\left(2\times 5\times 2\right)$
	${\alpha }^2=0.75$					${\alpha }^2=1$
MSE-1	0.396	0.195	0.154	0.167		0.353	0.243	0.173	0.184
MSE-2	6.311	0.241	0.194	0.204		5.807	0.335	0.259	0.246
PRIAL-1	0	0.469	0.589	0.557		0	0.218	0.469	0.440
PRIAL-2	0	0.959	0.966	0.966		0	0.934	0.948	0.953
size of LM	0.004	0.058	0.017	0.013		0.004	0.067	0.017	0.013
size of Wald	0.690	0.091	0.087	0.093		0.690	0.106	0.090	0.091

Table 6. Table 6: M T subscript 𝑀 𝑇 M_{T} , Σ ~ ~ Σ \tilde{\Sigma} , LW04 and LW17 stand for the sample covariance matrix, quadratic form estimator (factorisations given in parentheses), Ledoit and Wolf ( 2004 ) ’s linear shrinkage estimator, and Ledoit and Wolf ( 2017 ) ’s direct nonlinear shrinkage estimator, respectively. MSE-1 and MSE-2 are the MSE in terms of Σ Σ \Sigma and Σ − 1 superscript Σ 1 \Sigma^{-1} , respectively. PRIAL-1 and PRIAL-2 are the PRIAL in terms of Σ Σ \Sigma and Σ − 1 superscript Σ 1 \Sigma^{-1} , respectively. α 2 = 0.5 superscript 𝛼 2 0.5 \alpha^{2}=0.5 .

$n/T=0.32$	$M_T$	$\tilde{\Sigma }$	$\tilde{\Sigma }$	$\tilde{\Sigma }$	$\tilde{\Sigma }$	LW04	LW17
		$\left(2\times 2\times 2\times 2\right)$	$\left(4\times 4\right)$	$\left(4\times 2\times 2\right)$	$\left(2\times 8\right)$
MSE-1	0.292	0.118	0.122	0.120	0.145	0.098	0.109
MSE-2	1.491	0.134	0.142	0.137	0.190	0.110	0.118
PRIAL-1	0	0.580	0.571	0.576	0.492	0.655	0.618
PRIAL-2	0	0.907	0.902	0.905	0.870	0.924	0.919
size of LM	0.013	0.057	0.050	0.050	0.041	0.023	0.019
size of Wald	0.373	0.081	0.090	0.080	0.133	0.072	0.074
$n/T=2$	$M_T$	$\tilde{\Sigma }$	$\tilde{\Sigma }$	$\tilde{\Sigma }$	$\tilde{\Sigma }$	LW04	LW17
		$\left(5\times 2\times 2\times 2\right)$	$\left(5\times 2\times 4\right)$	$\left(5\times 8\right)$	$\left(10\times 2\times 2\right)$
MSE-1	1.684	0.168	0.175	0.216	0.234	0.159	0.196
MSE-2	NA	0.182	0.194	0.286	0.337	0.151	0.164
PRIAL-1	0	0.898	0.894	0.870	0.860	0.904	0.882
PRIAL-2	NA	NA	NA	NA	NA	NA	NA
size of LM	NA	0.051	0.054	0.050	0.049	0.051	0.070
size of Wald	NA	0.155	0.159	0.224	0.260	0.129	0.140

Equations

$$\Sigma ={\sigma }^2\times {\Sigma }_1\otimes \cdots \otimes {\Sigma }_v,$$

$$w_{i,j,t}={\mu }_{i,j}+{\gamma }_{i,t}{f}_{j,t},i=1,2,3,j=1,\dots ,50,t=1,\dots ,10$$

$$\mathbb{E}[(y_t-\mu )(y_t-\mu {)}^{⊺}]=\mathbb{E}[{\gamma }_t{\gamma }_t^{⊺}]\otimes \mathbb{E}[f_t{f}_t^{⊺}].$$

$$y_t\text{mathchar}58=(w_{1,1,\dots ,1,t},\dots ,w_{n_1,n_2,\dots ,n_v,t}{)}^{⊺}=\mu +{\epsilon }_t^1\otimes {\epsilon }_t^2\otimes \cdots \otimes {\epsilon }_t^v,$$

$$\mathbb{E}[(y_t-\mu )(y_t-\mu {)}^{⊺}]=\mathbb{E}[{\epsilon }_t^1{\epsilon }_t^{1⊺}]\otimes \mathbb{E}[{\epsilon }_t^2{\epsilon }_t^{2⊺}]\otimes \cdots \otimes \mathbb{E}[{\epsilon }_t^v{\epsilon }_t^{v⊺}].$$

$$\mathrm{tr}(\Sigma )={\sigma }^2\times \mathrm{tr}({\Sigma }_1\otimes \cdots \otimes {\Sigma }_v)={\sigma }^2\times \mathrm{tr}({\Sigma }_1)\times \cdots \times \mathrm{tr}({\Sigma }_v)={\sigma }^{2}n,$$

$$A=\left(\begin{array}[]{ccc}A_{-1;1,1}&\cdots&A_{-1;1,n_{1}}\\ &\ddots&\vdots\\ &&A_{-1;n_{1},n_{1}}\end{array}\right).$$

$$\mathrm{PTR}_{n_{1}}(A)=\left(\begin{array}[]{ccc}\mathrm{tr}(A_{-1;1,1})&\cdots&\mathrm{tr}(A_{-1;1,n_{1}})\\ &\ddots&\vdots\\ &&\mathrm{tr}(A_{-1;n_{1},n_{1}})\end{array}\right).$$

$$\Sigma_{-h}\mathrel{\mathop{\mathchar 58\relax}}=\left\{\begin{array}[]{ll}\Sigma_{h+1}\otimes\cdots\otimes\Sigma_{v}\otimes\Sigma_{1}\otimes\cdots\otimes\Sigma_{h-1}&\text{ for }h=2,\ldots,v-1\\ \Sigma_{1}\otimes\cdots\otimes\Sigma_{v-1}&\text{ for }h=v.\end{array}\right..$$

$${\Sigma }^{(h)}$$

$$=K_{n_h\times \cdots \times n_v,n_1\times \cdots \times n_{h-1}}({\sigma }^2\times {\Sigma }_1\otimes \cdots \otimes {\Sigma }_v)K_{n_1\times \cdots \times n_{h-1},n_h\times \cdots \times n_v}$$

$$={\sigma }^2\times {\Sigma }_h\otimes {\Sigma }_{h+1}\otimes \cdots \otimes {\Sigma }_v\otimes {\Sigma }_1\otimes \cdots \otimes {\Sigma }_{h-1}={\sigma }^2\times {\Sigma }_h\otimes {\Sigma }_{-h}.$$

$${\Sigma }_h=\frac{d^{(h)}}{\mathrm{tr}(d^{(h)})/n_h}.$$

$$M_T\text{mathchar}58=\frac{1}{T}\sum _{t=1}^T(y_t-\bar{y})(y_t-\bar{y}{)}^{⊺},$$

$$M_T^{(h)}\text{mathchar}58=K_{n_h\times \cdots \times n_v,n_1\times \cdots \times n_{h-1}}{M}_T{K}_{n_1\times \cdots \times n_{h-1},n_h\times \cdots \times n_v},$$

$${\tilde{\Sigma }}_h\text{mathchar}58=\frac{{\hat{d}}^{(h)}}{\mathrm{tr}({\hat{d}}^{(h)})/n_h},$$

$$\tilde{\Sigma }={\hat{\sigma }}^2\times {\tilde{\Sigma }}_1\otimes \cdots \otimes {\tilde{\Sigma }}_v,$$

$${\hat{\sigma }}^2\text{mathchar}58=\frac{\mathrm{tr}(M_T)}{n}.$$

$$\underset{1\le i\le n}{\max }\frac{1}{T}\sum _{t=1}^T\mathbb{E}|y_{t,i}{|}^m\le A^m,m=2,3,\dots ,$$

$$\operatorname*{var}\mathinner{\bigl{(}(y_{t}-\mu)\otimes(y_{t}-\mu)\bigr{)}}\leq C\operatorname*{var}\mathinner{\bigl{(}(z_{t}-\mu)\otimes(z_{t}-\mu)\bigr{)}},$$

$$\lim_{n\to\infty}\frac{1}{n^{\beta_{1}}}\|\Sigma\|_{F}^{2}=\lim_{n\to\infty}\frac{\sigma^{4}}{n^{\beta_{1}}}\mathinner{\biggl{(}\prod_{j=1}^{v}\|\Sigma_{j}\|_{F}^{2}\biggr{)}}=\omega<\infty.$$

$$\frac{1}{n}\Vert \Sigma {\Vert }_F^2$$

$$\frac{\mathinner{\!\bigl{\lVert}\tilde{\Sigma}-\Sigma\bigr{\rVert}}_{F}}{\|\Sigma\|_{F}}=O_{p}\mathinner{\biggl{(}\sqrt{\frac{\log^{3}n}{n^{2-\beta_{1}}T}}\biggr{)}}+O_{p}\mathinner{\biggl{(}\frac{\log^{2}n}{T}\biggr{)}}.$$

$$\frac{\mathinner{\!\bigl{\lVert}\tilde{\Sigma}^{-1}-\Sigma^{-1}\bigr{\rVert}}_{F}}{\|\Sigma^{-1}\|_{F}}=O_{p}\mathinner{\biggl{(}\sqrt{\frac{\log^{3}n}{n^{2-\beta_{1}}T}}\biggr{)}}+O_{p}\mathinner{\biggl{(}\frac{\log^{2}n}{T}\biggr{)}}.$$

$$\frac{\mathinner{\!\bigl{\lVert}\tilde{\Sigma}-\Sigma\bigr{\rVert}}_{1}}{\|\Sigma\|_{1}}=O_{p}\mathinner{\biggl{(}\sqrt{\frac{\log^{3}n}{n^{2-\beta_{1}}T}}\biggr{)}}+O_{p}\mathinner{\biggl{(}\frac{\log^{2}n}{T}\biggr{)}}.$$

$$\frac{\mathinner{\!\bigl{\lVert}\tilde{\Sigma}^{-1}-\Sigma^{-1}\bigr{\rVert}}_{1}}{\|\Sigma^{-1}\|_{1}}=O_{p}\mathinner{\biggl{(}\sqrt{\frac{\log^{3}n}{n^{2-\beta_{1}}T}}\biggr{)}}+O_{p}\mathinner{\biggl{(}\frac{\log^{2}n}{T}\biggr{)}}.$$

$$\frac{\mathinner{\!\bigl{\lVert}\tilde{\Sigma}-\Sigma\bigr{\rVert}}_{\ell_{2}}}{\|\Sigma\|_{\ell_{2}}}=O_{p}\mathinner{\biggl{(}\sqrt{\frac{\log^{3}n}{n^{2-\beta_{1}}T}}\biggr{)}}+O_{p}\mathinner{\biggl{(}\frac{\log^{2}n}{T}\biggr{)}}.$$

$$\frac{\mathinner{\!\bigl{\lVert}\tilde{\Sigma}^{-1}-\Sigma^{-1}\bigr{\rVert}}_{\ell_{2}}}{\|\Sigma^{-1}\|_{\ell_{2}}}=O_{p}\mathinner{\biggl{(}\sqrt{\frac{\log^{3}n}{n^{2-\beta_{1}}T}}\biggr{)}}+O_{p}\mathinner{\biggl{(}\frac{\log^{2}n}{T}\biggr{)}}.$$

$$L{M}_{n,T}=T(\bar{y}-{\mu }_0{)}^{⊺}{\tilde{\Sigma }}_{{\mu }_0}^{-1}(\bar{y}-{\mu }_0),$$

$$W_{n,T}=T(\bar{y}-{\mu }_0{)}^{⊺}{\tilde{\Sigma }}^{-1}(\bar{y}-{\mu }_0),$$

$$\lim_{n\to\infty}\frac{1}{n^{\beta_{2}}}\|\Sigma^{-1}\|_{1}=\lim_{n\to\infty}\frac{1}{n^{\beta_{2}}\sigma^{2}}\mathinner{\biggl{(}\prod_{j=1}^{v}\|\Sigma_{j}^{-1}\|_{1}\biggr{)}}=\omega^{\prime}<\infty.$$

Full text

Some *key words: *Covariance matrix; Kronecker product; Quadratic form; Lagrange multiplier test; Wald test

1 Introduction

Covariance matrices are of great importance in many fields. In finance, they are a key element in portfolio choice and risk management (Markowitz (1952)). In psychology, scholars have long assumed that some observed variables are related to certain latent traits through a factor model, and then use the covariance matrix of the observed variables to deduce properties of the latent traits. In econometrics, covariance matrices often appear in test statistics representing the sampling variability of a vector of parameter estimates. Anderson (1984) is a classic statistical reference that studies estimation of and hypothesis testing about covariance matrices in the low dimensional case (i.e., the dimension of the covariance matrix, $n$, is small compared with the sample size $T$).

There are many new methodological approaches to covariance and precision matrix estimation in the large dimensional case (i.e., $n$ is large compared with $T$);111Some studies have made a distinction between the large dimensional case and the high dimensional case (Hafner et al. (2020)). We no longer make this distinction in this article. As long as $n$ is large relative to $T$, regardless of $n$ exceeding $T$, we call it the large dimensional case. see, e.g., Ledoit and Wolf (2003), Bickel and Levina (2008), Fan et al. (2008), Ledoit and Wolf (2012), Fan et al. (2013), and Ledoit and Wolf (2015). Fan et al. (2016) gave an excellent account of the recent developments in theory and practice of estimating large dimensional covariance matrices. The usual approaches include: to impose some sparsity on the covariance matrix, meaning that many elements of the covariance matrix are assumed to be zero or small, thereby reducing the number of parameters to be estimated; or at least to ”shrink” towards a sparse matrix, or to use a factor model which reduces the dimensionality of the parameter space. Most of this literature assumes i.i.d. data.

We consider the problem of estimating a large covariance matrix $\Sigma.$ We impose a model structure that reduces the effective dimensionality. In particular, we consider the Kronecker product model. Let $n=n_{1}\times\cdots\times n_{v}$, where $n_{j}\in\mathbb{Z}$ and $n_{j}\geq 2$ for $j=1,\ldots,v$. We suppose that

[TABLE]

where $\Sigma_{j}$ is an $n_{j}\times n_{j}$ unknown covariance matrix satisfying $\mathrm{tr}(\Sigma_{j})=n_{j}$ for $j=1,\ldots,v,$ and $0<\sigma^{2}<\infty$ is a scalar parameter.

Kronecker product models arise naturally from multiway data (c.f. Kroonenberg (2008)). Multiway data are a generalization of two-way or three-way data that are widely encountered in social science. For example, the scores on 3 subjects (mathematics, English and music) of 50 students observed over 10 years are three-way data, the ”ways” being subjects, students and years. Let $w_{i,j,t}$ denote the score of subject $i$ of student $j$ in year $t$. To model $w_{i,j,t}$, one could use an interactive effects model similar to Bai (2009):

[TABLE]

where $\mu_{i,j}$ is the subject-student specific mean, while $\gamma_{i,t}$ and $f_{j,t}$ are the subject-time specific and student-time specific effects, respectively. Stacking all the observations $\{w_{i,j,t}\}$ of year $t$ into a $150\times 1$ column vector $y_{t}$, we have $y_{t}=\mu+\gamma_{t}\otimes f_{t}$, where $\mu$ is the $150\times 1$ mean vector containing stacked $\{\mu_{i,j}\}$, $\gamma_{t}=(\gamma_{1,t},\gamma_{2,t},\gamma_{3,t})^{\intercal}$, and $f_{t}=(f_{1,t},\ldots,f_{50,t})^{\intercal}$. Suppose that $\gamma_{t}$ is a random vector independent of $f_{t}$, and that both are mean-zero and stationary in time. Then,

[TABLE]

In this case the covariance matrix of $y_{t}$ is a Kronecker product of two sub-matrices, which describe the subject specific and individual specific dependencies.

Extending the idea to multiway data, one might think of a typical equity portfolio constructed by intersections of 5 size quintiles, 5 book-to-market ratio quintiles, and 10 industries, in the spirit of Fama and French (1993), over a number of years, as four-way data: sizes $\times$ B/P ratios $\times$ industries $\times$ years. Situations in which higher-way data are collected are also on the increase. For example, electroencephalography (EEG), a non-invasive way of detecting structural abnormalities such as brain tumors, also provide multiway data, such as EEG bands $\times$ patients $\times$ leads $\times$ doses $\times$ time $\times$ task conditions (Estienne et al. (2001)).

Consider $(v+1)$-way data $w_{i_{1},i_{2},\ldots,i_{v},t}$, where $i_{j}=1,\ldots,n_{j}$ for $j=1,\ldots,v$ and $t=1,\ldots,T$. We use subscript $t$ to denote the $(v+1)$th way of the data in the hope to broadly interpret the $(v+1)$th way as ”time”, $T$ as the sample size, all other ways as the ”cross-section”, and $n\mathrel{\mathop{\mathchar 58\relax}}=n_{1}\times\cdots\times n_{v}$ as the cross-sectional dimension. In other words, the $(v+1)$th way of the data need not correspond to the time dimension, should the multiway data contain such a dimension. In the rest of the article, we shall no longer stress this distinction. Suppose that $w_{i_{1},i_{2},\ldots,i_{v},t}=\mu_{i_{1},i_{2},\ldots,i_{v}}+\varepsilon_{i_{1},t}^{1}\varepsilon_{i_{2},t}^{2}\cdots\varepsilon_{i_{v},t}^{v},$ where $i_{j}=1,\ldots,n_{j}$ for $j=1,\ldots,v,$ and $t=1,\ldots,T$. Equivalently, in the stacked form

[TABLE]

where $\mu$ is the stacked mean vector, $\varepsilon_{t}^{j}\mathrel{\mathop{\mathchar 58\relax}}=(\varepsilon_{1,t}^{j},\ldots,\varepsilon_{n_{j},t}^{j})^{\intercal}$ is an $n_{j}\times 1$ mean-zero random vector with covariance matrix $\mathbb{E}[\varepsilon_{t}^{j}\varepsilon_{t}^{j\intercal}]$ for all $t$ for $j=1,\ldots,v$. If $\varepsilon_{t}^{1},\ldots,\varepsilon_{t}^{v}$ are mutually independent for all $t$, then

[TABLE]

We hence see that the covariance matrix of $y_{t}$ is a Kronecker product of $v$ sub-matrices.

Recent work on Kronecker product models for multiway data include Hoff (2011), Hoff (2015), Hoff (2016) etc. Kronecker product models have also been considered in the psychometric literature (Campbell and O’Connell (1967), Swain (1975), Cudeck (1988), Verhees and Wansbeek (1990) etc). In the spatial literature, there are a number of studies that consider Kronecker product models for the correlation matrix of a random field (Loh and Lam (2000)). Robinson (1998) and Hidalgo and Schafgans (2017) exploited separable error covariance matrix structures to develop inference methods without the need for smoothing.

These literatures have all focussed on the low dimensional case. Hafner et al. (2020) were the first to study Kronecker product models in the large dimensional case. The proper framework for studying the large dimensional case is the joint limit setting developed by Phillips and Moon (1999) in which $n$ and $T$ tend to infinity simultaneously.222Peter Phillips has made some fundamental contributions to large dimensional analysis. Phillips and Moon (1999) provided three asymptotic frameworks for analysing double-index ($n,T$) processes: sequential limit framework (e.g., $n\rightarrow\infty$ followed by $T\rightarrow\infty$), diagonal path limit framework (i.e., both $n$ and $T$ pass to infinity along some specific diagonal in the two dimensional array), and joint limit framework (i.e., $n,T\rightarrow\infty$ simultaneously). In particular, they provided a central limit theorem in joint limit framework for double-index processes (Theorem 2 of Phillips and Moon (1999)). However, the Lindeberg condition of that theorem is perhaps difficult to verify in practice. In Section B, we provide a variant (Theorem B.1), which relies on a Lyapounov’s condition. Moreover, the variant allows the central limit theorem to kick in from either the cross-sectional or time dimension. Since $n$ tends to infinity, there are two main cases when considering (1.1): (a) $\{n_{j}\}_{j=1}^{v}$ are all fixed while $v\rightarrow\infty$; (b) $n_{j}\rightarrow\infty$ for at least some $j$ while $v$ is fixed. Case (a) corresponds to practical situations where the data have a large number of ways but in each way the number of entities is small; case (b) often corresponds to, say, three-way or four-way data in which at least one way has a large number of entities. The methodologies developed in Hafner et al. (2020) and this article are perfectly geared for case (a) in the sense that (1.1) is correctly specified for the data.

We do not analyse case (b) theoretically, but our estimation and inference procedures can in principle be applied to case (b) also, but the theory will require more work and stronger restrictions on the relationship between $n$ and $T.$ For example, if $v=2$ and $n_{1}=n_{2}=\sqrt{n},$ then the sub-matrices $\Sigma_{1},\Sigma_{2}$ each contain order $n$ unknown quantities. If $n/T\rightarrow 0$ fast enough, then we may show some consistency of our estimators of the sub-matrices $\Sigma_{1},\Sigma_{2}.$ On the other hand, if this rate condition is not satisfied, one could combine the separable structure (i.e., the Kronecker product) with sparsity restrictions on the sub-matrices. This has been investigated in the literature. Other approaches have been considered in Akdemir and Gupta (2011), Hoff (2011), and Hoff (2015). Henceforth, when we say the Kronecker product model (1.1), we implicitly mean case (a).

The Kronecker product model leads to substantial dimension reduction even though it need not be sparse in the sense of (2.1) of Fan et al. (2016). Hafner et al. (2020) showed that the matrix logarithm of a Kronecker product covariance or correlation matrix is a sparse matrix (with $O(\log n)$ unknown quantities) and the logarithmic operator converts the multiplicative Kronecker product structure into an additive one. Therefore, the logarithm of a Kronecker product covariance or correlation matrix is a linear function of a much ”smaller” vector of unknown quantities. They used this to develop a closed-form estimator, they established its consistency and provided a central limit theorem (CLT). However, their results require strong, albeit sufficient but not necessary, conditions; in particular they obtained Frobenius norm consistency of the estimator under a condition that at least $n/T\rightarrow 0,$ which is very restrictive. On the contrary, other methodologies typically achieve average Frobenius norm consistency provided $s\log n/T\rightarrow 0,$ where $s$ is some sparsity index (e.g., see Bickel and Levina (2008) Theorem 2 with $q=0$).333Average Frobenius norm means dividing a Frobenius norm by $\sqrt{n}$, while relative Frobenius norm means dividing a Frobenius norm by the Frobenius norm of a target matrix, say, the unknown covariance matrix. These two concepts are similar, but not exactly the same.

In this article, we relax the rate restriction on $n$ imposed by Hafner et al. (2020) and allow $n$ to be possibly larger than $T$. We propose a new covariance matrix estimator called the quadratic form estimator based on the Kronecker product model. Our estimator averages elements of the sample covariance matrix, so we obtain a rate improvement by averaging. In particular, under a cross-sectional weak dependence condition, the quadratic form estimator achieves relative Frobenius norm consistency provided $\log^{3}n/T\to 0$. Moreover, this method automatically produces a symmetric and positive definite covariance matrix estimator, unlike some of the sparsifying methods considered by Fan et al. (2016).

We apply our methodology to a concrete testing problem; we consider the null hypothesis $H_{0}\mathrel{\mathop{\mathchar 58\relax}}\mu=\mu_{0}$, where $\mu$ is the mean of the large dimensional data $y_{t}$ and $\mu_{0}$ is some known vector. One practical example would be that $y_{t}$ corresponds to differences between treated and controlled groups and we want to test whether the mean cross-sectional differences are different from zero. We define the Lagrange multiplier (LM) and Wald test statistics based on our estimated precision matrix and establish their asymptotic distributions under both null and local alternatives of the form $H_{1}\mathrel{\mathop{\mathchar 58\relax}}\mu=\mu_{0}+\theta/\sqrt{T}$ for some vector $\theta$. We also provide two results regarding testing linear restrictions of $\mu$.

We compare our estimation and testing methods with Ledoit and Wolf (2004)’s linear shrinkage estimator and Ledoit and Wolf (2017)’s direct nonlinear shrinkage estimator in Monte Carlo simulations. Our methods perform very well in moderate-sized samples. In fact, they work well even in situations where a Kronecker product model is misspecified for a covariance matrix.

The rest of the article is structured as follows. In Section 2 we discuss the model and identification while in Section 3 we propose the quadratic form estimator. Section 4 gives the rate of convergence for the quadratic form estimator. In Section 5 we define the LM and Wald test statistics and establish their asymptotic distributions under both null and local alternatives. We also consider testing linear restrictions of $\mu$. Section 6 conducts Monte Carlo simulations comparing our approach with Ledoit and Wolf estimators. Section 7 concludes. All the major proofs are put in the Appendix while auxiliary lemmas and theorems are in Section B.

1.1 Notation

Let $A$ be an $m\times n$ matrix. Let $\operatorname*{vec}A$ denote the vector obtained by stacking the columns of $A$ one underneath the other. The commutation matrix $K_{m,n}$ is an $mn\times mn$ orthogonal matrix which translates $\operatorname*{vec}A$ to $\operatorname*{vec}(A^{\intercal})$, i.e., $\operatorname*{vec}(A^{\intercal})=K_{m,n}\operatorname*{vec}(A)$. If $A$ is a symmetric $n\times n$ matrix, its $n(n-1)/2$ supradiagonal elements are redundant in the sense that they can be deduced from symmetry. If we eliminate these redundant elements from $\operatorname*{vec}A$, we obtain a new $n(n+1)/2\times 1$ vector, denoted $\operatorname*{vech}A$. They are related by the full-column-rank, $n^{2}\times n(n+1)/2$ duplication matrix $D_{n}$: $\operatorname*{vec}A=D_{n}\operatorname*{vech}A$. Conversely, $\operatorname*{vech}A=D_{n}^{+}\operatorname*{vec}A$, where $D_{n}^{+}$ is $n(n+1)/2\times n^{2}$ and the Moore-Penrose generalized inverse of $D_{n}$. In particular, $D_{n}^{+}=(D_{n}^{\intercal}D_{n})^{-1}D_{n}^{\intercal}$ because $D_{n}$ is full-column-rank.

For $x\in\mathbb{R}^{n}$, let $\|x\|_{2}\mathrel{\mathop{\mathchar 58\relax}}=\sqrt{\sum_{i=1}^{n}x_{i}^{2}}$ and $\|x\|_{\infty}\mathrel{\mathop{\mathchar 58\relax}}=\max_{1\leq i\leq n}|x_{i}|$ denote the Euclidean ($\ell_{2}$) norm and the element-wise maximum ($\ell_{\infty}$) norm, respectively. Let $\lambda_{\max}(\cdot)$ and $\lambda_{\min}(\cdot)$ denote the maximum and minimum eigenvalues of some real symmetric matrix, respectively. For any real $m\times n$ matrix $A=(a_{i,j})_{1\leq i\leq m,1\leq j\leq n}$, let $\|A\|_{F}\mathrel{\mathop{\mathchar 58\relax}}=[\text{tr}(A^{\intercal}A)]^{1/2}\equiv[\text{tr}(AA^{\intercal})]^{1/2}\equiv\|\operatorname*{vec}A\|_{2}$, $\|A\|_{1}\mathrel{\mathop{\mathchar 58\relax}}=\sum_{i=1}^{m}\sum_{j=1}^{n}|a_{i,j}|$, $\|A\|_{\ell_{2}}\mathrel{\mathop{\mathchar 58\relax}}=\max_{\|x\|_{2}=1}\|Ax\|_{2}\equiv\sqrt{\lambda_{\max}(A^{\intercal}A)}$, $\|A\|_{\ell_{1}}\mathrel{\mathop{\mathchar 58\relax}}=\max_{1\leq j\leq n}\sum_{i=1}^{m}|a_{i,j}|$, and $\|A\|_{\ell_{\infty}}\mathrel{\mathop{\mathchar 58\relax}}=\max_{1\leq i\leq m}\sum_{j=1}^{n}|a_{i,j}|$ denote the Frobenius ($\ell_{2}$) norm, $\ell_{1}$ norm, and spectral norm ($\ell_{2}$-operator norm), maximum column sum matrix norm ($\ell_{1}$-operator norm), and maximum row sum matrix norm ($\ell_{\infty}$-operator norm) of $A$, respectively. Note that $\|\cdot\|_{\infty}$ can also be applied to matrix $A$, i.e., $\|A\|_{\infty}=\max_{1\leq i\leq m,1\leq j\leq n}|a_{i,j}|$; however $\|\cdot\|_{\infty}$ is not a matrix norm so it does not have the submultiplicative property of a matrix norm.

Landau (order) notation in this article, unless otherwise stated, should be interpreted in the sense that $n,T\rightarrow\infty$ simultaneously. An absolute positive constant refers to a constant independent of anything which is a function of $n$ and/or $T$. We write $a\asymp b$ if there exist absolute constants $0<c_{1}\leq c_{2}$ such that $c_{1}b\leq a\leq c_{2}b$. For real numbers $a,b$ let $a\vee b$ denote $\max(a,b)$.

2 The Model and Identification

We now directly work with the high-level $n$-dimensional random vector $y_{t}$ with $\mu\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}y_{t}$ and $\Sigma\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}[(y_{t}-\mu)(y_{t}-\mu)^{\intercal}]$ for every $t$. In particular, $\Sigma$ takes the form of (1.1). For each $j$, $\Sigma_{j}$ contains $n_{j}(n_{j}+1)/2-1$ (unrestricted) parameters. In total, model (1.1) contains $\sum_{j=1}^{v}n_{j}(n_{j}+1)/2-(v-1)$ unknown parameters. This model is the same as considered in Hafner et al. (2020) except that we make a different identifying restriction. The implied form for $\Sigma^{-1}$ is also Kronecker, i.e., $\Sigma^{-1}=\sigma^{-2}\times\Sigma_{1}^{-1}\otimes\cdots\otimes\Sigma_{v}^{-1}.$

We show that model (1.1) is indeed identified. First, the parameter $\sigma$ is identified because

[TABLE]

whence we have $\sigma^{2}=\mathrm{tr}(\Sigma)/n$. We next consider identification of the remaining parameters based on the partial trace operator (Filipiak et al. (2018)). Suppose that an $n\times n$ matrix $A$ can be written in terms of $n_{1}\times n_{1}$ blocks of $n_{-1}\times n_{-1}$ dimensional matrices $A_{-1;i,j}$, where $n_{-1}\mathrel{\mathop{\mathchar 58\relax}}=n/n_{1}$; that is,

[TABLE]

Then the partial trace operator $\mathrm{PTR}_{n_{1}}\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n\times n}\mathbb{\rightarrow R}^{n_{1}\times n_{1}}$ is defined as follows:

[TABLE]

Consider model (1.1), and let $\Sigma_{-1}\mathrel{\mathop{\mathchar 58\relax}}=\Sigma_{2}\otimes\cdots\otimes\Sigma_{v}$. Define the $n_{1}\times n_{1}$ matrix $d^{(1)}\mathrel{\mathop{\mathchar 58\relax}}=\mathrm{PTR}_{n_{1}}(\Sigma)=\sigma^{2}\mathrm{tr}(\Sigma_{-1})\times\Sigma_{1}.$ Then $\Sigma_{1}=d^{(1)}/(\mathrm{tr}(d^{(1)})/n_{1}).$ According to Definition 1.1(ii) of Filipiak et al. (2018), $\mathrm{PTR}_{n_{1}}(\Sigma)=\sum_{\ell=1}^{n_{-1}}(I_{n_{1}}\otimes e_{\ell,n_{-1}}^{\intercal})\Sigma(I_{n_{1}}\otimes e_{\ell,n_{-1}}),$ where $e_{\ell,n_{-1}}$ is the $n_{-1}\times 1$ elementary vector with one in position $\ell$ and zero elsewhere. In this sense, $d^{(1)}$ is a quadratic form of $\Sigma$.

We next consider the remaining components $\Sigma_{h},$ $h=2,\ldots,v$. Write

[TABLE]

Note that $\Sigma_{-h}$ is $n_{-h}\times n_{-h}$ dimensional, where $n_{-h}\mathrel{\mathop{\mathchar 58\relax}}=n/n_{h}$. Recalling the identity $B\otimes A=K_{p,m}(A\otimes B)K_{m,p}$ for $A$ ($m\times m$) and $B$ ($p\times p$) (Magnus and Neudecker (1986) Lemma 4), we write

[TABLE]

Define the $n_{h}\times n_{h}$ matrix $d^{(h)}\mathrel{\mathop{\mathchar 58\relax}}=\mathrm{PTR}_{n_{h}}(\Sigma^{(h)})=\sigma^{2}\mathrm{tr}(\Sigma_{-h})\times\Sigma_{h}$. Then

[TABLE]

3 Estimation

We observe an $n$-dimensional weakly stationary time series vector $\{y_{t}\}_{t=1}^{T}$ with mean $\mu$ and covariance matrix $\Sigma$. Define the sample covariance matrix

[TABLE]

where $\bar{y}\mathrel{\mathop{\mathchar 58\relax}}=\frac{1}{T}\sum_{t=1}^{T}y_{t}$. Define $\hat{d}^{(1)}\mathrel{\mathop{\mathchar 58\relax}}=\mathrm{PTR}_{n_{1}}(M_{T})$. Then let $\tilde{\Sigma}_{1}\mathrel{\mathop{\mathchar 58\relax}}=\hat{d}^{(1)}/(\mathrm{tr}(\hat{d}^{(1)})/n_{1}).$ Likewise, define the ”permuted” sample covariance matrix

[TABLE]

for $h=2,\ldots,v$. Define $\hat{d}^{(h)}\mathrel{\mathop{\mathchar 58\relax}}=\mathrm{PTR}_{n_{h}}(M_{T}^{(h)})$ for $h=2,\ldots,v$. Then

[TABLE]

for $h=1,\ldots,v$.

The quadratic form estimator $\tilde{\Sigma}$ for $\Sigma$ is

[TABLE]

[TABLE]

By Lemma 2.4 of Filipiak et al. (2018), if $M_{T}$ is symmetric and positive semidefinite, then so are $\{\tilde{\Sigma}_{j}\}_{j=1}^{v}$ and hence $\tilde{\Sigma}.$ Moreover, simulations show that even for positive semidefinite $M_{T}$, $\{\tilde{\Sigma}_{j}\}_{j=1}^{v},$ and hence $\tilde{\Sigma},$ are positive definite. As a result, the quadratic form estimator $\tilde{\Sigma}^{-1}$ for $\Sigma^{-1}$ is $\tilde{\Sigma}^{-1}=\hat{\sigma}^{-2}\times\tilde{\Sigma}_{1}^{-1}\otimes\cdots\otimes\tilde{\Sigma}_{v}^{-1}.$ We stress that $\tilde{\Sigma}^{-1}$ exists even if $n>T$. The quadratic form estimator is closely related to the quasi-maximum likelihood estimation (QMLE), but has the particular advantage in large dimensions in the sense that it is in closed form.444In the previous version of this article, we introduced a variant of the quadratic form estimator, which was derived by replacing the partial trace operator with a partial sum operator. Because of inferiority of that variant, we no longer include it in the current version.

In general we expect each element of $M_{T}$ to be $\sqrt{T}$-consistent, but here we are averaging over a large number of such elements. Under a cross-sectional weak dependence condition, like Assumption 4.3, we should have a rate improvement for the quadratic form estimator. We formally establish this in Section 4.

4 The Rate of Convergence

In this section, we shall derive the rate of convergence for the quadratic form estimator. We make the following assumptions:

Assumption 4.1**.**

1. (i)

The sample $\{y_{t}\}_{t=1}^{T}$ are independent over $t$. 2. (ii)

[TABLE]

for some absolute positive constant $A$. 3. (iii)

Consider a normal random vector $z_{t}$ which has the same mean vector and covariance matrix as those of $y_{t}$. The $n^{2}\times n^{2}$ kurtosis matrix of $y_{t}$ satisfies

[TABLE]

for some absolute positive constant $C$ for every $t$, where $\leq$ is to be interpreted componentwise.

Assumption 4.1(i) facilitates our technical analysis, but is perhaps not necessary. Assumption 4.1(ii) assumes the existence of an infinite number of moments of $y_{t}$, which allows one to invoke a concentration inequality such as the Bernstein’s inequality. Normal random vectors or random vectors that exhibit some exponential-type tail probability (e.g., subgaussianity, subexponentiality, semiexponentiality etc) satisfy this condition. Assumption 4.1(iii) supposes that the kurtosis matrix of $y_{t}$ is of the same order of magnitude as if it were a normal random vector. We impose this restriction on the kurtosis matrix of $y_{t}$ because not much research has touched on unrestricted kurtosis matrices in the large dimensional case.

Assumption 4.2**.**

1. (i)

$\{n_{j}\}_{j=1}^{v}$* are all fixed for any fixed $v$ and $v\rightarrow\infty$ (i.e., $n,T\to\infty$).* 2. (ii)

$\min_{1\leq j\leq v}\lambda_{\min}(\Sigma_{j})$* is bounded away from zero by an absolute positive constant as $n,T\to\infty$.*

Assumption 4.2(i) says that the dimensions of the sub-matrices are fixed when the number of sub-matrices tends to infinity. Note that Assumption 4.2(ii) does not necessarily imply that $\lambda_{\min}(\Sigma)$ is bounded away from zero by an absolute positive constant. This is because $\lambda_{\min}(\Sigma)=\sigma^{2}\times\prod_{j=1}^{v}\lambda_{\min}(\Sigma_{j})$ and $v\rightarrow\infty$.

Lemma 4.1**.**

Suppose Assumption 4.2(i) hold. We have

1. (i)

$v=O(\log n)$. 2. (ii)

$\max_{1\leq j\leq v}\lambda_{\max}(\Sigma_{j})$* is bounded from above by an absolute positive constant as $n,T\to\infty$.*

Note that Lemma 4.1(ii) does not necessarily imply that $\lambda_{\max}(\Sigma)$ is bounded from above by an absolute positive constant. This is because $\lambda_{\max}(\Sigma)=\sigma^{2}\times\prod_{j=1}^{v}\lambda_{\max}(\Sigma_{j})$ and $v\to\infty$.

Assumption 4.3**.**

Let $0\leq\beta_{1}\leq 2$.

[TABLE]

Assumption 4.3 characterises the cross-sectional dependence of $\{y_{t}\}_{t=1}^{T}$. According to Proposition 1 of Chudik and Pesaran (2013), $\{y_{t}\}_{t=1}^{T}$ is said to be cross-sectionally weakly dependent. The smaller $\beta_{1}$ is, the less cross-sectional dependence of $\{y_{t}\}_{t=1}^{T}$ is allowed and the stronger Assumption 4.3 is. When $\beta_{1}=2$, Assumption 4.3 is slack as we are not restricting cross-sectional dependence of $\{y_{t}\}_{t=1}^{T}$ at all ($\|\Sigma\|_{F}^{2}=O(n^{2})$ in general). On the one hand, we would like to assume $\beta_{1}$ as close to 2 as possible to make Assumption 4.3 as weak as possible. On the other hand, the smaller $\beta_{1}$ is, the weaker cross-sectional dependence $\{y_{t}\}_{t=1}^{T}$ exhibits, and the faster rate of convergence the quadratic form estimator will be able to achieve. There is a trade off.

One important case is $\beta_{1}=1$. In this case one sufficient condition for Assumption 4.3 is that $\Sigma$ has bounded maximum column sum matrix norm (i.e., $\|\Sigma\|_{\ell_{1}}=O(1)$) or bounded maximum row sum matrix norm (i.e., $\|\Sigma\|_{\ell_{\infty}}=O(1)$). To see this

[TABLE]

Note that for symmetric $\Sigma$, bounded maximum column sum matrix norm or bounded maximum row sum matrix norm implies that the maximum eigenvalue of $\Sigma$ is bounded from above by an absolute positive constant and the minimum eigenvalue of $\Sigma^{-1}$ is bounded away from zero by an absolute positive constant: $1/(\lambda_{\min}(\Sigma^{-1}))=\lambda_{\max}(\Sigma)=\|\Sigma\|_{\ell_{2}}\leq\|\Sigma\|_{\ell_{1}}=\|\Sigma\|_{\ell_{\infty}}=O(1)$. The assumption of bounded maximum column/row sum matrix norm has been used by Fan et al. (2015) (their Assumption 4.1(i)) and Pesaran and Yamagata (2012) (their Assumption 3).

Theorem 4.1**.**

Suppose Assumptions 4.1, 4.2 and 4.3 hold. If $\log^{3}n/T\to 0$ as $n,T\to\infty$, then we have

1. (i)

[TABLE] 2. (ii)

[TABLE] 3. (iii)

[TABLE] 4. (iv)

[TABLE] 5. (v)

[TABLE] 6. (vi)

[TABLE]

The reason that we divide the Frobenius norm of the estimation error, say, $\|\tilde{\Sigma}-\Sigma\|_{F}$, by the Frobenius norm of the target, i.e., $\|\Sigma\|_{F}$, is to define a proper notion of ”consistency”. This is necessary because the cross-sectional dimension $n$ is growing to infinity. In this case, even if every element of a matrix-valued estimator is converging in probability to the corresponding element of its target matrix, there is no guarantee that its overall estimation error will converge to zero in probability when $n,T\rightarrow\infty$. The rescaling of the Frobenius norm of the estimation error is standard in the large dimensional case, but in the literature scholars tend to divide the Frobenius norm of the estimation error by $\sqrt{n}$ (e.g., see Bickel and Levina (2008) Theorem 2, Fan et al. (2011) p3330, Ledoit and Wolf (2004) Definition 1 etc). The same reasoning applies to the $\ell_{1}$ and the spectral norm of the estimation error.

Note that there are two terms on the right side. The term O_{p}\mathinner{\bigl{(}\log^{2}n/T\bigr{)}} exists because we need to estimate the unknown $\mu$. If we knew $\mu$, this term would not be present.555If we knew $\mu$, the estimation procedure in Section 3 applies to $M_{T}^{0}\mathrel{\mathop{\mathchar 58\relax}}=T^{-1}\sum_{t=1}^{T}(y_{t}-\mu)(y_{t}-\mu)^{\intercal}$ instead of $M_{T}$. The rate of convergence, \mathinner{\bigl{(}\log^{3}n/(n^{2-\beta_{1}}T)\bigr{)}}^{1/2}, contains an additional, non-standard item $\sqrt{n^{2-\beta_{1}}}$ in the denominator. This non-standard item exists because of the cross-sectional weak dependence condition (Assumption 4.3). If $\beta_{1}=2$ (i.e., we are not restricting cross-sectional dependence of $\{y_{t}\}_{t=1}^{T}$ at all), this term vanishes. The rate of convergence of the quadratic form estimator then becomes $(\log^{3}n/T)^{1/2}$, which is comparable to the convergence rates of other existent estimators in the large dimensional case.

Take part (i) of the theorem as an illustration. If $\beta_{1}=2$ and we knew $\mu$, we have \mathinner{\!\bigl{\lVert}\tilde{\Sigma}-\Sigma\bigr{\rVert}}_{F}=O_{p}\mathinner{\bigl{(}\|\Sigma\|_{F}(\log^{3}n/T)^{1/2}\bigr{)}}. A typical threshold estimator $\hat{\Sigma}_{\text{thres}}$ has \|\hat{\Sigma}_{\text{thres}}-\Sigma^{\dagger}\|_{F}=O_{p}\mathinner{\bigl{(}(sn\log n/T)^{1/2}\bigr{)}}, where $\Sigma^{\dagger}$ is some sparse truth and $s$ is its sparsity index (see Bickel and Levina (2008) Theorem 2 with $q=0$). According to Bickel and Levina (2008), $s$ is the upper bound of non-zero elements for every row, so $\|\Sigma^{\dagger}\|_{F}=O(\sqrt{sn})$ under the sparsity model. If one assumes $\|\Sigma^{\dagger}\|_{F}\asymp\sqrt{sn}$, one can write \|\hat{\Sigma}_{\text{thres}}-\Sigma^{\dagger}\|_{F}=O_{p}\mathinner{\bigl{(}\|\Sigma^{\dagger}\|_{F}(\log n/T)^{1/2}\bigr{)}}. Then the two rates of convergence only differ by a logarithmic factor.

Because of the cross-sectional weak dependence condition (Assumption 4.3), the quadratic form estimator is able to achieve a faster rate of convergence than a typical estimator does.

5 Test Statistics

We apply our methodology to the testing issue. We consider the problem of testing the null hypothesis $H_{0}\mathrel{\mathop{\mathchar 58\relax}}\mu=\mu_{0}$ against the alternative $H_{1}\mathrel{\mathop{\mathchar 58\relax}}\mu\neq\mu_{0}$.

The classical Wald test statistic (based on the sample covariance matrix $M_{T}$) is not defined when $n\geq T$; there is a large literature that proposes alternative test statistics. Bai and Saranadasa (1996) proposed a statistic based on $\|\bar{y}\|_{2}^{2}$, thereby avoiding the inversion of the large sample covariance matrix, and established its asymptotic normality. Pesaran and Yamagata (2012) extended this approach to the Capital Asset Pricing Model (CAPM) regression setting and proposed several test statistics. One of the test statistics is based on $\|t\|_{2}^{2},$ where $t$ is a vector of individual $t$-statistics; Pesaran and Yamagata (2012) derived the limiting normal distribution of the centred and scaled version of this under cross-sectional weak dependence conditions. Fan et al. (2015) considered a Wald test statistic for testing the CAPM restrictions inside a linear regression in the large dimensional case. They regularized the estimated error covariance matrix by imposing a sparsity assumption, and used that to form a quadratic form. They established the null limiting distribution of their test statistic (they also proposed a novel power enhancement procedure, which we do not study here).

We now define the Lagrange multiplier (LM) test statistic

[TABLE]

where $\tilde{\Sigma}_{\mu_{0}}$ is the quadratic form estimator assuming that we know $\mu=\mu_{0}$. The Wald test statistic is

[TABLE]

which is the Hotelling $T^{2}$-statistic based on the quadratic form estimator. We next present the large sample properties of the binity $LM_{n,T}$ and $W_{n,T}$. We make one more cross-sectional dependence assumption.

Assumption 5.1**.**

Let $0\leq\beta_{2}\leq 2$.

[TABLE]

The bigger $\beta_{2}$ is, the weaker Assumption 5.1 is. This is because it is putting less restriction on the cross-sectional dependence of $\Sigma^{-1}$. When $\beta_{2}=2$, Assumption 5.1 is slack, as in essence we are not restricting anything. On the one hand, we wish to assume $\beta_{2}$ as close to 2 as possible to make Assumption 5.1 as weak as possible. On the other hand, we wish to assume that $\beta_{2}$ is as small as possible so that our methodology could accommodate an $n$ as large as possible.

One important case is $\beta_{2}=1$. In this case, a sufficient condition for Assumption 5.1 is that $\Sigma^{-1}$ has bounded maximum column sum matrix norm (i.e., $\|\Sigma^{-1}\|_{\ell_{1}}=O(1)$) or bounded maximum row sum matrix norm (i.e., $\|\Sigma^{-1}\|_{\ell_{\infty}}=O(1)$). To see this

[TABLE]

Note that for symmetric $\Sigma^{-1}$, bounded maximum column sum matrix norm or bounded maximum row sum matrix norm implies that the maximum eigenvalue of $\Sigma^{-1}$ is bounded from above by an absolute positive constant and the minimum eigenvalue of $\Sigma$ is bounded away from zero by an absolute positive constant: $1/(\lambda_{\min}(\Sigma))=\lambda_{\max}(\Sigma^{-1})=\|\Sigma^{-1}\|_{\ell_{2}}\leq\|\Sigma^{-1}\|_{\ell_{1}}=\|\Sigma^{-1}\|_{\ell_{\infty}}=O(1)$. The assumption of bounded maximum column/row sum matrix norm has been used by Fan et al. (2015) (their Assumption 4.1(i)) and Pesaran and Yamagata (2012) (their Assumption 3).

Theorem 5.1**.**

Suppose Assumptions 4.1, 4.2, 4.3, and 5.1 hold. We make the following assumptions:

1. (a)

[TABLE] 2. (b)

Consider the Cholesky decomposition of $\Sigma$, i.e., $\Sigma=LL^{\intercal}$, where $L$ is a nonsingular lower triangular matrix L with positive diagonal elements. Assume that $x_{t}\mathrel{\mathop{\mathchar 58\relax}}=L^{-1}(y_{t}-\mu)$ is cross-sectionally independent for any $t$, and for some $\delta>0$

[TABLE]

Then under $H_{0}\mathrel{\mathop{\mathchar 58\relax}}\mu=\mu_{0}$, as $n,T\rightarrow\infty$,

[TABLE]

If one additionally assumes

[TABLE]

then under $H_{0}\mathrel{\mathop{\mathchar 58\relax}}\mu=\mu_{0}$, as $n,T\rightarrow\infty$,

[TABLE]

For the LM test, if we want to allow the interesting case of $n/T\to\infty$, then assumption (a) necessarily implies that $2\beta_{2}+\beta_{1}<4$, which restricts both $\beta_{2}$ and $\beta_{1}$. In the special case of $\beta_{1}=\beta_{2}=1$, assumption (a) is reduced to $\log^{5}n/T=o(1)$, which is a weak condition.

Assumption (b) is standard in the literature. Fan et al. (2015) maintained normality (their Assumption 4.1(i)), which is a special case of assumption (b). Pesaran and Yamagata (2012) also maintained assumption (b) (their Assumption 2a). Assumption (b) implicitly assumes that $\lambda_{\min}(\Sigma)$ is bounded away from zero by an absolute positive constant, which strengthens Assumption 4.2(ii). Also note that $\operatorname*{var}(x_{t})=I_{n}$, so strengthening from cross-sectional uncorrelatedness to cross-sectional independence in assumption (b) is rather innocuous. In addition, we assume that the $(4+2\delta)$th moment of $x_{t,i}$ is (uniformly in $i$ and $t$) finite for $n,T$ sufficiently large, which is also a weak assumption. Under the more restricted sequential limit ($T\rightarrow\infty$ and then $n\rightarrow\infty),$ $\sqrt{T}(\bar{y}-\mu_{0})$ is approximately normal so the limiting properties could be calculated for the non-normal case as if normality held. However, in our framework of joint limits, such procedures break down, so we make assumption (b).

In the low-dimensional case ($n$ fixed, $T\rightarrow\infty$), the LM test statistic $LM_{n,T}$ and the Wald test statistic $W_{n,T}$ are asymptotically equivalent in the sense that they all converge in distribution to $\chi_{n}^{2}$.666The finite sample performance of these statistics is known to vary. Park and Phillips (1988) established higher order approximations for a Wald test of nonlinear restrictions in the finite dimensional case, and showed how to improve performance of the test statistic. It may be possible to apply their methodology to the large dimensional case. In the large dimensional case ($n,T\rightarrow\infty$), Theorem 5.1 shows that $LM_{n,T}$ and $W_{n,T}$ are, again, asymptotically equivalent. The Wald test requires an additional rate restriction (5.3), which is the price we pay for estimating $\Sigma^{-1}$ under the alternative $H_{1}\mathrel{\mathop{\mathchar 58\relax}}\mu\neq\mu_{0}$.

Recall that a typical threshold estimator $\hat{\Sigma}_{\text{thres}}$ has \|\hat{\Sigma}_{\text{thres}}^{-1}-(\Sigma^{\dagger})^{-1}\|_{\ell_{2}}=O_{p}\mathinner{\bigl{(}s(\log n/T)^{1/2}\bigr{)}}, where $\Sigma^{\dagger}$ is some sparse truth and $s$ is its sparsity index (see Bickel and Levina (2008) Theorem 1 with $q=0$). For this rate of convergence, a result like (5.4) requires, as both Pesaran and Yamagata (2012) and Fan et al. (2015) have pointed out, $n\log n/T=o(1)$, which is essentially a low-dimensional scenario. Pesaran and Yamagata (2012) and Fan et al. (2015) have hence come up with their own ingenious ways to relax the condition $n\log n/T=o(1)$ and established results similar to (5.4) for their Wald test statistics in the CAPM context.

In the case of our Wald test, if we also want to allow the interesting large dimension case of $n/T\rightarrow\infty$, then assumption (a) and (5.3) necessarily imply that $2\beta_{2}+\beta_{1}<4$ and $\beta_{2}<3/2$, respectively. For example, we can choose the special case $\beta_{1}=\beta_{2}=1$, so assumption (a) and (5.3) reduce to

[TABLE]

the latter of which is the binding rate condition and the same as the rate condition in Assumption 4.2 of Fan et al. (2015).

In the simulation study below we compare our tests with test statistics that use Ledoit and Wolf procedures to regularize the sample covariance matrix estimator.

5.1 Power Investigation

In this section, we analyse the asymptotic distributions of the proposed test statistics under the alternative hypothesis $H_{1}\mathrel{\mathop{\mathchar 58\relax}}\mu\neq\mu_{0}$. In particular, we shall focus on a sequence of local alternatives $H_{1}\mathrel{\mathop{\mathchar 58\relax}}\mu=\mu_{T}\mathrel{\mathop{\mathchar 58\relax}}=\mu_{0}+\theta/\sqrt{T}$, where $\max_{1\leq i\leq n}|\theta_{i}|=O(\sqrt{\log n})$. We focus on the Wald test without loss of generality.

Theorem 5.2**.**

Suppose Assumptions 4.1, 4.2, 4.3, and 5.1 hold. We make the following additional assumptions:

1. (a)

2. (i)

[TABLE] 2. (ii)

[TABLE] 2. (b)

Consider the Cholesky decomposition of $\Sigma$, i.e., $\Sigma=LL^{\intercal}$, where $L$ is an $n\times n$ nonsingular lower triangular matrix with positive diagonal elements. Assume that $x_{t}\mathrel{\mathop{\mathchar 58\relax}}=L^{-1}(y_{t}-\mu)$ is cross-sectionally independent for any $t$, and for some $\delta>0$

[TABLE]

Then under $H_{1}\mathrel{\mathop{\mathchar 58\relax}}\mu=\mu_{0}+\theta/\sqrt{T}$,

[TABLE]

The preceding theorem shows that the asymptotic distribution of $(W_{n,T}-n)/\sqrt{2n+4\theta^{\intercal}\Sigma^{-1}\theta}$ under $H_{1}$ has a center $\theta^{\intercal}\Sigma^{-1}\theta/\sqrt{2n+4\theta^{\intercal}\Sigma^{-1}\theta}$. Note that

[TABLE]

In the special case of $0<\lambda_{\min}(\Sigma)<\lambda_{\max}(\Sigma)<\infty$, we see that the test has power against local alternatives that satisfy $\max_{1\leq i\leq n}|\theta_{i}|=O(\sqrt{\log n})$ and $\theta^{\intercal}\theta=O(n^{\delta_{a}}),$ where $\delta_{a}\geq 1/2,$ and power tending to one in the case where $\delta_{a}>1/2.$ This specification requires that $\theta$ has a sufficiently large number of non-zero elements. It does not require that all the elements of $\theta$ are non-zero.

5.2 Testing Linear Restrictions of $\mu$

In this section, we consider testing linear restrictions of $\mu$ using two approaches. We first consider $H_{0}\mathrel{\mathop{\mathchar 58\relax}}R\mu=r$, where $R$ is a $q\times n$ matrix of rank $q$. We assume that $q$ is a fixed number; this case covers applications where a finite number of linear restrictions are coming from economic theory.

Theorem 5.3**.**

Suppose Assumptions 4.1, 4.2 and 4.3 hold. We also make the following assumptions:

1. (a)

$\lambda_{\min}(\Sigma)$* is bounded away from zero by an absolute positive constant.* 2. (b)

Consider $H_{0}\mathrel{\mathop{\mathchar 58\relax}}R\mu=r$, where $R$ is a $q\times n$ matrix of rank $q$ for any fixed $n$ and $n\rightarrow\infty$ ($q$* is a fixed number). Moreover, $R$ and $r$ are rescaled in such a way that $\lambda_{\min}(RR^{\intercal})$ is bounded away from zero by an absolute constant, and*

[TABLE]

Then under $H_{0}\mathrel{\mathop{\mathchar 58\relax}}R\mu=r$, if $\log^{3}n/T\to 0$ as $n,T\to\infty$,

[TABLE]

Assumption (a) strengthens Assumption 4.2(ii) slightly, which is a mild condition. A sufficient condition for (5.5) in assumption (b) is $\lambda_{\max}(RR^{\intercal})$ is bounded from above by an absolute positive constant and $\|\Sigma\|_{\ell_{2}}<\infty$. The requirement of $\lambda_{\min}(RR^{\intercal})$ and $\lambda_{\max}(RR^{\intercal})$ being bounded away from zero and from above by absolute positive constants, respectively, could be achieved by normalising each row of $R$ to have $\ell_{2}$ norm of 1.

We next take another approach to derive simultaneous confidence intervals for all linear combinations of $\mu$.

Lemma 5.1**.**

Suppose Assumptions 4.1, 4.2, 4.3, and 5.1 hold. Simultaneously for all $\phi\in\mathbb{R}^{n}$, the unknown $\mu$ satisfies the following inequalities with confidence $1-\alpha$:

[TABLE]

as $n,T\rightarrow\infty$, where $z_{\alpha}$ is the upper $\alpha$ percentile of a standard normal.

One disadvantage of this approach is that the confidence region for $\mu$ could be conservative.

6 Simulation Study

In this section, we provide some Monte Carlo simulations that evaluate performance of our procedures.

6.1 The Correctly Specified Case

We suppose that $y_{t}\sim N(\mu,\Sigma)$ with $\Sigma=\Sigma_{1}\otimes\cdots\otimes\Sigma_{v}$, where

[TABLE]

so in this case $\Sigma$ is also the correlation matrix. Given $\|\Sigma_{j}\|_{F}^{2}=2(1+\rho_{j}^{2})$, we have

[TABLE]

Since $\prod_{j=1}^{v}(1+\rho_{j}^{2})\geq 1$, Assumption 4.3 necessarily implies $\beta_{1}\geq 1$. When $\beta_{1}=1$, $\frac{1}{n^{\beta_{1}}}\|\Sigma\|_{F}^{2}=\prod_{j=1}^{v}(1+\rho_{j}^{2})$, which converges to a finite, non-zero limit as $v\rightarrow\infty$ if and only if $\sum_{j=1}^{v}\rho_{j}^{2}$ converges (Knopp (1947) Theorem 28.3). When $\beta_{1}>1$, Assumption 4.3 is satisfied if $\prod_{j=1}^{v}(1+\rho_{j}^{2})=O(n^{\beta_{1}-1})$.

Likewise

[TABLE]

so that via Lemma B.3 in Section B

[TABLE]

Since \prod_{j=1}^{v}\mathinner{\bigl{(}1-|\rho_{j}|\bigr{)}}\leq 1, Assumption 5.1 necessarily implies $\beta_{2}\geq 1$. When $\beta_{2}=1$, \frac{1}{n^{\beta_{2}}}\|\Sigma^{-1}\|_{1}=1/\prod_{j=1}^{v}\mathinner{\bigl{(}1-|\rho_{j}|\bigr{)}}, the denominator of which converges to a finite, non-zero limit as $v\rightarrow\infty$ if and only if $\sum_{j=1}^{v}|\rho_{j}|$ converges (Knopp (1947) Theorem 28.4). When $\beta_{2}>1$, Assumption 5.1 is satisfied if \mathinner{\bigl{[}\prod_{j=1}^{v}\mathinner{\bigl{(}1-|\rho_{j}|\bigr{)}}\bigr{]}}^{-1}=O(n^{\beta_{2}-1}).777Furthermore, the largest eigenvalue of $\Sigma$ is $\prod_{j=1}^{v}(1+|\rho_{j}|)$, which converges as $v\rightarrow\infty$ if and only if $\sum_{j=1}^{v}|\rho_{j}|$ converges (Knopp (1947) Theorem 28.3).

We consider $\mu=0$, $n=2^{v}$, and $\rho_{j}=\rho^{j}$ for $j=1,\ldots,v$. The number of Monte Carlo simulations is 1000. We compare the quadratic form estimator with Ledoit and Wolf (2004)’s linear shrinkage estimator (the LW04 estimator hereafter) and Ledoit and Wolf (2017)’s direct nonlinear shrinkage estimator (the LW17 estimator hereafter).888The Matlab code for the LW04 and LW17 estimators is downloaded from the website of Professor Michael Wolf from the Department of Economics at the University of Zurich. We are grateful for this.

The first evaluation criterion is the relative mean square error (MSE) in terms of $\Sigma$. Given a generic estimator $\hat{\Sigma}_{G}$ of the covariance matrix $\Sigma$, we compute

[TABLE]

where the expectation operator is taken with respect to all the simulations. Often the precision matrix $\Sigma^{-1}$ is of more interest than $\Sigma$, so we also compute the MSE of the estimator of $\Sigma^{-1}$:

[TABLE]

where the expectation operator is taken with respect to all the simulations. Note that this requires invertibility of the generic estimator $\hat{\Sigma}_{G}$ and therefore cannot be calculated for the sample covariance matrix $M_{T}$ when $n>T$.

We next calculate

[TABLE]

where the expectation operator is taken with respect to all the simulations. The preceding display is called the simulated percentage relative improvement in average loss (PRIAL) criterion in terms of $\Sigma$ by Ledoit and Wolf (2004). The PRIAL measures the performance of the generic estimator $\hat{\Sigma}_{G}$ with respect to the sample covariance estimator $M_{T}$. Note that PRIAL$\in(-\infty,1]$: A negative value means $\hat{\Sigma}_{G}$ performs worse than $M_{T}$ while a positive value means otherwise. Likewise we also compute

[TABLE]

Note that this requires invertibility of the sample covariance matrix $M_{T}$ and therefore can only be calculated for $n<T$.

Finally, we consider testing $H_{0}\mathrel{\mathop{\mathchar 58\relax}}\mu=0$ against $H_{1}\mathrel{\mathop{\mathchar 58\relax}}\mu\neq 0$. We compute sizes of the LM and Wald tests (Theorem 5.1). The significance level is 5%. To investigate power, we generate $\mu_{i}$ as $\mu_{i}\overset{i.i.d.}{\sim}N(0,1)/\sqrt{T}$ for $i=1,2,\ldots,\lfloor n^{0.7}\rfloor$, where $\lfloor x\rfloor$ is the largest integer less than or equal to $x$; $\mu_{i}=0$ for $i=\lfloor n^{0.7}\rfloor+1,\lfloor n^{0.7}\rfloor+2,\ldots,n$. These also require invertibility of $\hat{\Sigma}_{G}$.

The results are reported in Tables 1-3. In Table 1, we set $T=252$ and $v=10$ so that $n=2^{v}=1024$; we set $\rho=0.5,0.7,0.85$. First, consider the top panel ($\rho=0.5$). For the MSE in terms of $\Sigma$ (i.e., MSE-1), all the estimators beat the sample covariance matrix $M_{T}$ by a large margin. The quadratic form estimator $\tilde{\Sigma}$ also outperformed the LW04 and LW17 estimators considerably. For the MSE in terms of $\Sigma^{-1}$ (i.e., MSE-2), a similar pattern exists. Note that the MSE-2 cannot be computed for $M_{T}$ because $M_{T}$ is not invertible when $n>T$. For the PRIAL in terms of $\Sigma$ (i.e., PRIAL-1), again $\tilde{\Sigma}$ is better than the LW04 and LW17 estimators. The sample covariance matrix $M_{T}$ has zero PRIAL-1 by definition. The superiority of $\tilde{\Sigma}$ in this experiment is expected because the true covariance matrix is indeed a Kronecker product.

Considering the size of the Wald test, we realize that the quadratic form estimator $\tilde{\Sigma}$ has the correct size while the LW04 and LW17 estimators are over-sized. Note that the Wald test is not defined for $M_{T}$ because $M_{T}$ is not invertible. Size of the LM test is similar to that of the Wald test for $\tilde{\Sigma}$, but the LM test seems to perform poorly for both the LW04 and LW17 estimators. Undoubtedly, the quadratic form estimator $\tilde{\Sigma}$ is the best performing estimator.

As we increase the ”mother” correlation parameter $\rho$ from 0.5 to 0.85, performance of $\tilde{\Sigma}$ remains unchanged across all five criteria. In terms of MSE-1, performance of $M_{T}$ improves while performances of LW04 and LW17 estimators initially worsen and then improve. In terms of MSE-2, PRIAL-1, the size of the LM test, and the size of the Wald test, the performances of both the LW04 and LW17 estimators worsen. Again the quadratic form estimator $\tilde{\Sigma}$ is the best performing estimator.

Next, we fix $\rho$ at 0.7 and examine effects of $n$ and $T$; the results are reported in Table 2. If we fix $T$ at 252 and increase $v$ (and hence $n$), in terms of MSE-1, all the estimators except the quadratic form estimator $\tilde{\Sigma}$ worsen. The same pattern is observed when we use the MSE-2 criterion instead (the sample covariance matrix $M_{T}$ dropped out in this case). In terms of PRIAL-1, we see that all the candidate estimators are becoming increasingly superior to $M_{T}$. As $n$ increases, size of the Wald test worsens for all the estimators except $\tilde{\Sigma}$; a similar pattern is observed for the LM test. If we increase $T$ from 252 to 504, all the estimators improve in terms of both the MSE-1 and MSE-2 criteria. Also sizes of the Wald and LM tests in general improve for all the estimators.

The results of power investigation are reported in Table 3. We see that power of the quadratic form estimator $\tilde{\Sigma}$ is very good for the specified local alternative. Powers of the LW04 and LW17 estimators in terms of the Wald test are less good. Powers of the LW04 and LW17 estimators in terms of the LM test come at a price of their sizes.

6.2 The Misspecified Case

To gauge how well the Kronecker product model performs when the true covariance matrix does not have a Kronecker product form, we consider the Monte Carlo setting used by Ledoit and Wolf (2004). We still assume that $y_{t}\sim N(\mu,\Sigma)$. The true covariance matrix $\Sigma$ is diagonal without loss of generality. The diagonal entries $\Sigma_{ii}$ (i.e., the eigenvalues of $\Sigma$) are log normally distributed: $\log\Sigma_{ii}\sim N(\mu_{\text{LW}},\sigma_{\text{LW}}^{2})$. Ledoit and Wolf (2004) defined the grand mean $\mu_{\text{g}}$ and cross-sectional dispersion $\alpha^{2}$ of the eigenvalues of $\Sigma$ as, respectively,

[TABLE]

In the Monte Carlo simulations, we re-define $\mu_{\text{g}}$ and $\alpha^{2}$ as the corresponding population counterparts:

[TABLE]

Ledoit and Wolf (2004) set $\mu_{\text{g}}=1$, so we can solve $\mu_{\text{LW}}=-\log(1+\alpha^{2})/2$ and $\sigma_{\text{LW}}^{2}=\log(1+\alpha^{2}),$ whence we have

[TABLE]

Note that in this data generating process, there are two sources of randomness: one from the normal distribution of $y_{t}$ and the other from the log normal distribution of $\Sigma_{ii}$. Also note that a diagonal covariance matrix need not have a Kronecker product structure unless, say, the diagonal elements are all equal. The number of Monte Carlo simulations is again set at 1000. In the baseline setting of Ledoit and Wolf (2004), $\mu=0,n=20$, $T=40$, and $\alpha^{2}=0.5$.

There are a number of different Kronecker products that we can consider to approximate $\Sigma$ (see Hafner et al. (2020) for more discussions of model selection). The possible Kronecker factorizations are $5\times 2\times 2$, $4\times 5$, $2\times 10$. Within each Kronecker factorization, we can further permute the Kronecker sub-matrices to obtain different Kronecker models. We study all the Kronecker products and compare with the LW04 and LW17 estimators. All estimators do not use knowledge of $\mu=0$ and have to estimate it, except in the case of the LM test.

The results are reported in Table 4. The first observation is that the performance of the quadratic form estimator $\tilde{\Sigma}$ is relatively robust to the Kronecker product factorization; the best performing one is $2\times 5\times 2$. All the candidate estimators beat the sample covariance matrix $M_{T}$. In terms of MSE-1 and MSE-2, the LW04 and LW17 estimators are only slightly better than $\tilde{\Sigma}$ $(2\times 5\times 2)$. In terms of PRIAL-1 and PRIAL-2, $\tilde{\Sigma}$ $(2\times 5\times 2)$ is almost as good as the LW04 and LW17 estimators. In terms of size of the LM test, $\tilde{\Sigma}$ $(2\times 5\times 2)$ has the correct size while the LW04 and LW17 estimators are under-sized. In terms of size of the Wald test, all candidate estimators are slightly over-sized.

We next vary $\alpha^{2}$. We base the comparisons on the $2\times 5\times 2$ Kronecker product factorization. The results are reported in Table 5. As $\alpha^{2}$ increases, performance of $M_{T}$ actually improves in terms of MSE-1 and MSE-2. On the other hand, performances of $\tilde{\Sigma}$, the LW04 and LW17 estimators worsen in terms of MSE-1, MSE-2, PRIAL-1 and PRIAL-2. The worsening performance of $\tilde{\Sigma}$ is not surprising because $\alpha^{2}$ can be interpreted as the distance of $\Sigma$ from a Kronecker product model. The worsening performance of the LW04 estimator has also been documented by Ledoit and Wolf (2004). As $\alpha^{2}$ increases, $\tilde{\Sigma}$ has roughly correct size for the LM test while both the LW04 and LW17 estimators are under-sized. In terms of the Wald test, all the candidate estimators are slightly over-sized.

Finally, we vary the ratio $n/T$. In the baseline setting we have $n/T=0.5$. Here we consider two variations. The first variation is $n=16,T=50$ with a ratio of $n/T=0.32$. The second variation is $n=40,T=20$ with a ratio of $n/T=2$. For the first variation, we identify the Kronecker product factorizations: $2\times 2\times 2\times 2$, $4\times 4$, $4\times 2\times 2$ and $2\times 8$. For the second variation, we use the Kronecker product factorizations: $5\times 2\times 2\times 2$, $5\times 2\times 4$, $5\times 8$ and $10\times 2\times 2$. We also considered permutations of sub-matrices for each factorization, but the performances remained relatively unchanged, so we do not report them in the interest of space. The results are reported in Table 6.

Consider the top panel of Table 6 first. All the candidate estimators beat the sample covariance matrix $M_{T}$. Performance of the quadratic form estimator $\tilde{\Sigma}$ is relatively robust to the Kronecker product factorizations ($2\times 2\times 2\times 2$, $4\times 4$ and $4\times 2\times 2$); the best performing one is $4\times 2\times 2$. In terms of MSE-1, MSE-2, PRIAL-1 and PRIAL-2, the quadratic form estimator $\tilde{\Sigma}$ $(4\times 2\times 2)$ is only slightly worse than the LW04 and LW17 estimators. In terms of size of the LM test, $\tilde{\Sigma}$ $(4\times 2\times 2)$ has the correct size while both the LW04 and LW17 estimators are under-sized. In terms of size of the Wald test, all the candidate estimators are slightly over-sized.

Next consider the bottom panel of Table 6. All the candidate estimators beat the sample covariance matrix $M_{T}$ again. The best performing quadratic form estimator has a factorization ($5\times 2\times 2\times 2$). In terms of MSE-1, MSE-2 and PRIAL-1, $\tilde{\Sigma}$ $(5\times 2\times 2\times 2)$ is comparable to the LW04 and LW17 estimators. In terms of size of the LM test, $\tilde{\Sigma}$ $(5\times 2\times 2\times 2)$ and the LW04 estimator have correct size while the LW17 estimator is slightly over-sized. In terms of size of the Wald test, all the candidate estimators are slightly over-sized.

By looking at Tables 4 and 6 together, we observe that as $n/T$ increases, PRIAL-1 increases monotonically for the best performing quadratic form estimator as well as the LW04 and LW17 estimators. Such a pattern is consistent with Ledoit and Wolf (2004). In terms of MSE-1 and MSE-2, performances of the best performing quadratic form estimator as well as the LW04 and LW17 estimators worsen as $n/T$ increases. In terms of size of the LM test, the best performing quadratic form estimator always has the correct size, while sizes of the Wald tests increase monotonically with $n/T$.

7 Concluding Remarks

We have proposed a new estimator of the Kronecker product model for covariance matrices - the quadratic form estimator. We establish the rate of convergence and use the estimated precision matrix to form the LM and Wald test statistics. The asymptotic distributions of these test statistics are established under both null and local alternative hypotheses. Testing linear restrictions of the unknown mean vector is also investigated. In Monte Carlo simulations, the quadratic form estimator performs well both when the Kronecker product model is correctly specified and when it is misspecified.

We remark on a number of possible extensions. One can generalize to allow weakly time series dependent data (see Hafner et al. (2020) for some work in this direction), and perhaps to where the spectral density matrix is Kronecker product factored. We may also consider the two-sample case where $\Sigma_{1}\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}[(y_{1,t}-\mu_{1})(y_{1,t}-\mu_{1})^{\intercal}]$ $(n\times n)$, $\Sigma_{2}\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}[(y_{2,t}-\mu_{2})(y_{2,t}-\mu_{2})^{\intercal}]$ $(n\times n)$, $\mu_{1}\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}(y_{1,t})$, and $\mu_{2}\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}(y_{2,t}).$ Cho and Phillips (2018) showed that the hypothesis of $\Sigma_{1}=\Sigma_{2}$ can be tested based on $\mathrm{tr}(\Sigma_{1}\Sigma_{2}^{-1})=n$; if both the covariance matrices have a conformable Kronecker product structure, this simplifies to $\mathrm{tr}(\Sigma_{1,1}\Sigma_{2,1}^{-1})\times\cdots\times\mathrm{tr}(\Sigma_{1,v}\Sigma_{2,v}^{-1})=n.$

Appendix A Appendix

A.1 Proof of Lemma 4.1

Proof.

For part (i), since $\prod_{j=1}^{v}n_{j}=n$, we have \mathinner{\bigl{(}\min_{1\leq j\leq v}n_{j}\bigr{)}}^{v}\leq n. Thus

[TABLE]

For part (ii):

[TABLE]

∎

A.2 Proof of Theorem 4.1

We first give an auxiliary lemma and an auxiliary theorem leading to the proof of Theorem 4.1.

A.2.1 Lemma A.1

Lemma A.1**.**

Suppose Assumptions 4.1 and 4.2 hold. Then we have

1. (i)

Both $\max_{1\leq j\leq v}\|\Sigma_{j}\|_{F}$ and $\max_{1\leq j\leq v}\|\Sigma_{j}^{-1}\|_{F}$ are bounded from above by absolute positive constants. Moreover both $\min_{1\leq j\leq v}\|\Sigma_{j}\|_{F}$ and $\min_{1\leq j\leq v}\|\Sigma_{j}^{-1}\|_{F}$ are bounded away from zero by absolute positive constants. 2. (ii)

Both $\max_{1\leq j\leq v}\|\Sigma_{j}\|_{1}$ and $\max_{1\leq j\leq v}\|\Sigma_{j}^{-1}\|_{1}$ are bounded from above by absolute positive constants. Moreover both $\min_{1\leq j\leq v}\|\Sigma_{j}\|_{1}$ and $\min_{1\leq j\leq v}\|\Sigma_{j}^{-1}\|_{1}$ are bounded away from zero by absolute positive constants. 3. (iii)

Both $\max_{1\leq j\leq v}\|\Sigma_{j}\|_{\ell_{2}}$ and $\max_{1\leq j\leq v}\|\Sigma_{j}^{-1}\|_{\ell_{2}}$ are bounded from above by absolute positive constants. Moreover both $\min_{1\leq j\leq v}\|\Sigma_{j}\|_{\ell_{2}}$ and $\min_{1\leq j\leq v}\|\Sigma_{j}^{-1}\|_{\ell_{2}}$ are bounded away from zero by absolute positive constants.

Proof of Lemma A.1.

For part (i), note that

[TABLE]

whence we deduce that $\max_{1\leq j\leq v}\|\Sigma_{j}\|_{F}$ is bounded from above by an absolute positive constant and $\min_{1\leq j\leq v}\|\Sigma_{j}\|_{F}$ is bounded away from zero by an absolute positive constant via Assumption 4.2 and Lemma 4.1. Similarly, we have

[TABLE]

whence we deduce that $\max_{1\leq j\leq v}\|\Sigma_{j}^{-1}\|_{F}$ is bounded from above by an absolute positive constant and $\min_{1\leq j\leq v}\|\Sigma_{j}^{-1}\|_{F}$ is bounded away from zero by an absolute positive constant via Assumption 4.2 and Lemma 4.1.

For part (ii), note that

[TABLE]

whence we deduce that part (ii) holds via part (i).

For part (iii), we have

[TABLE]

whence we could deduce the first half of the statement via part (i). Next,

[TABLE]

which is bounded away from zero by an absolute positive constant via Assumption 4.2(ii). Finally,

[TABLE]

which is bounded away from zero by an absolute positive constant via Lemma 4.1(ii). ∎

A.2.2 Theorem A.1

Theorem A.1**.**

Suppose Assumptions 4.1, 4.2 and 4.3 hold. Then

1. (i)

[TABLE] 2. (ii)

We have $\operatorname*{tr}(d^{(h)})/(n_{h}n_{-h})=\sigma^{2}>0$ for $h=1,\ldots,v$. Also,

[TABLE]

As a result, $\min_{1\leq h\leq v}\operatorname*{tr}(\hat{d}^{(h)})/(n_{h}n_{-h})$ is bounded away from zero by an absolute positive constant in probability. 3. (iii)

[TABLE]

*where $[\tilde{\Sigma}_{h}]_{i,j}$ and $[\Sigma_{h}]_{i,j}$ are the *$(i,j)$th entry of $\tilde{\Sigma}_{h}$ and $\Sigma_{h}$, respectively. 4. (iv)

[TABLE] 5. (v)

[TABLE] 6. (vi)

[TABLE] 7. (vii)

[TABLE] 8. (viii)

[TABLE] 9. (ix)

[TABLE] 10. (x)

[TABLE]

Proof.

For part (i), note that $d_{i,j}^{(h)}=\operatorname*{tr}\Sigma^{(h)}_{\{[i,j]\}}$, where $\Sigma^{(h)}_{\{[i,j]\}}$ is the $[i,j]$th block of $\Sigma^{(h)}$ (each block is $n_{-h}\times n_{-h}$ dimensional) for $i,j=1,\ldots,n_{h}$. Similarly, $\hat{d}_{i,j}^{(h)}=\operatorname*{tr}M^{(h)}_{T,\{[i,j]\}}$, where $M^{(h)}_{T,\{[i,j]\}}$ is the $[i,j]$th block of $M_{T}^{(h)}$ (each block is $n_{-h}\times n_{-h}$ dimensional). Write

[TABLE]

where

[TABLE]

$M^{0,(h)}_{T,\{[i,j]\}}$ is the $[i,j]$th block of $M_{T}^{0,(h)}$ (each block is $n_{-h}\times n_{-h}$ dimensional), and \mathinner{\bigl{(}\mathinner{\bigl{[}(\bar{y}-\mu)(\bar{y}-\mu)^{\intercal}\bigr{]}}^{(h)}\bigr{)}}_{\{[i,j]\}} is the $[i,j]$th block of \mathinner{\bigl{[}(\bar{y}-\mu)(\bar{y}-\mu)^{\intercal}\bigr{]}}^{(h)} (each block is $n_{-h}\times n_{-h}$ dimensional). Thus we have

[TABLE]

We consider the first term of (A.1) first. Note that $\mathbb{E}[\hat{d}_{i,j}^{0,(h)}]=d_{i,j}^{(h)}$. Write for some $M>0$

[TABLE]

where the second inequality is due to Chebyshev’s inequality. We now show that

[TABLE]

For arbitrary $i,j=1,\ldots,n_{h}$,

[TABLE]

where $\dot{y}_{t}^{(h)}\mathrel{\mathop{\mathchar 58\relax}}=K_{n_{h}\times\cdots\times n_{v},n_{1}\times\cdots\times n_{h-1}}(y_{t}-\mu)$ such that $\mathbb{E}[\dot{y}_{t}^{(h)}\dot{y}_{t}^{(h)\intercal}]=\Sigma^{(h)}$ and $\dot{z}_{t}^{(h)}$ is to be interpreted similarly, the third equality is due to independence over $t$ of $y_{t}$ in Assumption 4.1(i), and the first inequality is due to Assumption 4.1(iii). Using Lemma 9 of Magnus and Neudecker (1986), we have

[TABLE]

where the last equality is due to (33) of Magnus and Neudecker (1986). Thus we recognise that the summand on the right side of (A.2) is some element of \mathinner{\bigl{(}I_{n^{2}}+K_{n,n}\bigr{)}}(\Sigma^{(h)}\otimes\Sigma^{(h)}). We need to determine the exact position of the summand on the right side of (A.2) in \mathinner{\bigl{(}I_{n^{2}}+K_{n,n}\bigr{)}}(\Sigma^{(h)}\otimes\Sigma^{(h)}). We consider $\Sigma^{(h)}\otimes\Sigma^{(h)}$ and $K_{n,n}(\Sigma^{(h)}\otimes\Sigma^{(h)})$ separately.

Consider $\Sigma^{(h)}\otimes\Sigma^{(h)}$ first. We now introduce a new way to locate an element in a matrix. Divide the $n^{2}\times n^{2}$ matrix $\Sigma^{(h)}\otimes\Sigma^{(h)}$ into $n\times n$ blocks of matrices, each of which is $n\times n$ dimensional. Then $(\Sigma^{(h)}\otimes\Sigma^{(h)})_{\{[x,w],[p,q]\}}$ refers the $[p,q]$th element of the $[x,w]$th block matrix of $\Sigma^{(h)}\otimes\Sigma^{(h)}$, where $x,w,p,q=1,\ldots,n$. It is not difficult to see that

[TABLE]

corresponds to

[TABLE]

We now consider $K_{n,n}(\Sigma^{(h)}\otimes\Sigma^{(h)})$. It is important to recognise that $K_{n,n}$ is a permutation matrix. Left multiplication of $\Sigma^{(h)}\otimes\Sigma^{(h)}$ by $K_{n,n}$ permutes the rows of $\Sigma^{(h)}\otimes\Sigma^{(h)}$. Since $K_{n,n}$ is $n\times n$, we can also divide $K_{n,n}$ into $n\times n$ blocks of matrices, each of which is $n\times n$ dimensional. Since $K_{n,n}$ is also a permutation matrix, its elements can only be either 0 or 1. It is not difficult to see that the $[q,p]$th element of the $[p,q]$th block matrix of $K_{n,n}$ is 1 for $p,q=1,\ldots,n$; all other elements of $K_{n,n}$ are 0. Switch back to the traditional way to locate an element in a matrix. For $p,q=1,\ldots,n$, $[K_{n,n}]_{(p-1)n+q,(q-1)n+p}=1$. This implies that the $((p-1)n+q)$th row of $K_{n,n}(\Sigma^{(h)}\otimes\Sigma^{(h)})$ is actually the $((q-1)n+p)$th row of $\Sigma^{(h)}\otimes\Sigma^{(h)}$. Switch back to the new way to locate an element in a matrix. This says that, for arbitrary $x,w=1,\ldots,n$, the $[q,x]$th element of the $[p,w]$th block matrix of $K_{n,n}(\Sigma^{(h)}\otimes\Sigma^{(h)})$ is the $[p,x]$th element of the $[q,w]$th block matrix of $\Sigma^{(h)}\otimes\Sigma^{(h)}$. Thus

[TABLE]

corresponds to

[TABLE]

Using (A.3) and (A.4), we have

[TABLE]

where the second inequality is due to Cauchy-Schwarz inequality $[\Sigma_{h}]_{i,j}\leq\sqrt{[\Sigma_{h}]_{i,i}}\sqrt{[\Sigma_{h}]_{j,j}}$ using the fact that $\Sigma_{h}$ is a covariance matrix, the fourth last equality uses the fact that $\max_{1\leq h\leq v}\max_{1\leq i\leq n_{h}}[\Sigma_{h}]_{i,i}\leq\max_{1\leq h\leq v}\lambda_{\max}(\Sigma_{h})<\infty$, the second last equality is due to Lemma A.1, and the last equality is due to Assumption 4.3. We hence have

[TABLE]

We now consider the second term of (A.1).

[TABLE]

where the last equality is due to Lemma B.1 in Section B. Inserting (A.5) and (A.6) into (A.1) delivers part (i).

For part (ii), note that for $h=1,\ldots,v$

[TABLE]

Now write

[TABLE]

where the last equality is due to part (i). The last part of part (ii) also follows.

For part (iii), write

[TABLE]

Consider the first term on the right side of (A.7).

[TABLE]

where the second equality is due to part (ii) and the last equality is due to part (i). Consider the second term on the right side of (A.7).

[TABLE]

where the second equality is due to part (ii), and the last equality is due to Lemma 4.1(ii). Part (iii) hence follows.

For part (iv), write

[TABLE]

where $\dot{y}_{t,i}\mathrel{\mathop{\mathchar 58\relax}}=y_{t,i}-\mu_{i}$, and the last equality is due to Lemma B.1 in Section B. We now establish a rate for the first term in (A.8). For some $M>0$,

[TABLE]

We now show \operatorname*{var}\mathinner{\bigl{(}\frac{1}{nT}\sum_{i=1}^{n}\sum_{t=1}^{T}y_{t,i}^{2}\bigr{)}}=O(1/(n^{2-\beta_{1}}T)).

[TABLE]

where the first equality is due to independence over $t$ of Assumption 4.1(i), the first inequality is due to Assumption 4.1(iii), the third and fourth equalities are due to the similar arguments which we used to prove part (i), and the last equality is due to Assumption 4.3. Thus we have

[TABLE]

Substituting this into (A.8) delivers part (iv).

For part (v), we have

[TABLE]

where the last equality is due to part (iii). For part (vi), invoke Lemma B.4 and use that $\max_{1\leq h\leq v}\|\Sigma_{h}^{-1}\|_{F}=O(1)$ in Lemma A.1.

For part (vii), we have

[TABLE]

For part (viii), we have

[TABLE]

For part (ix), we have

[TABLE]

For part (x), we have

[TABLE]

∎

A.2.3 Proof of Theorem 4.1

Proof of Theorem 4.1.

For part (i),

[TABLE]

We consider the first term in (A.9). By inserting terms like $\Sigma_{1}\otimes\tilde{\Sigma}_{2}\otimes\cdots\otimes\tilde{\Sigma}_{v}$ and the triangular inequality, we have

[TABLE]

We first divide the first term of (A.10) by $\prod_{\ell=1}^{v}\|\Sigma_{\ell}\|_{F}$. We have

[TABLE]

We next divide the summand of the second term of (A.10) by $\prod_{\ell=1}^{v}\|\Sigma_{\ell}\|_{F}$. We have for $j=2,\ldots,v-1$

[TABLE]

We finally divide the third term of (A.10) by $\prod_{\ell=1}^{v}\|\Sigma_{\ell}\|_{F}$. We have

[TABLE]

Thus we have

[TABLE]

where the first inequality is due to that $\|\Sigma\|_{F}=\sigma^{2}\prod_{j=1}^{v}\|\Sigma_{j}\|_{F}$ via Lemma B.3, (A.11), (A.12) and (A.13), the first equality is due to Lemma A.1 and Theorem A.1(v)999To see this:

\displaystyle\mathinner{\biggl{(}1+\frac{\max_{1\leq k\leq v}\|\tilde{\Sigma}_{k}-\Sigma_{k}\|_{F}}{\min_{1\leq k\leq v}\|\Sigma_{k}\|_{F}}\biggr{)}}^{v-1}=\mathinner{\biggl{(}1+O_{p}\mathinner{\biggl{(}\sqrt{\frac{\log n}{n^{2-\beta_{1}}T}}\vee\frac{\log n}{T}\biggr{)}}\biggr{)}}^{O(\log n)}=O_{p}(1)

where the last equality could be deduced from the fact that $\lim_{x\to\infty}(1+1/x)^{x}=e$ and $\log^{3}n/T\to 0$. , the second equality is due to Lemma A.1 and Theorem A.1(iv), and the third equality is due to Theorem A.1(v).

We now consider the second term in (A.9).

[TABLE]

where the last equality is due to Theorem A.1(iv). Part (ii)-(vi) of the theorem could be established in a similar manner, so we omit the details. ∎

A.3 Proof of Theorem 5.1

We first give an auxiliary theorem leading to the proof of Theorem 5.1.

A.3.1 Theorem A.2

The following theorem is adapted from Theorem 1 of Kelejian and Prucha (2001).

Theorem A.2**.**

Consider $\{\varepsilon_{T,i}\mathrel{\mathop{\mathchar 58\relax}}1\leq i\leq n,n\geq 1,T\geq 1\}$, an array of real numbers $\{b_{T,i}\mathrel{\mathop{\mathchar 58\relax}}1\leq i\leq n,n\geq 1,T\geq 1\}$ and $Q_{n,T}\mathrel{\mathop{\mathchar 58\relax}}=\sum_{i=1}^{n}\varepsilon_{T,i}^{2}+\sum_{i=1}^{n}b_{T,i}\varepsilon_{T,i}$. Suppose that

1. (i)

$\mathbb{E}[\varepsilon_{T,i}]=0$* for $1\leq i\leq n,n\geq 1,T\geq 1$. Furthermore, for each $n\geq 1,T\geq 1$, $\varepsilon_{T,1},\ldots,\varepsilon_{T,n}$ are (mutually) independent.* 2. (ii)

[TABLE]

for some $\delta>0$. 3. (iii)

[TABLE]

for some absolute positive constant $C$.

Then as $n,T\to\infty$,

[TABLE]

Proof.

We can calculate that

[TABLE]

where the last equality is due to independence of $\varepsilon_{T,i}^{2}$ across $i$. We now show that

[TABLE]

as $n,T\to\infty$. This reduces to verifying the Lyapounov’s condition in Theorem B.1 part (b); that is, for some $\delta>0$,

[TABLE]

We first find an upper bound for \mathbb{E}\mathinner{\!\bigl{\lvert}Y_{T,i}\bigr{\rvert}}^{2+\delta}.

[TABLE]

for absolute positive constants $K_{1}$ and $K_{2}$ for sufficiently large $T$, where the first inequality is due to Loeve’s $c_{r}$ inequality, and the last inequality is due to assumption (ii) of the theorem. Then we have

[TABLE]

as $n,T\to\infty$, where the convergence to 0 relies on assumptions (ii) and (iii) of the theorem. ∎

A.3.2 Proof of Theorem 5.1

Proof of Theorem 5.1.

Write

[TABLE]

We first show that under $H_{0}$ as $n,T\to\infty$,

[TABLE]

Write

[TABLE]

Note that for each $n\geq 1,T\geq 1$, $z_{T,1},\ldots,z_{T,n}$ are (mutually) independent under assumption (b) of the theorem and Assumption 4.1(i). Under $H_{0}$,

[TABLE]

where $\gamma_{z,i}$ is the excess kurtosis of $z_{T,i}$:

[TABLE]

We next calculate $\mathbb{E}[z_{T,i}^{4}]$ in terms of moments of $x_{t,i}$.

[TABLE]

Note that the summand in (A.14) is non-zero only if $t=s=k=\ell$, $t=s\neq k=\ell$, $t=k\neq s=\ell$, $t=\ell\neq k=s$. First, consider the case $t=s=k=\ell$. Collecting all the summands in (A.14) satisfying this, we have

[TABLE]

where $\gamma_{x,t,i}$ is the excess kurtosis of $x_{t,i}$:

[TABLE]

Second, consider the case $t=s\neq k=\ell$. Collecting all the summands in (A.14) satisfying this, we have

[TABLE]

Likewise, for cases $t=k\neq s=\ell$ and $t=\ell\neq k=s$, both sums are $1-1/T$. Substituting (A.15) and (A.16) into (A.14), we have

[TABLE]

whence we have $\gamma_{z,i}=\mathbb{E}[z_{T,i}^{4}]-3=\frac{1}{T^{2}}\sum_{t=1}^{T}\gamma_{x,t,i}$ and

[TABLE]

It remains to verify condition (ii)-(iii) of Theorem A.2. We have

[TABLE]

for large enough $T$ because $\gamma_{x,t,i}>-3$ for all $t$ and $i$ by definition of the excess kurtosis. Hence (iii) of Theorem A.2 is satisfied. Condition (ii) of Theorem A.2 is also satisfied: for some $\delta>0$

[TABLE]

by Theorem B.3 in Section B under assumption (b) of the theorem. Thus we have

[TABLE]

under $H_{0}$ as $n,T\to\infty$, where the second equality is due to

[TABLE]

under assumption (b) of the theorem, and the weak convergence is due to Theorem A.2.

The theorem would follow if we show that

[TABLE]

We now show this.

[TABLE]

where the fourth equality is due to Lemma B.1, the sixth equality is due to Assumption 5.1, and the last equality is due to assumption (a) of the theorem.

For the Wald statistic, write

[TABLE]

We have already shown in the proof of the LM test that under assumptions (a)-(b) of the theorem and under $H_{0}$ as $n,T\to\infty$,

[TABLE]

Display (5.4) would follow if we show that

[TABLE]

We now show this.

[TABLE]

where the fourth equality is due to Lemma B.1 and Theorem 4.1(iv), and the sixth equality is due to Assumption 5.1. ∎

A.4 Proof of Theorem 5.2

Proof.

Write

[TABLE]

whence we have

[TABLE]

We first consider the first term on the right side of (A.18).

[TABLE]

We show that the first term on the right side of (A.19) converges in distribution under the local alternatives.

[TABLE]

Note that for each $n\geq 1,T\geq 1$, $z_{T,1},\ldots,z_{T,n}$ are (mutually) independent under assumption (b) of the theorem and Assumption 4.1(i). Under $H_{1}$,

[TABLE]

We next calculate $\operatorname*{var}(Q_{n,T})$.

[TABLE]

In (A.17), we have already calculated that

[TABLE]

We now calculate $\mathbb{E}[z_{T,i}^{3}]$.

[TABLE]

Backing up, we have

[TABLE]

We now verify conditions (ii) and (iii) of Theorem A.2. For condition (ii), we have already verified in the proof of Theorem 5.1 that \limsup_{T\to\infty}\sup_{n\geq 1}\sup_{1\leq i\leq n}\mathbb{E}\mathinner{\!\bigl{\lvert}z_{T,i}\bigr{\rvert}}^{4+2\delta}<\infty. Next,

[TABLE]

via assumption (b) of the theorem. Thus condition (ii) of Theorem A.2 is met. Finally,

[TABLE]

for large enough $n$ and $T$ because $\gamma_{x,t,i}>-3$ for all $t$ and $i$ by definition of the excess kurtosis. Thus condition (iii) of Theorem A.2 is met.

Thus we have

[TABLE]

as $n,T\to\infty$.

We next show that the second term on the right side of (A.19) is $o_{p}(1)$ under $H_{1}$

[TABLE]

where the second equality is due to Lemma B.1 and Theorem 4.1(iv), and the fourth equality is due to Assumption 5.1.

Backing up, in (A.19), we hence have under $H_{1}$ as $n,T\to\infty$

[TABLE]

We finally consider the second term on the right side of (A.18), noting that

[TABLE]

using a similar device as that in (A.20). Backing up, in (A.18), we hence have under $H_{1}$ as $n,T\to\infty$

[TABLE]

∎

A.5 Proof of Theorem 5.3

Proof of Theorem 5.3.

[TABLE]

We now show that the first term of (A.21) is asymptotically chi square distributed under $H_{0}$. Since $R$ has full row rank $q$ and $\lambda_{\min}(\Sigma)$ is bounded away from zero by an absolute positive constant, $R\Sigma R^{\intercal}$ has full rank $q$. Consider the Cholesky decomposition of $R\Sigma R^{\intercal}=L_{R}L_{R}^{\intercal}$, where $L_{R}$ is a $q\times q$ nonsingular lower triangular matrix with positive diagonal elements. Write

[TABLE]

Note that $L_{R}^{-1}R(y_{1}-\mu),L_{R}^{-1}R(y_{2}-\mu),\ldots,L_{R}^{-1}R(y_{T}-\mu)$ are independent random vectors in $\mathbb{R}^{q}$ with mean zero and variance matrix $I_{q}$. Then we can invoke a version of the multivariate central limit theorem to show $T^{-1/2}\sum_{t=1}^{T}L_{R}^{-1}R(y_{t}-\mu)\xrightarrow{d}N(0,I_{q})$ as $n,T\to\infty$, whence we have $T(R\bar{y}-r)^{\intercal}(R\Sigma R^{\intercal})^{-1}(R\bar{y}-r)\xrightarrow{d}\chi^{2}_{q}$ as $n,T\to\infty$.

We now show that the second term of (A.21) is $o_{p}(1)$ under $H_{0}$.

[TABLE]

We need to find a rate for \mathinner{\!\bigl{\lVert}(R\tilde{\Sigma}R^{\intercal})^{-1}-(R\Sigma R^{\intercal})^{-1}\bigr{\rVert}}_{1}. First, note that

[TABLE]

where the second to last equality is due to Theorem 4.1(v), and the last equality is due to (5.5). Second,

[TABLE]

where the last inequality is due to Lemma B.2 in Section B and the last equality is due to the assumption of the theorem. Then via Lemma B.4 in Section B we have

[TABLE]

Backing up, we have proved that the second term of (A.21) is $o_{p}(1)$ under $H_{0}$. ∎

A.6 Proof of Lemma 5.1

Proof of Lemma 5.1.

The assumptions of the lemma allow us to invoke Theorem 5.1. Thus under $H_{0}\mathrel{\mathop{\mathchar 58\relax}}\mu=\mu_{0}$, as $n,T\rightarrow\infty$,

[TABLE]

This implies that for any unknown $\mu$

[TABLE]

as $n,T\to\infty$, where $z_{\alpha}$ is the upper $\alpha$ percentile of $N(0,1)$.

Invoking Lemma B.5 in Section B with $x=\bar{y}-\mu$ and $S=\tilde{\Sigma}$ yields: For any $\phi\in\mathbb{R}^{n}$

[TABLE]

whence we have

[TABLE]

Multiply both sides by $T$, minus $n$, and divide by $\sqrt{2n}$:

[TABLE]

Thus we assert with confidence $1-\alpha$ that the unknown $\mu$ satisfies simultaneously for all $\phi$ the inequalities:

[TABLE]

as $n,T\to\infty$. ∎

Appendix B Auxiliary Lemmas

Lemma B.1**.**

Suppose Assumption 4.1(i)-(ii) hold. Then we have

[TABLE]

Proof.

Under Assumption 4.1(ii), we have, for $i=1,\ldots,n$, $m=2,3,\ldots$,

[TABLE]

for some absolute positive constant $A$. Now invoke the Bernstein’s inequality in Section B with $\sigma_{0}^{2}=4A^{2}$: For all $\epsilon>0$

[TABLE]

Invoking Corollary B.1 in Section B, we have

[TABLE]

The lemma follows. ∎

We next give two central limit theorems for double-index ($n,T$) processes.

Theorem B.1**.**

1. (a)

Suppose $Y_{n,t}$ is a random variable independent across $1\leq t\leq T$ for $n\geq 1$ and $T\geq 1$. Assume that

[TABLE]

Define

[TABLE]

Assume that $s_{n,T}^{2}>0$ for large enough $n$ and $T$. Suppose the following Lyapounov’s condition holds: For some $\delta>0$,

[TABLE]

Then as $n,T\to\infty$

[TABLE] 2. (b)

Suppose $Y_{T,i}$ is a random variable independent across $1\leq i\leq n$ for $n\geq 1$ and $T\geq 1$. Assume that

[TABLE]

Define

[TABLE]

Assume that $s_{n,T}^{2}>0$ for large enough $n$ and $T$. Suppose the following Lyapounov’s condition holds: For some $\delta>0$,

[TABLE]

Then as $n,T\to\infty$

[TABLE]

Proof.

The proofs can be easily adapted from Lyapounov’s condition for triangular arrays (cf. page 362 Billingsley (1995)). ∎

Theorem B.2** (Bernstein’s inequality).**

We let $Z_{1},\ldots,Z_{T}$ be independent random variables, satisfying for absolute positive constants $A$ and $\sigma_{0}^{2}$

[TABLE]

Let $\epsilon>0$ be arbitrary. Then

[TABLE]

Proof.

Slightly adapted from Bühlmann and van de Geer (2011) p487. ∎

We can use Bernstein’s inequality to establish a rate for the maximum.

Corollary B.1**.**

Suppose via Bernstein’s inequality that we have for $1\leq i\leq n$,

[TABLE]

for some absolute positive constants $K$ and $\sigma_{0}^{2}$. Then

[TABLE]

Proof.

We need to use joint asymptotics $n,T\to\infty$. We shall use the preceding inequality with $\epsilon=(2\log n)/(T\sigma_{0}^{2})$. Fix $\varepsilon>0$. These exist $N_{\varepsilon}\mathrel{\mathop{\mathchar 58\relax}}=2/\varepsilon$, $T_{\varepsilon}$ and $M_{\varepsilon}\mathrel{\mathop{\mathchar 58\relax}}=\max(4K,4\sigma_{0})$ such that for all $n>N_{\varepsilon}$ and $T>T_{\varepsilon}$ we have

[TABLE]

∎

Lemma B.2**.**

Suppose matrix $A$ is real symmetric. Then for any comparable real matrix $B$

[TABLE]

Proof.

First, note that $BAB^{\intercal}$ is Hermitian. By Rayleigh-Ritz theorem, we have

[TABLE]

On the other hand,

[TABLE]

∎

Lemma B.3**.**

For any real matrices $A$ and $B$,

1. (i)

[TABLE] 2. (ii)

[TABLE] 3. (iii)

[TABLE]

Proof.

For part (i),

[TABLE]

For part (ii),

[TABLE]

where the fourth equality is due to the fact that both $A^{\intercal}A$ and $B^{\intercal}B$ are symmetric and positive semidefinite. For part (iii), suppose that $A$ is $m\times n$ and $B$ is $p\times q$.

[TABLE]

∎

Lemma B.4**.**

Let $\hat{\Omega}_{n,j}$ and $\Omega_{n,j}$ be invertible (both possibly stochastic) $n\times n$ square matrices for $j=1,\ldots,m$, where both $n$ and $m$ could be growing. Let $T$ be the sample size. For any matrix norm $\|\cdot\|$, suppose that $\max_{1\leq j\leq m}\|\Omega_{n,j}^{-1}\|=O_{p}(1)$ and $\max_{1\leq j\leq m}\|\hat{\Omega}_{n,j}-\Omega_{n,j}\|=O_{p}(a_{m,n,T})$ for some sequence $a_{m,n,T}$ with $a_{m,n,T}\to 0$ as $m,n,T\to\infty$ simultaneously. Then $\max_{1\leq j\leq m}\|\hat{\Omega}_{n,j}^{-1}-\Omega_{n,j}^{-1}\|=O_{p}(a_{m,n,T})$.

Proof.

The original proof could be found in Saikkonen and Lutkepohl (1996) Lemma A.2.

[TABLE]

Let $v_{j,n,T}$, $z_{j,n,T}$ and $x_{j,n,T}$ denote $\|\Omega_{j,n}^{-1}\|$, $\|\hat{\Omega}_{j,n}^{-1}-\Omega_{j,n}^{-1}\|$ and $\|\Omega_{j,n}-\hat{\Omega}_{j,n}\|$, respectively. From the preceding equation, we have

[TABLE]

whence we have $\max_{1\leq j\leq m}w_{j,n,T}\leq\max_{1\leq j\leq m}x_{j,n,T}=O_{p}(a_{m,n,T})=o_{p}(1)$. We now solve for $z_{j,n,T}$:

[TABLE]

Then we have

[TABLE]

where the second equality is due to the fact that $0\leq v_{j,n,T}w_{j,n,T}\leq 1$ for any $j$. ∎

Theorem B.3**.**

Let $\{x_{t,i}\}$ be a double-index process having zero mean and being independent across $1\leq t\leq T$ for $n\geq 1$ and $T\geq 1$. If there exists $k$, $k\geq 2$, such that

[TABLE]

then we have

[TABLE]

for some absolute positive constant $K$.

Proof.

Slightly adapted from Brillinger (1962). ∎

Lemma B.5** (Generalised Cauchy-Schwarz Inequality).**

For a positive definite matrix $S$ and any vectors $\phi$ and $x$

[TABLE]

Proof.

See Lemma 5.3.2 (p178) of Anderson (1984). ∎