Highly Efficient Geometric Brownian Motion Modeling with QMCPy

pip install qmcpy numpy matplotlib

import numpy as np
import matplotlib.pyplot as plt
import qmcpy as qp
sampler = qp.Lattice(252, seed=42)  # daily steps for 1 year
gbm = qp.GeometricBrownianMotion(sampler, t_final=1, initial_value=1, drift=0.05, diffusion=0.2)
paths = gbm.gen_samples(16)
t = np.linspace(0, 1.0, paths.shape[1])
plt.plot(t, paths[:5].T, alpha=0.8)
plt.xlabel('$t$'); plt.ylabel('$S_t$'); plt.title('GBM paths');
plt.show()
    

Introduction

In this blog, we demonstrate how to simulate and analyze a geometric Brownian motion (GBM) process using QMCPy in Python. GBM is widely used in finance to model stock prices and other assets. We will walk through key code snippets, plots, and insights. The numerical results can be reproduced using the Jupyter notebook available at https:/github.com/QMCSoftware/QMCSoftware/blob/master/demos/gbm_demo.ipynb.

GBM is a continuous stochastic process in which the natural logarithm of its values follows a Brownian motion (BM) [1]. Mathematically, it can be defined as follows:

where

• \(S_0\) is the initial value,
• \(\mu\) is a drift coefficient
• \(\sigma\) is the volatility (note: QMCPy uses diffusion = \(\sigma^2\))
• \(W_t\) is a (standard) BM.

At any time \(0 < t \le T\), where \(0\) and \(T\) represent the beginning and end time of the process, \(S_t\) follows a log-normal distribution with expected value and variance as follows (see Section 3.2 in \([1]\)):

• \(E[S_t] = S_0 e^{\mu t}\)
• \(\text{Var}[S_t] = S_0^2 e^{2\mu t}(e^{\sigma^2 t} – 1)\)
• \(
  \text{Cov}(S_{t_i}, S_{t_j}) = S_0^2 e^{\mu(t_i + t_j)} \left(e^{\sigma^2
  \min(t_i, t_j)} – 1\right).\)

GBM is commonly used to model stock prices driving option payoffs in derivatives pricing \([4, 5]\).

GBM Objects in QMCPy

GBM in QMCPy inherits from BrownianMotion class \([Choi2022,QMCPy2020a]\). We can instantiate a GBM class and generate sample paths to see the class in action; see Listing 1:

Listing 1: Generating 4 GBM sample paths evaluated at 2 time points using QMCPy.

# GBM/code/gbm_qmcpy.py
import qmcpy as qp
qp_gbm = qp.GeometricBrownianMotion(qp.Lattice(2, seed=42))
qp_gbm.gen_samples(n=4)

The output shows 4 sample paths evaluated at 2 time points, yielding a \((4\times 2)\) array where rows represent paths and columns represent time steps:

Log-Normality Property

The log-normal property is fundamental in financial modeling because it ensures asset prices remain strictly positive (as required in reality) while allowing for unlimited upside potential. This property makes GBM the cornerstone of the Black-Scholes model and many derivative pricing frameworks.

To validate theoretical properties, we generate \(2^{12} = 4096\) paths over 5 time steps and compare empirical moments with theoretical values. Note that the theoretical values match the last values captured in qp_gbm.mean_gbm and qp_gbm.covariance_gbm for the final time point. The results are
shown in Table 1.

Listing 2: Generating and validating GBM sample moments.

# GBM/code/gbm_moments_validation.py
# Generate GBM samples for theoretical validation
import qmcpy as qp
import numpy as np

S0, mu, sigma, T, n_samples = 100.0, 0.05, 0.20, 1.0, 2**12
diffusion = sigma**2
sampler = qp.Lattice(5, seed=42)
qp_gbm = qp.GeometricBrownianMotion(sampler, t_final=T, initial_value=S0, drift=mu, diffusion=diffusion)
paths = qp_gbm.gen_samples(n_samples)
S_T = paths[:, -1]  # Final values only

# Calculate theoretical vs empirical sample moments
theo_mean = S0 * np.exp(mu * T)
theo_var = S0**2 * np.exp(2*mu*T) * (np.exp(diffusion * T) - 1)
qp_emp_mean = np.mean(S_T)
qp_emp_var = np.var(S_T, ddof=1) 
print(f"Mean: {qp_emp_mean:.3f} (theoretical: {theo_mean:.3f})")
print(f"Variance: {qp_emp_var:.3f} (theoretical: {theo_var:.3f})")
qp_gbm

Table 1: Theoretical vs Empirical Validation of GBM Properties.

GBM vs BM

Below we compare BM and GBM using the same parameters: \(\texttt{drift} = 0\), \(\texttt{diffusion} = 1\), \(\texttt{initial_value} = 1\). The driftless BM paths should fluctuate symmetrically around the initial value (\(y = 1\)) and can take negative values, while those of the GBM remain strictly positive. We generate these samples in Listing 3.

