Another multi-dimensional integral with applications in Random Matrix Theory.

Question

This question generalizes my other question here.

Let $\vec{a}:= \left( a_j\right)_{j=1}^N \in {\mathbb N}^N$ and let $\vec{A} := \left(A_j\right)_{j=1}^N$ . We define a following multidimensional integral:

\begin{equation} {\mathfrak J}_N^{\vec{a}} \left( \vec{A} \right):= \int\limits_{\infty > \lambda_1 \ge \lambda_2 \ge \cdot \ge \lambda_N \ge 0} \prod\limits_{\xi=1}^N \lambda_\xi^{a_\xi} \cdot \prod\limits_{1 \le i < j \le N} (\lambda_i - \lambda_j) \cdot e^{-\vec{A} \cdot \vec{\lambda}} d^N\lambda \quad (1) \end{equation}

The quantity in $(1)$ is important in the Random Matrix Theory, in particular in computing spectral moments of the sample correlation matrix.

By using the technique from the the other post linked above we know that ${\mathfrak J}_N^{\vec{a}} \left( \vec{A} \right) = \prod\limits_{i=1}^N ( \partial_{{\bar A}_i} + \cdot \partial_{{\bar A}_N} )^{a_i} \cdot \prod\limits_{1 \le i < j \le N} ( \partial_{{\bar A}_i} + \cdots + \partial_{{\bar A}_{j-1}} ) \cdot \prod\limits_{j=1}^N 1/{\bar A}_j $ where ${\bar A}_j:= A_1+\cdots+A_j$ for $j=1,\cdots,N$. This fact allows us to find the functional form of the solution which reads as follows. We have:

\begin{equation} {\mathfrak J}_N^{\vec{a}} \left( \vec{A} \right) = \sum\limits_{\begin{array}{l} \lambda_1+\cdots \lambda_N = \left|\vec{a}\right| + \binom{N+1}{2} \\ \lambda_1 \ge 0, \cdots, \lambda_N \ge 0 \end{array}} {\mathfrak C}^{(N)}_{\lambda_1,\cdots,\lambda_N}(\vec{a}) \cdot \prod\limits_{j=1}^N \frac{1}{{\bar A}_j^{\lambda_j}} \quad (2) \end{equation}

And now , all what we need to do in order to compute our quantity is to find recurrence relations for our coefficients ${\mathfrak C}^{(N)}_{\vec{\lambda}}(\vec{a})$.

We will give those recurrences for particular values of $N=2,3,4$ and then we will verify the equations by using a Mathematica code snippet. Here we go:

