<< DiscreteMath`Combinatorica`

Unprotect[Times, Power];

(Subscript[S, k])_ (Subscript[S, j])_  :=  
      Module[{n1 = Length[k], n2 = Length[j], k1, k2, p, q, x, y, res = 0, 
          N = 0}, 
{k1, k2, n1, n2} = 
          If[n1 < n2, {j, k, n2, n1}, {k, j, n1, n2}]; 
p = 
          Sort[Table[
              NthSubset[i, k2], {i, 0, (2^n2) - 1}]]; 
q = 
          Reverse[p]; 
Do[
          x = Permutations[
              PadRight[p[[i]], n1, 
                N]]; 
y = Length[p[[i]]]!/
              Length[Permutations[p[[i]]]]; 
Do[
            res = res + 
                y (Subscript[S, Sort[Join[k1 + x[[m]], q[[i]]], Greater]]), {m, 1, Length[x]}], {i, 1, (2^n2)}];  res]

((((Subscript[S, k])_))^m)_  :=  
      Module[{res = (Subscript[S, k])}, Do[res = Expand[(Subscript[S, k]) res], {j, m - 1}]; 
        Expand[res]]

\(Expand[\((S\_{2, 1} + S\_{1})\) S\_{3, 1}]\)

(Subscript[S, k])_ (Subscript[S, j])_  :=  
      Module[{n1 = Length[k], n2 = Length[j], k1, k2, p, q, x, y, res = 0, 
          N = k[[1]] 0}, 
{k1, k2, n1, n2} = 
          If[n1 < n2, {j, k, n2, n1}, {k, j, n1, n2}]; 
p = 
          Sort[Table[
              NthSubset[i, k2], {i, 0, (2^n2) - 1}]]; 
q = 
          Reverse[p]; 
Do[
          x = Permutations[
              PadRight[p[[i]], n1, 
                N]]; 
y = Length[p[[i]]]!/
              Length[Permutations[p[[i]]]]; 
Do[
            res = res + 
                y (Subscript[S, Sort[Join[k1 + x[[m]], q[[i]]], lexi]]), {m, 
              1, Length[x]}], {i, 1, (2^n2)}];  res]

lexi[A_, B_] := 
      Module[{N = Length[A]}, 
Which[
          N == 0 || A[[1]] > B[[1]], 
          True, 
A[[1]] < B[[1]], 
          False, 
A[[1]] == B[[1]], 
          lexi[Drop[A, 1], Drop[B, 1]]]]

\(Expand[\(S\_{{2, 0}, {1, 0}}\) S\_{{2, 1}}^2]\)

\(Expand[S\_{{1, 1, 1}}^3]\)

EV[A_] := FixedPoint[Expand, A] /. (Subscript[S, k])_ -> esip[k] /. 
            esip -> ESIP // PowerExpand // Together

ESIP[k_] := 
      Module[{N = Length[k], res}, res = Binomial[n, N] N!; 
        Do[res = res MOMENT[k[[i]]], {i, N}]; res]

MOMENT[k_] := 
      D[Exp[λ (Exp[t] - 1)], {t, k[[1]]}] /. t -> 0 // 
        Simplify

\(EV[\((S\_{{1}}/n - λ)\)\^3]\)

MOMENT[k_] := 
      D[Exp[(Subscript[μ, 1]) (Subscript[t, 1]) + (Subscript[μ, 2]) 
                  (Subscript[t, 2]) + (((Subscript[σ, 1])^2) ((Subscript[t, 1])^2) + ((Subscript[σ, 2])^2) 
                        ((Subscript[t, 2])^2) + 
                      2  ρ (Subscript[σ, 1]) (Subscript[σ, 2]) (Subscript[t, 1]) 
                        (Subscript[t, 2]))/2], {(Subscript[t, 1]), k[[1]]}, {(Subscript[t, 2]), 
              k[[2]]}] /. {(Subscript[t, 1]) -> 0, (Subscript[t, 2]) -> 0} // 
        Simplify

\(EV[\((S\_{{1, 0}}/n - μ\_1)\) \((S\_{{0, 1}}/
              n - μ\_2)\)]\)

MOMENT[k_] := 
      D[p/(Exp[-t] - 1 + p), {t, k[[1]]}] /. t -> 0 // 
        Simplify

