# Generalized integrals and point interactions

Jan Dereziński1 Christian Gaß1 and Błażej Ruba2 1 Department of Mathematical
Methods in Physics, Faculty of Physics, University of Warsaw, Pasteura 5,
02-093 Warszawa, Poland  2 Department of Mathematics, University of
Copenhagen, Universitetsparken 5, DK-2100 Copenhagen Ø, Denmark
<EMAIL_ADDRESS>
###### Abstract

First we recall a method of computing scalar products of eigenfunctions of a
Sturm-Liouville operator. This method is then applied to Macdonald and
Gegenbauer functions, which are eigenfunctions of the Bessel, resp. Gegenbauer
operators. The computed scalar products are well defined only for a limited
range of parameters. To extend the obtained formulas to a much larger range of
parameters, we introduce the concept of a generalized integral. The (standard
as well as generalized) integrals of Macdonald and Gegenbauer functions have
important applications to operator theory. Macdonald functions can be used to
express the integral kernels of the resolvent (Green functions) of the
Laplacian on the Euclidean space in any dimension. Similarly, Gegenbauer
functions appear in Green functions of the Laplacian on the sphere and the
hyperbolic space. In dimensions 1,2,3 one can perturb these Laplacians with a
point potential, obtaining a well defined self-adjoint operator. Standard
integrals of Macdonald and Gegenbauer functions appear in the formulas for the
corresponding Green functions. In higher dimensions the Laplacian perturbed by
point potentials does not exist. However, the corresponding Green function can
be generalized to any dimension by using generalized integrals.

The following notes are a short version of our papers [7, 8].

##  1 Bilinear integrals of eigenfunctions of Sturm-Liouville operators

Consider a Sturm-Liouville operator

$\displaystyle{\mathcal{C}}:=-\rho(r)^{-1}\Big{(}\partial_{r}p(r)\partial_{r}+q(r)\Big{)}$