\begin{eqnarray} && \left. {\mathfrak C}^{(2)}_{\lambda_1,\lambda_2}(\vec{a}) = \binom{a_1}{\lambda_1 - 2} \cdot (\lambda_1-1)! \cdot (a_1+a_2+2-\lambda_1)! \cdot 1_{a_1+2 \ge \lambda_1 \ge 2} \cdot 1_{\lambda_1+\lambda_2 = \left|\vec{a}\right|+3} \right. \\ %%%%%%%%%%%%%%%%%%%%%% && \left. {\mathfrak C}^{(3)}_{\lambda_1,\lambda_2}(\vec{a}) = \sum\limits_{\begin{array}{l} l_1+l_2+l_3= \lambda_3-1 \\ l_1=0,\cdots,a_1 \\l_2=0,\cdots,a_2 \\ l_3 = a_3 \end{array}} \binom{a_1}{l_1} \binom{a_2}{l_2} (\lambda_3-1)! \right. \\ && \left[ \right. \\ && \left. (\lambda_1-1)_{(1)} (\lambda_2-1)_{(1)} {\mathfrak C}^{(2)}_{\lambda_1-1,\lambda_2-1}(a_1-l_1,a_2-l_2) + \right. \\ && \left. (\lambda_1-1)_{(0)} (\lambda_2-1)_{(2)} {\mathfrak C}^{(2)}_{\lambda_1-0,\lambda_2-2}(a_1-l_1,a_2-l_2) \right. \\ && \left. \right] \cdot 1_{\lambda_1+\lambda_2+\lambda_3 = \left|\vec{a}\right| +6} \cdot 1_{\lambda_1 \ge 2} \cdot 1_{\lambda_3 \ge 1} \\ %%%%%%%%%%%%%%%%%%%%%% && \left. {\mathfrak C}^{(4)}_{\lambda_1,\lambda_2,\lambda_3}(\vec{a}) = \sum\limits_{\begin{array}{l} l_1+l_2+l_3+l_4= \lambda_4-1 \\ l_1=0,\cdots,a_1 \\l_2=0,\cdots,a_2 \\ l_3=0,\cdots, a_3 \\ l_4 = a_4 \end{array}} \binom{a_1}{l_1} \binom{a_2}{l_2} \binom{a_3}{l_3} (\lambda_4-1)! \right. \\ && \left[ \right. \\ && \left. 1 \cdot (\lambda_1-1)_{(1)} (\lambda_2-1)_{(1)} (\lambda_3-1)_{(1)} {\mathfrak C}^{(3)}_{\lambda_1-1,\lambda_2-1,\lambda_3-1}(a_1-l_1,a_2-l_2,a_3-l_3) + \right. \\ && \left. 1 \cdot (\lambda_1-1)_{(0)} (\lambda_2-1)_{(2)} (\lambda_3-1)_{(1)} {\mathfrak C}^{(3)}_{\lambda_1-0,\lambda_2-2,\lambda_3-1}(a_1-l_1,a_2-l_2,a_3-l_3) + \right. \\ && \left. 1 \cdot (\lambda_1-1)_{(1)} (\lambda_2-1)_{(0)} (\lambda_3-1)_{(2)} {\mathfrak C}^{(3)}_{\lambda_1-1,\lambda_2-0,\lambda_3-2}(a_1-l_1,a_2-l_2,a_3-l_3) + \right. \\ && \left. 2 \cdot (\lambda_1-1)_{(0)} (\lambda_2-1)_{(1)} (\lambda_3-1)_{(2)} {\mathfrak C}^{(3)}_{\lambda_1-0,\lambda_2-1,\lambda_3-2}(a_1-l_1,a_2-l_2,a_3-l_3) + \right. \\ && \left. 1 \cdot (\lambda_1-1)_{(0)} (\lambda_2-1)_{(0)} (\lambda_3-1)_{(3)} {\mathfrak C}^{(3)}_{\lambda_1-0,\lambda_2-0,\lambda_3-3}(a_1-l_1,a_2-l_2,a_3-l_3) + \right. \\ && \left. \right] \cdot 1_{\lambda_1+\lambda_2+\lambda_3 +\lambda_4= \left|\vec{a}\right| +10} \cdot 1_{\lambda_1 \ge 2} \cdot 1_{\lambda_4 \ge 1} \\ \end{eqnarray}

where $a_{(n)}:= a(a-1) \cdots (a-n+1)$ is the lower Pochhammer symbol.

Now comes the verification.

CC[2, lmb_List, a_List] := 
  If[a[[1]] >= lmb[[1]] - 2 && lmb[[1]] >= 2, 
   Binomial[a[[1]], 
     lmb[[1]] - 2] (lmb[[1]] - 1)! (a[[1]] + a[[2]] + 2 - lmb[[1]])!, 
   0] ;
l2 =.;
CC[3, lmb_List, a_List] := Sum[With[{l2 = lmb[[3]] - 1 - a[[3]] - l1},
    Binomial[a[[1]], l1] Binomial[a[[2]], l2] (lmb[[3]] - 1)!
     (
         (lmb[[1]] - 1) ( lmb[[2]] - 1) CC[
         2, {lmb[[1]] - 1, lmb[[2]] - 1}, {a[[1]] - l1, 
          a[[2]] - l2}] + (lmb[[2]] - 1) (lmb[[2]] - 2) CC[
         2, {lmb[[1]] - 0, lmb[[2]] - 2}, {a[[1]] - l1, a[[2]] - l2}]
      )]
   , {l1, Max[0, lmb[[3]] - 1 - a[[2]] - a[[3]]], 
    Min[a[[1]], lmb[[3]] - 1 - a[[3]]]}];
