How can I numerically compute a stochastic integral?

I am trying simulate a to solve a 2-dimensional stochastic process and $$Y_t^1$$ is a mean-reverting square root process which I simulated on a time grid using its known conditional distribution. I want to compute the stochastic process $$Y_t^2$$ given below . $$\begin{eqnarray*} Y_t^2:=e^{\tilde{a}_{22}t}\tilde{s}_2+\int_0^te^{\tilde{a}_{22}(t-s)}\left[\left(\tilde{a}_{21}Y_s^1+ \tilde{b}_2\right)ds+\sqrt{ \alpha+\vert \beta\vert Y_s^1}d\tilde{W}^2_s\right] \end{eqnarray*}$$

That means that I want to compute the values of $$Y_t^2$$ along 100 paths at the time discretization grid (we know that the time discretization grid is $$[0 ,0.01, 0.02 ,\dots,0.99]$$. I want to compute the a $$100 \times 100$$ matrix where the columns vector represent the paths. $$\tilde{a}_{21}, \tilde{a}_{22},\alpha, \vert \beta\vert,\tilde{b}_2,\tilde{s}_2$$ are known constants and $$Y_t^1$$ is known (in other words $$100 \times 100$$ matrix with each column representing the paths of $$Y_t^1$$) is known.

Attempt: I tried to approximate the integral by using the following approximation with

$$t_1=0,t_2=0.01,t_3=0.02,\dots,t_{100}=0.99$$

$$Y_{t_k}^2(\omega)\approx e^{\tilde{a}_{22}t_k}\tilde{s}_2+\tilde{a}_{21} Y_0^1(\omega)\int_0^{t_2} e^{\tilde{a}_{22}(t_k-s)}ds+\tilde{b}_2\int_0^{t_2}e^{\tilde{a}_{22}(t_k-s)}ds+\sqrt{\alpha+\vert \beta \vert Y_0^1(\omega)}\int_0^{t_2}e^{\tilde{a}_{22}(t_k-s)}d\tilde{W}^2_s+\dots +\tilde{a}_{21} Y_{t_{k-1}}^1(\omega)\int_{t_{k-1}}^{t_k} e^{\tilde{a}_{22}(t_k-s)}ds+\tilde{b}_2\int_{t_{k-1}}^{t_k} e^{\tilde{a}_{22}(t_k-s)}ds+\sqrt{\alpha+\vert \beta \vert Y_{t_{k-1}}^1(\omega)}\int_{t_{k-1}}^{t_k}e^{\tilde{a}_{22}(t_k-s)}d\tilde{W}^2_s$$.

I compute the stochastic integral by assuming the integrand takes value at the left-endpoint of the interval which is $$t_1=0$$ in the case below and integrand comes out of the integral and we are just left with stochastic integral the Brownian motion which is simulated by simulating normally distributed random variables of mean zero and variance equal to the time step(in this case $$0.01$$)

$$\sqrt{\alpha+\vert \beta \vert Y_0^1(\omega)}\int_0^{t_2}e^{\tilde{a}_{22}(t_k-s)}d\tilde{W}^2_s=\sqrt{\alpha+\vert \beta \vert Y_0^1(\omega)}e^{\tilde{a}_{22}t_k}\int_0^{t_2}d\tilde{W}^2_s\\=\sqrt{\alpha+\vert \beta \vert Y_0^1(\omega)}e^{\tilde{a}_{22}t_k} *N(0,0.01) \text{sample}$$

I compute the lebesgue integral as follows

$$\tilde{a}_{21} Y_{t_{k-1}}^1(\omega)\int_{t_{k-1}}^{t_k} e^{\tilde{a}_{22}(t_k-s)}ds= Y_{t_{k-1}}^1(\omega)\frac{\tilde{a}_{21}}{\tilde{a}_{22}}\left(e^{\tilde{a}_{22}(t_k-t_{k-1})}-e^{\tilde{a}_{22}(t_k-t_k)}\right)$$

• I don’t have time to quote the entire idea right now but in Diffusions, Markov Processes, and Martingales volume 2 by L. C. G. Rogers and D. Williams, section 47 of chapter IV, there is the topic “Riemann-sum approximation to Itô and Stratonovich integrals; simulation” which may be of some help if you’re not already familiar. – Nap D. Lover May 1 at 1:37
• Thank you . I am going to have a look at it . But I am afraid this is more of a numerical issue – user3503589 May 1 at 9:46