Simulating processes

Prefer exact transitions where they exist: GBM, OU, Vasicek, CIR (non-central χ²), Hull-White, Ho-Lee and G2++ are all simulated without discretisation error, so a coarse grid is enough when you only need values on a few dates. Heston uses Andersen's quadratic-exponential scheme by default.

Interest rates fitted to a curve

Short-rate models return the short rate and the pathwise discount factor exp(-∫ r dt). Hull-White, Ho-Lee and G2++ are fitted to a ZeroCurve, so their Monte Carlo bond prices reproduce it:

"""Hull-White fitted to an upward-sloping curve: Monte Carlo bond prices reproduce the curve."""

import numpy as np

import stochPathSim as sps
from stochPathSim.analytics import hull_white_caplet

curve = sps.ZeroCurve([0.5, 1, 2, 5, 10], [0.030, 0.032, 0.035, 0.038, 0.040])
model = sps.HullWhite(a=0.05, sigma=0.01, curve=curve)
paths = model.simulate(T=10.0, n_steps=120, n_paths=50_000, seed=7)

for t in (1.0, 5.0, 10.0):
    mc = paths.at(t, "discount").mean()
    print(f"P(0,{t:>4}) curve {curve.discount(t):.5f}  Monte Carlo {mc:.5f}")

print(f"1y x 3m caplet at 4%: {hull_white_caplet(1.0, 1.25, 0.04, 0.05, 0.01, curve):.6f}")
print(f"Short-rate 99% quantile at 10y: {np.quantile(paths.terminal('rate'), 0.99):.4%}")

ZeroCurve interpolates zero rates with a monotone cubic by default, which keeps instantaneous forward rates continuous: curve-fitted models follow the forward curve, so a linear interpolation would show kinks at the pillars.

Several assets

MultiAssetGBM takes a correlation matrix and returns one state per asset (asset_0, asset_1, …), which Basket and Spread payoffs read directly.

Normal forwards and NIG equity returns

ArithmeticBrownianMotion(x0, mu, sigma) simulates dX = mu dt + sigma dW exactly, including negative values. Drift and volatility are in absolute units. To price a normal forward, set mu=0 under the forward measure and supply the discount rate separately to price_mc. Compare with bachelier_price; bachelier_greeks returns forward delta/gamma, vega per unit normal volatility and theta holding the forward and discount factor fixed. It is not a risk-neutral spot model obtained by setting mu=r.

NormalInverseGaussian(s0, mu, alpha, beta, delta, q=0) models positive equity prices with skewed, heavier-tailed log returns. alpha controls tail decay, beta asymmetry and delta activity per year. Compensation makes the expected spot s0*exp((mu-q)*t). Use mu=r for risk-neutral pricing. Exact inverse Gaussian clock increments work on arbitrary time grids.

The constructor requires alpha > abs(beta) and alpha > abs(beta+1). variance(t) additionally requires alpha > abs(beta+2) when activity is positive; a valid first moment does not guarantee a finite variance. Setting delta=0 gives deterministic growth. log_return_mean and log_return_variance describe log returns, not spot prices.

nig_price prices scalar European calls/puts using SciPy NIG probabilities under the original and stock measures. These are numerically evaluated probabilities, not Monte Carlo estimates. Parameters follow SciPy's NIG parameterisation with a=alpha*delta*t, b=beta*delta*t, scale=delta*t for the unshifted Levy increment.

With NIG, normals= fixes only the Gaussian shocks; reuse the seed to reproduce clock draws too. Antithetic sampling pairs only the Gaussian shocks and does not guarantee a variance reduction for every payoff.

Run the following example to compare both pricing methods and save an interactive chart to examples/output/normal-and-nig.html:

"""Normal forwards and heavy-tailed NIG equities, with a self-contained Plotly comparison.

Run from the repository root after installing the package:
    python examples/06_normal_and_nig.py
"""

from pathlib import Path

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

import stochPathSim as sps
from stochPathSim.analytics import bachelier_greeks, bachelier_price, nig_price
from stochPathSim.pricing import EuropeanCall, price_mc

