**I am specifically interested in computing:
$$\mathbb{E}[S_p S_q S_s S_t]$$**
where $S_t=\frac{dN_t}{dt}$ and $N_t$ is a Poisson process (so $S_t$ is a "Poisson pulse train"):
$$\mathbb{P}(N_t=n)=\frac{(\lambda t)^n}{n!}e^{-\lambda t} \; , \; t>0$$

**Here is my attempt:**

Assuming $t>s>p>q$, I tried to compute it from the process $N_t$ rewriting it with the independent random variables $N_t-N_s, N_s-N_p, N_p-N_q, N_p$, that we can rename $N_j$, $j=1,2,3,4$: $$\mathbb{E}[N_p N_q N_s N_t]=\mathbb{E}[\prod_{k=1}^4\sum_{j=1}^k N_j]$$

This unfolds a combinatorics which somewhat looks like the one we obtain for normally distributed variables using the Isserlis / Wick theorem.

Is there an analogous theorem about such correlation function which applies to the Poisson process?

This looks simpler for the pulse train since the only term which does not vanish when we take the derivative is: $$\mathbb{E}[S_p S_q S_s S_t]=\frac{d^4}{dpdqdsdt}\mathbb{E}[N_p N_q N_s N_t]=\frac{d^4}{dpdqdsdt}\big(\mathbb{E}[N_1]\mathbb{E}[N_2]\mathbb{E}[N_3]\mathbb{E}[N_4]\Big)=\lambda^4$$ However, this is without considering coincidences among $t,s,p,q$. From here, I see (p.10): $$\mathbb{E}[S_s S_t]=\lambda^2+\lambda \delta(s-t)$$

I deduce that we should have something like $$\mathbb{E}[S_p S_q S_s S_t]=\lambda^4+\lambda^3\sum_{i\neq j}\delta(t_i-t_j)+\lambda^2\sum_{i\neq j,\; k\neq l }\delta(t_i-t_j)\delta(t_k-t_l)+\lambda\sum_{i\neq j,\; k\neq l,\; m\neq n }\delta(t_i-t_j)\delta(t_k-t_l)\delta(t_m-t_n)$$

Is this correct (except the sloppily specified summation indices)? Do you know any related general result?