Listing 3: Generate 16 BM and 16 GBM sample paths for comparison.

# GBM/code/bm_gbm_16.py
n = 16
sampler = qp.Lattice(2**7, seed=42)
plot_paths('BM', sampler, t_final=1, initial_value=1, drift=0, diffusion=1, n=n)
plot_paths('GBM', sampler, t_final=1, initial_value=1, drift=0, diffusion=1, n=n)

Next, we demonstrate how easily one can swap samplers or change parameters in QMCPy. For example, to model a stock price with initial value \(S_0=50\), drift \(\mu=0.1\), and volatility \(\sigma=\sqrt{0.2}\) over a 5‐year horizon using IID sampling; see Listing 4.

Listing 4: Generate 32 GBM sample paths with IID sampling.

# GBM/code/gbm_iid_32.py
gbm_iid = plot_paths('GBM', qp.IIDStdUniform(2**8, seed=42), t_final=5, initial_value=50, drift=0.1, diffusion=0.2, n=32)

Also, we can use a low‐discrepancy lattice sampler with the same parameters (Listing 5):

Listing 5: Generate 32 GBM sample paths with lattice sampling.

# GBM/code/gbm_lattice_32.py
gbm_lattice = plot_paths('GBM', qp.Lattice(2**8, seed=42), t_final=5, initial_value=50, drift=0.1, diffusion=0.2, n=32)

The generated sample paths are plotted in Figure 1. Subplots (\(a\)), (\(b\)), (\(c\)), and (\(d\)) show, respectively:

• Figure 1: Comparison of sample paths for BM and GBM using different samplers..
• BM with lattice sampler (\(T=1\), \(S_0=1\), \(\mu=0\), \(\sigma^2=1\), 16
  paths)
• GBM with lattice sampler (\(T=1\), \(S_0=1\), \(\mu=0\), \(\sigma^2=1\), 16
  paths)
• GBM with IID sampler (\(T=5\), \(S_0=50\), \(\mu=0.1\), \(\sigma^2=0.2\), 32
  paths)
• GBM with lattice sampler (\(T=5\), \(S_0=50\), \(\mu=0.1\), \(\sigma^2=0.2\),
  32 paths).

(a)

(b)

(c)

(d)

QuantLib vs QMCPy Comparison

In this section, we compare QMCPy’s GeometricBrownianMotion implementation with the industry-standard QuantLib library \([6]\) to validate its accuracy and performance. The numerical results are summarized in Table 2.

Both libraries produce statistically equivalent GBM simulations that match theoretical values. QMCPy typically runs 1.5 to 2 times faster due to vectorized operations, lazy loading, and optimized memory management. More importantly, it demonstrates superior numerical accuracy (lower mean absolute errors) with
Sobol, lattice and Halton samplers, making it excellent for research and high-performance applications. QuantLib remains the industry standard for production systems that require comprehensive support for financial modeling and risk management.

Listing 6: Function to generate GBM paths with QuantLib.

# GBM/code/quantlib_util.py
import QuantLib as ql
import numpy as np

def generate_quantlib_paths(initial_value, mu, sigma, maturity, n_steps, n_paths, sampler_type='IIDStdUniform', seed=7):
    """Generate GBM paths using QuantLib with configurable sampler type and seed."""
    gbm = ql.GeometricBrownianMotionProcess(initial_value, mu, sigma)
    times = ql.TimeGrid(maturity, n_steps)
    dimension = n_steps
    if sampler_type == 'IIDStdUniform':
        uniform_rng = ql.UniformRandomGenerator(seed)
        sequence_gen = ql.GaussianRandomSequenceGenerator(ql.UniformRandomSequenceGenerator(n_steps, uniform_rng))
        path_gen = ql.GaussianPathGenerator(gbm, maturity, n_steps, sequence_gen, False)
        paths = np.zeros((n_paths, n_steps + 1))
        for i in range(n_paths):
            sample_path = path_gen.next().value()
            paths[i, :] = np.array([sample_path[j] for j in range(n_steps + 1)])
        return paths, gbm
    elif sampler_type == 'Sobol':
        uniform_rsg = ql.UniformLowDiscrepancySequenceGenerator(dimension, seed)
        sequence_gen = ql.GaussianLowDiscrepancySequenceGenerator(uniform_rsg)
        path_gen = ql.GaussianSobolMultiPathGenerator(gbm, list(times), sequence_gen, False)
        paths = np.zeros((n_paths, n_steps + 1))
        paths[:, 0] = initial_value  # Set initial value
        for i in range(n_paths):
            sample_path = path_gen.next().value()
            path_values = sample_path[0]  # For 1D process, get the first (and only) path
            paths[i, :] = np.array([path_values[j] for j in range(n_steps + 1)])
        return paths, gbm
    else:
        raise ValueError(f"Unsupported sampler type: {sampler_type}. Use 'IIDStdUniform' or 'Sobol'")

