As derived in PRB 73, 035120 (2006) in appendix, the Green's function moments are,
\begin{eqnarray} G = \frac{M_0}{i\omega} + \frac{M_1}{(i\omega)^2} + \frac{M_2}{(i\omega)^3}+\cdots \end{eqnarray}and can be computed from
\begin{eqnarray} M_k = (-1)^{k+1} \left[\frac{d^k}{d\tau^k}G(0^+)-\frac{d^k}{d\tau^k}G(0^-)\right] \end{eqnarray}This follows from integration by parts of $G(i\omega)=\int_0^\beta G(\tau)e^{i\omega\tau}d\tau$.
We thus have
\begin{eqnarray} M_0 &=& G(O^-)-G(0^+)=\delta_{\alpha\beta}\\ M_1 &=& -\langle [H,\psi_\alpha]\psi^\dagger_\beta\rangle-\langle \psi^\dagger_\beta[H,\psi_\alpha]\rangle =-\langle \{[H,\psi_\alpha],\psi^\dagger_\beta\}\rangle\\ M_2 &=& \langle [H,[H,\psi_\alpha]]\psi^\dagger_\beta\rangle+\langle \psi^\dagger_\beta[H,[H,\psi_\alpha]]\rangle= \langle \{[H,\psi_\alpha],[\psi^\dagger_\beta,H]\}\rangle\\ M_3 &=& -\langle \{[H,[H,\psi_\alpha]],[\psi_\beta^\dagger,H]\}\rangle \end{eqnarray}where we took into account that $H$ commutes with $e^{-\beta H}$ and hence $H$ can be moved from the first to the last place in the brackets.
First we compute
\begin{eqnarray} [H,\psi_\alpha] = -\varepsilon_{\alpha i}\;\psi_i - U_{\alpha i j k}\;\psi^\dagger_i \psi_j \psi_k - V_{\alpha k}c_k\\ [\psi^\dagger_\alpha,H]=-\psi_i^\dagger\; \varepsilon_{i\alpha} - U_{\alpha i j k}\;\psi_k^\dagger \psi^\dagger_j \psi_i - c^\dagger_k V_{k\alpha} \end{eqnarray}where $H$ is the general Anderson impurity Hamiltonian.
$M_1$ is thus
\begin{eqnarray} {M_1}_{\alpha\beta} = \varepsilon_{\alpha\beta} + (U_{\alpha i j \beta}-U_{\alpha i\beta j}) n_{ij} \end{eqnarray}$M_2$ becomes
\begin{eqnarray} {M_2}_{\alpha\beta}=\varepsilon_{\alpha i}\varepsilon_{i\beta} + V_{\alpha k}V_{k\beta}+ (U_{\alpha i j k}-U_{\alpha i k j})\varepsilon_{k \beta}\langle\psi_i^\dagger \psi_j\rangle+\varepsilon_{\alpha k}(U_{\beta j i k}-U_{\beta j k i})\langle\psi_i^\dagger \psi_j\rangle \\ + (U_{\alpha i j k}-U_{\alpha i k j})(U_{\beta l p k}-U_{\beta l k p}) \langle \psi_p^\dagger \psi_i^\dagger \psi_j \psi_l\rangle +U_{\alpha k j l} U_{\beta k i p}\langle \psi_p^\dagger \psi_i^\dagger \psi_j \psi_l\rangle \\ +(U_{\alpha i j k}-U_{\alpha i k j})U_{\beta l j k} \langle \psi^\dagger_i \psi_l\rangle \end{eqnarray}To calculate the moments of the self-energy, we need to use
\begin{eqnarray} G = \frac{1}{i\omega-\varepsilon-\Sigma_\infty-\frac{S}{i\omega+\cdots}-\Delta}=\frac{1}{i\omega}+\frac{M_1}{(i\omega)^2}+\frac{M_2}{(i\omega)^3}+\cdots=\frac{1}{i\omega}+\frac{\varepsilon+\Sigma_{\infty}}{(i\omega)^2}+ \frac{S+(\varepsilon+\Sigma_\infty)^2+(\Delta*(i\omega))}{(i\omega)^3}+\cdots \end{eqnarray}We thus see
\begin{eqnarray} \Sigma_{\infty} &=& M_1-\varepsilon\\ S &=& M_2 - (M_1)^2 -V^2 \end{eqnarray}Finally, we obtain
\begin{eqnarray} \Sigma_{\alpha\beta}(\infty) = (U_{\alpha i j \beta}-U_{\alpha i \beta j})n_{ij}\\ \end{eqnarray}and
\begin{eqnarray} S_{\alpha \beta} = (U_{\alpha i j k}-U_{\alpha i k j})(U_{\beta l p k}-U_{\beta l k p}) \left(\langle \psi_p^\dagger \psi_i^\dagger \psi_j \psi_l\rangle - \langle \psi_p^\dagger\psi_l\rangle\langle \psi_i^\dagger \psi_j \rangle\right) \\ +U_{\alpha k j l} U_{\beta k i p}\langle \psi_p^\dagger \psi_i^\dagger \psi_j \psi_l\rangle \\ +(U_{\alpha i j k}-U_{\alpha i k j})U_{\beta l j k} \langle \psi^\dagger_i \psi_l\rangle \end{eqnarray}