l3 =.; Clear[lPochh];
lPochh[a_, n_] := Pochhammer[a - n + 1, n];
CC[4, lmb_List, a_List] := 
  Sum[With[{l3 = lmb[[4]] - 1 - a[[4]] - l1 - l2},
    Binomial[a[[1]], l1] Binomial[a[[2]], l2] Binomial[a[[3]], 
      l3] (lmb[[4]] - 1)!
     (
        1 lPochh[lmb[[1]] - 1, 1] lPochh[lmb[[2]] - 1, 1] lPochh[
         lmb[[3]] - 1, 1] CC[3, Drop[lmb, -1] - {1, 1, 1}, 
         Drop[a, -1] - {l1, l2, l3}] + 
       1 lPochh[lmb[[1]] - 1, 0] lPochh[lmb[[2]] - 1, 2] lPochh[
         lmb[[3]] - 1, 1] CC[3, Drop[lmb, -1] - {0, 2, 1}, 
         Drop[a, -1] - {l1, l2, l3}] +
       1 lPochh[lmb[[1]] - 1, 1] lPochh[lmb[[2]] - 1, 0] lPochh[
         lmb[[3]] - 1, 2] CC[3, Drop[lmb, -1] - {1, 0, 2}, 
         Drop[a, -1] - {l1, l2, l3}] +
       2 lPochh[lmb[[1]] - 1, 0] lPochh[lmb[[2]] - 1, 1] lPochh[
         lmb[[3]] - 1, 2] CC[3, Drop[lmb, -1] - {0, 1, 2}, 
         Drop[a, -1] - {l1, l2, l3}] +
       1 lPochh[lmb[[1]] - 1, 0] lPochh[lmb[[2]] - 1, 0] lPochh[
         lmb[[3]] - 1, 3] CC[3, Drop[lmb, -1] - {0, 0, 3}, 
         Drop[a, -1] - {l1, l2, l3}]
      )]
   , {l1, 0, a[[1]]}, {l2, 
    Max[0, lmb[[4]] - 1 - a[[3]] - a[[4]] - l1], 
    Min[a[[2]], lmb[[4]] - 1 - a[[4]] - l1]}];
(*******)
NN = 2;
a = RandomInteger[{1, 5}, NN];
A = RandomInteger[{1, 5}, NN];
Ab = Accumulate[A];
NIntegrate[
 Product[l[xi]^a[[xi]], {xi, 1, NN}] (l[1] - 
    l[2]) Exp[-Sum[A[[xi]] l[xi], {xi, 1, NN}]], 
 Evaluate[Sequence @@ 
   Table[{l[xi], 0, If[xi == 1, Infinity, l[xi - 1]]}, {xi, 1, NN}]]]
Clear[lmb]; 
Table[CC[NN, {lmb[1], lmb[2]}, a] Product[ 1/Ab[[xi]]^
      lmb[xi], {xi, 1, NN}] /. 
    lmb[NN] :> Total[a] + 3 - lmb[1], {lmb[1], 2, a[[1]] + 3}] // 
  Total // N
(***********)
NN = 3;
a = RandomInteger[{1, 5}, NN];
A = RandomInteger[{1, 5}, NN];
Ab = Accumulate[A];
NIntegrate[
 Product[l[xi]^a[[xi]], {xi, 1, NN}] Product[
   l[i] - l[j], {i, 1, NN}, {j, i + 1, 
    NN}] Exp[-Sum[A[[xi]] l[xi], {xi, 1, NN}]], 
 Evaluate[Sequence @@ 
   Table[{l[xi], 0, If[xi == 1, Infinity, l[xi - 1]]}, {xi, 1, NN}]]]
