Finding covariance between internally studentized residuals

Question

Consider a linear regression model $\boldsymbol y=X\boldsymbol\beta+\boldsymbol\varepsilon$, where $\boldsymbol y$ is an $n\times 1$ response vector, $X$ is an $n\times p$ matrix of covariates (fixed), and the error vector $\boldsymbol\varepsilon$ is multivariate normal $N_n(\boldsymbol 0,\sigma^2I)$.

The $i$th internally studentized residual is

$$r_i=\frac{e_i}{\hat\sigma\sqrt{1-h_{ii}}}\,,$$

where $e_i=y_i-\hat y_i$ is the $i$th residual, $h_{ii}$ is the $i$th diagonal entry of the hat matrix $H=X(X^TX)^{-1}X^T$, and $\hat\sigma^2=\frac{\boldsymbol e^T\boldsymbol e}{n-p}$ is the usual unbiased estimator of $\sigma^2$.

I am trying to find $\operatorname{Cov}(r_i,r_j)$ without explicitly finding the joint distribution of $(r_i,r_j)$.

I can see that $\operatorname E(r_i)=0$, since $\frac{e_i}{\lVert \boldsymbol e\rVert}$ is symmetric about $0$.

Now, for $\operatorname{Cov}(r_i,r_j)=\operatorname E(r_ir_j)$, I need $$\operatorname{Cov}\left(\frac{e_i}{\lVert \boldsymbol e\rVert},\frac{e_j}{\lVert \boldsymbol e\rVert}\right)=\operatorname E\left[\frac{e_ie_j}{\boldsymbol e^T\boldsymbol e}\right]$$

Any suggestion in finding this quantity is welcome.

I know that $\boldsymbol e\sim N_n(\boldsymbol 0,\sigma^2(I-H))$ and $\frac{r_i^2}{n-p}\sim \text{Beta}\left(\frac12,\frac{n-p-1}{2}\right)$.

---

The answer is supposed to be the same as the covariance if $\sigma$ was not estimated:

$$\operatorname{Cov}(r_i,r_j)=\operatorname{Cov}\left(\frac{e_i}{\sigma\sqrt{1-h_{ii}}},\frac{e_j}{\sigma\sqrt{1-h_{jj}}}\right)=\frac{-h_{ij}}{\sqrt{1-h_{ii}}\sqrt{1-h_{jj}}}\,,$$

where $h_{ij}$ is the $(i,j)$th entry of $H$.

Answer 1

Preparations

Lemma 1. If $X = (X_1, X_2, \ldots, X_d) \sim N_d(0, I_{(d)})$, then \begin{align} & E\left[\frac{X_j}{\|X\|}\right] = 0, 1 \leq j \leq d, \\ & E\left[\frac{X_iX_j}{\|X\|^2}\right] = \begin{cases} 0, & 1 \leq i \neq j \leq d, \\ d^{-1}, & 1 \leq i = j \leq d. \end{cases} \end{align}

Lemma 2. Suppose an order $n$ matrix $A$ is idempotent and symmetric with $\operatorname{rank}(A) = r$, then there exists an order $n$ orthogonal matrix $P$ such that \begin{align} A = P\operatorname{diag}(I_{(r)}, 0)P'. \end{align}

Proof

Without loss of generality, assume $\sigma^2 \equiv 1$. I am also assuming you are only interested in the case $i \neq j$ (the case $i = j$ is similar).

Denote $I_{(n)} - H$ by $\bar{H}$, then $\bar{H}$ is symmetric, idempotent, and $\operatorname{rank}(\bar{H}) = n - p =: r$. By Lemma 2, there exists an order $n$ orthogonal matrix $P$ such that $\bar{H} = P\operatorname{diag}(I_{(r)}, 0)P'$. Denote $P'\varepsilon$ by $X := (\xi, \eta)$, where $\xi \in \mathbb{R}^r$ and $\eta \in \mathbb{R}^{n - r}$. Since $P$ is orthogonal and $\varepsilon\sim N_n(0, I_{(n)})$, it follows that $X \sim N_n(0, I_{(n)})$ and $\xi \sim N_r(0, I_{(r)})$. Denote the matrix formed by the first $r$ columns of $P$ by $Q$ and the matrix formed by the last $n - r$ columns of $P$ by $R$, then \begin{align} & e = \bar{H}\varepsilon = P\operatorname{diag}(I_{(r)}, 0)P'\varepsilon = \begin{bmatrix} Q & R \end{bmatrix}\begin{bmatrix} \xi \\ 0 \end{bmatrix} = Q\xi. \tag{1} \\ & e'e = \|e\|^2 = \xi'Q'Q\xi = \|\xi\|^2. \tag{2} \end{align}

Denote the $i$-th row and $j$-th row of $Q$ by $\begin{bmatrix}q_{i1} & \cdots & q_{ir} \end{bmatrix}$ and $\begin{bmatrix}q_{j1} & \cdots & q_{jr} \end{bmatrix}$ respectively, it then follows by $(1)$ and $(2)$ that \begin{align} E\left[\frac{e_ie_j}{\|e\|^2}\right] = E\left[\frac{\sum_{k = 1}^r q_{ik}\xi_k \cdot \sum_{l = 1}^r q_{jl}\xi_l}{\|\xi\|^2}\right], \tag{3} \end{align} which by Lemma 1 and $\xi \sim N_r(0, I_{(r)})$ becomes to \begin{align} \sum_{k = 1}^r q_{ik}q_{jk}E\left[\frac{\xi_k^2}{\|\xi\|^2} \right] = r^{-1}\sum_{k = 1}^r q_{ik}q_{jk}. \tag{4} \end{align}

Finally, $\bar{H} = I_{(n)} - H = QQ'$ implies that \begin{align} \sum_{k = 1}^r q_{ik}q_{jk} = \bar{H}(i, j) = -h_{ij}. \tag{5} \end{align} Combining $(3), (4), (5)$ yields: \begin{align} E\left[\frac{e_ie_j}{\|e\|^2}\right] = -r^{-1}h_{ij}, \end{align} which immediately gives \begin{align} \operatorname{Cov}(r_i, r_j) = \frac{-h_{ij}}{\sqrt{1 - h_{ii}}\sqrt{1 - h_{jj}}}. \end{align}

---

Proof of Lemmas

Lemma 2 is a standard linear algebra result. So I will only discuss how to prove Lemma 1, which is a quite interesting result. Lemma 1 is a corollary of the following well-known theorem of spherical distributions (cf. Theorem 1.5.6 in Aspects of Multivariate Statistical Theory by R. Muirhead):

Theorem. If $X$ has an $m$-variate spherical distribution with $P(X = 0) = 0$ and $r = \|X\| = (X'X)^{1/2}, T(X) = \|X\|^{-1}X$, then $T(X)$ is uniformly distributed on $S_m$ and $T(X)$ and $r$ are independent.

The proof to this theorem is not hard (but might be tedious) by laying out the spherical coordinates of $X$.

To apply this theorem to the proof of Lemma 1, write \begin{align} & X_j = \frac{X_j}{\|X\|} \times r, \tag{6} \\ & X_iX_j = \frac{X_iX_j}{\|X\|^2} \times r^2. \tag{7} \end{align}

Then taking expectations on both sides of $(6), (7)$ and applying the independence of multipliers (by Theorem) finishes the proof.