What is the best way of evaluating time-ordered integrals numerically? In time-dependent perturbation theories, one encounters multi-dimensional time-ordered integrals
$$
(-i)^n\int_0^tdt_1\int_0^{t_1} dt_2 \cdots \int_0^{t_{n-1}}dt_n f(t_1,t_2,\cdots,t_n)
$$
What is the best way of numerically evaluating such multi-dimensional time-ordered integrals? Monte Carlo integration is the first way that comes into my mind. How does Monte Carlo deal with the time-ordering?
 A: The time-ordered integral can be transformed to a normal multi-dimensional integral over a rectangular volume by changing variables. Setting $t=\beta$, and using the following change of variables:
$$
\begin{cases}
t_1 = y_1\\
t_2 = y_2\frac{y_1}{\beta}\\
t_3 = y_3\frac{y_2y_1}{\beta^2}\\
\vdots\\
t_n = y_n\frac{y_{n-1}\cdots y_1}{\beta^{n-1}}
\end{cases}
$$
and the corresponding Jacobian
$$
\left(\frac{y_1}{\beta}\right)^{n-1}\left(\frac{y_2}{\beta}\right)^{n-2}\cdots \left(\frac{y_{n-1}}{\beta}\right)^{1},
$$
the integral $(-i)^n\int_0^tdt_1\int_0^{t_1} dt_2 \cdots \int_0^{t_{n-1}}dt_n f(t_1,t_2,\cdots,t_n)$
becomes
$$
(-i)^n\int_0^\beta dy_1 \int_0^{\beta} dy_2 \cdots \int_0^{\beta} dy_n \, f(y_1, \frac{y_2y_1}{\beta},\cdots,\frac{y_ny_{n-1}\cdots y_1}{\beta^{n-1}}) \left(\frac{y_1}{\beta}\right)^{n-1}\left(\frac{y_2}{\beta}\right)^{n-2}\cdots \left(\frac{y_{n-1}}{\beta}\right)^{1}.
$$
It is then easy to use a regular Monte Carlo integrator to complete the integral. For example, one may use the vegas algorithm https://vegas.readthedocs.io/en/latest/tutorial.html .