Table 2: GBM Final Value Statistics and Performance Comparison

• Figure 2: QMCPy vs QuantLib comparison. Top: sample paths from QMCPy (left) and
  QuantLib (right). Bottom left: marginal distribution at \(t=1\) showing close
  agreement between libraries and with Table 2. Bottom right: QMCPy
  covariance heatmap (time 0 at top), consistent with the numerical matrix in
  Table 1.
  

The evaluation of computational efficiency was done by creating comprehensive performance benchmarks comparing QMCPy and QuantLib across two key scaling dimensions. The benchmarks were performed using the perfplot library, which automatically handles warm-up, multiple runs, and statistical analysis to ensure reliable timing measurements.

Figure 3 presents the results of our performance analysis. The left panel (a) shows how execution time scales with the number of time steps while keeping the number of paths fixed at 4,096. Both libraries exhibit approximately linear scaling, but QMCPy demonstrates superior performance at smaller time step counts, with QuantLib becoming more competitive as the number of time steps increases. The right panel (b) examines scaling
behavior with respect to the number of paths while fixing the time steps at 252 (representing a typical trading year). Here, QMCPy maintains a consistent performance advantage across all path counts, with the gap becoming more pronounced at higher path numbers. This performance difference is particularly relevant for Monte Carlo applications that require large numbers of simulation paths for accurate estimation.

• Figure 3: GBM Performance Comparison: QuantLib vs QMCPy. Left plot shows performance scaling with number of time steps for fixed paths. Right plot shows performance scaling with number of paths for fixed number of time steps.
  

In order to further validate these performance findings, we conducted an extended parameter sweep analysis across a broader range of configurations. Figure 4 presents the comprehensive results of this analysis, systematically examining performance across varying time steps and path counts.

Regarding accuracy, QMCPy’s Halton sampler generally achieves the lowest mean absolute error (MAE), particularly when a larger number of paths are used, reaching errors below \(10^{-4}\). QMCPy’s Sobol sampler demonstrates remarkable consistency, maintaining relatively stable MAEs across different time steps and showing clear convergence improvements with more paths. The lattice sampler also provides competitive accuracy, especially with higher path counts. While QuantLib’s standard uniform IID sampler yields slightly lower MAE for a larger number of paths compared to QMCPy’s IID sampler, it is notably slower.

Regarding speed, low-discrepancy samplers like Sobol, lattice, and Halton demonstrate superior convergence rates with increasing path counts compared to IID methods. With QMCPy, Sobol and lattice samplers offer the best
speed-accuracy trade-off. QuantLib achieves comparable runtime to QMCPy’s faster samplers, but without the accuracy benefits. The Halton sampler, while yielding the most accurate results, incurs significantly higher computational costs. These results highlight QMCPy’s quasi-Monte Carlo methods as particularly well-suited for applications requiring high accuracy, with Sobol and lattice samplers providing an optimal balance of speed and precision for most practical scenarios.

• Figure 4: Comprehensive parameter sweep performance analysis comparing QuantLib and QMCPy across varying time steps and path counts. The results demonstrate QMCPy’s consistent performance advantages and superior scaling characteristics for both dimensions of the parameter space.
  

Internals

The GeometricBrownianMotion class in QMCPy is engineered for speed, robustness, and mathematical correctness. Its design leverages object-oriented inheritance and vectorized operations, resulting in both flexibility and high performance. GeometricBrownianMotion inherits from BrownianMotion which itself inherits from Gaussian. This layered design allows the GBM class to reuse and extend efficient implementations for Gaussian random vectors and BM increments. The constructor rigorously checks input parameters (e.g., positivity of initial value and diffusion, valid decomposition type), ensuring mathematical integrity and preventing run-time errors.

The class uses vectorized NumPy operations to generate entire arrays of GBM paths in a single call, minimizing Python loops and maximizing computational throughput. Sample generation proceeds in two stages:

1. The parent class BrownianMotion generates standard BM sample paths using the specified sampler (e.g., low-discrepancy lattice, IID uniform), with drift and diffusion handled in the mean and covariance structure.
2. The GBM class transforms the BM samples via the exponential mapping, performed in a fully vectorized fashion, ensuring that thousands of paths can be efficiently simulated.

The class computes and stores the theoretical mean and covariance matrices for GBM at initialization, which can be used for validation and theoretical comparisons. Both mean and covariance are calculated using analytical formulas, leveraging NumPy’s broadcasting for efficient computations. (Briefly, broadcasting in NumPy allows
arithmetic operations between arrays of different shapes by automatically expanding the smaller array to match the shape of the larger array.)

