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

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?

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

## 1 Answer

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 <- 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))

• 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. – Sextus Empiricus Mar 29 at 15:55
• 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.. – runr Mar 29 at 16:18
• 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. – runr Mar 29 at 16:31