Limiting distribution of $\sum_{j=1}^{p}\lambda_j U_j$

Question

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?

Answer 1

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}$$

---

Simulations

Below is a simulation compared with the normal curve.

library(MASS)

### 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)
  return(out)
}

### settings
set.seed(1)
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))