First note, I had asked a similar question here, but the thread seems to have died, so I'll revive it here with more details. As a simplification of my real problem, I want to compute $\mathbb{E}[\int_o^T N_{t-}dS_t]$, where $S_t$ is a $(\mu,\sigma)$-geom. Brownian motion, and $N_t$ an independent Poisson process with const. intensity $\lambda$. Hence $\mathbb{E}[\int_o^T N_{t-}dS_t] = \mathbb{E}[\int_o^T N_{t-}S_t(\mu dt + \sigma dW_t)]$. Ito formula shows that the latter part of the integral should be square integrable, with integrable quadratic variation, so it is zero in expectation. So we should get

$\mathbb{E}[\int_o^T N_{t-}dS_t] = \mu\int_0^T\underbrace{\mathbb{E}[N_{t-}]}_{=\mathbb{E}[N_t-\mathbb{1}_{\Delta N_t \neq 0}]=\lambda t}\mathbb{E}[S_t]dt = \int_0^T \mu\lambda t e^{\mu t}dt$

$N_{t-}$ is a simple, left-continuous process, so I tried to confirm that numerically by calculating (MC)

$\mathbb{E}\left[ \sum_{i=0}^{N_t} i(S_{\tau_{i+1}}-S_{\tau_i})\right]$,

with $\tau_0=0$, $\tau_{N_t+1}=T$ and $\tau_i$ the exp-$(\lambda)$-distributed jump times. But values don't seem to converge to the analytic solution above. Is anything wrong?