I am using Monte Carlo methods to calculate the following $$ f (t, \lambda) = 1 - {\mathbb E}_{t,\lambda} \left[ \int^T_t f (t, \lambda(s)) \Gamma (\lambda(s), f (s,\lambda(s))) d s \right] , $$ where $\{ \lambda (t) \}_{t \in [0, T]}$ is a one- or multi-dimensional positive-valued Markov process with an explicit SDE (e.g. CIR process), $f$ is scalar-valued, and $\Gamma (\lambda, f)$ is a nonlinear function in both $\lambda$ and $f$.
The terminal condition $f (T, \lambda) = 1$, $\forall \lambda \in \mathbb R^+$ is given. My task is to calculate $f (0, \lambda_0)$, where $\lambda_0$ is the initial value of the state process, i.e., $\lambda (0) = \lambda_0$.
The difficulty in the numerical analysis is that $f (t, \lambda)$ is defined by the expectation recursively. That is, it enters into the integral of the expectation formula itself.
Although I can still use the idea of Monte Carlo simulations for a one-dimensional state process $\{ \lambda (t) \}_{t\in[0,T]}$, when it comes to a two-dimensional state process the computation time and accuracy becomes an issue.
I wondered if there are any good algorithms and references to implement Monte Carlo simulations for such a recursive expectation formula.