**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](http://robotics.caltech.edu/%7Ezoran/Research/poisson.pdf), 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?



General result for the N-point correlation of the Poisson process (and its derivative)?
This question (and its generalization) is conventiently addressed by considering the moment-generating functional. Let $N\_t$ be a counting process with (possibly stochastic and time-dependent) intensity $\lambda\_t$ (you recover your case by setting $\lambda\_t=\lambda$, where $\lambda$ is a constant).


Suppose that under $\mathbb{P}$, the compensated counting process $N\_t-\int\_0^t\lambda\_s ds$ is independent of $\mathcal{F}^{\lambda}\_t$, the natural filtration of $\lambda\_t$.
Then we have



> 
> $$Z[a]\doteq \mathbb{E}\Bigg[\exp\Big[\int\_0^Ta(t)dN\_t\Big]\Bigg]=\mathbb{E}\Bigg[\exp\Big[\int\_0^T\big(e^{a(t)}-1\big)\lambda\_t dt\Big]\Bigg], \quad (\star)$$
> 
> 
> 


for e.g. continuous $a$ for which the RHS exists.
Indeed, under the above assumption there is an equivalent measure $\mathbb{P}^{\ast}$ under which $N\_t$ has unit intensity and the law of $(\lambda\_t)\_{t\geq 0}$ is the same as under the original measure $\mathbb{P}$. By using a Girsanov-type change-of-measure formula (see e.g. Theorems 26.8 and 26.9 in *Kallenberg, O. (2006). Foundations of Modern Probability.*) for the measure $\mathbb{P}^{\ast}$, we obtain 


$$
\begin{split}
\mathbb{E}\Bigg[\exp\Big[\int\_0^Ta(t)dN\_t\Big]\Bigg]&=\mathbb{E}^{\ast}\Bigg[\frac{d\mathbb{P}}{d\mathbb{P}^{\ast}}\exp\Big[\int\_0^Ta(t)dN\_t\Big]\Bigg] \\
&=\mathbb{E}^{\ast}\Bigg[\exp\Big[\int\_0^T\log \lambda\_tdN\_t+\int\_0^T\big(1-\lambda\_t\big)dt\Big]\exp\Big[\int\_0^Ta(t)dN\_t\Big]\Bigg]\\
&=\mathbb{E}^{\ast}\Bigg[\exp\Big[\int\_0^T\log \big( e^{a(t)}\lambda\_t\big)dN\_t+\int\_0^T\big(1-e^{a(t)}\lambda\_t+(e^{a(t)}-1)\lambda\_t\big)dt\Big]\Bigg] \\
&=\mathbb{E}^{\ast}\Bigg[\exp\Big[\int\_0^T\log \big( e^{a(t)}\lambda\_t\big)dN\_t+\int\_0^T\big(1-e^{a(t)}\lambda\_t\big)dt\Big] \\ &\qquad\quad\times\exp\Big[\int\_0^T(e^{a(t)}-1)\lambda\_tdt\Big] \Bigg]\\
&=\mathbb{E}'\Bigg[\exp\Big[\int\_0^T\big(e^{a(t)}-1\big)\lambda\_t dt\Big]\Bigg].
\end{split}$$
In the last line, the measure $\mathbb{P}'$ is such that $N\_t$ has rate $a(t)\lambda\_t$, but the law of $(\lambda\_t)\_{t\geq 0}$ is the same as under both $\mathbb{P}$ and $\mathbb{P}^{\ast}$. Since the integrand only depends on $(\lambda\_t)\_{t\geq 0}$, $(\star)$ follows. 




---


From $(\star)$, we may extract arbitrary correlation functions using functional derivatives. We have 


$$\frac{\delta Z[a]}{\delta a(t)}=\mathbb{E}\Bigg[e^{a(t)}\lambda\_t\exp\Big[\int\_0^T\big(e^{a(s)}-1\big)\lambda\_s ds\Big]\Bigg],$$


on which we can iteratively apply additional functional derivatives. E.g.


$$\mathbb{E}[S\_t]=\frac{\delta Z[a]}{\delta a(t)}\Bigg|\_{a=0}=\mathbb{E}[\lambda\_t],$$


$$\mathbb{E}[S\_tS\_{t'}]=\frac{\delta^2 Z[a]}{\delta a(t)\delta a(t')}\Bigg|\_{a=0}=\delta(t-t')\mathbb{E}[\lambda\_t]+\mathbb{E}[\lambda\_t\lambda\_{t'}].$$


By specializing to the case of fixed deterministic rate, you can check your guesses for higher-order moments of a Poisson pulse train.


EDIT: I had to add an assumption for $(\star)$ to hold.


