Let $T,N$ be positive integers subject to $T\ge N+1$. A sample correlation matrix(to be termed a Wishart matrix) is defined as $c:= 1/T {\bf X}\cdot {\bf X}^T$ where ${\bf X}:= {\tilde C} \cdot {\bf Y}$ with ${\bf Y}:= \left( Y_{i,t} \right)_{i=1,t=1}^{N,T}$ being a random matrix with iid entries conforming to a standard Gaussian distribution. Here ${\tilde C}$ has a size $N \times N$ and is such that the matrix $C:={\tilde C} \cdot {\tilde C}^T$ (to be termed the underlying correlation matrix) is positively definite .
Having all this in mind and following the lead from my previous question let us define a normalization constant ${\mathfrak P}_{N,T} = \frac{\pi ^{\frac{N^2}{2}-\frac{1}{2} (N-1) N} 2^{-\frac{N T}{2}}}{\left(\prod _{j=0}^{N-1} \Gamma \left(\frac{N-j}{2}\right)\right) \prod _{j=0}^{N-1} \Gamma \left(\frac{T-j}{2}\right)}$ and a parameter ${\mathfrak a} := (T-N-1)/2 $ . Now braced with this notation we can write down the joint probability density $\rho_{\Lambda_1,\cdots,\Lambda_N}(\lambda_1,\cdots,\lambda_N)$ of eigenvalues $\left( \lambda_j \right)_{j=1}^N $, subject to $\lambda_ \ge \lambda_2 \ge \cdots \lambda_N \ge 0 $, of our sample correlation matrix (see Theorem 7.1.2 page 445 in:
Harnad, John (ed.), Random matrices, random processes and integrable systems, CRM Series in Mathematical Physics. Berlin: Springer (ISBN 978-1-4419-9513-1/hbk; 978-1-4419-9514-8/ebook). xviii, 524 p. (2011). ZBL1215.15002. ) as follows:
\begin{eqnarray} &&\rho_{\Lambda_1,\cdots,\Lambda_N}(\lambda_1,\cdots,\lambda_N) = \\ && \frac{{\mathfrak P}_{N,T}}{\det(C)^{\frac{T}{2}}} \cdot \left(\prod\limits_{j=1}^N \lambda_j^{\mathfrak a}\right) \cdot \left( \prod\limits_{1 \le \lambda_j < \lambda_k \le N} (\lambda_j - \lambda_k) \right) \cdot \int\limits_{{\mathbb O}(N)} \exp\left( -\frac{1}{2} Tr[C^{-1} \cdot H \cdot diag(\lambda_j)_{j=1}^N \cdot H^T] \right) d H \tag{1} \end{eqnarray}
where the integral is over the de Haar measure of the orthogonal group ${\mathbb O}(N)$ of $N\times N $ matrices.
Now, the spectral moments $m_p:= 1/N \cdot Tr[c^p] = \int\limits_0^\infty \lambda^p \rho(\lambda) d\lambda $ for $p \in {\mathbb Z} $ of the sample correlation matrix are defined as the moments of the spectral density $\rho(\lambda) := 1/N \int\limits_{\lambda_1 \ge \lambda_2 \ge \cdots \lambda_N \ge 0} \left( \sum\limits_{j=1}^N \delta(\lambda- \lambda_j) \right) \rho_{\Lambda_1,\cdots,\Lambda_N}(\lambda_1,\cdots,\lambda_N) \prod\limits_{j=1}^N d \lambda_j $. It is our objective to find closed form expressions for those moments. We denote by $A_j:= @[\lambda_j]\left( Tr[ C^{-1} \cdot H \cdot diag(\lambda_j)_{j=1}^N \cdot H^{'}] \right)$ for $j=1,\cdots,N$. Then we can integrate over the eigenvalues by using the trick from my other question and then the general result reads:
\begin{eqnarray} &&m_p = \frac{1}{N T^p} \frac{{\mathfrak P}_{N,T}}{(\det(C))^{\frac{T}{2}}} \cdot 2^{N {\mathfrak a}+ p+ \binom{N+1}{2}} \sum\limits_{J=1}^N \\ && \int\limits_{O(N)} \prod\limits_{j=1}^N \left( \partial_{{\bar A}_j} + \cdots +\partial_{{\bar A}_N}\right)^{{\mathfrak a} + \delta_{j,J} p} \cdot \prod\limits_{1 \le i < j \le N} \left( \partial_{{\bar A}_i} + \cdots \partial_{{\bar A}_{j-1}}\right) \left. \prod\limits_{j=1}^N \frac{1}{{\bar A}_j} \right|_{{\bar A}_j = A_1+\cdots+A_j} dH \tag{1a} \end{eqnarray}
Now, starting from $(1a)$ and using the parametrization $H = \left( \begin{array}{lll} \cos(\phi) & -\sin(\phi) \\ \sin(\phi) & \cos(\phi) \end{array} \right) $ with $\phi \in (0, 2\pi) $ of the two-dimensional orthogonal group we managed to compute the negative moments in question in the case $N=2$ and $C= \left( \begin{array}{lll} 1 & \rho \\ \rho & 1 \end{array} \right) $. We denote by $(s_1,s_2):= (1+(N-1) \rho, 1-\rho)$ the eigenvalues of the underlying correlation matrix and then (for the unnormalized moments) we have:
\begin{eqnarray} m_p &=& \frac{ {\mathfrak P}_{N,T} }{ (\prod\limits_{\xi=1}^N s_\xi)^{\frac{T}{2}} } \cdot \frac{1}{N \cdot T^p} \cdot 2^{\binom{N+1}{2} + N {\mathfrak a} +p} \cdot (\frac{ s_1 s_2}{s_1+s_2})^{2 {\mathfrak a}+p+3} \cdot \\ && \sum\limits_{J=1}^N \sum\limits_{k_1=0}^{{\mathfrak a} + p \cdot 1_{J=1}} \binom{{\mathfrak a}+ p \cdot 1_{J=1}}{k_1} (k_1+1)! (2 {\mathfrak a}+p-k_1)! 2^{k_1+2} F_{2,1} \left[ \begin{array}{lll} \frac{k_1+2}{2} & \frac{k_1+3}{2} \\ 1 \end{array}; \left( \frac{s_1-s_2}{s_1+s_2} \right)^2 \right] \tag{2a} \\ \end{eqnarray}
The final result (for the normalized moments) reads: \begin{eqnarray} \frac{m_{p}}{m_0} &=& \frac{T^{-p}}{\left(\prod\limits_{j=1}^{-p} (T-j)\right)\left( \prod\limits_{j=1}^{-p} (\frac{T-1}{2} - j) \right)} \cdot \frac{1}{(1-\rho^2)^{-p}} \cdot \\ && \sum\limits_{q=0}^{\lfloor - \frac{p}{2} \rfloor} \binom{-p}{2 q} \sum\limits_{l=0}^{-{\tilde p}-{\bar p}-q} \sum\limits_{l_2=0}^q \binom{q}{l_2} \frac{(l+q)!}{(l+l_2)!} \binom{-{\tilde p}-{\bar p}-q}{l} \cdot \\ && \left( \prod\limits_{j=p+q+l}^{-1} (\frac{T+1}{2} + j) \right) \cdot \left( \prod\limits_{j={\tilde p}+{\bar p}}^{{\tilde p}+{\bar p}+l+l_2-1} (\frac{T}{2} + j) \right) (\rho^2)^{l+l_2} \cdot (1-\rho^2)^{-{\tilde p}-{\bar p}-l-l_2} \tag{2b} \end{eqnarray}
where ${\tilde p}:= \lfloor p/2 \rfloor$, ${\bar p}:= p\%2$ and $p=0,-1,-2,\cdots$.
As always, the code snippet below verifies the result $(2b)$ by comparing them with the generic result, that we managed to guess for $p\ge -3$. Here we go:
(*Here we obtain all spectral moments for small values of N=2,3,..*)
Clear[Pfct, M];
mGamma[NN_, x_] :=
Pi^(NN (NN - 1)/4) Product[Gamma[(2 x - j)/2], {j, 0, NN - 1}];
Pfct[NN_, T_] :=
Pi^(NN^2/2) 2^(-NN T/2)/(mGamma[NN, NN/2] mGamma[NN, T/2]);
(*Firstly N\[Equal]2. Here all can be done analytically.*)
Clear[k1, k2, m, M, p, s, c1, c2, r]; NN = 2; Xi =
RandomInteger[{3, 10}];
T = NN + 1 + 2 Xi; a = (T - NN - 1)/2; rr = NN/T;
{NN, T};
rho =.;
rho = RandomReal[{0, 1}, WorkingPrecision -> 50];
s[1] = 1 + (NN - 1) rho;
s[2] = 1 - rho;
M[p_] := 1/2 (s[1]^p + s[2]^p);
m[p_] :=
1/ Product[s[xi], {xi, 1, NN}]^(T/2) Pfct[NN,
T]/(NN T^p) 2^(Binomial[NN + 1, 2] + NN a +
p) ((s[1] s[2])/(s[1] + s[2]))^(2 a + p + 3) Sum[
Binomial[a + If[J == 1, 1, 0] p,
k1] (k1 + 1)! (2 a + p - k1)! 2^(k1 + 2) Hypergeometric2F1[(
k1 + 2)/2, (3 + k1)/2, 1, ((s[1] - s[2])/(s[1] + s[2]))^2], {J,
1, NN}, {k1, 0, a + If[J == 1, 1, 0] p}];
(Print["NN,T,rho=",{NN,T,rho}];)
SetOptions[NIntegrate, WorkingPrecision -> 20, PrecisionGoal -> 15];
ll = Table[m[p], {p, 0, -3, -1}]; lmax = 20;
ll /= ll[[1]];
res1 = ll // Simplify
res2 = Table[((1 - rho) (1 + rho))^(p + T/2) T^-p Gamma[p + T] /(
4 Gamma[((T - 1))])
NIntegrate[(1 - t)^(
1/2 (-3 + T)) ((1 - Sqrt[t])^p + (1 + Sqrt[t])^p) (1 -
rho^2 t)^(1/2 (-p - T))
LegendreP[-1 + p + T, 1/Sqrt[1 - rho^2 t]], {t, 0, 1}], {p,
0, -3, -1}];
res2a = Table[
T^-p Pochhammer[T, p] Pochhammer[(T - 1)/2,
p] ((1 - rho) (1 + rho))^(p + T/2)
Sum[Binomial[-p, 2 q] 1/
Pochhammer[(T + 1)/
2 , +p + q] (D[
x^q HypergeometricPFQ[{p/2 + T/2,
1/2 + p/2 + T/2}, {1/2 + p + q + T/2}, x rho^2], {x,
q}] /. x :> 1), {q, 0, Floor[-p/2]}], {p, 0, -3, -1}];
res2b = Table[ T^-p/((!(
*UnderoverscriptBox[([Product]), (j = 1), (-p)]((T -
j)))) (!(
*UnderoverscriptBox[([Product]), (j = 1), (-p)]((
*FractionBox[(T - 1), (2)] - j)))))
1/(1 - rho^2)^-p Sum[
Binomial[-p, 2 q] (With[{pt = Floor[p/2], pm = Mod[p, 2]},
Sum[(Binomial[q, l2] Gamma[1 + l + q]/
Gamma[1 + l + l2] Binomial[-pt - q - pm, l]) (!(
*UnderoverscriptBox[([Product]), (j = p + q + l), (-1)]((
*FractionBox[(T + 1), (2)] + j)))) (!(
*UnderoverscriptBox[([Product]), (j = pm + pt), (pm + pt + l +
l2 - 1)]((
*FractionBox[(T), (2)] + j)))) (rho^2)^(
l + l2) (1 - rho^2)^(-pm - pt - l - l2), {l,
0, -pt - q - pm}, {l2, 0, q}]]), {q, 0, Floor[-p/2]}], {p,
0, -3, -1}];
res3 = {1,
T /(T - NN - 1) M[-1], (
T^2 (T - 1))/((T - NN) Product[T - NN - 2 j - 1, {j, 0, 1}]) (
M[-2] - NN/(T - 1) (M[-2] - M[-1]^2)),
T^5/( Product[T - NN + j, {j, -1, 1}] Product[
T - NN - 2 j - 1, {j, 1, 2}]) (M[-3] (1 - 1/T)^2 +
NN/T (1 - 1/T) ( -2 M[-3] +
3 M[-2] M[-1]) - (NN/(T))^2 (-M[-3] + 3 M[-2] M[-1] -
2 M[-1]^3))} // Simplify;
res2/res1
res2a/res1
res2b/res1
res3/res1
$N=$." />In view of all this my question is two fold. First of all have those results been derived before, if yes then please provide some reference from the literature. My second question would be how would one derive the results for higher values of $N$.
Update:
From $(1a)$ it is clear that computing the spectral moments in general is equivalent to averaging over the orthogonal group. In view of this my question would be as follows. Let $k_j \in {\mathbb N}_+$ for $j=1,\cdots, N-1$. Define an following integral over the orthogonal group $H \in O(N)$. We have:
\begin{eqnarray} {\mathfrak I}_{k_1,\cdots,k_{N-1}}\left(C^{-1}\right) &:=& \int\limits_{O(N)} \prod\limits_{j=1}^{N-1} \frac{1}{\left(Tr\left[C^{-1} \cdot H \cdot diag(\underbrace{1,\cdots,1}_{j},\underbrace{0,\cdots,0}_{N-j} )\cdot H^T\right]\right)^{k_j}} dH \end{eqnarray}
How do we evaluate this integral?
