# Limiting Spectral Distribution of Large-Dimensional Sample Covariance Matrices Generated by the Periodic Autoregressive Model

Authors: Jin Zou; Dong Han
Journal: Journal of Mathematics
 (2021)

---

## Abstract

The explicit representation for the limiting spectral moments of sample covariance matrices generated by the periodic autoregressive model (PAR) is established. We propose to use the moment-constrained maximum entropy method to estimate the spectral density function. The experiments show that the maximum entropy spectral density function curve obtained based on the fourth-order limiting spectral moment can match histograms of the eigenvalues of the covariance matrices very well.


---

## Body

## 1. Introduction

Periodic time-series models are used for modeling stationary time series with periodic characteristics. Periodicity refers to wave or oscillating movements around long-term trends presented in the time series, such as the water sector, where the supply of water varies with seasonal rainfall. In the economic time series, it is often accompanied by some cycle or quasi-periodic fluctuations, such as day and night changes of urban electricity consumption, seasonal changes of railway passenger transport, and repeated acceleration, interruption, stagnation, and recovery of economic growth; see Franses and Paap [1], Bell et al. [2], and Aliat and Hamdi [3]. The periodic autoregressive (PAR) model is a generalisation of the classical autoregressive (AR) model by allowing the parameters to vary with the variables over a period:(1)Xip+v=∑k=1qvαkvXip+v−k+εip+v.For a givenv and a predetermined value p, the random variable Xip+v denotes the value in the vth quarter of the ith cycle, where v∈1,…,p and i∈1,…,n. In the ith cycle, Xip+v is an AR process. The model order is given by qv, and αkv,k=1,…,qv, represents the AR model coefficients for the vth quarter. The error process ε=εip+v corresponds to a periodic white noise, with Eεip+v=0 and varεip+v=σ2v>0. Let Xi=Xip+1,…,Xip+p⊤, where the transpose (indicated by ⊤) of a row vector is a column vector, and we define the large-dimensional sample covariance matrices as Bn=1/n∑i=1nXiXi⊤, with the periodic length p tends to infinity in proportion to the sample size n, namely, p/n⟶c∈0,∞. Sample covariance matrices are very important in multivariate statistical inference since many test statistics are defined by their eigenvalues or functions.Random matrix theory (RMT) originated from the study of the energy levels of a large number of particles in quantum mechanics. Many of the laws of mathematical physics were discovered through numerical studies. In the late 1950s, E. P. Wigner formulated the problem in terms of the empirical distribution of random matrices [4, 5], which led to the study of the semicircular law of Gaussian matrices. Since then, RMT has formed an active branch in modern probability theory. The LSD of the sample covariance matrix of a large-dimensional random matrix composed of independent random variables, the MP law, was first established by Marcenko and Pastur [6]. Further research efforts were conducted to estimate the LSD of a product of two random matrices. Yin and Krishnaiah [7] proved the existence of LSD for the matrix sequence SnTn when the two matrix sequences Sn and Tn are standard Wishart matrices and positive definite matrices, respectively. In particular, Bai et al. [8] obtained a display representation of the LSD of the matrix sequence SnTn when the two matrix sequences Sn and Tn are the sample covariance matrix and the Wigner matrix, respectively. To relax the independence of the entries of Xi, Silverstein [9] considered the case of Xi=Mxi, where M is a nonnegative definite matrix and xi consists of independent and identically distribution (i.i.d.) entries, and the sample covariance matrices take the form of Bn=1/n∑i=1nMxixi⊤M⊤. Random matrices of the form Bn were extensively investigated by many researchers, including Yin [10], Silverstein [9], Bai and Zhou [11], Jin et al. [12], Zhang et al. [13] and Yao [14]. Silverstein’s important work [9] aimed to relax the independence structure between Xi coordinates and considered random vectors of the form Mxi, assuming that the spectral norm of the matrix sequence MM⊤ is bounded and that the ESD of the matrix sequence MMT is convergent. Silverstein proved that the LSD of the sample covariance matrix Bn exists, given by the LSD of equation representation of the Stieltjes transformation.This paper investigates the LSD of large-dimensional sample covariance matrices generated by the PAR model, and the sample covariance matrices take the form ofBn. A similar ARMA-type processes problem has been considered by Bai and Zhou [11], Jin et al. [12], Yao [14], etc. It follows that, in the ARMA-type processes, the covariance matrices MM⊤ equals to the Toeplitz matrix, according to the fundamental eigenvalue distribution theorem for Toeplitz-type matrices [15], the ESD of MM⊤ weakly converges to the nonrandom distribution H. For some situations, the ESD of MM⊤ have explicit form; then, the LSD for the sample covariance matrix Bn can be solved explicitly. For the random matrix generated by the first-order vector autoregressive (VAR(1)) model and the first-order vector moving average (VMA(1)) model, Jin et al. [12] and Bai and Zhou [11] obtained a display representation of this LSD. However, in our first-order variable-coefficient PAR process, we consider the process be interrupted by certain factors that the odd terms and the even terms of the coefficient αv converge to α and β, respectively, and that α∈0,1andβ∈0,1. The mixing matrix Q is a lower triangular matrix whose main diagonal elements are all 1 and subdiagonal elements differ from each other. When α=β, the limiting form of the matrices QQ⊤ are Toeplitz matrices, and the eigenvalues of Toeplitz matrices have been proved by Gregory [16]; when α≠β, the ESD of QQ⊤ has no explicit form, see [17, 18], and the structure of the limiting matrices Q−1Q−1⊤ are very complex. Our main contributions in this paper include the following. (i) For the first-order variable-coefficient PAR model, the explicit forms of limiting spectral moments for Q−1Q−1⊤ are given; through verifying the Carleman condition, we derive the explicit spectral moments of 1/n∑i=1nXiXi⊤. (ii) A framework of the maximum entropy-based method is provided for estimate the LSD.The paper is organized as follows. In Section2, we employ the moment method to deduce the spectral density of sample covariance matrices. In Section 3, we use the maximum entropy method to find the approximate solution and analyze the convergence of the maximum entropy solution. We simulate the PAR model to obtain the periodic sample covariance matrices, and the experiment shows that the maximum entropy spectral density function curve can match histograms of the eigenvalues of the covariance matrices very well.


