Assume $U_j$ are $\chi^2(1)$ random variables and $\lambda_1, \ldots, \lambda_p$ are the eigenvalues of a covariance matrix $\Sigma = (r^{|i-j|})_{ij}$ with a Toeplitz-type structure (for some fixed $|r|<1$).

Assuming $S_p = \sum_{j=1}^{p}\lambda_j U_j$, what is the limiting distribution of $S_p/p$? Knowing that $\sum_{j=1}^{p}\lambda_j = p$ and that $\lambda_j U_j \sim \Gamma(1/2, 1/2\lambda_j)$, can we claim the result straightforwardly through LLN or CLT?

  • $\begingroup$ How does $\Sigma$ change as $p$ is increased? Without that information, there's no way to answer this question. $\endgroup$ – whuber Mar 29 at 12:45
  • $\begingroup$ @whuber edited. $\endgroup$ – runr Mar 29 at 12:59
  • $\begingroup$ How are $n$ and $p$ related?? $\endgroup$ – whuber Mar 29 at 13:46
  • 1
    $\begingroup$ @whuber typo! Will fix, $n=p$ $\endgroup$ – runr Mar 29 at 13:47

An example in which case we get this sum is when we consider the distribution of the squared Euclidian distance (from the origin) of a multivariate normal distributed variable $X$ with covariance $\Sigma$.

Instead of considering the sum $\sum_{j=1}^p \lambda_j U_j$ we can also consider the sum of $\sum_{j=1}^p X_j^2$. Thus a sum of correlated variables. The correlation is $Cor(X_j^2,X_i^2) = Cor(X_j,X_i)^2$

Next we apply the CLT for correlated processes to obtain

$$\frac{\sum_{j=1}^p \lambda_jU_j^2 - \sum_{j=1}^p\lambda_j}{\sigma_p \sqrt{p}} = \frac{\sum_{j=1}^p (X_j^2 -1)}{\sigma_p \sqrt{p}} \to N(0,1)$$

with $$\begin{array}{} \sigma_p^2 &=& Var(X_1^2) + 2 \sum_{k=2}^{p} Cov(X_1^2,X_k^2) \\ &=& 2 \sum_{k=1}^p \rho^{k-1} = 2 \frac{\rho-\rho^p}{1-\rho} \end{array}$$


Below is a simulation compared with the normal curve.

simulation compared with normal curve


### function for simulation
simu <- function(p, rho) {
#  The mvrnorm requires lots of computation effort
#  instead we use the ar(1) process below, which is the same
#  Sigma = toeplitz(rho^(0:(p-1)))
#  mu = rep(0,p)
#  X <- mvrnorm(n,mu,Sigma)
  Y <- rnorm(p,0,sqrt(1-rho^2))
  X <- rep(0,p)
  X[1] <- rnorm(1,0,1)
  for (i in 2:p) {
    X[i] <- rho*X[i-1] + Y[i]
  out <- sum(X^2)

### settings
p = 10000
rho = 0.9
n = 10^4

### plot histogram with normal approximation

### simulate n times
xn <- replicate(n,simu(p, rho)) 
## compute the scaled variable
sig2 = (rho-rho^p)/(1-rho)*2
v <- (xn-p)/sqrt(p)/sqrt(sig)

### plot
hist(v, breaks = seq(floor(min(v)),ceiling(max(v)),0.1), xlim = c(-5,5), freq = 0,
     main = "when p = 10^4 and rho = 0.9",
     xlab = "scaled and shifted variable")
### normal distribution as comparison
xs <- seq(-5,5,0.1)
lines(xs, dnorm(xs))
  • $\begingroup$ Alternatively, you might figure out the eigenvalues. I started trying to figure it out and stranded with this answer here math.stackexchange.com/a/2816081 . That answer refers to a reference that I could not directly obtain, so then I thought it might be easier to tackle it as an ar(1) process or a sum of correlated variables. $\endgroup$ – Sextus Empiricus Mar 29 at 15:55
  • $\begingroup$ Thanks for the great writeup. The link to the eigenvalues is also very useful, didn't know that the form of $\Sigma$ was called a KMS matrix. Would be interesting to see an equivalent result directly through the eigenvalues $\lambda_j$, but judging by the linked paper I can see why the straightforward derivation might be more practical.. $\endgroup$ – runr Mar 29 at 16:18
  • $\begingroup$ On a slightly related note, regarding my other question, could I apply these results in a similar manner to a product $p^{-2}(X'AX)(X'BX)$ and, via Slutsky, arrive at a limiting product normal distribution? Or is it not that simple? In the previous question, I went to the eigenvalues in order to avoid products of correlated squared normals, however, I'm not entirely sure I can go around them. $\endgroup$ – runr Mar 29 at 16:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.