The Gaussian and BM classes both implement Cholesky and PCA factorization of the covariance matrix

respectively. In the Cholesky method, one computes the lower‐triangular \(L=\operatorname{chol}(\Sigma)\) and obtains correlated increments via \(X=L\,Z\), where \(Z\sim\mathcal{N}(0,I)\). In the PCA approach, one first diagonalizes \(\Sigma=PDP^{\!\top}\), then forms \(X = P\,D^{1/2}\,Z\). By default PCA is used for its superior numerical stability in high dimensions and slightly lower cost when many eigenvalues are near zero. These correlated normals feed directly into the GBM update

ensuring that the simulated paths respect the intended covariance structure and
remain strictly positive (with strict‐positivity checks raising warnings or
errors if violated).

Conclusions and Future Work

This blog demonstrates that QMCPy’s quasi-Monte Carlo implementations provide significant advantages over traditional Monte Carlo methods for modeling geometric Brownian motion. QMCPy’s approach combines superior numerical accuracy with enhanced computational efficiency, making it particularly well-suited for high-performance financial modeling applications.

In the future, we aim to investigate ways to improve the runtime of the Halton sampler. For example, we may consider starting with Halton sampling points for high accuracy in early iterations, then switch to Sobol for faster convergence as sample size increases. It would also be interesting to experiment with ensemble sampling by running multiple samplers in parallel and combine the results using weighted averaging based on their relative accuracies.

Acknowledgments

The authors thank Joshua Jay Herman and Jiangrui Kang for their insightful feedback and help with the blog post.

References

1. Glasserman, P. (2003). Monte Carlo Methods in Financial Engineering. Springer.
2. Choi, S.-C. T., Hickernell, F. J., Jagadeeswaran, R., McCourt, M. J., & Sorokin, A. G. (2022). Quasi-Monte Carlo Software. In A. Keller (Ed.), Monte Carlo and Quasi-Monte Carlo Methods. Springer International Publishing. DOI:
   10.1007/978-3-030-98319-2_2.
3. Choi, S.-C. T., Hickernell, F. J., Jagadeeswaran, R., McCourt, M. J., & Sorokin, A. G. (2020–2025). QMCPy: A quasi-Monte Carlo Python Library (Version 2.1). DOI:
   10.5281/zenodo.16822646. Retrieved from
   www.qmcpy.org.
4. Hull, J. C. (2017). Options, Futures, and Other Derivatives (10th ed.). Pearson.
5. Ross, S. M. (2014). Introduction to Probability Models (11th ed.). Academic Press.
6. The QuantLib contributors (2003–2025). QuantLib: A free/open-source library for quantitative finance (Version 1.38). DOI:
   10.5281/zenodo.1440997. Retrieved from
   www.quantlib.org.

Larysa Matiukha

+ posts Bio

• Larysa Matiukha

  Highly Efficient Geometric Brownian Motion Modeling with QMCPy

Aleksei Sorokin

Website |  + posts Bio

Aleksei Sorokin is an Applied Mathematics and Data Science student at the Illinois Institute of Technology.

• Aleksei Sorokin

  A QMCPy Quick Start
• Aleksei Sorokin

  A Collection of QMCPy Posters
• Aleksei Sorokin

  A Seminar at PyData Chicago
• Aleksei Sorokin

  QMCPy Version 1.0
• Aleksei Sorokin

  Speeding up QMCPy with Distributable C Code
• Aleksei Sorokin

  Visualizing the Internals of Object Classes in QMCPy
• Aleksei Sorokin

  Quasi-Monte Carlo Software Article
• Aleksei Sorokin

  A Talk at the Chicago Area SIAM Student Conference
• Aleksei Sorokin

  A Presentation at IIT’s Computational Mathematics Seminar
• Aleksei Sorokin

  A Talk at the Great Lakes SIAM Conference
• Aleksei Sorokin

  QMCPy Events Coming Soon
• Aleksei Sorokin

  Highly Efficient Geometric Brownian Motion Modeling with QMCPy

Sou-Cheng Choi

Website |  + posts Bio

Dr. Sou-Cheng T. Choi is Research Associate Professor of Applied Mathematics at the Illinois Institute of Technology and Founder of SouLab LLC. She formerly served as Chief Data Scientist at the Kamakura Corporation, acquired by SAS Institute Inc. in 2022.

• Sou-Cheng Choi

  A QMCPy Quick Start
• Sou-Cheng Choi

  QMCPy Version 1.0
• Sou-Cheng Choi

  Visualizing the Internals of Object Classes in QMCPy
• Sou-Cheng Choi

  Quasi-Monte Carlo Software Article
• Sou-Cheng Choi

  Highly Efficient Geometric Brownian Motion Modeling with QMCPy
• Sou-Cheng Choi

  Parsl Accelerated QMCPy Notebook Tests