## 2. Methodology and Main Results

In this section, we first derive the explicit spectral moments form ofQ−1Q−1⊤; subsequently, we verify the spectral moments statistics satisfying the Carleman condition, and then, we employ the moment method to deduce the spectral density of sample covariance.We consider first-order PAR model:(2)Xip+v=αvXip+v−1+εip+v.Let(3)Q=10⋯⋯0−α110⋯00−α21⋯0⋮⋱⋱⋱⋮0⋯⋯−αp−11,Xi=Xip+1,…,Xip+p⊤,εi=εip+1,…,εip+p⊤.Then, the first-order PAR model can be defined as(4)QXi=εi⇒Xi⊤Q⊤=εi⊤,i=1,…,n,where ⇒ denotes the equation on the left can push the equation on the right. From equation (3), Q is a lower triangular matrix whose main diagonal elements are all 1, and therefore, the matrix Q is invertible. Hence,(5)QXiXi⊤Q⊤=εiεi⊤⇒XiXi⊤=Q−1εiεi⊤Q−1⊤⇒1n∑i=1nXiXiT=Q−11n∑i=1nεiεi⊤Q−1⊤.The sample covariance matrices1/n∑i=1nXiXi⊤ and Q−11/n∑i=1nεiεi⊤Q−1⊤ have the same ESD. Let S=1/n∑i=1nεiεi⊤, and the ESD of Q−11/n∑i=1nεiεi⊤Q−1⊤ is the same as that of Q−1Q−1⊤S. Next, we establish the existence of the LSD of the covariance matrix sequence 1/n∑i=1nXiXi⊤.Theorem 1.
LetXip+v be generated by a first-order PAR model in equation (2), and the odd terms and the even terms of the coefficient αv converge to α and β, respectively, and α∈0,1,β∈0,1. The error εip+v is independent random variables satisfying for any δ, 1/δ2np∑i,tEεi,t2Iεi,t≥δn⟶0; furthermore, we assume that p/n⟶c∈0,∞; then, the spectral moments of sample covariance matrix 1/n∑i=1nXiXi⊤ tend to(6)β1n∑i=1nXiXi⊤k=∫λkdFλ⟶a.s.∑s=1kck−s∑i1+⋯+is=k−s+1i1+⋯+sis=kk!s!∏m=1sμmimim!,where μm is the limit of spectral moments of Q−1Q−1⊤, which have an explicit form:(7)μm=α+β24αβ1−αβ2m−1∑j=0m−1Cm−1j2αβj−α−β24αβ1+αβ2m−1∑j=0m−1Cm−1j2−αβj.Proof.
The proof is divided into two steps. We first show the existence of the limit of spectral moments which is defined aslimp⟶∞1/ptrQ−1mQ−1m⊤; then, we verify the Carleman condition. From equation (3), the inverse of matrix Q can be written as(8)Q−1=10⋯⋯0α110⋯0α1α2α21⋯0⋮⋱⋱⋱⋮∏k=1p−1αk∏k=2p−1αk⋯∏k=p−1p−1αk1=10⋯⋯0a2110⋯0a31a321⋯0⋮⋱⋱⋱⋮ap1⋯⋯ap,p−11,where aij can be written as(9)aij=∏k=ji−1αk,j<i,aii=1,i,j=1,2,…,p−1.
Let(10)cmi=Cm−1+i−1i−1=m−1+i−1!i−1!m−1!.
Then,(11)Q−1m=10⋯⋯0cm2a2110⋯0cm3a31cm2a321⋯0⋮⋱⋱⋱⋮cmpap1⋯⋯cm2ap,p−11.
Furthermore, we obtain that(12)trQ−1Q−1⊤m=trQ−1mQ−1⊤m=trQ−1mQ−1m⊤.
Applying the trace operator, we have(13)1ptrQ−1Q−1⊤m=1+∑N=2p1p∑k=NpcmN2akk−N−12.
Based on the assumption thatα2v+1⟶α,α2v+1⟶β,asv⟶∞,α∈0,1,β∈0,1, we obtain(14)1p∑k=NpcmN2akk−N−12⟶cmN2−1N+1+12αN−1βN−1+−1N+12αNβN−2+αN−2βN2,asp⟶∞,and therefore,(15)1ptrQ−1mQ−1m⊤⟶α+β24αβ∑k=0∞cmk+12αβk−α−β24αβ∑k=0∞cmk+12−αβk.
Let(16)hx=∑k=0∞cmk+12xk.
Then,(17)1ptrQ−1mQ−1m⊤⟶α+β24αβhαβ−α−β24αβh−αβ.
The following section wishes to obtain a simplified form of the functionhx. From equations (10) and (16), we have(18)hx=∑k=0∞m−1+k!k!m−1!2xk.
Define(19)fx=1m−1!2∑k=0∞k+1k+2,…,m−1+kxm−1−k.
Then,hx can be rewritten as(20)hx=dm−1fxdxm−1.
Let(21)gx=1m−1!2∑k=0∞xm−1+k=1m−1!2xm−11−x.
Therefore,fx can have a simplified representation as(22)fx=dm−1gxdxm−1⋅xm−1.
By the derivation operations,(23)dm−1gxdxm−1=1m−1!2∑i=0m−1Cm−1i∏k=1im−kxm−1−im−1−i!1−x−m−i=1m−1!∑i=0m−1Cm−1ixm−1−i1−x−m−i=1m−1!1−xm,we obtain that(24)fx=xm−1m−1!1−xm,dm−1fxdxm−1=1m−1!∑i=0m−1Cm−1i∏k=1im−kxm−1−i∏r=0m−i−2m+r1−x−2m−1−i.
By partially replacingx in equation (24) with y, we define a new function(25)f˜x,y=1m−1!xm−11−y2m−1∑i=0m−1Cm−1i∏k=m−i2m−i−2kx2m−2−i1−yi.
Define(26)Hx,y=∑i=0m−1Cm−1im−im−i+1,…,2m−i−2x2m−2−i1−yi,Mx,y=∑i=0m−1Cm−1ix2m−2−i1−yi=xm−11−y+xm−1,by which(27)Hx,y=dm−1Mx,ydxm−1⋅xm−1,dm−1Mx,ydxm−1=∑i=0m−1Cm−1i∏k=1im−kxm−1−i∏k=1m−1−im−k1−y+xi=∑i=0m−1Cm−1i2xm−1−im−1!1−y+xi,so(28)Hx,y=xm−1∑i=0m−1Cm−1i2xm−1−im−1!1−y+xi.
Then,(29)f˜x,y=xm−1m−1!xm−11−y2m−1∑i=0m−1Cm−1i2xm−1−im−1!1−y+xi=xm−11−y2m−1∑i=0m−1Cm−1i2x−i1−y+xi.
We usex to take place of y and obtain that(30)hx=xm−11−x2m−1∑i=0m−1Cm−1i2x−i=11−x2m−1∑j=0m−1Cm−1j2xj.
Therefore,(31)μm=limp⟶∞1ptrQ−1mQ−1m⊤=α+β24αβ1−αβ2m−1∑j=0m−1Cm−1j2αβj−α−β24αβ1+αβ2m−1∑j=0m−1Cm−1j2−αβj.
Next, we verify the Carleman condition∑k=1∞μ2k−1/2k=∞.
From equation (30), it follows that(32)μ2k−1/2k=α+β24αβ1−αβ2m−1∑j=0m−1Cm−1j2αβj−α−β24αβ1+αβ2m−1∑j=0m−1Cm−1j2−αβj−1/2k.
Then, applying the equalityC2nn=∑j=0nCnj2 and α∈0,1,β∈0,1, we obtain(33)μ2k−1/2k≥α+β22αβ1−αβ4k−1C4k−22k−1−1/2k.
From Stirling’s formulan!=1+o12πnnne−n, we get n!≥nne−n. Then, C4k−22k−1=∏j=02k−24k−2−j/2k−1!≤4k−22k−1/2k−1!≤4k−2e/2k−12k−1; applying this inequality, μ2k−1/2k has a finite lower bound:(34)μ2k−1/2k≥1−αβ22e.
Then,(35)∑k=1∞μ2k−1/2k≥∑k=1∞1−αβ22e=∞.
With the moment convergence theorem, the spectral distribution ofQ−1Q−1⊤ tends to a nonrandom probability distribution H. Let S=1/n∑i=1nεiεi⊤ and p/n⟶c∈0,∞; by the important Theorem 2.10 of Bai [19], the ESD of 1/n∑i=1nXX⊤ tends to a nonrandom limit distribution in probability (or almost surely). The spectral moments of 1/n∑i=1nXX⊤ satisfy the following equation of Yin [10]:(36)β1n∑i=1nXX⊤k=∫λkdFλ⟶a.s.∑s=1kck−s∑i1+⋯+is=k−s+1i1+⋯+sis=kk!s!∏m=1sμmimim!.
From the moment convergence theorem, we can obtain that the ESD of1/n∑k=1nXX⊤ converges with probability close to 1 to a nonrandom distribution, determined by the spectral moment statistic.Remark 1.
By the proof of Theorem1, the limiting spectral moments of Q−1Q−1⊤ depend only on the limiting coefficient α and β. When α=β, the limiting form of the matrices QQ⊤ are Toeplitz matrices, and the eigenvalues of Toeplitz matrices have been proved, see [16]. When α≠β, the limiting matrices of QQ⊤ are tridiagonal matrices and take the form like the following:(37)−α1+α2−β0000⋯00−β1+β2−α000⋯000−α1+α2−β00⋯0000−β1+β2−α0⋯0⋮⋱⋱⋱⋱⋱⋱⋱⋮0⋯⋯⋯⋯⋯⋯⋯⋯.
The eigenvalues of this limiting matrix have no explicit form, see [18] and [17]. The limiting form of the matrices Q−1Q−1⊤ are very complex, and we obtain the explicit spectral moments of Q−1Q−1⊤.