# Normal volatility is in absolute forward units. A forward has zero drift
# under its forward measure, even though the discount rate is nonzero.
forward, strike, T, rate, normal_vol = -0.005, 0.0, 1.0, 0.015, 0.01
normal = sps.ArithmeticBrownianMotion(forward, mu=0, sigma=normal_vol)
normal_mc = price_mc(normal, EuropeanCall(strike), T, discount=rate,
                     n_steps=1, n_paths=100_000, seed=42, antithetic=True)
normal_exact = bachelier_price(forward, strike, T, np.exp(-rate*T), normal_vol)
print(f"Normal call: MC {normal_mc.price:.7f} +/- {normal_mc.stderr:.7f}; formula {normal_exact:.7f}")
print("Forward Greeks:", bachelier_greeks(forward, strike, T, np.exp(-rate*T), normal_vol))

# NIG delta is an activity parameter, not a Gaussian volatility or an option Greek.
params = dict(alpha=8.0, beta=-3.0, delta=0.25)
nig = sps.NormalInverseGaussian(s0=100, mu=rate, q=0.01, **params)
nig_mc = price_mc(nig, EuropeanCall(105), T, discount=rate, n_steps=1,
                  n_paths=100_000, seed=42, antithetic=True)
nig_exact = nig_price(100, 105, T, rate, **params, q=0.01)
print(f"NIG call: MC {nig_mc.price:.5f} +/- {nig_mc.stderr:.5f}; NIG CDF price {nig_exact:.5f}")

# Equal log-return variance isolates the distribution-shape comparison.
sigma_gbm = np.sqrt(nig.log_return_variance(1.0))
gbm = sps.GBM(100, rate, sigma_gbm, q=0.01)
nig_result = nig.simulate(T=T, n_steps=64, n_paths=5_000, seed=13)
gbm_result = gbm.simulate(T=T, n_steps=64, n_paths=5_000, seed=13)
normal_result = normal.simulate(T=T, n_steps=64, n_paths=40, seed=13)
fig = make_subplots(rows=1, cols=2, subplot_titles=("Normal forward paths (negative values allowed)",
                                                  "Log returns: matched variance, different shape"))
for i, path in enumerate(normal_result["spot"][:15]):
    fig.add_trace(go.Scatter(x=normal_result.times, y=path, mode="lines", opacity=0.45,
                             name="Normal forward", legendgroup="normal", showlegend=i == 0,
                             line=dict(color="#0A66C2", width=1)), row=1, col=1)
for label, result, color in [("NIG", nig_result, "#d97706"), ("GBM", gbm_result, "#0f766e")]:
    fig.add_trace(go.Histogram(x=np.log(result.terminal()/100), histnorm="probability density",
                               nbinsx=80, opacity=0.55, name=label, marker_color=color), row=1, col=2)
fig.update_layout(template="plotly_white", barmode="overlay", height=480,
                  title="stochPathSim: normal forwards and NIG returns (simulated data)")
fig.update_xaxes(title_text="Years", row=1, col=1)
fig.update_yaxes(title_text="Forward", row=1, col=1)
fig.update_xaxes(title_text="log(S_T / S_0)", row=1, col=2)
fig.update_yaxes(title_text="Density", row=1, col=2)
output = Path(__file__).parent / "output" / "normal-and-nig.html"
output.parent.mkdir(exist_ok=True)
fig.write_html(output, include_plotlyjs=True)
print(f"Interactive chart: {output}")

Displaced diffusion and CIR++

Use DisplacedDiffusion(s0, mu, sigma, shift, q=0) for a forward whose value can be negative, but whose shifted value S + shift is positive. The shifted value is simulated exactly as a GBM. For a forward under its forward measure set mu=q=0, so its expectation is the input forward. Call displaced_diffusion_price(F, K, T, df, sigma, shift) to price options using shifted Black-76, including negative strikes; the matching displaced_diffusion_greeks gives forward delta/gamma, sigma-vega and theta at fixed forward and discount factor. Drift acts on S + shift, not on S itself.