\(\((Series[
            EV[Series[\((n - 1)\)/\((ϵ \((S\_{{1}} - n/p)\) + n/p - 
                          1)\), {ϵ, 0, 6}] // 
                  Normal] /. {ϵ -> 1, n -> 1/v}, {v, 0, 
              3}] // Simplify)\) /. v -> 1/n\)

\(\((Series[
            EV[Series[\((\((n - 1)\)/\((ϵ \((S\_{{1}} - n/p)\) + n/p - 1)\) - p)\)\^2, {ϵ, 0, 6}] // Normal] /. {ϵ -> 1, 
                n -> 1/v}, {v, 0, 3}] // Simplify)\) /. 
      v -> 1/n\)

\(\((Series[
            EV[Series[\((\((n - 1)\)/\((ϵ \((S\_{{1}} - n/p)\) + n/p - 1)\) - p)\)\^3, {ϵ, 0, 6}] // Normal] /. {ϵ -> 1, 
                n -> 1/v}, {v, 0, 3}] // Simplify)\) /. 
      v -> 1/n\)

\(\((Series[
            EV[Series[\((\((n - 1)\)/\((ϵ \((S\_{{1}} - n/p)\) + n/p - 1)\) - p)\)\^4, {ϵ, 0, 6}] // Normal] /. {ϵ -> 1, 
                n -> 1/v}, {v, 0, 3}] // Simplify)\) /. 
      v -> 1/n\)

MOMENT[k_] := 
      D[Exp[(((Subscript[t, 1])^2) + ((Subscript[t, 2])^2) + 2  ρ  (Subscript[t, 1]) (Subscript[t, 2]))/2], {(Subscript[t, 1]), 
              k[[1]]}, {(Subscript[t, 2]), k[[2]]}] /. {(Subscript[t, 1]) -> 0, 
            (Subscript[t, 2]) -> 0} // Simplify

\(\(r = \((ϵ \((\ 
                  S\_{{1, 1}}/
                      n - ρ)\) + ρ - \(ϵ\^2\) \(S\_{{0, 
                        1}}\) S\_{{1, 0}}/
                  n\^2)\)/\@\(\((ϵ \((S\_{{2, 0}}/n - 1)\) + 1 - \(ϵ\^2\) S\_{{1, 0}}\^2/n\^2)\) \((ϵ \((S\_{{0, 2}}/n - 1)\) + 1 - \(ϵ\^2\) S\_{{0, 1}}\^2/n\^2)\)\);\)\)

μ = 
      Collect[Collect[Series[r, {ϵ, 0, 6}] // Normal, ϵ, 
              EV] /. ϵ -> 1, n, Factor] /. 
        (n^(i_ /; i < -3)) -> 0

(Subscript[μ, 2]) = 
        Collect[Collect[
                Series[(r^2), {ϵ, 0, 6}] // Normal, ϵ, 
                EV] /. ϵ -> 1, n, Factor] /. 
          (n^(i_ /; i < -3)) -> 0;

Var = 
      Collect[(Subscript[μ, 2]) - (μ^2), n, Factor] /. 
        (n^(i_ /; i < -3)) -> 0

(Subscript[μ, 3]) = 
        Collect[Collect[
                Series[(r^3), {ϵ, 0, 6}] // Normal, ϵ, 
                EV] /. ϵ -> 1, n, Factor] /. 
          (n^(i_ /; i < -3)) -> 0;

Collect[
        Series[(((Subscript[μ, 3]) - 3  (Subscript[μ, 2]) μ + 2  (μ^3))/(Var^(3/2))) /. n -> 1/w, {w, 0, 2}] // Normal, w, Simplify[#, (ρ^2) < 1] &] /. 
      w -> 1/n

(Subscript[μ, 4]) = 
        Collect[Collect[
                Series[(r^4), {ϵ, 0, 6}] // Normal, ϵ, 
                EV] /. ϵ -> 1, n, Factor] /. 
          (n^(i_ /; i < -3)) -> 0;

Collect[
        Series[(((Subscript[μ, 4]) - 4  (Subscript[μ, 3]) μ + 6  (Subscript[μ, 2]) (μ^2) - 3  (μ^4))/(Var^2)) /. n -> 1/w, {w, 0, 1}] // Normal, w, 
        Simplify[#, (ρ^2) < 1] &] /. w -> 1/n

μ = 
      Collect[Collect[
              Series[ArcTanh[r], {ϵ, 0, 6}] // Normal, ϵ, 
              EV] /. ϵ -> 1, n, Factor] /. 
        (n^(i_ /; i < -3)) -> 0

(Subscript[μ, 2]) = 
        Collect[Collect[
                Series[ArcTanh[r]^2, {ϵ, 0, 6}] // 
                  Normal, ϵ, EV] /. ϵ -> 1, n, 
            Factor] /. (n^(i_ /; i < -3)) -> 0;

Var = 
      Collect[(Subscript[μ, 2]) - (μ^2), n, Factor] /. 
        (n^(i_ /; i < -3)) -> 0

(Subscript[μ, 3]) = 
        Collect[Collect[
                Series[ArcTanh[r]^3, {ϵ, 0, 6}] // 
                  Normal, ϵ, EV] /. ϵ -> 1, n, 
            Factor] /. (n^(i_ /; i < -3)) -> 0;

Collect[
        Series[(((Subscript[μ, 3]) - 3  (Subscript[μ, 2]) μ + 2  (μ^3))/(Var^(3/2))) /. n -> 1/w, {w, 0, 2}] // Normal, w, Simplify[#, (ρ^2) < 1] &] /. 
      w -> 1/n

(Subscript[μ, 4]) = 
        Collect[Collect[
                Series[ArcTanh[r]^4, {ϵ, 0, 6}] // 
                  Normal, ϵ, EV] /. ϵ -> 1, n, 
            Factor] /. (n^(i_ /; i < -3)) -> 0;

Collect[
        Series[(((Subscript[μ, 4]) - 4  (Subscript[μ, 3]) μ + 6  (Subscript[μ, 2]) (μ^2) - 3  (μ^4))/(Var^2)) /. n -> 1/w, {w, 0, 1}] // Normal, w, 
        Simplify[#, (ρ^2) < 1] &] /. w -> 1/n