How to estimate $\int_0^\tau f(X_t)\:{\rm d}t$ when $X$ is a diffusion process?

Question

Say we have Markov processes $\left(X^{(i)}_t\right)_{t\ge0}$ with lifetime $\tau_i$ such that $\left(\left(X^{(i)}_t\right)_{t\ge0},\tau_i\right)_{i\in\mathbb N}$ is independent and identically distributed. The lifetimes can be chosen such that (see p. $$\frac{\sum_{i=0}^{n-1}Z_i}{\sum_{i=0}^{n-1}\tau_i}\xrightarrow{n\to\infty}\pi f,\tag1$$ where $$Z_i:=\int_0^{\tau_i}f(X_t)\:{\rm d}t$$ and $\pi$ is a given probability distribution. That's very nice, but how do we make use of $(1)$ in practice? The "problem" is clearly the question of how we can compute $Z_i$. Each $\tau_i$ can be computed via Poisson thinning as introduced in Non-Uniform Random Variate Generation (p. 264). An algorithm for simulating $Z_i$ could this way be

1. Choose $x_0$
2. Set $j=0$ and $\tau=0$
3. Set $j=j+1$
4. Sample $\Delta\tau\sim\operatorname{Exp}(\lambda)$, where $\lambda$ is an upper bound of the killing rate associated to $\tau_i$
5. Compute a sample $x_j$ from $X^{(i)}_{\Delta\tau}$ given $X_0^{(i)}=x_{j-1}$
6. Set $\tau^j=\tau^{j-1}+\Delta\tau$
7. Compute the killing rate $c$ at that sample
8. Sample a uniform $u$ on $[0,1)$
9. If $u<1-c/\lambda$ go to 3.
10. Otherwise $\tau^j\sim\tau_i$

Now, this only gives a "skeleton" $(x_0,\tau^0,\ldots,x_j,\tau^j)$ of $X^{(i)}$ on $[0,\tau_i)$. We may now estimate $Z_i$ by $$\frac1{\tau^j}\sum_{k=0}^jf(x_j)\tag2,$$ but is this a sensible estimate? The algorithm above actually cares only about obtaining a sample from $\tau_i$ (if I got it right) and it is only by its construction that we obtain some intermediate steps $\tau^j$ before reaching the final sample. But for computing $Z_i$ properly, we should need a "finer" discretization of $[0,\tau_i)$. But even if we would compute at a finer scale, I guess we still get some error due to the discretization.

So, what should we do here? You can assume that $X^{(i)}$ is a $d$-dimensional diffusion $${\rm d}X_t=\frac{\sigma^2}2\nabla\ln p(X_t){\rm d}t+\sigma{\rm d}W_t\tag3.$$ In that case, step 5 above is also not so easy to implement. Should we use Euler-Maruyama discretization or try "exact sampling". See Exact Simulation of Diffusions