CIRPlusPlus(x0, kappa, theta, sigma, curve) fits any ZeroCurve by adding the deterministic difference between market and CIR instantaneous forwards to an exact CIR factor. It returns rate, factor and discount; the positive CIR factor does not guarantee positive shifted short rates. The deterministic discount correction is exact, but the factor's integrated discount is trapezoidal, so simulated bond prices have time-step bias. Refine n_steps when accuracy matters. scheme='euler' allows common Gaussian draws, while the default exact factor draws additional random variates from the seed and ignores the supplied normals.

A runnable example with interactive Plotly charts:

"""Negative-rate forward option and CIR++ yield-curve fitting demonstration.

Run: python examples/07_displaced_and_cirpp.py
Writes a self-contained interactive Plotly figure in examples/output/.
"""

from pathlib import Path

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

import stochPathSim as sps
from stochPathSim.analytics import displaced_diffusion_price
from stochPathSim.pricing import EuropeanCall, price_mc

# A forward can be negative, but F + displacement remains strictly positive.
F, K, shift, sigma, T, r = -0.015, -0.01, 0.06, 0.4, 1.0, 0.03
forward = sps.DisplacedDiffusion(s0=F, mu=0.0, sigma=sigma, shift=shift)
mc = price_mc(forward, EuropeanCall(K), T, discount=r,
              n_paths=100_000, n_steps=1, seed=42, antithetic=True)
closed = displaced_diffusion_price(F, K, T, np.exp(-r * T), sigma, shift)
print(f"Displaced forward call: MC {mc.price:.6f} +/- {mc.stderr:.6f}; closed {closed:.6f}")

# CIR++ shifts a positive CIR factor to reproduce the initial market curve.
curve = sps.ZeroCurve([0.25, 1.0, 3.0], [0.02, 0.026, 0.035])
rates = sps.CIRPlusPlus(x0=0.03, kappa=1.2, theta=0.035, sigma=0.14, curve=curve)
result = rates.simulate(T=3, n_steps=120, n_paths=10_000, seed=42)
print(f"CIR++ 3y bond: MC {result['discount'][:, -1].mean():.6f}; market {curve.discount(3):.6f}")
print("Pathwise CIR++ discounts use a trapezoidal integral of the CIR factor.")

short = forward.simulate(T=1, n_steps=60, n_paths=30, seed=9)
fig = make_subplots(rows=1, cols=2, subplot_titles=(
    "Displaced forward paths (negative values allowed)",
    "CIR++ zero-curve check (Monte Carlo vs market)",
))
for idx, path in enumerate(short["spot"][:20]):
    fig.add_trace(go.Scatter(x=short.times, y=path, mode="lines", opacity=0.35,
                             line=dict(color="#3465A4", width=1.2),
                             showlegend=idx == 0, name="Forward paths"), row=1, col=1)
fig.add_trace(go.Scatter(x=result.times[1:], y=result["discount"][:, 1:].mean(axis=0),
                         name="CIR++ Monte Carlo", mode="lines",
                         line=dict(color="#0F766E")), row=1, col=2)
fig.add_trace(go.Scatter(x=result.times[1:], y=curve.discount(result.times[1:]),
                         name="Market discount curve", mode="lines",
                         line=dict(color="#D97706", dash="dash")), row=1, col=2)
fig.update_layout(template="plotly_white", height=490,
                  title="stochPathSim: displaced forwards and curve-fitted CIR++")
fig.update_xaxes(title_text="Years")
fig.update_yaxes(title_text="Forward", row=1, col=1)
fig.update_yaxes(title_text="Discount factor", row=1, col=2)
output = Path(__file__).parent / "output" / "displaced-and-cirpp.html"
output.parent.mkdir(exist_ok=True)
fig.write_html(output, include_plotlyjs=True)
print(f"Interactive Plotly demo: {output}")

Advanced stochastic models and variance reduction

GARCH(1,1) and GJR-GARCH(1,1) represent volatility clustering and leverage effects as per-observation conditional variance recursions. They are discrete-time statistical models, not risk-neutral pricing models. Do not interpret them as continuous-time diffusions or compare prices across different numbers of time steps without recalibrating the parameters.

import stochPathSim as sps
model = sps.GJRGARCH11(100, 0.03, omega=0.002, alpha=0.05, beta=0.85, gamma=0.08)
result = model.simulate(T=1, n_steps=252, n_paths=1024, seed=42)
result.plot(name="variance").show()