## 3. Estimating Spectral Density Function with Maximum Entropy

Unfortunately, there are no explicit LSD forms of this sample covariance matricesBn. In this section, with the spectral moments of β1/n∑i=1nXX⊤k, we employ the maximum entropy method to estimate the LSD. Mead and Papanicolaou [20] studied the maximum entropy method to solve the uncertainty problem. Granziol et al. [21] demonstrated the equivalence between the MaxEnt principle and the optimization of variables with constraints. Gao and Han [22] used the maximum entropy method to solve nonlinear expectation. Every possible event has a probability of occurring, information entropy is related to all possibilities, and a system with the highest information entropy is most likely, when all known information, taken into account. In mathematics, the variable is generally estimated by maximising the Shannon entropy of the information of the moments of the random variable of the density function:(38)S=−∫Ωpxlnpxdx,where Ω is the support of the distribution. In this paper, the sample covariance matrices Q−11/n∑i=1nεiεi⊤Q−1⊤ is nonnegative definite, so the estimated limiting distribution should be supported on Ω=ℝ+∪0.Given the firstN + 1 moments β0,β1,…,βN, the idea of determining the density function of a probability distribution based on the principle of maximum entropy is to find the density function px such that(39)∫Ωxkpxdx=βk,k=0,1,…,N,where β0=1.Lagrange multiplier is a method for finding the extremes of a multivariate function under a set of constraints. With the help of the Lagrange multiplier method, the entropy function is defined as(40)H≡∫Ω−pxlnpx+∑k=0Nλkxkpxdx−∑k=0Nλkβk.When the first-order derivative of the functionH with respect to λk and px are equal to zero, the constraint defined in equation (39) is returned; then, the first-order derivative of the function do not provide any additional information. The second of these derivatives, however, evaluate to(41)lnpNx=−1+∑k=0Nλkxk.The general solution of equation (41) is(42)pNx=e−1+∑k=0Nλkxk=e∑k=0Nλkxk,where we have set λ0−1⟶λ0.We obtain the maximum entropy solution by solving for the Lagrange multiplierλk. To do this, we must solve the following nonlinear equation:(43)∫Ωeλ0+λ1x+⋯+λkxNdx=β0∫Ωxeλ0+λ1x+⋯+λkxNdx=β1⋮∫Ωxkeλ0+λ1x+⋯+λkxNdx=βk.WhenN≤2, we can obtain a display representation of the Lagrange multiplier. Taking N = 2 as an example, from equation (36), we have(44)μ1=2+α2+β221−α2β2,μ2=α+β21+αβ4−α−β21−αβ44αβ1−α2β23,β1=μ1,β2=yμ12+μ2.Based on equation (40), we obtain λ0,λ1,λ2 as(45)λ0=−β12+ln2πβ2−β12+1,λ1=β1β2−β12,λ2=−12β2−β12.Then,(46)px=expλ0+λ1x+λ2x2.WhenN≥3, we cannot obtain a display solution for the Lagrange multiplier λk, but at this point, the numerical solution for λk can be calculated using numerical optimization algorithms. From equation (40), the analytic expressions of the gradient is(47)∂H∂λj=βj−∫Ωxjexp∑i=0mλixidx.If we denote the maximum entropy solution obtained by the method based on moments as(48)μ˜k≡∫Ωxkeλ0+λ1x+⋯+λkxkdx,then the Hessian matrix is given by(49)G=μ˜0μ˜1⋯μ˜kμ˜1μ˜2⋯μ˜k+1μ˜2μ˜3⋯μ˜k+2⋮⋮⋮⋮μ˜kμ˜k+1⋯μ˜k+k.μ˜k is the kth moment of the simulated distribution, and the Newton iteration is done by searching λk until the termination criterion μ˜k−μk<ε holds. We obtain the spectral density px through Algorithm 1.Algorithm 1: The proposed algorithm.
Input:β=β0,…,βN, Ending accuracy ε, Diitter variance in Hessian η=10−4;Output: Coefficientsλi;initialλ0=−ln2π,0,…,0, k = 0;repeatCompute gradientg=g0,…,gN: gj=βj−∫Ωxjexp∑i=0Nλkixidx;Compute HessianG;Direction of Newton’s law:dk=−G+ηI−1⋅g;Next iteration point:λk+1=λk+dk;k = k+1;untilβ−μ˜2<ε.We can obtain that the maximum entropy solutionpNx is weakly convergent, i.e.,(50)limN⟶+∞∫ΩpNxlnpNxdx=∫Ωpxlnpxdx.In our simulations, we use the Newton method from the Python-based SciPy package to solve for minimisationμ˜k−μk. It computes the gradient within a ending accuracy ε as well as the Hessian matrix. To give the Hessian matrix a faster rate of convergence, we symmetrize it and add a disturbance of intensity η=10−4 along the diagonal. In the simulation of the algorithm, we observe that the algorithm is not too sensitive to these choices.From equations (31) and (36), we obtain the four limiting spectral moments 1/n∑i=1nXiXi⊤ of the sample covariance matrices as follows:(51)β1=μ1,β2=yμ12+μ2,β3=y2μ13+3yμ1μ2+μ3,β4=y3μ14+3y2μ12μ2+2yμ22+2μ1μ3+μ4.The authors of [23, 24] state that, as the number of moments increases (especially when n≥5), the maximum entropy method becomes sensitive to numerical differences between the order moments, the gradient, and the state of the Hessian matrix, as well as improper arithmetic accuracy becomes very sensitive and can be improved by introducing different basis functions, such as a drift with Chebyshev polynomial basis functions with drift to overcome these effects. In this section, to simplify the operations, we estimate the maximum entropy density function in terms of 4th-order spectral moments. The value of (d; n; a; b) determines the value of each spectral moment. Figures 1–6 show histograms of eigenvalues and maximum entropy spectral density curves (MaxEntPDF) for different values of (d; n; a; b). From these plots, it can be seen that the maximum entropy spectral density function curve can match histograms of the eigenvalues of the covariance matrices very well.Figure 1
d = 200, n = 100, alpha = 0.5, and beta = 0.3.Figure 2
d = 200, n = 200, alpha = 0.5, and beta = 0.3.Figure 3
d = 200, n = 400, alpha = 0.5, and beta = 0.3.Figure 4
d = 1000, n = 500, alpha = 0.5, and beta = 0.3.Figure 5
d = 1000, n = 1000, alpha = 0.5, and beta = 0.3.Figure 6
d = 1000, n = 2000, alpha = 0.5, and beta = 0.3.


## 4. Conclusion

In this paper, we investigate the LSD of large-dimensional sample covariance matrix generated by the first-order variable-coefficient PAR process. We consider the process be interrupted by certain factors that the odd terms and the even terms of the coefficientαv converge to α and β, respectively, and that α∈0,1andβ∈0,1. The main contributions in this paper include the following. (i) For the first-order variable-coefficient PAR model, the explicit form of limiting spectral moments for Q−1Q−1⊤ is given; through verifying the Carleman condition, we derive the explicit spectral moments of 1/n∑i=1nXiXi⊤. (ii) A framework of the maximum entropy-based method is provided for estimate the LSD. We have some ideas for the next step of the research. On the one hand, when generalizing the convergence of the coefficients of the PAR model to a more general structure, how about obtaining an expression for the spectral moments of the limiting covariance matrix? On the other hand, it is of importance to extend the method presented in this paper to PAR with larger order.



---

*Source: 9526991-2021-08-21.xml*