This notebook contains an illustration of the use of Monte Carlo methods for numerical integration. The associated lecture slides on risk-engineering.org provide an introduction to the use of stochastic simulation methods. Numerical integration is often used to evaluate risk measures in the finance industry.
Task: resolve the integral
$\int_1^5 x^2 dx$
We start by resolving this integral using the standard analytical method, assisted by the SymPy symbolic mathematics library.
import sympy from sympy.interactive import printing printing.init_printing(use_latex='mathjax') x = sympy.Symbol('x') integral = sympy.integrate(x**2) integral.subs(x, 5) - integral.subs(x, 1)
float(integral.subs(x, 5) - integral.subs(x, 1))
The expression we wish to integrate is very simple here, so calculating its integral analytically is easy (even without resorting to Python’s symbolic mathematics package!). In many cases, however, analytical approaches to integration are not feasible:
In such situations, stochastic simulation (“Monte Carlo”) methods allow us to generate an approximation of the integral, simply by evaluating the expression a large number of times at randomly selected points in the input space and counting the proportion that are less than the integrand at that point. The larger the number of simulations we run, the better the approximation (and note that computers are very very fast today, so on simple problems the simulation will run in a blink of an eye!).
N = 100000 accum = 0 for i in range(N): x = numpy.random.uniform(1, 5) accum += x**2 measure = 5 - 1 measure * accum/float(N)
x = sympy.Symbol('x') d1 = sympy.integrate(x**2 + 4 * x * sympy.sin(x)) sol = d1.subs(x, 3) - d1.subs(x, 2) float(sol)
N = 100000 accum = 0 for i in range(N): x = numpy.random.uniform(2, 3) accum += x**2 + 4*x*numpy.sin(x) measure = 3 - 2 measure * accum/float(N)
Exercise: now undertake the same comparison of analytical and stochastic simulation methods to evaluate the integral
(Hint: the result should be approximately 1464.2)
Now we move to a 2D integral, $\int\int cos(x^4) + 3y^2 dx dy$ over $x ∈ [4,6]$ and $y ∈ [0, 1]$.
Monte Carlo integration methods become increasingly attractive when compared to other numerical integration methods (such as quadrature) as the number of dimensions increases.
x = sympy.Symbol('x') y = sympy.Symbol('y') d1 = sympy.integrate(sympy.cos(x**4) + 3 * y**2, x) d2 = sympy.integrate(d1.subs(x, 6) - d1.subs(x, 4), y) sol = d2.subs(y, 1) - d2.subs(y, 0) float(sol)
N = 100000 accum = 0 for i in range(N): x = numpy.random.uniform(4, 6) y = numpy.random.uniform(0, 1) accum += numpy.cos(x**4) + 3 * y * y measure = 2 * 1 measure * accum/float(N)