Clear[lmb]; 
Table[CC[NN, {lmb[1], lmb[2] - lmb[1], Total[a] + 6 - lmb[2]}, 
       a] Product[ 1/Ab[[xi]]^(
       lmb[xi] - lmb[xi - 1]), {xi, 1, NN}] /. {lmb[0] :> 0, 
      lmb[NN] ->  Total[a] + 6}, {lmb[1], 2, Total[a] + 5}, {lmb[2], 
     lmb[1], Total[a] + 5}] // Flatten // Total // N
(*******)
NN = 4;
a = RandomInteger[{1, 5}, NN];
A = RandomInteger[{1, 5}, NN];
Ab = Accumulate[A];
NIntegrate[
 Product[l[xi]^a[[xi]], {xi, 1, NN}] Product[
   l[i] - l[j], {i, 1, NN}, {j, i + 1, 
    NN}] Exp[-Sum[A[[xi]] l[xi], {xi, 1, NN}]], 
 Evaluate[Sequence @@ 
   Table[{l[xi], 0, If[xi == 1, Infinity, l[xi - 1]]}, {xi, 1, NN}]]]
Clear[lmb]; 
Flatten[Table[
    CC[NN, {lmb[1], lmb[2] - lmb[1], lmb[3] - lmb[2], 
        Total[a] + 10 - lmb[3]}, a] Product[ 1/Ab[[xi]]^(
       lmb[xi] - lmb[xi - 1]), {xi, 1, NN}] /. {lmb[0] :> 0, 
      lmb[NN] ->  Total[a] + 10}, {lmb[1], 2, Total[a] + 9}, {lmb[2], 
     lmb[1], Total[a] + 9}, {lmb[3], lmb[2], Total[a] + 9}]] // 
  Total // N

---

Now, after having said all this my questions would be the following. Firstly, can we write down a recurrence relation for the coefficients for an arbitrary values of $N \ge 2$. Secondly, how does the computational expense, needed for computing the coefficients in question, increase with $N$.

Answer 1

For arbitrary $N \ge 2$ the recurrence relation for the coefficients reads as follows:

\begin{eqnarray} &&\left.{\mathfrak C}^{(N)}_{\lambda_1,\cdots,\lambda_N} ( \vec{a} )= \sum\limits_{\begin{array}{lll} l_1+l_2+ \cdots +l_N &=& \lambda_N-1 \\ l_1&=&0,\cdots, a_1 \\ \vdots \\ l_{N-1}&=&0, \cdots, a_{N-1} \\ l_N &=& a_N\end{array}} \prod\limits_{\xi=1}^{N-1} \binom{a_\xi}{l_\xi} \cdot (\lambda_N-1)! \right. \\ && \left[ \right. \\ % &&\sum\limits_{\begin{array}{lll} d_1+\cdots + d_{N-1} &=& N-1 \\ d_1 &=& 0,1 \\ d_2 &=& 0,\cdots, 2\\ \vdots \\ d_{N-1} &=& 1,\cdots, N-1 \end{array}} % {\mathfrak N}^{(N-1)}_{d_1,\cdots,d_{N-1}} \cdot \prod\limits_{\xi=1}^{N-1} (\lambda_\xi-1)_{(d_\xi)} \cdot {\mathfrak C}^{(N-1)}_{\lambda_1-d_1,\cdots,\lambda_{N-1}-d_{N-1}} (\vec{a} - \vec{l}) % \\ % && \left. \right] \cdot 1_{N+a_1 \ge \lambda_1 \ge 2} \cdot \prod\limits_{\xi=2}^{N-1} 1_{\lambda_\xi \ge 2} \cdot 1_{\lambda_N \ge 1} \quad (1) \end{eqnarray}

subject to $\sum\limits_{\xi=1}^N \lambda_\xi = |\vec{a}| + \binom{N+1}{2}$.

