Monte Carlo Methods: Computing via Random Sampling¶

Monte Carlo (MC) methods are a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results. The underlying concept is to use randomness to solve problems that might be deterministic in principle, such as multidimensional integration or optimization.

1. Core Mathematical Principle¶

The foundation of Monte Carlo methods is the Law of Large Numbers. If we want to find the expected value of a function \(f(x)\) where \(x\) is a random variable with probability density \(p(x)\), we can approximate it by taking the empirical mean of \(N\) independent samples \(x_1, x_2, \dots, x_N\):

As \(N \to \infty\), the sample mean converges to the true expected value.

1.1 Error Convergence¶

Unlike many numerical integration methods (like the Trapezoidal rule), the error in Monte Carlo methods is independent of the dimensionality of the problem. The standard error decreases at a rate of:

This makes MC methods exceptionally powerful for high-dimensional physics and financial modeling.

Examples¶

Classic Example: Estimating \(\pi\)¶

A common way to visualize the Monte Carlo method is by estimating the value of \(\pi\) using a square and an inscribed circle.

1. Consider a unit square with area \(A_{square} = 4\) (from \(x,y \in [-1, 1]\)).
2. Inscribe a unit circle with area \(A_{circle} = \pi(1)^2 = \pi\).
3. The ratio of the areas is:

1. By randomly "dropping" points in the square and counting how many land inside the circle, we can estimate \(\pi\):

The following code implements the \(\pi\) estimation algorithm and a basic Monte Carlo integration.

import numpy as np

def estimate_pi(n_samples=10000):
    """
    Estimates Pi using random sampling in a 2D plane.
    """
    # Generate random (x, y) coordinates between -1 and 1
    points = np.random.uniform(-1, 1, (n_samples, 2))

    # Calculate distance from origin (x^2 + y^2)
    distances = np.sum(points**2, axis=1)

    # Check if points are inside the unit circle (distance <= 1)
    inside_circle = np.sum(distances <= 1)

    pi_estimate = 4 * (inside_circle / n_samples)
    return pi_estimate

def mc_integrate(func, a, b, n_samples=100000):
    """
    Performs 1D Monte Carlo integration of func from a to b.
    """
    # Sample x values uniformly in [a, b]
    x = np.random.uniform(a, b, n_samples)

    # Compute the average value of the function
    f_mean = np.mean(func(x))

    # Integral = (b - a) * E[f(x)]
    return (b - a) * f_mean

# Example: Integrate x^2 from 0 to 1 (Exact answer is 1/3)
result = mc_integrate(lambda x: x**2, 0, 1)
print(f"Estimated Pi: {estimate_pi()}")
print(f"Estimated Integral of x^2: {result:.5f}")

Importance Sampling¶

When the function \(f(x)\) has high variance or most of its "mass" is in a small region, uniform sampling is inefficient. Importance Sampling uses a different distribution \(q(x)\) that focuses on the important regions:

Where \(w_i = \frac{p(x_i)}{q(x_i)}\) is known as the likelihood ratio or importance weight.

In this example, we'll integrate a Gaussian-like function \(f(x)\) that is very "skinny."

import numpy as np
import matplotlib.pyplot as plt

# 1. Define the function we want to integrate: f(x) * p(x)
# Let's say p(x) is uniform [0, 5], but f(x) has a sharp peak at x=4
def target_function(x):
    return np.exp(-(x - 4)**2 / (2 * 0.1**2))

# Parameters
a, b = 0, 5  # Integration limits
n_samples = 1000

# --- METHOD 1: STANDARD MONTE CARLO (Uniform Sampling) ---
x_uniform = np.random.uniform(a, b, n_samples)
f_values_uniform = target_function(x_uniform)
# Estimate = (b-a) * average(f(x))
mc_estimate = (b - a) * np.mean(f_values_uniform)

# --- METHOD 2: IMPORTANCE SAMPLING ---
# We choose a proposal distribution q(x) that "looks like" our target.
# Let's use a Normal distribution centered at 4.
mu_q, sigma_q = 4, 0.5
x_importance = np.random.normal(mu_q, sigma_q, n_samples)

# We must calculate the weights: w = p(x) / q(x)
# p(x) is our original uniform density: 1 / (b-a)
p_x = 1.0 / (b - a)

# q(x) is the Normal PDF we sampled from
from scipy.stats import norm
q_x = norm.pdf(x_importance, mu_q, sigma_q)

weights = p_x / q_x
is_estimate = np.mean(target_function(x_importance) * weights) * (b - a)

print(f"Standard MC Estimate: {mc_estimate:.6f}")
print(f"Importance Sampling Estimate: {is_estimate:.6f}")
print(f"Analytical (Approx) Integral: {0.25066:.6f}") # Constant for this specific bell curve

Value at Risk (VaR): Measuring Portfolio Risk with Monte Carlo¶

Value at Risk (VaR) is the financial industry's primary metric for quantifying potential losses. It answers a single, critical question: "What is the maximum amount I could lose over a specific time period with a given level of confidence?"

While traditional methods assume markets follow a simple, predictable path, Monte Carlo Simulation is the preferred approach for complex portfolios because it models thousands of random market scenarios to capture extreme risks.

How the Simulation Works¶

The process transforms uncertain market variables into a clear probability distribution:

1. Define the Engine: Analysts use Geometric Brownian Motion (GBM) to model price paths. This formula combines the "Drift" (expected return) with "Volatility" (random market shocks).
2. Generate Scenarios: The computer simulates 10,000+ "possible futures" for the portfolio. Each path represents a unique combination of interest rate shifts, stock price movements, and asset correlations.
3. Calculate Outcomes: For every simulated path, the final portfolio value is recorded at the end of the time horizon (e.g., 1 day or 30 days).
4. Identify the Cut-off: The results are sorted from worst to best. If calculating a 95% VaR, the analyst identifies the 5th percentile—the threshold where only 5% of outcomes resulted in a larger loss.

Key Advantages¶

• Captures "Fat Tails": Unlike standard models, Monte Carlo can simulate rare, high-impact market crashes (Black Swan events).
• Non-Linear Risk: It accurately prices complex instruments like options and derivatives, where risk changes rapidly as prices move.
• Correlation Modeling: It accounts for how different assets might fail simultaneously during a systemic crisis.

Summary of Results¶

A typical Monte Carlo VaR report provides a spectrum of risk: