Often it is necessary to compute discrete-time approximation of continuous-time stochastic processes and to resort to Monte Carlo replications to have an approximation of the distribution of certain functionals of the processes. I better explain my point with a simple example. Consider the standard Ito semi-martingale model \begin{eqnarray} dY_t &=& \mu_t\,dt+\sigma_t\,dW_t\\ d\sigma_t &=& \alpha_t\,dt+\beta_t\,dW_t+\gamma_t\,dW^{\prime}_t\quad (1) \end{eqnarray} where $W_t$ and $W^{\prime}_t$ are independent Brownian motion and $\mu$, $\alpha$, $\beta$ and $\gamma$ are process satisfying some regularity conditions.
Suppose that I want to compute $$ \mathbb{E}\left[\sum_{j=1}^{n}\left(Y_{t_{j,n}}-Y_{t_{j-1,n}}\right)^2\right]\quad (2) $$ where $0=t_{0,n}<t_{1,n}<t_{2,n}<\cdots<t_{n,n}=1$ is a discrete partition of the time interval $\left[0,1\right]$ over which the process $Y$ is assumed to be observed. To do that, one typically generates $R$ discrete-time approximated replica of the system of equations in $(1)$ via an Euler-like scheme
\begin{eqnarray} Y_{j/n}\left(\omega_r\right)-Y_{(j-1)/n}\left(\omega_r\right) &=& \mu_0\,\Delta_n+\sigma_{(j-1)/n}\left(\omega_r\right)\,\varepsilon_{j}\left(\omega_r\right)\,n^{-1/2}\\ \sigma_{j/n}\left(\omega_r\right)-\sigma_{(j-1)/n}\left(\omega_r\right) &=& \alpha_0\,\Delta_n+\beta_0\,\varepsilon_{j}\left(\omega_r\right)\,n^{-1/2}+\gamma_0\,\varepsilon_{j}^{\prime}\left(\omega_r\right)\,n^{-1/2} \end{eqnarray} having assumed, for simplicity, that $\mu$, $\alpha$, $beta$ and $\gamma$ are constant, that the partition is equi-spaced (with spacing $\Delta_n=1/n$) and where $\varepsilon_j$ and $\varepsilon_{j}^{\prime}$ are both iid $\textrm{N}\left(0,1\right)$ and where, with $\left(\omega_r\right)$ I am indicating the $r-$th replication of the model, with $r=1,...,R$. The moment in equation (2) is then approximated assuming that, as $R\rightarrow +\infty$, the following LNN holds (in some sense) $$ \frac{1}{R}\sum_{j=1}^R\left(Y_{j/n}\left(\omega_r\right)-Y_{(j-1)/n}\left(\omega_r\right)\right)^2\longrightarrow \mathbb{E}\left[\sum_{j=1}^{n}\left(Y_{t_{j,n}}-Y_{t_{j-1,n}}\right)^2\right] $$ This approach is, however, quite heuristics, since the computation of the $R$ replica of the discrete time approximation requires to summon some pseudo-random number generator. My question is: is there any reference textbook or any contribution in the literature that deals with a formalization of Monte Carlo simulations and pseudo-random number generations?