In here the coefficients ${\mathfrak N}_{\cdots}^{(N-1)}$ are such that $\prod\limits_{i=1}^{N-1} \left( \partial_{{\bar A}_i} + \cdots + \partial_{{\bar A}_{N-1}} \right) = \sum\limits_{d_1=0}^1 \sum\limits_{d_2=0}^2 \cdots \sum\limits_{d_{N-2}=0}^{N-2} \sum\limits_{d_{N-1}=1}^{N-1} {\mathfrak N}^{(N-1)}_{d_1,\cdots,d_{N-1}} \partial^{(d_1)}_{{\bar A}_1} \partial^{(d_2)}_{{\bar A}_2}\cdots \partial^{(d_{N-1})}_{{\bar A}_{N-1}} $ where $d_1+\cdots+d_{N-1}=N-1$.

---

Having said all this how do we make sure that all this is correct? I will explain that in remainder of this answer. As I was saying the spectral moments of the sample correlation matrix, i.e. the quantities $m_p := 1/T \cdot Tr[ X \cdot X^T ]$ where $X$ is an $N\times T $ random matrix composed of iid entries Gaussian distributed with mean zero and variance one, are easily expressed through the integral in question. We define a normalization constant as ${\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 $\vec{e}_j = ( \delta_{\xi,j} )_{\xi=1}^N $ and ${\mathfrak a} := (T-N-1)/2 $ and now we are ready to write down the moments as follows:

\begin{equation} m_p = \frac{{\mathfrak P}_{N,T}}{N \cdot T^p} \cdot 2^{\binom{N+1}{2} + N {\mathfrak a} +p} \cdot \sum\limits_{j=1}^N {\mathfrak J}_N^{({\mathfrak a} \vec{1}_N + p \vec{e}_j )} (\vec{1}_N ) \quad (2) \end{equation}

where $p \in {\mathbb Z}$.

In particular the negative moments admit the following representation:

\begin{eqnarray} %\left( \begin{array}{lll} m_{-1} &=& \frac{T}{-N+T-1} \\ m_{-2} &=& \frac{(T-1) T^2}{(-N+T-3) (-N+T-1) (T-N)} \\ m_{-3} &=& \frac{(T-1) T^3 (N+T-1)}{(-N+T-5) (-N+T-3) (-N+T-1) (T-N) (-N+T+1)} \\ m_{-4} &=& \frac{(T-1) T^4 \left(-\left(2 N^2+8 N-7\right) T+(N-2) \left(-N^2-N+2\right)+2 (N-2) T^2+T^3\right)}{(-N+T-7) (-N+T-5) (-N+T-3) (-N+T-2) (-N+T-1) (T-N) (-N+T+1) (-N+T+2)} \\ \vdots \end{array}\quad (3) %\right) \end{eqnarray}

Now we are going to present a code snippet that uses equation $(2)$ along with the recurrence $(1)$ and confirms the result $(3)$. Here we took $N=5$, $p=-1$, $T= N+1-2 p+ 2 (0,1)$ and here we go:

(*Definitions *)
Clear[l]; Clear[lmb]; Clear[a];
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]);
(Precompute mll[[]])
mll =.; Clear[Op];
mll = Table[{Exponent[#, 
       Table[Op[xi], {xi, 1, mNN - 1}]], (# /. 
        Op[j_] :> 1)} & /@ (List @@ 
      Expand[Product[
        Sum[Op[j], {j, i, mNN - 1}], {i, 1, mNN - 1}]]), {mNN, 2, 5}];
CC[2, lmb_List, a_List] := 
  If[a[[1]] >= lmb[[1]] - 2 && lmb[[1]] >= 2, 
   Binomial[a[[1]], 
     lmb[[1]] - 2] (lmb[[1]] - 1)! (a[[1]] + a[[2]] + 2 - lmb[[1]])!, 
   0] ;