FractionalBrownianMotion uses a Cholesky factorisation of the exact increment covariance for any grid; hurst=0.5 is standard Brownian motion. For other Hurst exponents it is not a semimartingale: do not use it as a risk-neutral traded underlying. The O(n_steps^3) covariance factorisation is intended for grids with a few hundred observations.

For Gaussian-factor Monte Carlo use scrambled Sobol quasi-random normals, optionally ordered by a Brownian bridge so early coordinates determine coarse-scale moves. Use independent scrambled replications to estimate QMC uncertainty; ordinary iid standard errors computed across Sobol points are not valid uncertainty estimates.

from stochPathSim.core import sobol_normals
z = sobol_normals(4096, 32, seed=7, brownian_bridge=True)
paths = sps.GBM(100, 0.03, 0.2).simulate(T=1, n_steps=32, normals=z)
paths.fan().show()

The source example demonstrates all the new capabilities:

"""Advanced models: GARCH/GJR, fractional Brownian motion, Sobol, portfolio attribution.

Run from repository root:
    python examples/08_advanced_models.py

Creates four self-contained interactive Plotly HTML figures in examples/output/.
All data in the example are synthetic.
"""

from pathlib import Path

import numpy as np

import stochPathSim as sps
from stochPathSim import viz
from stochPathSim.core import sobol_normals
from stochPathSim.risk import normal_portfolio_risk

output = Path(__file__).parent / "output"
output.mkdir(exist_ok=True)

# GJR-GARCH reacts more strongly to a negative innovation than GARCH.
common_shocks = np.random.default_rng(12).normal(size=(1000, 252, 1))
garch = sps.GARCH11(100, 0.03, omega=0.0008, alpha=0.07, beta=0.9, v0=0.04)
gjr = sps.GJRGARCH11(100, 0.03, omega=0.0008, alpha=0.05,
                      beta=0.88, gamma=0.1, v0=0.04)
plain = garch.simulate(n_steps=252, normals=common_shocks)
leverage = gjr.simulate(n_steps=252, normals=common_shocks)
fig = viz.paths(leverage, name="variance", n_paths=15)
fig.update_layout(title="GJR-GARCH: volatility clustering (synthetic paths)")
fig.write_html(output / "gjr-volatility.html", include_plotlyjs=True)
print("Terminal mean variance, GARCH / GJR:", plain["variance"][:, -1].mean(),
      leverage["variance"][:, -1].mean())

# Hurst exponent controls persistence and roughness, not Gaussian variance alone.
fbm = sps.FractionalBrownianMotion(x0=0, hurst=0.25, sigma=1.0)
fpaths = fbm.simulate(T=1, n_steps=80, n_paths=150, seed=7)
fig = viz.paths(fpaths, name="value", n_paths=25)
fig.update_layout(title="Fractional Brownian motion (H = 0.25)")
fig.write_html(output / "fractional-brownian.html", include_plotlyjs=True)

# Scrambled Sobol with a Brownian bridge can feed any Gaussian-factor simulator.
z = sobol_normals(4096, 32, seed=5, brownian_bridge=True)
qmc_result = sps.GBM(100, 0.03, 0.20).simulate(T=1, n_steps=32, normals=z)
fig = viz.fan(qmc_result)
fig.update_layout(title="GBM: Sobol / Brownian-bridge scenario fan")
fig.write_html(output / "sobol-fan.html", include_plotlyjs=True)

# Signed Euler contributions correctly handle hedging positions.
covariance = np.array([[0.04, 0.012, 0.004],
                       [0.012, 0.09, 0.008],
                       [0.004, 0.008, 0.025]])
portfolio = normal_portfolio_risk(weights=[1.0, -0.2, 0.7],
                                  covariance=covariance, level=0.99)
fig = viz.portfolio_risk_chart(portfolio, labels=["Equity", "Rates hedge", "Credit"])
fig.write_html(output / "portfolio-risk-contributions.html", include_plotlyjs=True)
print(f"99% VaR: {portfolio.var:.5f}, ES: {portfolio.es:.5f}")
print("VaR contributions:", portfolio.component_var)
print("ES contributions:", portfolio.component_es)
print("Generated HTML charts:", output)