acting on functions on an interval $]a,b[$. ${\mathcal{C}}$ is formally
symmetric for the bilinear scalar product with the density $\rho$:

$\displaystyle\langle f|g\rangle:=\int_{a}^{b}f(r)g(r)\rho(r)\mathrm{d}r.$

Let us describe a method to compute the scalar product of two eigenfunctions
of ${\mathcal{C}}$.

First, consider eigenfunctions $f_{i}$ corresponding to two distinct
eigenvalues $E_{i}$, $i=1,2$. That is, suppose that
${\mathcal{C}}f_{i}=E_{i}f_{i}.$ Then the following is true:

$\displaystyle\int_{a}^{b}f_{1}(r)f_{2}(r)\rho(r)\mathrm{d}r$
$\displaystyle=\frac{{\mathcal{W}}(b)-{\mathcal{W}}(a)}{E_{1}-E_{2}}$ (1)
$\displaystyle\text{where }{\mathcal{W}}(r):=f_{1}(r)p(r)f_{2}^{\prime}(r)$
$\displaystyle-f_{1}^{\prime}(r)p(r)f_{2}(r)\text{ is the \em Wronskian}.$

(1) is sometimes called Green’s identity, or the integrated Lagrange identity.
Note that if $a,b$ are singular points of the corresponding differential
equation, then the right hand side of (1) can often be easily evaluated.

Using an appropriate limiting procedure we can also often evaluate (1) for
$f=f_{1}=f_{2}$:

$\displaystyle\langle f|f\rangle=\int_{a}^{b}f(r)^{2}\rho(r)\mathrm{d}r.$ (2)

##  2 Bilinear integrals of Macdonald and Gegenbauer functions

The following two families of Sturm-Liouville operators are especially
important for applications [17, 14, 13]:

$\displaystyle\text{the \em Bessel operator:}\qquad{\mathcal{B}}_{\alpha}$
$\displaystyle:=-\frac{1}{r}\partial_{r}r\partial_{r}+\frac{\alpha^{2}}{r^{2}},$
$\displaystyle\text{and the \em Gegenbauer
operator}:\quad{\mathcal{G}}_{\alpha}$
$\displaystyle:=-(1-w^{2})^{-\alpha}\partial_{w}(1-w^{2})^{\alpha+1}\partial_{w}.$

The modified Bessel equation is the eigenequation of ${\mathcal{B}}_{\alpha}$
with eigenvalue $-1$. The (standard) Bessel equation is its eigenequation for
eigenvalue $1$. The separation of variables in the Laplacian on the Euclidean
space $\mathbb{R}^{d}$ leads to the Bessel operator on $[0,\infty[$ with
density $\rho=2r$ and $\alpha=\frac{d-2}{2}$.

The Gegenbauer equation is the eigenequation of ${\mathcal{G}}_{\alpha}$ with
eigenvalue $\lambda^{2}-\big{(}\alpha+\frac{1}{2}\big{)}^{2}$. The Gegenbauer
operator on $[-1,1]$ with density $\rho(w)=(1-w^{2})^{\frac{d-2}{2}}$ arises
when we separate variables in the Laplacian on the sphere $\mathbb{S}^{d}$. We
obtain the Gegenbauer operator on $[1,\infty[$ with density
$\rho(w)=(w^{2}-1)^{\frac{d-2}{2}}$ when we separate variables in the
Laplacian on the hyperbolic space $\mathbb{H}^{d}$.

Scaled Macdonald functions $K_{\alpha}(br)$ are exponentially decaying
eigenfunctions of the Bessel operator with eigenvalue $-b^{2}$. Applying the
above method for $a>0$, $b>0$ we find

$\displaystyle\int_{0}^{\infty}K_{\alpha}(ar)K_{\alpha}(br)2r\mathrm{d}r$
$\displaystyle=$
$\displaystyle\frac{\pi\big{(}(a/b)^{\alpha}-(b/a)^{\alpha}\big{)}}{\sin(\pi\alpha)(a^{2}-b^{2})},\quad\begin{array}[]{c}|\operatorname{Re}(\alpha)|<1,\\\
\alpha\neq 0;\end{array}$ (5)
$\displaystyle\int_{0}^{\infty}K_{0}(ar)K_{0}(br)2r\mathrm{d}r$
$\displaystyle=$ $\displaystyle\frac{2\ln\frac{a}{b}}{a^{2}-b^{2}},$ (6)
$\displaystyle\int_{0}^{\infty}K_{\alpha}(br)^{2}2r\mathrm{d}r$
$\displaystyle=$
$\displaystyle\frac{\pi\alpha}{b^{2}\sin(\pi\alpha)},\qquad\begin{array}[]{c}|\operatorname{Re}(\alpha)|<1,\\\
\alpha\neq 0;\end{array}$ (9)
$\displaystyle\int_{0}^{\infty}K_{0}(br)^{2}2r\mathrm{d}r$ $\displaystyle=$
$\displaystyle\frac{1}{b^{2}}.$ (10)

The first formula follows directly by Green’s identity. The next three
identities are obtained by applying the de l’Hôpital rule to $\alpha=0$ and
$a=b$.

The identities (5)–(10) can be found in standard collections of integrals,
such as [9].

The Gegenbauer equation is the special case of the hypergeometric equation
with the symmetry $w\to-w$ and the singular points at $-1,1,\infty$:

$\displaystyle\Bigg{(}(1-w^{2})\partial_{w}^{2}-2(1+\alpha)w\partial_{w}+\lambda^{2}-\Big{(}\alpha+\frac{1}{2}\Big{)}^{2}\Bigg{)}f(w)=0.$

In the present context, the Gegenbauer equation is arguably more convenient
than the equivalent, but more frequently encountered associated Legendre
equation.

We will use two kinds of Gegenbauer functions: one is characterized by its
asymptotics $\sim\frac{1}{\Gamma(\alpha+1)}$ at $1$:

$\displaystyle{\bf S}_{\alpha,\pm\lambda}(w)$
$\displaystyle:=\sum_{j=0}^{\infty}\frac{\big{(}\frac{1}{2}+\alpha+\lambda\big{)}_{j}\big{(}\frac{1}{2}+\alpha-\lambda)\big{)}_{j}}{\Gamma(\alpha+1+j)j!}\Big{(}\frac{1-w}{2}\Big{)}^{j},$

where $(z)_{n}$ is the _Pochhammer symbol._ The other has the asymptotics
$\sim\frac{1}{w^{\frac{1}{2}+\alpha+\lambda}\Gamma(\lambda+1)}$ at $\infty$:

$\displaystyle\mathbf{Z}_{\alpha,\lambda}(w):$ $\displaystyle=\frac{1}{(w\pm
1)^{\frac{1}{2}+\alpha+\lambda}}\sum_{j=0}^{\infty}\frac{\big{(}\frac{1}{2}+\lambda\big{)}_{j}\big{(}\frac{1}{2}+\lambda+\alpha\big{)}_{j}}{\Gamma(\lambda+1)(1+2\lambda)_{j}j!}\Big{(}\frac{2}{1\pm
w}\Big{)}^{j}.$

We note the identities

$\displaystyle{\bf S}_{\alpha,\lambda}(w)={\bf S}_{\alpha,-\lambda}(w),$
$\displaystyle\quad{\bf Z}_{\alpha,\lambda}(w)=\frac{{\bf
Z}_{-\alpha,\lambda}(w)}{(w^{2}-1)_{\bullet}^{\alpha}},$

as well as the slightly more subtle Whipple identity:

$\displaystyle{\bf Z}_{\alpha,\lambda}(w)$
$\displaystyle:=(w^{2}-1)_{\bullet}^{-\frac{1}{4}-\frac{\alpha}{2}-\frac{\lambda}{2}}{\bf
S}_{\lambda,\alpha}\left(\frac{w}{(w^{2}-1)_{\bullet}^{\frac{1}{2}}}\right),$
$\displaystyle{\bf S}_{\alpha,\lambda}(w)$
$\displaystyle:=(w^{2}-1)_{\bullet}^{-\frac{1}{4}-\frac{\alpha}{2}-\frac{\lambda}{2}}{\bf
Z}_{\lambda,\alpha}\left(\frac{w}{(w^{2}-1)_{\bullet}^{\frac{1}{2}}}\right),\qquad\operatorname{Re}(w)>0,$

where $(w^{2}-1)_{\bullet}^{\alpha}:=(w-1)^{\alpha}(w+1)^{\alpha}$ and we use
the principal branch of the power function.

Here are the basic bilinear integrals of Gegenbauer functions. We assume
$|\operatorname{Re}(\alpha)|<1$, $\alpha\neq 0$ and
$\operatorname{Re}(\lambda)>0$:

$\displaystyle\int_{-2}^{2}{\bf S}_{\alpha,\mathrm{i}\beta_{1}}(w){\bf
S}_{\alpha,\mathrm{i}\beta_{2}}(w)(1-w^{2})^{\alpha}\mathrm{d}2w$ (11)
$\displaystyle=$
$\displaystyle\frac{2^{2\alpha+2}}{(\beta_{1}^{2}-\beta_{2}^{2})\sin\pi\alpha}\Big{(}\frac{\cosh(\pi\beta_{1})}{\Gamma(\frac{1}{2}+\alpha-\mathrm{i}\beta_{2})\Gamma(\frac{1}{2}+\alpha+\mathrm{i}\beta_{2})}-(\beta_{1}\leftrightarrow\beta_{2})\Big{)}$
$\displaystyle\int_{2}^{\infty}{\bf Z}_{\alpha,\lambda_{1}}(w){\bf
Z}_{\alpha,\lambda_{2}}(w)(w^{2}-1)^{\alpha}\mathrm{d}2w$ (12)
$\displaystyle=$
$\displaystyle\frac{2^{\lambda_{1}+\lambda_{2}+1}}{(\lambda_{1}^{2}-\lambda_{2}^{2})\sin\pi\alpha}\Big{(}\frac{1}{\Gamma(\frac{1}{2}-\alpha+\lambda_{1})\Gamma(\frac{1}{2}+\alpha+\lambda_{2})}-(\lambda_{1}\leftrightarrow\lambda_{2})\Big{)}.$

Applying the de l’Hôpital rule we extend these identities to $\alpha=0$,
$\lambda_{1}=\lambda_{2}$, and $\beta_{1}=\beta_{2}$, see [7]. In contrast to
the well-known integrals over Macdonald functions, the integrals (11) and (12)
seem to be new.

##  3 Generalized integral

The integrals (5)–(10) and (11)–(12) are divergent for
$|\operatorname{Re}(\alpha)|\geq 1$. Using the generalized integral, (5)–(10)
and (12) can be extended to all $\alpha\in\mathbb{C}$, and (11) can be
extended to $\operatorname{Re}(\alpha)>-1$. The concept of the generalized
integral can be traced back to Hadamard [10, 11] and Riesz [16]; see [12, 15]
for modern expositions.

Note that several variations of the generalized integral are possible. In our
note we restrict ourselves to functions which are non-integrable near a finite
point and have a finite number of homogeneous singularities at this point. In
some of the above references one can also find other variations of the
generalized integral applicable to functions which are non-integrable near
infinity and/or exhibit _almost homogeneous_ singularities, i.e.,
singularities that are homogeneous up to powers of logarithms.

Without loss of generality, we can put the singular point to 0. We then say
that a function $f$ on $]0,\infty[$ is integrable in the generalized sense if
it is integrable on $]1,\infty[$ and if there exists a finite set
$\Omega\subset\mathbb{C}$ and $f_{k}\in\mathbb{C}$, $k\in\Omega$, such that
$f-\sum\limits_{k\in\Omega}f_{k}r^{k}$ is integrable on $]0,1[$. We define the
generalized integral as

$\displaystyle\mathrm{gen}\int_{0}^{\infty}f(r)\mathrm{d}r:=$
$\displaystyle\sum_{k\in\Omega\backslash\\{-1\\}}\frac{f_{k}}{k+1}+\int_{0}^{1}\Big{(}f(r)-\sum_{k\in\Omega}f_{k}r^{k}\Big{)}\mathrm{d}r+\int_{1}^{\infty}f(r)\mathrm{d}r.$

For $f\in L^{1}[0,\infty[$ the generalized and standard integrals coincide:

$\displaystyle\mathrm{gen}\int_{0}^{\infty}f(r)\mathrm{d}r=\int_{0}^{\infty}f(r)\mathrm{d}r.$

If $f_{-1}\neq 0$, the generalized integral has a scaling anomaly but it is
always invariant with respect to a power transformation:

$\displaystyle\mathrm{gen}\int_{0}^{\infty}f(r)\mathrm{d}r$
$\displaystyle=\mathrm{gen}\int_{0}^{\infty}f(\alpha
u)\,\alpha\mathrm{d}u+f_{-1}\ln(\alpha),\quad$
$\displaystyle\qquad\mathrm{gen}\int_{0}^{\infty}f(r)\mathrm{d}r$
$\displaystyle=\mathrm{gen}\int_{0}^{\infty}f(u^{\alpha})\,\alpha
u^{\alpha-1}\mathrm{d}u.$

Under more general coordinate transformations $g(u)$, which are smooth and
preserve the leading scaling behavior of the integrand (i.e., $g(0)=0$ and
$g^{\prime}(0)\neq 0$), the generalized integral transforms as

$\displaystyle\mathrm{gen}\int_{0}^{\infty}f(g(u))g^{\prime}(u)\mathrm{d}u-\mathrm{gen}\int_{0}^{\infty}f(r)\mathrm{d}r$
(13) $\displaystyle=$ $\displaystyle-f_{-1}\ln
g^{\prime}(0)+\sum_{l\in(\mathbb{N}+1)\cap\Omega}\frac{f_{-l}}{(l-1)(l-1)!}\left.\frac{\mathrm{d}^{l-1}}{\mathrm{d}u^{l-1}}\left(\frac{u}{g(u)}\right)^{l-1}\right|_{u=0}.$

That is to say, the generalized integral transforms non-trivially under such
transformations if there is an $n\in\mathbb{N}$ such that $f_{-n}\neq 0$. For
this reason, we call the generalized integral anomalous if $f_{-n}\neq 0$ for
some $n\in\mathbb{N}$.

Non-anomalous generalized integrals have much better properties than anomalous
ones. They are often easy to compute: one just applies analytic continuation.

A systematic analysis of properties of the generalized integral and a proof of
(13) can be found in [7]. To our knowledge, the formula (13) is new.

##  4 Bilinear generalized integrals of Macdonald and Gegenbauer functions

Let $a,b>0$. For $\alpha\not\in\mathbb{Z}$ generalized integrals of Macdonald
functions are analytic continuations of standard integrals:

$\displaystyle\mathrm{gen}\int_{0}^{\infty}K_{\alpha}(ar)K_{\alpha}(br)2r\mathrm{d}r$
$\displaystyle=\frac{\pi}{\sin(\pi\alpha)}\frac{\big{(}\frac{a}{b}\big{)}^{\alpha}-\big{(}\frac{b}{a}\big{)}^{\alpha}}{a^{2}-b^{2}},$
$\displaystyle\mathrm{gen}\int_{0}^{\infty}K_{\alpha}(br)^{2}2r\mathrm{d}r$
$\displaystyle=\frac{\pi\alpha}{b^{2}\sin(\pi\alpha)}.$

They have poles at $\alpha\in\mathbb{Z}$.

For $\alpha\in\mathbb{Z}$ the generalized integrals are anomalous and more
complicated to compute. In particular, they do not coincide with the finite
parts of the above expressions:

$\displaystyle\mathrm{gen}\int_{0}^{\infty}K_{\alpha}(ar)K_{\alpha}(br)$
$\displaystyle
2r\mathrm{d}r=(-1)^{\alpha}2\frac{\big{(}\frac{a}{b}\big{)}^{\alpha}\ln\big{(}\frac{a}{2}\big{)}-\big{(}\frac{b}{a}\big{)}^{\alpha}\ln\big{(}\frac{b}{2}\big{)}}{a^{2}-b^{2}}$
$\displaystyle-\frac{(-1)^{\alpha}}{ab}\sum_{k=0}^{|\alpha|-1}\Big{(}\frac{a}{b}\Big{)}^{2k-|\alpha|+1}\big{(}\psi(1+k)+\psi(|\alpha|-k)\big{)};$
$\displaystyle\mathrm{gen}\int_{0}^{\infty}K_{\alpha}(br)^{2}2r\mathrm{d}r=$
$\displaystyle\frac{(-1)^{\alpha}}{b^{2}}\Big{(}|\alpha|\ln\big{(}\tfrac{b^{2}}{4}\big{)}+1+2|\alpha|\big{(}1-\psi(1+|\alpha|)\big{)}\Big{)}.$

Similarly we can be compute generalized bilinear integrals of ${\bf
S}_{\alpha,\mathrm{i}\beta}$, ${\bf Z}_{\alpha,\lambda}$. All these formulas
together with their derivations can be found in [7].

##  5 Convergence of Gegenbauer functions and their integrals

As $\beta,\lambda\to\infty$, Gegenbauer functions converge to Macdonald
functions in the following sense:

$\displaystyle\frac{\pi\mathrm{e}^{-\pi\beta}(\sin\theta)^{\alpha+\frac{1}{2}}}{2^{\alpha}\theta^{\alpha+\frac{1}{2}}}{\bf
S}_{\alpha,\pm\mathrm{i}\beta}(-\cos\theta)$
$\displaystyle=(\theta\beta)^{-\alpha}K_{\alpha}(\beta\theta)\big{(}1+O(\beta^{-1})\big{)};$
$\displaystyle\frac{\sqrt{\pi}\Gamma(\tfrac{1}{2}-\alpha+\lambda)(\sinh\theta)^{\alpha+\frac{1}{2}}}{2^{\lambda+\frac{1}{2}}\theta^{\alpha+\frac{1}{2}}}{\bf
Z}_{\alpha,\lambda}(\cosh\theta)$
$\displaystyle=(\lambda\theta)^{-\alpha}K_{\alpha}(\lambda\theta)\big{(}1+O(\lambda^{-1})\big{)}.$

The generalized integrals of Gegenbauer functions converge to the
corresponding generalized integrals of Macdonald functions:

$\displaystyle\frac{\pi^{2}{\mathrm{e}}^{-2\pi\beta}\beta^{2\alpha}}{2^{2\alpha}}\;\mathrm{gen}\int_{-2}^{2}{\bf
S}_{\alpha,\mathrm{i}\beta}(w)^{2}(1-w^{2})^{\alpha}\mathrm{d}2w=\Big{(}1+\mathcal{O}\big{(}\tfrac{1}{\beta}\big{)}\Big{)}\;\mathrm{gen}\int_{0}^{\infty}K_{\alpha}(\beta
r)^{2}2r\mathrm{d}r;$
$\displaystyle\frac{\pi\Gamma\big{(}\tfrac{1}{2}+\alpha+\lambda\big{)}^{2}}{2^{2\lambda+1}\lambda^{2\alpha}}\;\mathrm{gen}\int_{2}^{\infty}{\bf
Z}_{\alpha,\lambda}(w)^{2}(w^{2}-1)^{\alpha}\mathrm{d}2w=\Big{(}1+\mathcal{O}\big{(}\tfrac{1}{\lambda}\big{)}\Big{)}\;\mathrm{gen}\int_{0}^{\infty}K_{\alpha}(\lambda
r)^{2}2r\mathrm{d}r.$

The convergence of these generalized integrals is straightforward in the non-
anomalous case. In the anomalous case one has to choose the variables
carefully, which we did:

$\displaystyle 2r\mathrm{d}r=\mathrm{d}r^{2},\quad 2(\cosh r-1)\simeq
r^{2},\quad 2(1-\cos r)\simeq r^{2}.$

Note that the generalized integral is invariant with respect to the change of
variables $r\to r^{2}$, but not with respect to scaling. Proofs of these
convergence statements can be found in [7].

##  6 Laplacian on the Euclidean space, the hyperbolic space and the sphere

In the remaining part of our manuscript we describe an application of
generalized integrals to operator theory. These applications will involve
point interactions of the Laplacian. We start by recalling some basic
information about Green functions of Laplacians. The following formulas are
well-known, confer for example [5, 8].

Consider the Laplacian $\Delta_{d}$ on the Euclidean space $\mathbb{R}^{d}$.
Let $G_{d}(z;x,x^{\prime})$ be the Euclidean Green function, that is the
integral kernel of the resolvent $(-z-\Delta_{d})^{-1}$. For
$\operatorname{Re}\beta>0$, we have

$G_{d}(-\beta^{2};x,x^{\prime})=\frac{1}{(2\pi)^{\frac{d}{2}}}\Big{(}\frac{\beta}{|x-x^{\prime}|}\Big{)}^{\frac{d}{2}-1}K_{\frac{d}{2}-1}\big{(}\beta|x-x^{\prime}|\big{)}.$

The hyperbolic space is

$\mathbb{H}^{d}:=\\{x\in\mathbb{R}^{1,d}\mkern 5.0mu|\mkern 5.0mu[x|x]=1\\}$

where $[x|y]=x^{0}y^{0}-x^{1}y^{1}-\dots-x^{d}y^{d}$ is the Minkowskian
pseudoscalar product. The hyperbolic distance between
$x,x^{\prime}\in\mathbb{H}^{d}$ is given by $\cosh(r)=[x|x^{\prime}]$.

Let $\Delta_{d}^{\mathrm{h}}$ denote the Laplace-Beltrami operator on
$\mathbb{H}^{d}$ and let $G_{d}^{\mathrm{h}}(z;x,x^{\prime})$ be the
hyperbolic Green function, that is, the integral kernel of
$\big{(}-z-\Delta_{d}^{\mathrm{h}}-\big{(}\tfrac{d-1}{2}\big{)}^{2}\big{)}^{-1}$.
Then

$G_{d}^{\mathrm{h}}\Big{(}-\beta^{2};x,x^{\prime}\Big{)}=\frac{\sqrt{\pi}\Gamma(\frac{d-1}{2}+\beta)}{\sqrt{2}(2\pi)^{\frac{d}{2}}2^{\beta}}{\bf
Z}_{\frac{d}{2}-1,\beta}\big{(}[x|x^{\prime}]\big{)}.$

The unit sphere is

$\mathbb{S}^{d}:=\\{x\in\mathbb{R}^{1+d}\mkern 5.0mu|\mkern 5.0mu(x|x)=1\\},$

where $(x|y)=x^{0}y^{0}+x^{1}y^{1}+\dots+x^{d}y^{d}$ is the Euclidean scalar
product. The spherical distance between $x,x^{\prime}\in\mathbb{S}^{d}$ is
given by $\cos(r)=(x|x^{\prime})$. Let $\Delta_{d}^{\mathrm{s}}$ denote the
Laplace-Beltrami operator on $\mathbb{S}^{d}$. Let
$G_{d}^{\mathrm{s}}(z;x,x^{\prime})$ be the spherical Green function, that is,
the integral kernel of
$\big{(}-z-\Delta_{d}^{\mathrm{s}}+\big{(}\tfrac{d-1}{2}\big{)}^{2}\big{)}^{-1}$.
Then

$G_{d}^{\mathrm{s}}(-\beta^{2};x,x^{\prime})=\frac{\Gamma(\frac{d}{2}-\frac{1}{2}+\mathrm{i}\beta)\Gamma(\frac{d}{2}-\frac{1}{2}-\mathrm{i}\beta)}{2^{d}\pi^{\frac{d}{2}}}{\bf
S}_{\frac{d}{2}-1,\mathrm{i}\beta}\big{(}-(x|x^{\prime})\big{)}.$

We remark that the integral kernels of spectral projections of the discussed
operators may be expressed explicitly in terms of special functions of the
same type as for resolvents, see e.g. [8].

Let us also mention some basic properties of operators. Let $H$ be a self-
adjoint operator and let $G(-\rho):=(\rho+H)^{-1}$ be its resolvent. Then
$G(-\rho)$ satisfies

$\displaystyle(H+\rho)G(-\rho)$ $\displaystyle={\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}},$
(14) $\displaystyle G(-\rho)^{*}$ $\displaystyle=G(-\overline{\rho}),$ (15)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\rho}G(-\rho)$
$\displaystyle=-G(-\rho)^{2}.$ (16)

(16) is called the resolvent formula in the differential form.

##  7 Laplacian with point interactions

Let us now try to define the Laplacian with a perturbation localized in a
single point. For simplicity, we display details only for the case of the
Euclidean space. We look for an operator $-\Delta_{d}^{\gamma}$, which is a
self-adjoint extension of $-\Delta_{d}$ restricted to
$C_{\mathrm{c}}^{\infty}(\mathbb{R}^{d}\backslash\\{0\\})$. The parameter
$\gamma\in\mathbb{R}\cup\\{\infty\\}$ will parametrize these self-adjoint
extensions.

Actually, instead of $-\Delta_{d}^{\gamma}$ it is more convenient to look for
its resolvent

$\displaystyle G_{d}^{\gamma}(-\rho)=(-\Delta_{d}^{\gamma}+\rho)^{-1}.$

By the conditions (14), (15) and (16) its integral kernel
$G_{d}^{\gamma}(-\rho,x,x^{\prime})$ should satisfy

$\displaystyle(-\Delta_{x}+\rho)G_{d}^{\gamma}(-\rho,x,x^{\prime})$
$\displaystyle=\delta(x-x^{\prime}),\quad x\neq 0,$ $\displaystyle
G_{d}^{\gamma}(-\rho,x,x^{\prime})$
$\displaystyle=G_{d}^{\gamma}(-\rho,x^{\prime},x),$
$\displaystyle\partial_{\rho}G_{d}^{\gamma}(-\rho,x,x^{\prime})$
$\displaystyle=-\int
G_{d}^{\gamma}(-\rho,x,y)G_{d}^{\gamma}(-\rho,y,x^{\prime})\mathrm{d}y.$

These conditions are solved by a Krein-type resolvent

$\displaystyle G_{d}^{\gamma}(-\rho,x,x^{\prime})$
$\displaystyle=G_{d}(-\rho,x,x^{\prime})+\frac{G_{d}(-\rho,x,0)G_{d}(-\rho,0,x^{\prime})}{\gamma+\Sigma_{d}(\rho)},$

where

$\partial_{\rho}\Sigma_{d}(\rho)=\int_{\mathbb{R}^{d}}G_{d}(-\rho,0,y)^{2}\mathrm{d}y,$

and $\gamma$ is an arbitrary constant. In dimensions $d=1,2,3$ the above
integral is finite and we obtain

$\Sigma_{d}(\beta^{2})=\begin{cases}-\frac{1}{2\beta}&d=1;\\\
\frac{\ln(\beta^{2})}{4\pi}&d=2;\\\ \frac{\beta}{4\pi}&d=3.\end{cases}$

Thus we obtain formulas for the Green functions with a point potential in
dimensions $d=1,2,3$ parametrized by a real parameter $\gamma\in\mathbb{R}$:

$G_{d}^{\gamma}(-\beta^{2};x,x^{\prime})=\begin{cases}\frac{\mathrm{e}^{-\beta|x-x^{\prime}|}}{2\beta}+\frac{\mathrm{e}^{-\beta|x|}\mathrm{e}^{-\beta|x^{\prime}|}}{(2\beta)^{2}\big{(}\gamma-\frac{1}{2\beta}\big{)}},&d=1;\\\\[8.61108pt]
\frac{K_{0}(\beta|x-x^{\prime}|)}{2\pi}+\frac{K_{0}(\beta|x|)K_{0}(\beta|x^{\prime}|)}{(2\pi)^{2}(\gamma+\frac{\ln\beta^{2}}{4\pi})},&d=2;\\\\[8.61108pt]
\frac{\mathrm{e}^{-\beta|x-x^{\prime}|}}{4\pi|x-x^{\prime}|}+\frac{\mathrm{e}^{-\beta|x|}\mathrm{e}^{-\beta|x^{\prime}|}}{(4\pi)^{2}|x||x^{\prime}|(\gamma+\frac{\beta}{4\pi})},&d=3.\end{cases}$

For dimensions $d=1,3$ we used the fact that for half-integer parameters the
Macdonald function reduces to elementary functions. Of course, the above
construction is well known from the literature [1, 2, 3] and often used in the
physics literature.

The operators $-\Delta_{d}^{\gamma}$, strictly speaking, have no analogs for
$d\geq 4$. However, the functions $G_{d}^{\gamma}(-\rho;x,x^{\prime})$ can be
generalized to $d\geq 4$ using generalized integrals:

$\partial_{\rho}\Sigma_{d}(\rho)=\frac{(\beta^{2})^{\frac{d}{2}-1}\pi^{\frac{d}{2}}}{(2\pi)^{d}\Gamma(\frac{d}{2})}\mathrm{gen}\int_{0}^{\infty}K_{\frac{d}{2}-1}(\sqrt{\rho}r)^{2}2r\mathrm{d}r.$

We obtain

$\Sigma_{d}(\beta^{2})=\begin{cases}\frac{(-1)^{\frac{d+1}{2}}\beta^{d-2}}{(4\pi)^{\frac{d-1}{2}}2(\frac{1}{2})_{\frac{d-1}{2}}}&d\text{
odd};\\\
\frac{(-1)^{\frac{d}{2}+1}\beta^{d-2}}{(4\pi)^{\frac{d}{2}}\big{(}\frac{d}{2}-1\big{)}!}\left(2-2\psi\big{(}\tfrac{d}{2}\big{)}+\ln\tfrac{\beta^{2}}{4}\right)&d\text{
even}.\end{cases}$

Thus for each dimension $d$ we obtain a family of Green functions

$\displaystyle G_{d}^{\gamma}(-\beta^{2};x,x^{\prime})$
$\displaystyle=\frac{1}{(2\pi)^{\frac{d}{2}}}\Big{(}\frac{\beta}{|x-x^{\prime}|}\Big{)}^{\frac{d}{2}-1}K_{\frac{d}{2}-1}\big{(}\beta|x-x^{\prime}|\big{)}$
$\displaystyle+$
$\displaystyle\frac{1}{(2\pi)^{d}}\Big{(}\frac{\beta^{2}}{|x||x^{\prime}|}\Big{)}^{\frac{d}{2}-1}\frac{K_{\frac{d}{2}-1}(\beta|x|)K_{\frac{d}{2}-1}(\beta|x^{\prime}|)}{\gamma+\Sigma_{d}(\beta^{2})}.$
(17)

describing point interaction of strength controlled by the parameter $\gamma$.

A similar analysis can be performed for the hyperbolic and spherical Green
functions in all dimensions [8].

One can ask what is the meaning of $G_{d}^{\gamma}(-\beta^{2};x,x^{\prime})$
in dimensions $d\geq 4$, for wich it is not the kernel of a bounded operator
(and in particular, not a kernel of a resolvent of a self-adjoint operator).
Our expectation is as follows. Suppose that $V$ is a potential on
$\mathbb{R}^{d}$, $\mathbb{H}^{d}$ or $\mathbb{S}^{d}$, possibly strong but
with a small support. Consider the the Schrödinger operator $-\Delta_{d}+V$.
Let

$G^{V}(-\beta^{2})=(\beta^{2}-\Delta_{d}+V)^{-1}$

be its resolvent with the integral kernel $G^{V}(-\beta^{2};x,x^{\prime})$.
Then far from the support of $V$ we can approximate
$G^{V}(-\beta^{2};x,x^{\prime})$ by $G^{\gamma}(-\beta^{2};x,x^{\prime})$ as
in (17), possibly adding to $\Sigma_{d}(\beta^{2})$ a polynomial in the energy
of degree $<\frac{d-1}{2}$ if $d$ is odd and $<\frac{d-2}{2}$ if $d$ is even.
Thus possible infrared behaviors of $G^{V}(-\beta^{2};x,x^{\prime})$ are
controlled by a finite number of parameters (coefficients of the above
mentioned polynomials).

The choice given by the zero polynomial that is defined by the above
generalized integrals can be viewed as a “standard reference point”. The
generalized integral computed in different coordinates, which are related to
$r^{2}$ by a well-behaved coordinate transformation in the sense of (13), will
yield a different choice of polynomial.

The observation that a well-behaved change of variables will, by (13), only
change this polynomial is particularly important in the curved cases, where we
chose $2(\cosh r-1)$, resp. $2(1-\cos r)$ as integration variables instead of
$r^{2}$.

This is analogous to the renormalization of quantum quantum field theory using
dimensional regularization, see e.g. [4].

The above situation resembles the idea often expressed in the context of
quantum field theory and of the theory of critical phenomena, attributed to
Keneth Wilson: for large distances correlation functions have a universal
behavior independent of the details of the interaction, described by few
parameters.

The work of J.D. and G.G. was supported by National Science Center (Poland)
under the Grant UMO-2019/35/B/ST1/01651.

## References

## References

  * [1] Albeverio S, Gesztesy F, Høegh-Krohn R and Holden H 1988 _Solvable Models in Quantum Mechanics_ (Berlin, Heidelberg: Springer)
  * [2] Albeverio S and Kurasov P 2000 _Singular Perturbations of Differential Operators: Solvable Schrödinger-type Operators_ (Cambridge: Cambridge Univ. Press)
  * [3] Berezin F A and Faddeev L D 1961 Remark on the Schrödinger equation with singular potential Dokl. Akad. Nauk. SSSR 137 1011
  * [4] Bollini C G and Giambiagi J J 1996 Dimensional regularization in configuration space Phys. Rev. D 53 5761
  * [5] Cohl H S, Dang T H and Dunster T M 2018 Fundamental Solutions and Gegenbauer Expansions of Helmholtz Operators in Riemannian Spaces of Constant Curvature SIGMA 14 136
  * [6] Dereziński J 2014 Hypergeometric type functions and their symmetries Annales Henri Poincaré 15 1569
  * [7] Dereziński J, Ruba B and Gaß C 2023 Generalized integrals of Macdonald and Gegenbauer functions Preprint arXiv:2304.06515
  * [8] Dereziński J, Ruba B and Gaß C Point potentials on Euclidean space, hyperbolic space and sphere in any dimension in preparation
  * [9] Gradshteyn I S and Ryzhik I M 1980 Table of Integrals, Series, and Products (New York: Academic Press)
  * [10] Hadamard J 1923 _Lectures on Cauchy’s problem in linear partial differential equations_ (New York: Dover Phoenix editions, Dover Publications)
  * [11] Hadamard J 1932 _Le problème de Cauchy et les équations aux dérivées partielles linéaires hyperboliques_ (Paris: Hermann & Cie)
  * [12] Lesch M 1997 Differential operators of Fuchs type, conical singularities, and asymptotic methods (Vieweg+Teubner Verlag)
  * [13] NIST Digital Library of Mathematical Functions Release 1.1.6 of 2022-06-30 ed F W J Olver et al
  * [14] Olver F W J 1974 _Asymptotics and special functions_ (New York: Academic Press)
  * [15] Paycha S 2012 _University Lecture Series: Regularised Integrals, Sums and Traces: An Analytic Point of View_ (Providence, Rhode Island: American Mathematical Society) volume 59
  * [16] Riesz M 1949 L’intégrale de Riemann-Liouville et le problème de Cauchy Acta Mathematica 81 1
  * [17] Whittaker E T and Watson G N 1996 _A Course of Modern Analysis: An Introduction to the General Theory of Infinite Processes and of Analytic Functions; with an Account of the Principal Transcendental Functions_ (Cambridge: Cambridge Mathematical Library Series) (1st ed 1902 Cambridge: Cambridge University Press)


# Adaptive Testing for Alphas in High-dimensional Factor Pricing Models

Qiang Xia
College of Mathematics and Informatics, South China Agricultural University
and
Xianyang Zhang
Department of Statistics, Texas A&M University  Address correspondence to
Xianyang Zhang<EMAIL_ADDRESS>
###### Abstract

This paper proposes a new procedure to validate the multi-factor pricing
theory by testing the presence of alpha in linear factor pricing models with a
large number of assets. Because the market’s inefficient pricing is likely to
occur to a small fraction of exceptional assets, we develop a testing
procedure that is particularly powerful against sparse signals. Based on the
high-dimensional Gaussian approximation theory, we propose a simulation-based
approach to approximate the limiting null distribution of the test. Our
numerical studies show that the new procedure can deliver a reasonable size
and achieve substantial power improvement compared to the existing tests under
sparse alternatives, and especially for weak signals.

Keywords: Cross-sectional dependence, Factor pricing model, High-dimensional
intercepts, Sparse alternatives

##  1 Introduction

As one of the most fundamental results in empirical finance, the factor
pricing model postulates how financial returns are related to market risks and
has many critical applications such as portfolio selection, fund performance
evaluation, and corporate budgeting. There are several well-known factor
models, including the capital asset pricing model (CAPM) proposed by Sharpe
(1964) and Lintner (1965), the arbitrage pricing theory (APT) model by Ross
(1976), the Fama-French three-factor model by Fama and French (1993), and the
Fama-French five-factor model by Fama and French (2015).

Let $y_{it}$ be the excess return of the $i$th asset at time $t$ and
$X_{t}=(x_{1t},\dots,x_{rt})^{\top}\in\mathbb{R}^{r\times 1}$ be a vector of
observable factors such as tradable market risk factors. The form of the
factor pricing model is given by

$\displaystyle y_{it}=\alpha_{i}+\beta^{\top}_{i}X_{t}+u_{it},\quad
i=1,2,\dots,N,\quad t=1,2,\dots,T,$ (1)

where $\beta_{i}\in\mathbb{R}^{r\times 1}$ is a vector of factor loadings and
$u_{it}$ represents the idiosyncratic error which is assumed to be independent
of the factors $X_{t}$ with $\mathbb{E}[u_{it}]=0$ and
$\text{cov}(u_{it},u_{jt})=\sigma_{ij}$ for $i,j=1,\ldots,N$. The multi-factor
pricing theory implies that the vector of intercepts
$\bm{\alpha}=(\alpha_{1},\dots,\alpha_{N})^{\top}$ should be zero. Therefore,
to empirically validate the factor pricing theory, researchers often consider
the problem of testing

$\displaystyle H_{0}:\bm{\alpha}=\mathbf{0}\quad\text{versus}\quad
H_{a}:\bm{\alpha}\neq\mathbf{0}.$

An example of a commonly used test for this problem is the exact multivariate
$F$-test proposed by Gibbons et al. (1989), which assumes Gaussian errors and
$N<T$. However, its application has been confined to a relatively small number
of portfolios using monthly returns observed over a relatively long period.
With the advent of modern financial markets, where thousands of securities are
traded daily, there is a need for tests that can handle a large or diverging
number of assets $N$. Consequently, recent efforts have focused on developing
tests that allow $N$ to be larger than the time dimension $T$, given the
popularity of high-dimensional factor models. Some notable attempts along this
direction include Pesaran and Yamagata (2012), Gungor and Luger (2013),
Pesaran and Yamagata (2017), Ma et al. (2020), and Giglio et al. (2021).

The above methods are all sum-of-squares types, generally having good power
performance against dense alternatives. For instance, Pesaran and Yamagata
(2012, 2017) proposed Wald-type tests for validating the capital asset pricing
model based on the diagonal elements of the sample covariance matrix. However,
according to Fan et al. (2015), inefficient market pricing is more likely to
occur in exceptional assets rather than systematic mispricing of the entire
market. Therefore, it is desirable to develop a test that has good power
against sparse alternatives. The Wald-type test can suffer from low power when
$\bm{\alpha}$ is a sparse vector (i.e., sparse alternatives). To address this
issue, Fan et al. (2015) proposed adding a power enhancement component to the
Wald test (referred to as the FLY test) to improve its power under sparse
alternatives. The power enhancement component approaches zero under the null
hypothesis (so it does not affect the null distribution of the test
asymptotically) but diverges and dominates the Wald statistic under some
specific regions of sparse alternatives. Another effective way to deal with
sparse alternatives is the maximum type test proposed by Gungor and Luger
(2016), although they did not provide an in-depth analysis of the test.
Recently, Feng et al. (2022) re-examined the maximum type statistic and proved
its asymptotic null distribution to be the extreme type-I value distribution.

In general, the maximum type statistic only considers the largest entry of
$\bm{\alpha}$ and is not adaptive to the unknown sparsity level. This paper
proposes a novel testing procedure that overcomes the limitations of existing
methods and achieves considerable power improvement. The new test is designed
to

  1. 1.

be adaptive to the unknown sparsity level;

  2. 2.

borrow the cross-sectional dependence information to boost the statistical
power.

These features enable the test to outperform existing solutions over a wide
range of settings. To motivate our approach, we derive a likelihood ratio test
(LRT) under the assumption of joint normality for the errors against a sparse
alternative that assumes precisely $k$ nonzero components in $\bm{\alpha}$.
The resulting test statistic captures both the sparsity level of $\bm{\alpha}$
and the cross-sectional dependence via the error precision matrix.
Unfortunately, computing the test statistic involves solving a set
optimization problem with complexity $O(N^{k})$, which is computationally
infeasible for moderate $k$. To mitigate this issue, we propose replacing the
error precision matrix with a diagonal matrix in the test statistic, which
reduces computational complexity to $O(N\log N)$ for any fixed $k$.

An important consideration in this framework is how to select the appropriate
number of non-zero components $k$ that are adaptive to the signals. We propose
choosing $k$ to maximize the standardized test statistic, which draws
inspiration from the literature on adaptive testing (Ingster (1993); Fan
(1996); Donoho and Jin (2004)). To implement the testing procedure, we develop
a simulation approach to approximate the empirical null distributions of the
test statistics, based on recent high-dimensional Gaussian approximation
theory.

For a high-dimensional Central Limit Theorem (CLT) and bootstrap theory for
independent and identically distributed (i.i.d) sequences, we refer readers to
Chernozhukov et al. (2017), while for a CLT for high-dimensional time series,
we refer readers to Chang et al. (2023) (also see Zhang and Cheng (2018);
Zhang and Wu (2018)). In addition to the problem considered in this paper,
Gaussian approximation theory has been applied to other large-scale testing
problems, such as inference for high-dimensional mean (Chang et al. (2017)),
high-dimensional covariance matrices (Chang et al. (2017)), high-dimensional
precision matrices (Chang et al. (2018)), high-dimensional white noise testing
(Chang et al. (2017)), and inference for high-dimensional linear regression
(Zhang and Cheng (2017); Dezeure et al. (2017)).

We demonstrate theoretically that our simulation approach accurately
approximates the null distribution and that the proposed test is consistent
against a class of sparse alternatives. While our procedure is motivated by
the LR test derived under the joint normality assumption of the errors, the
validity of our approach does not depend on the Gaussian assumption. Our
empirical results show that the proposed procedure maintains reasonable size
and achieves considerable power improvement compared to existing tests under
weak and sparse alternatives. Furthermore, unlike other tests, the number of
factors in our proposed test can diverge as $N$ and $T$ increase.

The paper is organized as follows. Section 2 introduces the test statistic and
its adaptive version, and proposes a simulation-based procedure to approximate
the limiting null distributions. Section 3 investigates the asymptotic
properties of the tests. Sections 4 and 5 are dedicated to simulation studies
and real data analysis, respectively. Technical arguments are presented in the
Appendix.

Finally we introduce some notation that will be used throughout the paper. For
$a\in\mathbb{R}$, let $\lfloor a\rfloor$ be the integer part of $a$. For a
vector $z=(z_{1},\ldots,z_{N})\in\mathbb{R}^{N}$ and $q>0$, define
$|z|_{q}=(\sum_{i=1}^{N}|z_{i}|^{q})^{1/q}$ with $|z|=|z|_{2}$. Set
$|z|_{0}=\sum^{N}_{i=1}\mathbf{1}\\{z_{i}\neq 0\\}$ which is the number of
nonzero elements in $z$, and $|z|_{\infty}=\max_{1\leq i\leq N}|z_{i}|$. For
$V=(v_{ij})\in\mathbb{R}^{N\times N}$, define the matrix norms:
$\|V\|_{1}=\max_{1\leq j\leq N}\sum_{i=1}^{N}|v_{ij}|$,
$\|V\|_{2}=\max_{|z|=1}|Vz|$, $\|V\|_{\infty}=\max_{1\leq i\leq
N}\sum_{j=1}^{N}|v_{ij}|$ and $|V|_{\infty}=\max_{1\leq i,j\leq N}|v_{ij}|$.
Write $\text{diag}(V)=\text{diag}(v_{11},\ldots,v_{NN})$. Let $\|\cdot\|_{0}$
be the cardinality of a set. For $S\subseteq\\{1,\ldots,N\\}$, let $V_{S,S}$
be the submatrix of $V$ that contains the rows and columns in $S$. Similarly,
define $V_{S,-S}$ as the submatrix of $V$ with the rows in $S$ and the columns
in $\\{1,\ldots,N\\}\setminus S$.

##  2 Methodology

###  2.1 Test statistics

Let $Y_{t}=(y_{1t},\ldots,y_{Nt})^{\top}$ and
$\bm{\beta}=(\beta_{1},\ldots,\beta_{N})^{\top}$. Then the factor pricing
model (1) can be written in the matrix form as

$\displaystyle Y_{t}=\bm{\alpha}+\bm{\beta}X_{t}+U_{t},\quad t=1,2,\dots,T,$
(2)

where $U_{t}=(u_{1t},\ldots,u_{Nt})^{\top}$ is a martingale difference
sequence satisfying $\mathbb{E}(U_{t}|\mathfrak{F}_{t-1})=0$ with
$\mathfrak{F}_{s}$ denoting the $\sigma$-field generated by
$\\{U_{t}\\}_{t\leq s}$. Note that
$\\{\bm{\alpha}\in\mathbb{R}^{N}:\bm{\alpha}\neq\mathbf{0}\\}=\cup_{k=1}^{N}\mathcal{A}_{k},$
where
$\mathcal{A}_{k}=\\{\bm{\alpha}\in\mathbb{R}^{N}:|\bm{\alpha}|_{0}=k\\}$. To
motivate the subsequent derivations, we consider the testing problem:

$\displaystyle H_{0}:\bm{\alpha}=\mathbf{0}\quad\text{versus}\quad
H_{a,k}:\bm{\alpha}\in\mathcal{A}_{k}$

where the sparsity level $k$ under the alternative is assumed to be known as a
priori. In Section 2.3, we shall introduce a data-driven procedure to select
$k$.

Under the Gaussian assumption $U_{t}\sim\mathbb{N}(\mathbf{0},\Sigma)$, the
negative log-likelihood (up to a constant term) is equal to

$\mathcal{L}(\bm{\alpha};\bm{\beta};\Gamma)=\frac{1}{2}\sum_{t=1}^{T}(Y_{t}-\bm{\alpha}-\bm{\beta}X_{t})^{\top}\Gamma(Y_{t}-\bm{\alpha}-\bm{\beta}X_{t}),$

where $\Gamma=\Sigma^{-1}=(\gamma_{i,j})$ is the precision matrix. Define
$Z(\bm{\beta})=\bar{Y}-\bm{\beta}\bar{X}$,
$\mathbf{X}=(X_{1},\ldots,X_{T})\in\mathbb{R}^{r\times T}$,
$\mathbf{Y}=(Y_{1},\ldots,Y_{T})\in\mathbb{R}^{N\times T}$ and
$\mathbf{1}=(1,\ldots,1)^{\top}\in\mathbb{R}^{T\times 1}$, where
$\bar{Y}=T^{-1}\sum_{t=1}^{T}Y_{t}$ and $\bar{X}=T^{-1}\sum_{t=1}^{T}X_{t}.$
Fixing $\Gamma$ and using Lemma 2 in the appendix, we derive the likelihood
ratio test (LRT) for testing $H_{0}$ against $H_{a,k}$, which takes the form
of

$\displaystyle T\underset{\|S\|_{0}=k}{\max}\\{\Gamma
Z(\widehat{\bm{\beta}})\\}^{\top}_{S}\Gamma_{S,S}^{-1}\\{\Gamma
Z(\widehat{\bm{\beta}})\\}_{S}.$

Here $\widehat{\bm{\beta}}$ is the least squares estimator of $\bm{\beta}$
given by

$\widehat{\bm{\beta}}=\big{(}\mathbf{Y}-\bar{Y}\mathbf{1}^{\top}\big{)}\big{(}\mathbf{X}-\bar{X}\mathbf{1}^{\top}\big{)}^{\top}\left\\{\big{(}\mathbf{X}-\bar{X}\mathbf{1}^{\top}\big{)}\big{(}\mathbf{X}-\bar{X}\mathbf{1}^{\top}\big{)}^{\top}\right\\}^{-1}.$

Note that the MLE of $\bm{\alpha}$ can be expressed as

$\widehat{\bm{\alpha}}=(\widehat{\alpha}_{1},\ldots,\widehat{\alpha}_{N})^{\top}=\bar{Y}-\widehat{\bm{\beta}}\bar{X}=Z(\widehat{\bm{\beta}}),$

where
$\widehat{\alpha}_{i}=T^{-1}\sum_{t=1}^{T}y_{it}\\{1-(X_{t}-\bar{X})^{\top}\mathbf{W}\\}$
with
$\mathbf{W}=\\{T^{-1}\sum_{t=1}^{T}(X_{t}-\bar{X})(X_{t}-\bar{X})^{\top}\\}^{-1}\bar{X}\in\mathbb{R}^{r\times
1}$. Under the independence between $U_{t}$ and $X_{t}$, it can be shown that
the variance of $\widehat{\alpha}_{i}$ converges in probability to
$\sigma^{2}_{i}=T^{-1}(1+\mathbb{E}X_{t}^{\top}\Sigma_{X}^{-1}\mathbb{E}X_{t})\sigma_{ii},$
which can be estimated consistently by
$\widehat{\sigma}^{2}_{i}=T^{-2}\big{(}1+\bar{X}^{\top}\mathbf{W}\big{)}\sum_{t=1}^{T}\widehat{u}^{2}_{it},$
where
$\widehat{u}_{it}=y_{it}-\widehat{\alpha}_{i}-\widehat{\bm{\beta}}^{\top}_{i}X_{t}$.

Next, we introduce a procedure to estimate the unknown precision matrix.
Following Fan et al. (2015), we define a screening set:

$\hat{\textbf{s}}=\\{i:|\widehat{\alpha}_{i}|>\widehat{\sigma}_{i}\delta_{N,T},i=1,\ldots,N\\},$

where $\delta_{N,T}=C\log(\log T)\sqrt{\log N}$ with $C$ being some positive
constant. Under a minimal signal assumption, we can prove the sure screening
property (see Lemma 5 in the appendix) that $\hat{\textbf{s}}$ contains all
the nonzero components of $\bm{\alpha}$ with probability converging to one.
Set

$\widetilde{\bm{\alpha}}=\widehat{\bm{\alpha}}I(\hat{\textbf{s}})=\Big{(}\widehat{\alpha}_{1}I(|\widehat{\alpha}_{1}|>\widehat{\sigma}_{1}\delta_{N,T}),\ldots,\widehat{\alpha}_{N}I(|\widehat{\alpha}_{N}|>\widehat{\sigma}_{N}\delta_{N,T})\Big{)}^{\top}$

and let
$\widehat{U}=\mathbf{Y}-\widetilde{\bm{\alpha}}\mathbf{1}^{\top}-\widehat{\bm{\beta}}\mathbf{X}$
and
$\widehat{\Sigma}_{U}=\widehat{U}\widehat{U}^{\top}/T=(\widehat{\sigma}_{uij})^{N}_{i,j=1}$.
Further write $V=\textmd{diag}^{1/2}(\Sigma)$ and $R$ as the correlation
matrix of $\Sigma$. Accordingly, we define
$\widehat{V}=\textmd{diag}^{1/2}(\widehat{\Sigma}_{U})$ and $\widehat{R}$ as
the sample correlation matrix. Following Chen et al. (2013) and Shu and Nan
(2019), we use a slightly modified version of the graphical lasso estimator:

$\displaystyle\widehat{\textmd{K}}_{\rho}=\arg\min_{\Psi>0}\\{\textmd{tr}(\Psi\widehat{R})-\log\det(\Psi)+\rho|\Psi^{-}|_{1}\\},$
(3)

where the minimization is over a positive-definite matrix $\Psi=(\psi_{i,j})$,
$\rho$ is a non-negative tuning parameter, and $|\Psi^{-}|_{1}=\sum_{1\leq
i\neq j\leq N}|\psi_{i,j}|$. Then we estimate $\Gamma$ by
$\widehat{\Gamma}=\widehat{V}^{-1}\widehat{\textmd{K}}_{\rho}\widehat{V}^{-1}$.

######  Remark 1.

To estimate the sparse precision matrix, we can also use the nodewise Lasso
[Meinshausen and Bühlmann (2006)] or the Constrained $l_{1}$-Minimization for
Inverse Matrix Estimation (CLIME) [Cai et al. (2011)].

Combining the above derivations, we obtain the test statistic in the form of

$\displaystyle
G_{T}(k)=T\underset{\|S\|_{0}=k}{\max}\\{\widehat{\Gamma}Z(\widehat{\bm{\beta}})\\}^{\top}_{S}\widehat{\Gamma}_{S,S}^{-1}\\{\widehat{\Gamma}Z(\widehat{\bm{\beta}})\\}_{S}.$
(4)

However, computing the value of $G_{T}(k)$ is challenging as it involves a set
optimization problem with computational complexity of order $O(N^{k})$. This
complexity makes computation prohibitive even for moderate values of $k$. In
the following, we propose an approach to address this issue.

###  2.2 Modified tests

The computational cost of (4) is high, especially when combined with the
simulation-based approach in Section 2.4. To reduce the computational burden,
we propose a modified test by replacing $\widehat{\Gamma}_{S,S}$ with
$\text{diag}(\widehat{\Gamma}_{S,S})$, which includes only the diagonal
elements of $\widehat{\Gamma}$. This eliminates the need for best subset
selection, significantly improving computational efficiency. The resulting
test statistic is given by:

$\displaystyle\widetilde{G}_{T}(k)=T\underset{\|S\|_{0}=k}{\max}\\{\widehat{\Gamma}Z(\widehat{\bm{\beta}})\\}^{\top}_{S}\text{diag}^{-1}(\widehat{\Gamma}_{S,S})\\{\widehat{\Gamma}Z(\widehat{\bm{\beta}})\\}_{S}.$
(5)

Let $\widehat{\Gamma}Z(\widehat{\bm{\beta}})=(\widehat{z}_{i})^{N}_{i=1}$ and
$\widehat{\Gamma}=(\widehat{\gamma}_{i,j})_{i,j=1}^{N}$. Sort the values
$|\widehat{z}_{j}|^{2}/\widehat{\gamma}_{j,j}$ in descending order, say
$|\widehat{z}_{j_{1}^{*}}|^{2}/\widehat{\gamma}_{j_{1}^{*},j_{1}^{*}}\geq|\widehat{z}_{j_{2}^{*}}|^{2}/\widehat{\gamma}_{j_{2}^{*},j_{2}^{*}}\geq\cdots\geq|\widehat{z}_{j_{N}^{*}}|^{2}/\widehat{\gamma}_{j_{N}^{*},j_{N}^{*}}$.
Direct calculation yields that

$\displaystyle\widetilde{G}_{T}(k)=T\sum_{i=1}^{k}\frac{\widehat{z}_{j_{i}^{*}}^{2}}{\widehat{\gamma}_{j_{i}^{*},j_{i}^{*}}}.$

Therefore, the maximization in (5) can be solved using sorting algorithms,
which leads to the time complexity of the order $O(N\log N)$ for any fixed
$k$.

###  2.3 Selection of $k$ and adaptive testing

In this section, we propose a data-dependent approach to select $k$, inspired
by the literature on adaptive testing. For clarity, we focus on
$\widetilde{G}_{T}(k)$ but emphasize that the same idea applies equally to
$G_{T}(k)$. Let $\widehat{\mathbb{E}}[\widetilde{G}_{T}(k)]$ and
$\widehat{\text{var}}[\widetilde{G}_{T}(k)]$ be the estimates of the mean and
variance of $\widetilde{G}_{T}(k)$ under the null, which can be constructed
using the simulation-based approach in Section 2.4 below. Our idea is to find
$k$ which maximizes the standardized statistic, i.e.,

$\displaystyle k_{0}=\arg\max_{1\leq k\leq
K}\frac{\widetilde{G}_{T}(k)-\widehat{\mathbb{E}}[\widetilde{G}_{T}(k)]}{\sqrt{\widehat{\text{var}}[\widetilde{G}_{T}(k)]}},$

where $K$ is a user-specified upper bound. The corresponding adaptive test
statistic is given by

$\displaystyle\widetilde{\mathcal{G}}_{T}(K)=\frac{\widetilde{G}_{T}(k_{0})-\widehat{\mathbb{E}}[\widetilde{G}_{T}(k_{0})]}{\sqrt{\widehat{\text{var}}[\widetilde{G}_{T}(k_{0})]}}=\max_{1\leq
k\leq
K}\frac{\widetilde{G}_{T}(k)-\widehat{\mathbb{E}}[\widetilde{G}_{T}(k)]}{\sqrt{\widehat{\text{var}}[\widetilde{G}_{T}(k)]}}.$

###  2.4 A simulation-based approach

To approximate the sampling distributions of $G_{T}(k)$ and
$\widetilde{G}_{T}(k)$ for $k\geq 1$ and $\widetilde{\mathcal{G}}_{T}(K)$
under sparsity assumption, we propose a simulation-based approach, which is
closely related to the multiplier bootstrap in Chernozhukov et al. (2017),
Belloni and Oliveira (2018) and Chang et al. (2023). A key step in deriving
the simulation-based procedure is to understand the estimation effect of
$\widehat{\bm{\beta}}$ in $Z(\widehat{\bm{\beta}})$. We note that

$\displaystyle Z(\widehat{\bm{\beta}})=$
$\displaystyle\,\bm{\alpha}-T^{-1}\big{(}\mathbf{U}-\bar{U}\mathbf{1}^{\top}\big{)}\big{(}\mathbf{X}-\bar{X}\mathbf{1}^{\top}\big{)}^{\top}\mathbf{W}+\bar{U}$
$\displaystyle\approx$
$\displaystyle\,\bm{\alpha}+T^{-1}\sum^{T}_{t=1}U_{t}[1-(X_{t}-\bar{X})^{\top}\Sigma_{X}^{-1}\mathbb{E}X_{t}]:=\bm{\alpha}+\widetilde{Z},$

where $\mathbf{U}=(U_{1},\ldots,U_{T})$. Under the independence between
$\mathbf{U}$ and $\mathbf{X}$, it can be shown that the covariance matrix of
$\sqrt{T}\widetilde{Z}$ converges in probability to $(1+\zeta)\Sigma$ for
$\zeta=\mathbb{E}X_{t}^{\top}\Sigma_{X}^{-1}\mathbb{E}X_{t}$. The core of the
simulation approach is to construct a quantity using simulated normal
variables that approximates the distribution of $\widetilde{Z}$. The details
of our procedure can be described as follow:

  * 1.

Generate a sequence of independent standard normal random variables
$\\{\varepsilon_{t}\\}^{T}_{t=1}$ that are independent of the sample, and let

$\widetilde{Z}^{*}=T^{-1}\sqrt{1+\bar{X}^{\top}\widehat{\Sigma}_{X}^{-1}\bar{X}}\sum_{t=1}^{T}\left\\{\big{(}Y_{t}-\bar{Y}\big{)}-\widehat{\bm{\beta}}\big{(}X_{t}-\bar{X}\big{)}\right\\}\varepsilon_{t},$

where
$\widehat{\Sigma}_{X}=T^{-1}\sum_{t=1}^{T}(X_{t}-\bar{X})(X_{t}-\bar{X})^{\top}$.

  * 2.

Compute the simulation-based statistics for $\widetilde{G}_{T}(k)$ as

$\displaystyle\widetilde{G}^{*}_{T}(k)=T\underset{\|S\|_{0}=k}{\max}\\{\widehat{\Gamma}\widetilde{Z}^{*}\\}^{\top}_{S}\text{diag}^{-1}(\widehat{\Gamma}_{S,S})\\{\widehat{\Gamma}\widetilde{Z}^{*}\\}_{S}.$

  * 3.

Repeat Steps 1-2 $B$ times to get the $1-\alpha$ quantiles of
$\\{\widetilde{G}^{*}_{T,j}(k)\\}^{B}_{j=1}$, which are the simulation-based
critical values of $\widetilde{G}_{T}(k)$.

  * $3^{\prime}$.

We estimate the mean and variance of $\widetilde{G}_{T}(k)$ by

$\displaystyle\widehat{\mathbb{E}}[\widetilde{G}_{T}(k)]=B^{-1}\sum^{B}_{j=1}\widetilde{G}_{T,j}^{*}(k),$
$\displaystyle\widehat{\text{var}}[\widetilde{G}_{T}(k)]=B^{-1}\sum^{B}_{j=1}\left\\{\widetilde{G}_{T,j}^{*}(k)-B^{-1}\sum^{B}_{j=1}\widetilde{G}_{T,j}^{*}(k)\right\\}^{2},$

respectively. To approximate the distribution of the adaptive test, we
consider

$\displaystyle\widetilde{\mathcal{G}}_{T,j}^{*}(K)=\max_{1\leq k\leq
K}\frac{\widetilde{G}_{T,j}^{*}(k)-\widehat{\mathbb{E}}[\widetilde{G}_{T}(k)]}{\sqrt{\widehat{\text{var}}[\widetilde{G}_{T}(k)]}}.$

The $1-\alpha$ quantiles of
$\\{\widetilde{\mathcal{G}}^{*}_{T,j}(K)\\}^{B}_{j=1}$ will be the critical
value of our adaptive test.

##  3 Theory

###  3.1 Consistency under the null

In this subsection, we study the theoretical properties of the proposed tests
and justify the validity of the simulation-based approach. Define
$\alpha(T)=\sup_{A\in\mathcal{F}^{0}_{-\infty},B\in\mathcal{F}_{T}^{\infty}}|P(A)P(B)-P(AB)|,$
where $\mathcal{F}^{0}_{-\infty}$ and $\mathcal{F}_{T}^{\infty}$ denote the
$\sigma$-algebras generated by $\\{(X_{t},U_{t}):t\leq 0\\}$ and
$\\{(X_{t},U_{t}):T\leq t\\}$ respectively. Let $m=\underset{1\leq i\leq
N}{\max}\|\\{\gamma_{ij},j\neq i,1\leq j\leq N\\}\|_{0}$. For a random
variable $\xi$, we define its Orlicz norm to be
$\|\xi\|_{\psi_{\gamma_{1}}}=\inf\\{q\geq
0:\mathbb{E}[\psi_{\gamma_{1}}(|\xi|/q)]\leq 1\\}$ with
$\psi_{\gamma_{1}}(x):=\exp(x^{\gamma_{1}})-1$. To facilitate the derivations,
we make the following assumptions. Below we let $c_{i}$ with $0\leq i\leq 6$
be some positive constants independent of $N$ and $T$.

######  Assumption 1.

  * (a)

Suppose $\\{U_{t}\\}^{T}_{t=1}$ is a sequence of martingale difference
sequence with $\mathbb{E}(U_{t})=0$ and $\text{cov}(U_{t})=\Sigma$, where
$\lambda_{\min}(\Sigma)\geq c_{1}>0$. Assume that there exists a sequence of
constants $B_{T}>1$ and a universal constant $\gamma_{1}\geq 1$ such that
$\|u_{it}\|_{\psi_{\gamma_{1}}}\leq B_{T}$ for all $i$ and $t$. Moreover,
assume that $\sup_{1\leq i,j\leq N}\text{var}(u_{it}u_{jt})<\infty$,
$\|\Sigma\|_{1}<\infty$ and $\|\Gamma\|_{1}<\infty$.

  * (b)

Suppose $\\{X_{t}\\}^{T}_{t=1}$ is strictly stationary and is independent of
$\\{U_{t}\\}^{T}_{t=1}$. There exist $\omega,c_{2}>0$ such that

$P(|x_{it}|>s)\leq\exp\left\\{-(s/w)^{c_{2}}\right\\}.$

  * (c)

Suppose $\\{(X_{t},U_{t})\\}^{T}_{t=1}$ is $\alpha$-mixing with the mixing
coefficient $\alpha(T)$. There exist $c_{3}>0$ with
$c^{-1}_{2}+c^{-1}_{3}>0.5$ and $c_{4}>0$ such that for all
$T\in\mathbb{Z}^{+}$,

$\alpha(T)\leq\exp(-c_{4}T^{c_{3}}).$

  * (d)

Suppose $\lambda_{\min}(\Sigma_{X})\geq c_{5}>0$ and $\underset{1\leq i\leq
N}{\max}\|\beta_{i}\|<c_{6}$ for some $c_{5},c_{6}>0$.

######  Remark 2.

In contrast to the i.i.d assumption on $\\{U_{t}\\}^{T}_{t=1}$ in Assumption
4.1 (i) of Fan et al. (2015) and (A1) of Feng et al. (2022), we relax this
assumption in our paper because we utilize the high-dimensional CLT from Chang
et al. (2023) [see Theorem 7 therein]. Assumption 1(a) implies that
$\mathbb{E}\\{\exp(|u_{it}|^{\gamma_{1}}B_{T}^{-\gamma_{1}})\leq 2$, which is
equivalent to the condition that the tail of $u_{it}$ satisfy
$P(|u_{it}|>u)\leq 2\exp\left\\{-u^{\gamma_{1}}B_{T}^{-\gamma_{1}}\right\\}.$
Assumptions 1 (b) (d) are identical to those in Fan et al. (2015).
Specifically, conditions (b) and (c) require the tail decay and strong mixing
assumptions, respectively, which allow us to use Bernstein-type inequalities
for weakly dependent sequences in our technical proofs.

Denote
$\mathcal{C}_{N,T,k}=\max\big{\\{}B^{2}_{T}k^{\frac{2+6c_{3}}{c_{3}}}(\log
N)^{\frac{1+2c_{3}}{c_{3}}},B^{3}_{T}k^{10}(\log
N)^{\frac{7}{2}},k^{\frac{6-2c_{3}}{c_{3}}}(\log
N)^{\frac{3-c_{3}}{c_{3}}}\big{\\}}$.

######  Theorem 1.

Assume $\mathcal{C}_{N,T,k}=o(T^{1/3})$ and $k^{2}(r^{3}+r^{2}m)(\log
NT)^{5/2}/\sqrt{T}=o(1)$. Then under Assumption 1 (a)-(d) and the null
hypothesis $H_{0}$, we have

$\displaystyle\underset{\upsilon\geq
0}{\sup}\bigg{|}P\Big{(}\widetilde{G}^{*}_{T}(k)\leq\upsilon\Big{|}X_{1}^{T},Y_{1}^{T}\Big{)}-P\Big{(}\widetilde{G}_{T}(k)\leq\upsilon\Big{)}\bigg{|}=o_{p}(1),$
(6)

where $X_{1}^{T}=\\{X_{1},\ldots,X_{T}\\}$ and
$Y_{1}^{T}=\\{Y_{1},\ldots,Y_{T}\\}$.

######  Theorem 2.

Assume $\mathcal{C}_{N,T,k}=o(T^{1/3})$ and $K^{4}(r^{3}+r^{2}m)(\log
NT)^{5/2}/\sqrt{T}=o(1)$. Under Assumption 1 (a)-(d) and the null hypothesis
$H_{0}$, we have

$\displaystyle\underset{\upsilon\geq
0}{\sup}\bigg{|}P\Big{(}\widetilde{\mathcal{G}}^{*}_{T}(K)\leq\upsilon\Big{|}X_{1}^{T},Y_{1}^{T}\Big{)}-P\Big{(}\widetilde{\mathcal{G}}_{T}(K)\leq\upsilon\Big{)}\bigg{|}=o_{p}(1),$
(7)

where $X_{1}^{T}=\\{X_{1},\ldots,X_{T}\\}$ and
$Y_{1}^{T}=\\{Y_{1},\ldots,Y_{T}\\}$.

######  Remark 3.

If the $U_{t}$’s are assumed to be i.i.d, then the condition
$\mathcal{C}_{N,T,k}=o(T^{1/3})$ in Theorems 1 and 2 can be relaxed to
${k(\log NT/k)}^{7/3}=o(T^{1/3})$. However, if we allow the $U_{t}$’s to be
alpha-mixing in Theorems 1 and 2, a more restricted condition on the order of
$k$ relative to $T$ is required.

It is important to note that our theoretical results allow for the number of
factors $r$ to grow slowly with $N$ and $T$. However, the restrictions on $k$
and $K$ in Theorems 1-2 cannot be easily relaxed due to the use of the
Gaussian approximation theory in Chang et al. (2023). Nevertheless, our
numerical results demonstrate that the simulation-based approach performs
reasonably well even for large values of $K$.

###  3.2 Power analysis

Below we study the power properties of the proposed testing procedures. To
proceed, we impose the following conditions.

######  Assumption 2.

  * (a)

$\max_{j_{1}\leq\dots\leq
j_{k}}\sum_{i=1}^{k}\gamma_{j_{i},j_{i}}\alpha^{2}_{j_{i}}\geq(2k+\delta)\log(N)$
for some $\delta>0$.

  * (b)

$|\bm{\alpha}|_{0}=\lceil N^{\kappa}\rceil$ for some $0\leq\kappa<1/4$ and the
non-zero locations are randomly and uniformly drawn from $\\{1,2,\ldots,N\\}$.

  * (c)

Let
$\text{diag}^{-1/2}(\Gamma)\Gamma\text{diag}^{-1/2}(\Gamma)=(\phi_{i,j})^{N}_{i,j=1}$.
Assume that $\max_{1\leq i<j\leq N}|\phi_{i,j}|\leq c_{0}<1$.

Define
$\widetilde{C}^{*}_{\alpha}=\inf\\{\upsilon>0:P(\widetilde{G}^{*}_{T}(k)\leq\upsilon|X_{1}^{T},Y_{1}^{T})\geq
1-\alpha\\}$ as the simulation-based critical value. Let
$\widetilde{\mathcal{C}}_{\alpha}^{*}$ be the analogous quantity defined based
on $\widetilde{\mathcal{G}}^{*}_{T}(K)$. The following theorem establishes the
consistency of the testing procedures.

######  Theorem 3.

Assume that $\mathcal{C}_{N,T,k}=o(T^{1/3})$ and
$k^{2}(r^{3}+r^{2}m)\big{(}\log NT\big{)}^{5/2}/\sqrt{T}=o(1)$. Then under
Assumptions 1 and 2 (a) (b) (c), we have

$\displaystyle P(\widetilde{G}_{T}(k)>\widetilde{C}^{*}_{\alpha})\rightarrow
1.$ (8)

Moreover, suppose Assumption (a) holds by replacing $k$ with $K$ and
$K^{4}(r^{2}+rm)\big{(}\log NT\big{)}^{5/2}/\sqrt{T}=o(1)$. Then we have

$\displaystyle
P(\widetilde{\mathcal{G}}_{T}(K)>\widetilde{\mathcal{C}}^{*}_{\alpha})\rightarrow
1.$ (9)

##  4 Numerical studies

###  4.1 Competing methods

We carry out simulation experiments to evaluate the finite sample performance
of our tests. As a comparison, we also consider the FLY test of Fan et al.
(2015), the PY test of Pesaran and Yamagata (2017), and the testing procedures
recently proposed by Feng et al. (2022). Specifically, we consider the
adjusted Wald-statistic defined as

$\mathcal{T}_{\text{FLY}}=\frac{T(1-\bar{X}^{\top}\mathbf{\widetilde{W}})\widehat{\bm{\alpha}}^{\top}{\widehat{\Sigma^{\mathrm{T}}}}^{-1}\widehat{\bm{\alpha}}-\frac{N(T-r-1)}{T-r-3}}{\frac{T-r-1}{T-r-3}\sqrt{\frac{2N(T-r-2)}{T-r-5}\big{(}1+(N-1)\tilde{\rho}^{2}_{N,T}\big{)}}},$

where
$\mathbf{\widetilde{W}}=\big{(}\frac{1}{T}\sum_{t=1}^{T}X_{t}X^{\top}_{t}\big{)}^{-1}\bar{X}\in\mathbb{R}^{r\times
1}$, and $\widehat{\Sigma^{\mathrm{T}}}$ is a thresholding covariance
estimator [see e.g., Antoniadis and Fan (2001), Rothman et al. (2009)].

Pesaran and Yamagata (2012, 2017) proposed the modified Wald-type test:

$\displaystyle\mathcal{T}_{\text{PY}}=\frac{N^{-1/2}\sum_{i=1}^{N}\Big{(}t^{2}_{i}-\frac{T-r-1}{T-r-3}\Big{)}}{\frac{T-r-1}{T-r-3}\sqrt{\frac{2(T-r-2)}{T-r-5}\big{(}1+(N-1)\tilde{\rho}^{2}_{N,T}\big{)}}},$

where

$\displaystyle
t^{2}_{i}=\frac{\widehat{\alpha}^{2}_{i}\big{(}\mathbf{1}^{\top}\mathbf{M}_{\mathbf{X}}\mathbf{1}\big{)}}{\widehat{U}^{\top}_{\cdot}\widehat{U}_{\cdot}/(T-r-1)},\quad\widehat{V}_{U}=\sum_{t=1}^{T}\widehat{U}_{t}\widehat{U}_{t}^{\top}/(T-r-1),$
$\displaystyle\tilde{\rho}^{2}_{N,T}=\frac{2}{N(N-1)}\sum_{i=2}^{N}\sum_{j=1}^{i-1}\hat{\rho}^{2}_{ij}I\Big{(}(T-r-1)\hat{\rho}^{2}_{ij}\geq\varrho_{N,\alpha}\Big{)},$

with $\mathbf{Y}_{i\cdot}=(y_{i1},\ldots,y_{iT})^{\top}$,
$\widehat{U}_{\cdot}=(\widehat{u}_{i1},\ldots,\widehat{u}_{iT})^{\top}=\mathbf{M}_{\mathbf{X}}\big{(}\mathbf{Y}_{i\cdot}-\widehat{\alpha}_{i}\mathbf{1}\big{)}$,
$\hat{\rho}_{ij}=\frac{\widehat{U}^{\top}_{i}\widehat{U}_{j}}{\sqrt{\widehat{U}^{\top}_{i}\widehat{U}_{i}}\sqrt{\widehat{U}^{\top}_{j}\widehat{U}_{j}}}$,
$\widehat{U}_{t}=Y_{t}-\widehat{\bm{\alpha}}^{\top}-\widehat{\bm{\beta}}X_{t}$
and $\sqrt{\varrho_{N,\alpha}}=\Phi^{-1}(1-p_{N}/2)$. Here $\Phi$ is the
standard normal distribution function and $p_{N}=0.1/(N-1)$.

Feng et al. (2022) found that the maximum type statistic [Gungor and Luger
(2016)] defined as

$\displaystyle\mathcal{T}_{\text{MAX}}=\max_{1\leq i\leq N}t^{2}_{i}.$

follows an extreme value distribution. Moreover, Feng et al. (2022) combined
the $\mathcal{T}_{\text{PY}}$ and $\mathcal{T}_{\text{MAX}}$ tests to obtain a
combination test $\mathcal{T}_{\text{COM}}$ that is capable of detecting both
sparse and dense alternatives. See Feng et al. (2022) for the details.

###  4.2 Simulation results

Following the simulation studies of Pesaran and Yamagata (2017) and Feng et
al. (2022), the experiments are designed to mimic the commonly used Fama-
French three-factor model, where the factors $X_{t}$ have strong serial
correlation and heterogeneous variance. In particular, we consider the model

$\displaystyle y_{it}=\alpha_{i}+\beta^{\top}_{i}X_{t}+u_{it},\quad
i=1,2,\dots,N,\quad t=1,2,\dots,T.$

Here $\beta_{i}=(\beta_{i1},\beta_{i2},\beta_{i3})^{\top}$ for $1\leq i\leq N$
and $X_{t}=(x_{1t},x_{2t},x_{3t})^{\top}$ for $1\leq t\leq T$, which
corresponds to the three factors in the Fama-French model, namely the market
factor, the small minus big (SMB) and high minus low (HML). Each factor
follows a GARCH(1, 1) process, and all the coefficients are the same as those
in Pesaran and Yamagata (2017). Specifically,

$\displaystyle\textmd{Market
factor:}~{}~{}x_{1t}=0.53+0.06x_{1,t-1}+h_{1t}^{1/2}\zeta_{1t},$
$\displaystyle\textmd{SMB
factor:}~{}~{}x_{2t}=0.19+0.19x_{2,t-1}+h_{2t}^{1/2}\zeta_{2t},$
$\displaystyle\textmd{HML
factor:}~{}~{}x_{3t}=0.19+0.05x_{3,t-1}+h_{3t}^{1/2}\zeta_{3t},$

where $\zeta_{jt}^{\prime}s$ are simulated from a standard normal
distribution, and the variance terms $h_{jt}$ are generated from

$\displaystyle\textmd{Market
factor:}~{}~{}h_{1t}=0.89+0.85h_{1,t-1}+0.11h_{1t-1}\zeta^{2}_{1t-1},$
$\displaystyle\textmd{SMB
factor:}~{}~{}h_{2t}=0.62+0.74h_{2,t-1}+0.19h_{2t-1}\zeta^{2}_{2t-1},$
$\displaystyle\textmd{HML
factor:}~{}~{}h_{3t}=0.80+0.76h_{3,t-1}+0.15h_{3t-1}\zeta^{2}_{3t-1}.$

Similar to Pesaran and Yamagata (2017) and Feng et al. (2022), the above
process is generated over the periods $t=-49,\ldots,0,1,\ldots,T$ with
$x_{j,-50}=0$ and $h_{j,-50}=1$ for any $j=1,2$ and $3$. We discard the first
50 observations and use the data from time 1 to time $T$ to evaluate the
performance of the tests.

Following Pesaran and Yamagata (2017) and Feng et al. (2022), to capture the
main features of the individual asset returns and their cross-correlations, we
set the idiosyncratic errors to be

$U_{t}=\Sigma^{1/2}_{U}\varepsilon_{t},$

and consider three types of distribution for $\varepsilon_{t}$ as follows

  * (1)

$\varepsilon_{t}\sim N(\textbf{0},\mathrm{I}_{N})$, $t=1,2,\ldots,T$;

  * (2)

$\varepsilon_{it}\overset{iid}{\sim}t(3)/\sqrt{3}$, $i=1,2,\ldots,N$,
$t=1,2,\ldots,T$;

  * (3)

$\varepsilon_{t}$ consists of $N$ independent ARCH processes, i.e., each
component process is of the form $\varepsilon_{i,t}=\sigma_{i,t}e_{i,t}$,
where the $e_{i,t}$ are independent random variables following $N(0,1)$ or
$t(3)/\sqrt{3}$, and
$\sigma_{i,t}=\gamma_{0}+\gamma_{1}\varepsilon^{2}_{i,t-1}$ with $\gamma_{0}$
and $\gamma_{1}$ generated from, respectively, the uniform distributions on
$(0.25,0.5)$ and $(0,0.5)$ independently for different component processes.

Meanwhile, the factor loadings $\beta_{i1}$, $\beta_{i2}$ and $\beta_{i3}$ are
generated independently from the uniform distributions on $(0.2,2),(-1,1.5)$
and $(-1.5,1.5)$, respectively.

To better understand the power performance of all the methods under different
sparsity levels, we first present a small numerical study to investigate the
power behavior when $k=1,5,10$, and $15$, and $T=100,N=100$. We focus on the
AR(1) covariance structure $\Sigma_{U}=\big{(}\sigma_{ij}\big{)}_{i,j=1}^{N}$
with $\sigma_{ij}=0.6^{|i-j|}$ and $\Gamma=\Sigma^{-1}_{U}$. The non-zero
entries of $\bm{\alpha}$ are equal to
$\alpha_{i}=\omega_{i}\sqrt{a\log(N)/T}$, where $\omega_{i}$ are i.i.d. random
variables with $P(\omega_{i}=\pm 1)=1/2$. Figure 1 shows the power curves as a
function of the signal strength parameter $a$. The results clearly indicate
that the proposed test outperforms the other tests in all cases considered. In
particular, the COM, PY, and FLY tests are more sensitive to dense signals,
while the MAX test performs better than these tests for sparse signals.
Notably, the power of $\widetilde{\mathcal{G}}_{T}(K)$ with $K>1$ is higher
than that of $\widetilde{\mathcal{G}}_{T}(1)$, suggesting that it is
beneficial to choose a larger value of $K$ in practice. It should be noted
that while the restrictions on $k$ and $K$ in our theoretical results cannot
be relaxed, the simulation-based approach performs reasonably well even when
$K$ takes large values, as suggested by our numerical results.

Next, we consider the two cases for the error covariance $\Sigma_{U}$:

  * (1)

In the first case, we assume $\Sigma_{U}=\Lambda^{1/2}R\Lambda^{1/2}$. Similar
to Pesaran and Yamagata (2017) and Feng et al. (2022), we let
$\Lambda=\textmd{diag}(\sigma_{11},\cdots,\sigma_{NN})$ with the diagonal
entries $\sigma_{ii}$ drawn rom the uniform distribution on (1, 2). We set the
correlation matrix of $U_{t}$ to be

$R=\mathrm{I}_{N}+\textbf{b}\textbf{b}^{\top}-\textmd{diag}(\textbf{b})^{2},$

where $\textbf{b}=(b_{1},\ldots,b_{N})^{\top}$. The first and the last
$\lfloor N^{\delta_{\gamma}}\rfloor$ with $\delta_{\gamma}(=1/4,1/2,3/5)$
elements of b are drawn from the uniform distribution on $(0.7,0.9)$ to
generate different degrees of error cross-sectional dependence, where
$\delta_{\gamma}=1/4$ corresponds to the case of weak correlation and
$\delta_{\gamma}=3/5$ represents strong correlation. The remaining elements
are set to 0.

  * (2)

In the second case, we modify the example in Section 5.3 of Pesaran and
Yamagata (2017). In particular,

$\Sigma_{U}=\mathfrak{L}\mathfrak{L}^{\top}+(\mathrm{I}_{N}-\rho_{\epsilon}\mathcal{W})^{-1}(\mathrm{I}_{N}-\rho_{\epsilon}\mathcal{W}^{\top})^{-1},$

with $\mathfrak{L}=(\mathfrak{l}_{1},\ldots,\mathfrak{l}_{\lfloor
N^{\delta_{\gamma}}\rfloor},0,\ldots,0)^{\top}$. Here $\mathfrak{l}_{i}$ is
generated independently from the uniform distribution on (0.7, 0.9) for
$i=1,\ldots,\lfloor N^{\delta_{\gamma}}\rfloor$. Let $\rho_{\epsilon}=0.5$,
and $\mathcal{W}=(w_{ij})_{N\times N}$. All elements in $\mathcal{W}$ are zero
except that $w_{i+1,i}=w_{j-1,j}=0.5$ for $i=1,\ldots,N-2$ and $j=3,\ldots,N$,
and $w_{12}=w_{N,N-1}=1$.

We set $\bm{\alpha}=0$ under the null hypothesis. To evaluate the power of our
tests and the other competing tests under the alternative hypothesis, we
consider the two following scenarios.

  * (1)

Scenario 1: $k=\lfloor N^{1/4}\rfloor$ and the non-zero entries of
$\bm{\alpha}$ are equal to $\alpha_{i}=\omega_{i}\sqrt{2a\log(N)/T}$, where
$\omega_{i}$ are i.i.d. random variables with the Uniform distribution on (0,
1) and $a=1,2,\ldots,5$;

  * (2)

Scenario 2: $k=\lfloor N^{1/3}\rfloor$ and the non-zero entries of
$\bm{\alpha}$ are equal to $\alpha_{i}=\omega_{i}\sqrt{2a\log(N)/T}$, where
$\omega_{i}$ are i.i.d. random variables with $P(\omega_{i}=\pm 1)=1/2$ and
$a=0.4,0.8,1.2,1.6,2.0$.

We present the empirical rejection rates for $\Sigma_{U}$ under the null in
Case (1) in Table 1. Overall, all tests have reasonable empirical size, with
the PY, MAX, and COM tests performing better than the other tests. When
$\varepsilon_{t}$ follows the normal distribution and $N=100,200$, the
proposed adaptive test slightly over-rejects, while the other tests provide
quite accurate control of the size. However, when $\delta_{\gamma}=1/2$ and
$\delta_{\gamma}=3/5$, the FLY test is quite conservative.

Figures 2-10 provide a summary of the rejection rates of different procedures
under the alternatives for $\Sigma_{U}$ in Case (1), with $k=\lfloor
N^{1/4}\rfloor$. The proposed test, as well as the MAX and COM tests,
outperform the PY and FLY tests in all cases, regardless of the values of
$\delta_{\gamma}$. In the case of weak signals (i.e., $N=100,200$), the
proposed test delivers more power than the MAX and COM tests. As discussed in
Feng et al. (2022), the PY and FLY tests are sum-of-squares type tests that
are sensitive to dense signals, and are likely to perform poorly in the case
of sparse signals. Not surprisingly, strong cross-sectional dependence
controlled by the parameter $\delta_{\gamma}$ harms the power of the PY and
FLY tests. As $\delta_{\gamma}$ increases, the power of the PY and FLY tests
decreases. Similar observations are found in the other cases; see Figures 1-9
for $\Sigma_{U}$ in Case(1) in the Supplementary material.

We evaluate the performance of different procedures under the alternative
hypothesis for $\Sigma_{U}$ in Case (2) when $k=\lfloor N^{1/3}\rfloor$. The
results are summarized in Figures 11-19. We observe that the proposed test
consistently outperforms the other tests, indicating its robustness to strong
cross-sectional dependence. We also find that the power of the PY and FLY
tests decreases as $\delta_{\gamma}$ gets larger, which is a similar
observation to the case of $\Sigma_{U}$ in Case (1). This trend is further
highlighted in Figures 10-18 in Supplementary material. Additional simulation
results are provided in the Supplementary material, where we consider the case
where $\varepsilon_{t}$ follows an ARCH process. Table 1 and Figures 19-24
present the corresponding empirical rejection rates and power curves of the
different tests.

It is important to highlight that the power of the PY and FLY tests is
sensitive to the strength of the cross-sectional dependence, which is
controlled by the parameter $\delta_{\gamma}$. Additionally, we found that the
proposed test is relatively robust to the choice of $K$. In summary, our test
provides reasonable control of the type I error and delivers better power than
the other tests under weak signals.

###  4.3 Real data analysis

We examine the monthly excess returns on all the constituents of the S&P 500
from the CRSP database from January 1980 to December 2012, as analyzed by Fan
et al. (2015). One of the key empirical findings reported in Fan et al. (2015)
is that market inefficiency is predominantly caused by a small fraction of
stocks with positive intercepts, rather than a large proportion of slightly
mispriced assets. This finding provides empirical evidence of sparse
alternatives. Our goal here is to test market efficiency on a rolling window
basis. For each month from December 1987 to December 2012, we calculate test
statistics using the preceding 96 months’ returns ($T=96$). A total of 400
stocks were considered for this study, and we only include stocks without
missing observations in the past eight years in the data for each testing
month. We fit the Fama and French three-factor model to the data in each
rolling window ($l=12/1987,12/1988,\ldots,12/2012$):

$\displaystyle
r^{(l)}_{it}-f^{(l)}_{t}=\alpha^{(l)}_{i}+\beta^{(l)}_{i,\text{MKT}}(\text{MKT}^{(l)}_{t}-f^{(l)}_{t})+\beta^{(l)}_{i,\text{SMB}}\text{SMB}^{(l)}_{t}+\beta^{(l)}_{i,\text{HML}}\text{HML}^{(l)}_{t}+u_{it}^{(l)}$
(10)

for $i=1,\ldots,N$ and $t=l-95,\dots,l$, where $r_{it}$ is the return for
stock $i$ at month $t$, $f_{t}$ is the risk free rate, and MKT, SMB and HML
represent the market, size and value factors, respectively.

We apply the proposed approach and the other four completing procedures to
examine the market efficiency hypothesis for the data in each rolling window.
The screening set suggested by Fan et al. (2015) is computed, which shows that
only a few significant nonzero “alpha” components exist in seven rolling
windows, which directly implies the presence of sparse or weak signals. For
instance, with $l=12/1992$, we select $284$ stocks (Window 1) without missing
observations and fit the Fama-French three-factor model, and record the
resulting $p$-values from all the tests in Table 2. At the 5% significant
level, only the proposed test shows significant evidence against the null,
while the other tests fail to reject it, indicating that $\bm{\alpha}$ is
possibly sparse in this window. Similarly, for $l=12/1996$ (Window 2) based on
$279$ stocks, the proposed adaptive test provides the strongest evidence
against the null, followed by the MAX test. For $l=12/1999$ (Window 3) based
on $241$ stocks, both the PY and the proposed adaptive tests reject the null
hypothesis, while the other tests fail to reject it. Lastly, with $l=12/2011$
(Window 4) based on $170$ stocks, all the tests find significant evidence
against the null.

To better understand these results and conduct a fair comparison among
different procedures, we performed a simulation study. We generated data from
the fitted factor models for each of the rolling windows. Specifically, for
the $l$th rolling window, we generated data from either the null or
alternative hypothesis.

  * S1.

Let
$(\hat{\alpha}^{(l)}_{i},\hat{\beta}^{(l)}_{i,\text{MKT}},\hat{\beta}^{(l)}_{i,\text{SMB}},\hat{\beta}^{(l)}_{i,\text{HML}})$
be the estimates of the coefficients and $\hat{u}^{(l)}_{it}$ be the
corresponding residuals. Under the null, we generate the excess returns as
follows

$\displaystyle
r^{\ast}_{it}=\hat{\beta}^{(l)}_{i,\text{MKT}}(\text{MKT}^{(l)}_{t}-f^{(l)}_{t})+\hat{\beta}^{(l)}_{i,\text{SMB}}\text{SMB}^{(l)}_{t}+\hat{\beta}^{(l)}_{i,\text{HML}}\text{HML}^{(l)}_{t}+\epsilon_{t}\hat{u}^{(l)}_{it},$

for $i=1,2,\ldots,N$ and $t=l-95,\ldots,l$, where $\epsilon_{t}$ follows the
Rademacher distribution.

  * S2.

We let

$\tilde{\alpha}^{(l)}_{i}=\hat{\alpha}^{(l)}_{i}I\\{|\hat{\alpha}_{i}^{(l)}|>\hat{\sigma}_{i}\log(\log
T)\sqrt{\log N}\\}.$

Under the alternative, we generate the data according to the following model

$\displaystyle
r^{\ast}_{it}=\tilde{\alpha}^{(l)}_{i}+\hat{\beta}^{(l)}_{i,\text{MKT}}(\text{MKT}^{(l)}_{t}-f^{(l)}_{t})+\hat{\beta}^{(l)}_{i,\text{SMB}}\text{SMB}^{(l)}_{t}+\hat{\beta}^{(l)}_{i,\text{HML}}\text{HML}^{(l)}_{t}+\epsilon_{t}\hat{u}^{(l)}_{it},$

for $i=1,2,\ldots,N$ and $t=l-95,\ldots,l$, where $\epsilon_{t}$ follows the
Rademacher distribution.

We present the results of a simulation study in Table 2, which allows us to
compare the power and empirical size of the different tests. We conducted 1000
replications for each scenario and calculated the average power and empirical
size. All tests showed reasonable control over the empirical size, although
the FLY test was conservative. The proposed adaptive test was consistently the
most powerful in all scenarios, outperforming the MAX and COM tests by
approximately 10% in scenarios D3 and D4. These findings are consistent with
the results obtained from the real data analysis, indicating that the proposed
test may be a useful addition to the existing methods for testing the market
efficiency hypothesis.

SUPPLEMENTARY MATERIAL

Title:

The supplement contains some intermediate theoretical results and the proofs
of the theorems.

## References

  * Antoniadis and Fan (2001) Antoniadis, A. and J. Fan (2001).  Regularized wavelet approximations.  Journal of the American Statistical Association 96(455), 939–967.
  * Belloni and Oliveira (2018) Belloni and Oliveira (2018).  Approximate group context tree.  The Annals of Statistics 45(1), 355–385.
  * Cai et al. (2011) Cai, T., W. D. Liu, and X. Luo (2011).  A constrained l1 minimization approach to sparse precision matrix estimation.  Journal of the American Statistical Association 1064, 594–607.
  * Chang et al. (2023) Chang, J., X. Chen, and M. Wu (2023).  Central limit theorems for high dimensional dependent data.  Bernoulli.
  * Chang et al. (2017) Chang, J., Q. Yao, and W. Zhou (2017).  Testing for high-dimensional white noise using maximum cross-correlations.  Biometrika 104(1), 31–41.
  * Chang et al. (2018) Chang, J., Q. Yao, and T. Zou (2018).  Confidence regions for entries of a large precision matrix.  Journal of Econometrics 206(1), 57–82.
  * Chang et al. (2017) Chang, J., C. Zheng, W. Zhou, and W. Zhou (2017).  Simulation-based hypothesis testing of high dimensional means under covariance heterogeneity.  Biometrics 73(4), 1300–1310.
  * Chang et al. (2017) Chang, J., W. Zhou, W. Zhou, and L. Wang (2017).  Comparing large covariance matrices under weak conditions on the dependence structure and its application to gene clustering.  Biometrics 73(1), 31–41.
  * Chen et al. (2013) Chen, X., M. Xu, and W. Wu (2013).  Covariance and precision matrix estimation for high-dimensional time series.  The Annals of Statistics 41, 2994–3021.
  * Chernozhukov et al. (2017) Chernozhukov, V., D. Chetverikov, and K. Kato (2017).  Central limit theorems and bootstrap in high dimensions.  The Annals of Probability 45(4), 2309–2352.
  * Dezeure et al. (2017) Dezeure, R., P. Bühlmann, and C.-H. Zhang (2017).  High-dimensional simultaneous inference with the bootstrap.  Test 26, 685–719.
  * Donoho and Jin (2004) Donoho, D. and J. Jin (2004).  Higher criticism for detecting sparse heterogeneous mixtures.  The Annals of Statistics 32(3), 85–114.
  * Fama and French (1993) Fama, E. F. and K. R. French (1993).  Common risk factors in the returns on stocks and bonds.  Journal of Financial Economics 33(1), 3–56.
  * Fama and French (2015) Fama, E. F. and K. R. French (2015).  A five-factor asset pricing model.  Journal of Financial Economics 116(1), 1–22.
  * Fan (1996) Fan, J. (1996).  Test of significance based on wavelet thresholding and neyman’s truncation.  Journal of the American Statistical Association 91(434), 674–688.
  * Fan et al. (2015) Fan, J., Y. Liao, and J. Yao (2015).  Power enhancement in high dimensional cross-sectional tests.  Econometrica 83(4), 1497–1541.
  * Feng et al. (2022) Feng, L., W. Lan, B. Liu, and Y. Ma (2022).  High-dimensional test for alpha in linear factor pricing models with sparse alternatives.  Journal of Econometrics 229, 152–175.
  * Gibbons et al. (1989) Gibbons, M. R., S. A. Ross, and J. Shanken (1989).  A test of the efficiency of a given portfolio.  Econometrica 57(5), 1121–1152.
  * Giglio et al. (2021) Giglio, S., Y. Liao, and D. Xiu (2021).  Thousands of alpha tests.  The Review of Financial Studies 34(7), 3456–3496.
  * Gungor and Luger (2013) Gungor, S. and R. Luger (2013).  Testing linear factor pricing models with large cross sections: a distribution-free approach.  Journal of Business and Economic Statistics 31(1), 66–77.
  * Gungor and Luger (2016) Gungor, S. and R. Luger (2016).  Multivariate tests of mean-variance efficiency and spanning with a large number of assets and time-varying covariances.  Journal of Business and Economic Statistics 34(2), 161–175.
  * Ingster (1993) Ingster, Y. I. (1993).  Asymptotically minimax hypothesis testing for nonparametric alternatives. i, ii, iii.  Math. Methods Statist 2(2), 85–114.
  * Lintner (1965) Lintner, J. (1965).  The valuation of risky assets and the selection of risky investments in portfolio selection and capital budgets.  Journal of Finance 47(1), 13–37.
  * Ma et al. (2020) Ma, S., W. Lan, L. Su, and C. L. Tsai (2020).  Testing alphas in conditional time-varying factor models with high-dimensional assets.  Journal of Business and Economic Statistics 38(1), 1–34.
  * Meinshausen and Bühlmann (2006) Meinshausen, N. and P. Bühlmann (2006).  High-dimensional graphs and variable selection with the lasso.  The Annals of Statistics 34, 1436–1462.
  * Pesaran and Yamagata (2012) Pesaran, M. H. and T. Yamagata (2012).  Testing capm with a large number of assets.  SSRN Discussion Paper No. 6469.
  * Pesaran and Yamagata (2017) Pesaran, M. H. and T. Yamagata (2017).  Testing for alpha in linear factor pricing models with a large number of securities.  Working Paper.
  * Ross (1976) Ross, S. (1976).  The arbitrage theory of capital asset pricing.  Journal of Finance 13(3), 341–360.
  * Rothman et al. (2009) Rothman, A., E. Levina, and J. Zhu (2009).  Generalized thresholding of large covariance matrices.  Journal of the American Statistical Association 104(485), 177–186.
  * Sharpe (1964) Sharpe, W. (1964).  Capital asset prices: a theory of capital market equilibrium under conditions of risk.  Journal of Finance 19(3), 425–442.
  * Shu and Nan (2019) Shu, H. and B. Nan (2019).  Estimation of large covariance and precision matrices from temporally dependent observations.  The Annals of Statistics 47, 1321–1350.
  * Zhang and Wu (2018) Zhang, D. and W. B. Wu (2018).  Gaussian approximation for high-dimensional time series.  The Annals of Statistics 45, 1895–1919.
  * Zhang and Cheng (2017) Zhang, X. and G. Cheng (2017).  simultaneous inference for high-dimensional linear models.  Journal of the American Statistical Association 112, 757–768.
  * Zhang and Cheng (2018) Zhang, X. and G. Cheng (2018).  Gaussian approximation for high dimensional vector under physical dependence.  Bernoulli 24, 2640 – 2675.

Table 1: Size (%) of different tests for $\Sigma_{U}$ in Case (1) at the 5% significance level |  | $\delta_{\gamma}=1/4$ |  | $\delta_{\gamma}=1/2$ |  | $\delta_{\gamma}=3/5$
---|---|---|---|---|---|---
| $N=$ | $100$ | $200$ | $500$ |  | $100$ | $200$ | $500$ |  | $100$ | $200$ | $500$
Method | $T=100$ |  |  | Normal errors |  |  |  |  |
$\mathcal{T}_{\text{FLY}}$ |  | 5.1 | 5.9 | 5.3 |  | 2.8 | 2.1 | 1.7 |  | 0.6 | 0.4 | 0.3
$\mathcal{T}_{\text{PY}}$ |  | 5.5 | 6.3 | 5.6 |  | 5.4 | 4.4 | 5.4 |  | 6.0 | 5.9 | 5.0
$\mathcal{T}_{\text{MAX}}$ |  | 4.9 | 5.6 | 7.1 |  | 3.6 | 4.7 | 6.1 |  | 3.9 | 4.4 | 6.4
$\mathcal{T}_{\text{COM}}$ |  | 4.6 | 5.5 | 6.3 |  | 4.9 | 4.5 | 5.6 |  | 5.8 | 5.9 | 6.0
$\widetilde{\mathcal{G}}_{T}(5)$ |  | 7.2 | 7.5 | 6.5 |  | 7.2 | 7.2 | 5.9 |  | 6.4 | 6.6 | 6.2
$\widetilde{\mathcal{G}}_{T}(10)$ |  | 6.9 | 7.0 | 6.5 |  | 6.8 | 6.9 | 5.8 |  | 6.0 | 6.1 | 6.3
$\widetilde{\mathcal{G}}_{T}(15)$ |  | 6.7 | 6.7 | 6.4 |  | 6.5 | 6.8 | 5.8 |  | 5.8 | 6.0 | 6.2
$\widetilde{\mathcal{G}}_{T}(30)$ |  | 6.7 | 6.6 | 6.3 |  | 6.3 | 6.8 | 5.8 |  | 5.6 | 5.9 | 6.1
| $T=100$ |  |  | Student-t errors |  |  |  |  |
$\mathcal{T}_{\text{FLY}}$ |  | 4.0 | 3.6 | 2.9 |  | 1.9 | 1.0 | 0.8 |  | 0.2 | 0.2 | 0.3
$\mathcal{T}_{\text{PY}}$ |  | 5.2 | 3.7 | 4.1 |  | 5.4 | 5.0 | 3.9 |  | 5.1 | 6.8 | 6.2
$\mathcal{T}_{\text{MAX}}$ |  | 3.0 | 3.8 | 3.3 |  | 3.3 | 3.4 | 4.1 |  | 2.4 | 3.0 | 4.4
$\mathcal{T}_{\text{COM}}$ |  | 3.8 | 4.2 | 3.5 |  | 5.0 | 4.1 | 4.4 |  | 4.2 | 5.6 | 6.2
$\widetilde{\mathcal{G}}_{T}(5)$ |  | 5.3 | 5.0 | 3.1 |  | 5.2 | 5.1 | 3.8 |  | 4.1 | 3.7 | 3.8
$\widetilde{\mathcal{G}}_{T}(10)$ |  | 4.6 | 4.7 | 3.1 |  | 5.0 | 4.9 | 3.8 |  | 4.0 | 3.8 | 3.8
$\widetilde{\mathcal{G}}_{T}(15)$ |  | 4.4 | 4.5 | 3.1 |  | 4.8 | 4.9 | 3.8 |  | 4.0 | 3.8 | 3.8
$\widetilde{\mathcal{G}}_{T}(30)$ |  | 4.4 | 4.2 | 3.1 |  | 4.5 | 4.7 | 3.7 |  | 4.1 | 3.7 | 3.8

Figure 1: Power curves for different tests, where $k=1,5,10,15$ with
$T=100,N=100$, and $\Sigma_{U}=\big{(}\sigma_{ij}\big{)}_{i,j=1}^{N}$ with
$\sigma_{ij}=0.6^{|i-j|}$. “AT” represents the proposed test statistics.

Figure 2: Power (%) of different tests for $\Sigma_{U}$ in Case (1) with
$\delta_{\gamma}=1/4$ and $k=\lfloor N^{1/4}\rfloor$ Figure 3: Power (%) of
different tests for $\Sigma_{U}$ in Case (1) with $\delta_{\gamma}=1/2$ and
$k=\lfloor N^{1/4}\rfloor$ Figure 4: Power (%) of different tests for
$\Sigma_{U}$ in Case (1) with $\delta_{\gamma}=3/5$ and $k=\lfloor
N^{1/4}\rfloor$

Figure 5: Power (%) of different tests for $\Sigma_{U}$ in Case (1) with
$\delta_{\gamma}=1/4$ and $k=\lfloor N^{1/4}\rfloor$ Figure 6: Power (%) of
different tests for $\Sigma_{U}$ in Case (1) with $\delta_{\gamma}=1/2$ and
$k=\lfloor N^{1/4}\rfloor$ Figure 7: Power (%) of different tests for
$\Sigma_{U}$ in Case (1) with $\delta_{\gamma}=3/5$ and $k=\lfloor
N^{1/4}\rfloor$

Figure 8: Power (%) of different tests for $\Sigma_{U}$ in Case (1) with
$\delta_{\gamma}=1/4$ and $k=\lfloor N^{1/4}\rfloor$ Figure 9: Power (%) of
different tests for $\Sigma_{U}$ in Case (1) with $\delta_{\gamma}=1/2$ and
$k=\lfloor N^{1/4}\rfloor$ Figure 10: Power (%) of different tests for
$\Sigma_{U}$ in Case (1) with $\delta_{\gamma}=3/5$ and $k=\lfloor
N^{1/4}\rfloor$

Figure 11: Power (%) of different tests for $\Sigma_{U}$ in Case (2) with
$\delta_{\gamma}=1/4$ and $k=\lfloor N^{1/3}\rfloor$ Figure 12: Power (%) of
different tests for $\Sigma_{U}$ in Case (2) with $\delta_{\gamma}=1/2$ and
$k=\lfloor N^{1/3}\rfloor$ Figure 13: Power (%) of different tests for
$\Sigma_{U}$ in Case (2) with $\delta_{\gamma}=3/5$ and $k=\lfloor
N^{1/3}\rfloor$

Figure 14: Power (%) of different tests for $\Sigma_{U}$ in Case (2) with
$\delta_{\gamma}=1/4$ and $k=\lfloor N^{1/3}\rfloor$ Figure 15: Power (%) of
different tests for $\Sigma_{U}$ in Case (2) with $\delta_{\gamma}=1/2$ and
$k=\lfloor N^{1/3}\rfloor$ Figure 16: Power (%) of different tests for
$\Sigma_{U}$ in Case (2) with $\delta_{\gamma}=3/5$ and $k=\lfloor
N^{1/3}\rfloor$

Figure 17: Power (%) of different tests for $\Sigma_{U}$ in Case (2) with $\delta_{\gamma}=1/4$ and $k=\lfloor N^{1/3}\rfloor$ Figure 18: Power (%) of different tests for $\Sigma_{U}$ in Case (2) with $\delta_{\gamma}=1/2$ and $k=\lfloor N^{1/3}\rfloor$ Figure 19: Power (%) of different tests for $\Sigma_{U}$ in Case (2) with $\delta_{\gamma}=3/5$ and $k=\lfloor N^{1/3}\rfloor$ Table 2: Testing results of different procedures for the real data and simulated data mimicking the real data structure |  Data $\backslash$ Test | $\mathcal{T}_{\text{PY}}$ | $\mathcal{T}_{\text{MAX}}$ | $\mathcal{T}_{\text{COM}}$ | $\widetilde{\mathcal{G}}_{T}(5)$ | $\widetilde{\mathcal{G}}_{T}(10)$ | $\widetilde{\mathcal{G}}_{T}(30)$ | $\mathcal{T}_{\text{FLY}}$
---|---|---|---|---|---|---|---|---
| Window 1 |  |  |  |  $p$-value |  |  |
|  | 0.0942 | 0.0721 | 0.0832 | 0.0240 | 0.0240 | 0.0260 | 0.1267
$H_{0}$ | Bootstrap data |  |  |  | Size (%) |  |  |
|  | 8.6 | 4.5 | 7.8 | 5.8 | 5.9 | 5.7 | 2.6
$H_{a}$ | Bootstrap data |  |  |  | Power (%) |  |  |
|  | 18.2 | 69.6 | 65.0 | 74.4 | 74.5 | 74.4 | 9.2
| Window 2 |  |  |  |  $p$-value |  |  |
|  | 0.1189 | 0.0432 | 0.0861 | 0.0140 | 0.0140 | 0.0140 | 0.1368
$H_{0}$ | Bootstrap data |  |  |  | Size (%) |  |  |
|  | 6.4 | 5.3 | 6.4 | 5.6 | 5.7 | 5.6 | 0.16
$H_{a}$ | Bootstrap data |  |  |  | Power (%) |  |  |
|  | 12.3 | 59.8 | 53.4 | 63.3 | 63.2 | 63.1 | 4.4
| Window 3 |  |  |  |  $p$-value |  |  |
|  | 0.0364 | 0.0936 | 0.0650 | 0.0160 | 0.0160 | 0.0180 | 0.1005
$H_{0}$ | Bootstrap data |  |  |  | Size (%) |  |  |
|  | 8.7 | 5.3 | 8.2 | 5.5 | 5.5 | 5.4 | 2.0
$H_{a}$ | Bootstrap data |  |  |  | Power (%) |  |  |
|  | 14.2 | 55.3 | 49.2 | 60.8 | 60.8 | 60.6 | 3.2
| Window 4 |  |  |  |  $p$-value |  |  |
|  | 0.0035 | 0.0168 | 0.0102 | 0.0000 | 0.0000 | 0.0000 | 0.0482
$H_{0}$ | Bootstrap data |  |  |  | Size (%) |  |  |
|  | 8.1 | 4.2 | 6.8 | 4.4 | 4.4 | 4.4 | 2.6
$H_{a}$ | Bootstrap data |  |  |  | Power (%) |  |  |
|  | 21.4 | 73.3 | 68.1 | 79.8 | 79.6 | 79.6 | 10.2