l2 =.;
CC[3, lmb_List, a_List] := Sum[With[{l2 = lmb[[3]] - 1 - a[[3]] - l1},
    Binomial[a[[1]], l1] Binomial[a[[2]], l2] (lmb[[3]] - 1)!
     (
         (lmb[[1]] - 1) ( lmb[[2]] - 1) CC[
         2, {lmb[[1]] - 1, lmb[[2]] - 1}, {a[[1]] - l1, 
          a[[2]] - l2}] + (lmb[[2]] - 1) (lmb[[2]] - 2) CC[
         2, {lmb[[1]] - 0, lmb[[2]] - 2}, {a[[1]] - l1, a[[2]] - l2}]
      )]
   , {l1, Max[0, lmb[[3]] - 1 - a[[2]] - a[[3]]], 
    Min[a[[1]], lmb[[3]] - 1 - a[[3]]]}];
l3 =.; Clear[lPochh];
lPochh[a_, n_] := Pochhammer[a - n + 1, n];
CC[4, lmb_List, a_List] := 
  Sum[With[{l3 = lmb[[4]] - 1 - a[[4]] - l1 - l2},
    Binomial[a[[1]], l1] Binomial[a[[2]], l2] Binomial[a[[3]], 
      l3] (lmb[[4]] - 1)!
     (
        1 lPochh[lmb[[1]] - 1, 1] lPochh[lmb[[2]] - 1, 1] lPochh[
         lmb[[3]] - 1, 1] CC[3, Drop[lmb, -1] - {1, 1, 1}, 
         Drop[a, -1] - {l1, l2, l3}] + 
       1 lPochh[lmb[[1]] - 1, 0] lPochh[lmb[[2]] - 1, 2] lPochh[
         lmb[[3]] - 1, 1] CC[3, Drop[lmb, -1] - {0, 2, 1}, 
         Drop[a, -1] - {l1, l2, l3}] +
       1 lPochh[lmb[[1]] - 1, 1] lPochh[lmb[[2]] - 1, 0] lPochh[
         lmb[[3]] - 1, 2] CC[3, Drop[lmb, -1] - {1, 0, 2}, 
         Drop[a, -1] - {l1, l2, l3}] +
       2 lPochh[lmb[[1]] - 1, 0] lPochh[lmb[[2]] - 1, 1] lPochh[
         lmb[[3]] - 1, 2] CC[3, Drop[lmb, -1] - {0, 1, 2}, 
         Drop[a, -1] - {l1, l2, l3}] +
       1 lPochh[lmb[[1]] - 1, 0] lPochh[lmb[[2]] - 1, 0] lPochh[
         lmb[[3]] - 1, 3] CC[3, Drop[lmb, -1] - {0, 0, 3}, 
         Drop[a, -1] - {l1, l2, l3}]
      )]
   , {l1, 0, a[[1]]}, {l2, 
    Max[0, lmb[[4]] - 1 - a[[3]] - a[[4]] - l1], 
    Min[a[[2]], lmb[[4]] - 1 - a[[4]] - l1]}];
CC[NN_Integer, lmb_List, a_List] := Sum[
    With[{lastl = 
       lmb[[NN]] - 1 - a[[NN]] - Sum[l[xi], {xi, 1, NN - 2}]},
     Product[ 
       Binomial[a[[xi]], If[xi == NN - 1, lastl, l[xi]]], {xi, 1, 
        NN - 1}] (lmb[[NN]] - 1)!
      Sum[
       If[
        And @@ Table[
          lmb[[xi]] >= mll[[NN - 1]][[which, 1, xi]], {xi, 1, NN - 1}],
        mll[[NN - 1]][[which, 2]] Product[
          lPochh[lmb[[xi]] - 1, mll[[NN - 1]][[which, 1, xi]]], {xi, 
           1, NN - 1}] CC[NN - 1, 
          Drop[lmb, -1] - mll[[NN - 1]][[which, 1]], 
          Drop[a, -1] - 
           Table[If[xi == NN - 1, lastl, l[xi]], {xi, 1, NN - 1}]], 0]
       , {which, 1, Length[mll[[NN - 1]]]}]
     ]
    , Evaluate[
     Sequence @@ 
      Join[Table[{l[xi], 0, a[[xi]]}, {xi, 1, NN - 3}], 
       With[{tS = Sum[l[xi], {xi, 1, NN - 3}]}, {{l[NN - 2], 
          Max[0, lmb[[NN]] - 1 - a[[NN - 1]] - a[[NN]] - tS], 
          Min[a[[NN - 2]], lmb[[NN]] - 1 - a[[NN]] - tS]}}]]]
    ] /; NN >= 5;
NN = 5; p = -1;
Do[(Loop over T values.)
  T = NN + 1 - 2 p + 2 xi;
  aa = (T - NN - 1)/2;
  S = 0; Clear[lmb];
  Do[
   a = Table[aa + If[xi == j, p, 0], {xi, 1, NN}];
   A = ConstantArray[1, NN];
   Ab = Accumulate[A];
   t0 = TimeUsed[];
   res = Flatten[
     Table[{Table[
        If[xi == NN, Total[a] + Binomial[NN + 1, 2], lmb[xi]] - 
         If[xi == 1, 0, lmb[xi - 1]], {xi, 1, NN}], 
       CC[NN, Table[
          If[xi == NN, Total[a] + Binomial[NN + 1, 2], lmb[xi]] - 
           If[xi == 1, 0, lmb[xi - 1]], {xi, 1, NN}], a] Product[ 1/
         Ab[[xi]]^(
         If[xi == NN, Total[a] + Binomial[NN + 1, 2], lmb[xi]] - 
          If[xi == 1, 0, lmb[xi - 1]]), {xi, 1, NN}]}, 
      Evaluate[
       Sequence @@ 
        Table[{lmb[xi], If[xi == 1, 2, lmb[xi - 1] + 2], 
          If[xi == 1, NN + a[[1]], 
           Total[a] + Binomial[NN + 1, 2] - 1]}, {xi, 1, NN - 1}]]], 
     NN - 2];
   mres = 
    Pfct[NN, T]/(NN T^p) 2^(Binomial[NN + 1, 2] + NN aa + p)
      Total[res[[All, 2]]];
   S += mres;
   Print["j=", j, 
    " Analytical calculation done ", {#, N[#]} & /@ {mres}, 
    " .Done in t0=", TimeUsed[] - t0, " secs."];
   , {j, 1, NN}];
  S1 = Which[
    p == -1, T/(T - NN - 1),
    p == -2, (T^2 (T - 1))/(
    Product[T - NN + j, {j, -0, 0}] Product[
      T - NN - 2 j - 1, {j, 0, 1}]),
    p == -3, (T^3 (T - 1))/( 
     Product[T - NN + j, {j, -1, 1}] Product[
       T - NN - 2 j - 1, {j, 1, 2}]) (T + NN - 1),
    p == -4, (T^4 (T - 1))/( 
     Product[T - NN + j, {j, -3, 2}] Product[
       T - NN - 2 j - 1, {j, 2, 
        3}]) ((-2 + NN) (2 - NN - NN^2) - (-7 + 8 NN + 2 NN^2) T + 
       2 (NN - 2) T^2 + T^3), 
    p == -5, (T^5 (T - 1))/(
     Product[T - NN + j, {j, -5, 3}] Product[
       T - NN - 2 j - 1, {j, 3, 4}]) (-NN - 4 + T) (NN - 1 + 
       T) (-(-3 + NN) (-2 + NN) (2 + NN) + (15 - 18 NN - 4 NN^2) T + 
       4 (NN - 1) T^2 + T^3)];
  Print["N,T = ", {NN, T}, " m[", p, "]=", {S, S1}];
  , {xi, 0, 1}];

As we can see the match is perfect. Unfortunately the calculations are pretty slow because for every $N,T,p$ we have to change the vector $\vec{a}$ in ${\mathfrak J}_N^{\vec{a}}(\vec{A})$. As an overall conclusion we can say that this method is best suited for the case $\vec{a}$ being fixed and $\vec{A}$ being arbitrary because in that case we just pre-compute the coefficients once and then use them for all possible $\vec{A}$'s.