financialkit.monte_carlo#
Crude Monte Carlo with standard errors and confidence intervals,
Box-Muller normal variates and option prices by simulation through
financialkit.integrators; variance reduction by stratification
(Neyman allocation), importance sampling, antithetic variates, Latin
hypercube sampling, regression control variates and moment matching;
quasi-Monte Carlo with van der Corput, Halton and Sobol points
(scipy.stats.qmc), exact and \(L^2\) star discrepancies, the
Koksma-Hlawka bound, Owen-scrambled randomized QMC and Brownian-bridge
path construction; likelihood-ratio and pathwise Greeks; and Giles’s
adaptive multilevel Monte Carlo.
financialkit.monte_carlo: Monte Carlo estimation, variance reduction, quasi-Monte Carlo and Greeks.
Crude Monte Carlo after Metropolis and Ulam, with standard errors and
confidence intervals on every estimate; Box-Muller normal variates; and
Boyle’s option prices under Black-Scholes, simulated through
financialkit.integrators. Variance reduction: Neyman’s stratified
sampling, importance sampling by a mean shift, antithetic variates,
Latin hypercube sampling (scipy.stats.qmc), regression control
variates and Barraquand’s moment matching. Quasi-Monte Carlo: van der
Corput, Halton and Sobol sequences, exact and \(L^2\) star
discrepancies with the Koksma-Hlawka bound, Owen-scrambled randomized
QMC with error bars, and high-dimensional path integrals with and
without Brownian-bridge ordering. Greeks by the likelihood-ratio and
pathwise methods, and Giles’s adaptive multilevel Monte Carlo.
- class financialkit.monte_carlo.ControlVariateResult(estimate, std_error, n_samples, confidence=0.95, method='crude', beta=<factory>, r_squared=0.0, crude_estimate=nan, crude_std_error=nan)[source]#
Bases:
MonteCarloResultA control-variate estimate \(\bar Y - \beta^\top(\bar X - \mu_X)\) (Lavenberg-Welch 1981).
- Parameters:
- class financialkit.monte_carlo.MLMCResult(estimate, std_error, n_samples, confidence=0.95, method='crude', level_means=<factory>, level_variances=<factory>, level_samples=<factory>, level_costs=<factory>, cost=0.0, target_rmse=nan, M=2)[source]#
Bases:
MonteCarloResultA multilevel Monte Carlo estimate \(\sum_l E[P_l - P_{l-1}]\) (Giles 2008).
Level
lusesM**ltime steps; level 0 estimates \(E[P_0]\) itself.costcounts fine-grid time steps, the work unit of Giles’s complexity theorem, andn_samplesthe paths over all levels.- Parameters:
- level_costs: ndarray#
Time steps per sample on each level (fine plus coarse).
- Type:
ndarray, shape (L + 1,)
- level_variances: ndarray#
Sample variance of \(P_l - P_{l-1}\) on each level.
- Type:
ndarray, shape (L + 1,)
- property standard_mc_cost: float#
Time steps standard Monte Carlo on the finest level would need for the same variance.
That is \(V_L^{\text{fine}} C_L / (\varepsilon^2 / 2)\), with the finest level’s payoff variance approximated by the level-0 variance (the payoff variance barely depends on the step size).
- Type:
- class financialkit.monte_carlo.MonteCarloGreeks(price, delta, vega, gamma, price_se, delta_se, vega_se, gamma_se, n_paths, method='')[source]#
Bases:
objectMonte Carlo estimates of an option’s price, delta, vega and gamma, each with a standard error.
A Greek an estimator cannot produce (for example a pathwise gamma of a payoff with a kink) is
nan.- Parameters:
- class financialkit.monte_carlo.MonteCarloResult(estimate, std_error, n_samples, confidence=0.95, method='crude')[source]#
Bases:
objectA Monte Carlo estimate with its standard error and confidence interval.
- property half_width: float#
Half-width \(z_{(1+c)/2}\,\mathrm{SE}\) of the confidence interval.
- Type:
- property variance: float#
Variance per evaluation, \(n\,\mathrm{SE}^2\); the cost-adjusted figure of merit.
- Type:
- variance_ratio(baseline)[source]#
Variance reduction factor relative to baseline, per function evaluation.
Above 1 means this estimator needs that many times fewer evaluations than baseline for the same standard error.
- Return type:
- Parameters:
baseline (MonteCarloResult)
- class financialkit.monte_carlo.StratifiedResult(estimate, std_error, n_samples, confidence=0.95, method='crude', allocation=<factory>, stratum_means=<factory>, stratum_stds=<factory>)[source]#
Bases:
MonteCarloResultA stratified estimate over equal-probability strata of \([0, 1)\) (Neyman 1934).
- Parameters:
- financialkit.monte_carlo.antithetic_estimate(f, n_pairs, dim=None, rng=None, confidence=0.95)[source]#
Antithetic-variates estimate of \(E[f(Z)]\), \(Z \sim N(0, I)\) (Hammersley-Morton 1956).
Each normal draw is paired with its mirror image and the estimator averages the pair means \(\tfrac12[f(Z) + f(-Z)]\), which are i.i.d. Per evaluation the variance falls by the factor \(1 + \rho\), with \(\rho = \operatorname{corr}(f(Z), f(-Z))\), so it helps exactly when \(f\) is monotone.
- Parameters:
- Return type:
- Returns:
MonteCarloResult
Examples
For a linear integrand the pairs cancel exactly:
>>> antithetic_estimate(lambda z: 3.0 + z, 1000, rng=0).std_error < 1e-12 True
- financialkit.monte_carlo.box_muller(u1, u2)[source]#
Box-Muller transform of two independent uniforms into two independent standard normals.
\[R = \sqrt{-2 \ln u_1}, \quad z_1 = R \cos(2\pi u_2), \quad z_2 = R \sin(2\pi u_2).\]\(R^2\) is exponential with mean 2, the law of the squared radius of a standard bivariate normal, and the angle is uniform.
- financialkit.monte_carlo.box_muller_normals(n, rng=None)[source]#
nstandard normals by the Box-Muller transform of fresh uniforms.- Parameters:
- Return type:
- Returns:
ndarray, shape (n,)
Examples
>>> z = box_muller_normals(200_000, rng=0) >>> round(float(z.mean()), 2) == 0.0, round(float(z.std()), 2) (True, 1.0)
- financialkit.monte_carlo.control_variate_estimate(y, controls, control_means, confidence=0.95)[source]#
Control-variate estimate of \(E[Y]\) from controls with known means (Lavenberg-Welch 1981).
The coefficients \(\beta\) are the least-squares slopes of \(Y\) on the controls \(X\), and
\[\hat\mu = \bar Y - \beta^\top (\bar X - \mu_X), \qquad \operatorname{Var}\hat\mu \approx (1 - R^2) \frac{s_Y^2}{n}.\]The standard error uses the regression residuals with
n - k - 1degrees of freedom. Discounted \(S_T\) (mean \(S e^{-qT}\)) is the classic control for an option payoff.- Parameters:
- Return type:
- Returns:
ControlVariateResult
Examples
>>> import numpy as np >>> x = np.random.default_rng(0).standard_normal(10_000) >>> res = control_variate_estimate(np.exp(x), x, 0.0) >>> res.covers(np.exp(0.5)), res.std_error < res.crude_std_error (True, True)
- financialkit.monte_carlo.european_option_mc(S, K, T, r, sigma, q=0.0, option_type='call', n_paths=100000, rng=None, z=None, confidence=0.95)[source]#
European option price by crude Monte Carlo under Black-Scholes (Boyle 1977).
\[V = e^{-rT} E^Q[\text{payoff}(S_T)] \approx \frac{e^{-rT}}{n} \sum_i \text{payoff}(S_T^{(i)})\]- Parameters:
S (
float) – Spot, strike, maturity, risk-free rate and volatility.K (
float) – Spot, strike, maturity, risk-free rate and volatility.T (
float) – Spot, strike, maturity, risk-free rate and volatility.r (
float) – Spot, strike, maturity, risk-free rate and volatility.sigma (
float) – Spot, strike, maturity, risk-free rate and volatility.q (
float) – Continuous dividend yield.option_type (
str)n_paths (
int) – Number of terminal prices; ignored if z is given.z (
ArrayLike|None) – Standard normals to use instead of fresh ones.confidence (
float) – Coverage of the confidence interval.
- Return type:
- Returns:
MonteCarloResult
Examples
Hull, 10th ed., Example 15.6 (Black-Scholes value 4.76):
>>> res = european_option_mc(42.0, 40.0, 0.5, 0.10, 0.20, n_paths=200_000, rng=0) >>> res.covers(4.7594) True
- financialkit.monte_carlo.gbm_path_estimate(payoff, S, T, r, sigma, n_steps, n_paths, sampler='pseudo', q=0.0, n_replications=8, rng=None, confidence=0.95)[source]#
Discounted path-dependent payoff under Black-Scholes, by pseudo- or quasi-random paths.
Paths are exact GBM on
n_stepsdates fromfinancialkit.integrators.gbm_paths(), so the integral has dimensionn_steps. The increments come frombrownian_increments()("pseudo"), from scrambled Sobol points with coordinateidriving stepi("sobol", Paskov-Traub 1995), or from the same points assembled by Brownian bridge so the leading coordinates set the coarse shape ("bridge", Caflisch-Morokoff-Owen 1997). Every sampler is replicated n_replications times with fresh seeds and the standard error is taken across replications, so the three are compared on one footing.- Parameters:
payoff (
Callable[[ndarray],ArrayLike]) – Undiscounted payoff of the price paths of shape(n_paths, n_steps + 1)(column 0 isS).S (
float) – Spot, maturity, risk-free rate and volatility.T (
float) – Spot, maturity, risk-free rate and volatility.r (
float) – Spot, maturity, risk-free rate and volatility.sigma (
float) – Spot, maturity, risk-free rate and volatility.n_steps (
int) – Monitoring dates, the dimension of the integral.n_paths (
int) – Paths per replication (a power of two for the Sobol samplers).sampler (
str)q (
float) – Continuous dividend yield.n_replications (
int) – Independent replications.confidence (
float) – Coverage of the confidence interval.
- Return type:
- Returns:
MonteCarloResult
- financialkit.monte_carlo.halton(n, dim, scramble=False, rng=None)[source]#
npoints of thedim-dimensional Halton sequence (Halton 1960).Coordinate \(i\) is the van der Corput sequence in the \(i\)-th prime. Unscrambled points start at index 0 (the origin), as
scipy.stats.qmc.Haltondoes; neighbouring large primes give visibly correlated projections, which scrambling removes.
- financialkit.monte_carlo.importance_sampling_normal(f, n, shift, rng=None, confidence=0.95)[source]#
Importance-sampling estimate of \(E[f(Z)]\), \(Z \sim N(0, I)\), by a mean shift.
Sampling \(Z \sim N(\theta, I)\) instead and weighting by the likelihood ratio
\[E[f(Z)] = E_\theta\left[f(Z)\,e^{-\theta^\top Z + |\theta|^2/2}\right]\]concentrates samples where \(f\) matters (Kahn-Marshall 1953; Glasserman 2004, Sec. 4.6). A shift towards the strike turns a deep out-of-the-money option or a rare loss from a few lucky hits into a routine estimate.
- Parameters:
- Return type:
- Returns:
MonteCarloResult
Examples
\(P(Z > 4) \approx 3.17 \times 10^{-5}\), from 20,000 samples:
>>> res = importance_sampling_normal(lambda z: (z > 4.0).astype(float), 20_000, shift=4.0, rng=0) >>> abs(res.estimate - 3.1671e-5) < 4 * res.std_error, res.std_error < 1e-6 (True, True)
- financialkit.monte_carlo.koksma_hlawka_bound(points, variation)[source]#
Koksma-Hlawka bound \(V_{HK}(f)\,D_N^*\) on the error of the equal-weight rule at points.
\[\left| \frac1N \sum_i f(x_i) - \int_{[0,1]^d} f \right| \le V_{HK}(f)\,D_N^*(x_1, \dots, x_N).\]
- financialkit.monte_carlo.l2_star_discrepancy(points)[source]#
The \(L^2\)-star discrepancy, in any dimension (Warnock’s formula,
scipy.stats.qmc.discrepancy()).The root-mean-square of the local discrepancy over anchored boxes; it never exceeds \(D_N^*\) and is cheap in high dimension.
- financialkit.monte_carlo.latin_hypercube_estimate(f, n, dim, n_replications=10, rng=None, confidence=0.95)[source]#
Latin hypercube estimate of \(\int_{[0,1)^d} f(u)\,du\) (McKay-Beckman-Conover 1979).
Each coordinate of the
npoints takes exactly one value in each of thenintervals \([k/n, (k+1)/n)\), in an independent random order per coordinate, so every one-dimensional margin is stratified at once. The variance falls most for nearly additive integrands. The points of one design are dependent, so the standard error comes from n_replications independent designs.- Parameters:
- Return type:
- Returns:
MonteCarloResult
- financialkit.monte_carlo.likelihood_ratio_greeks(S, K, T, r, sigma, q=0.0, option_type='call', digital=False, n_paths=100000, rng=None, z=None)[source]#
Price, delta, vega and gamma by the likelihood-ratio (score-function) method (Glynn 1987).
Each Greek is \(E[e^{-rT}\,\text{payoff}(S_T)\,w(Z)]\) with the score weight \(w\) from the module docstring, estimated on the same paths as the price.
- Parameters:
S (
float) – Spot, strike, maturity, risk-free rate and volatility.K (
float) – Spot, strike, maturity, risk-free rate and volatility.T (
float) – Spot, strike, maturity, risk-free rate and volatility.r (
float) – Spot, strike, maturity, risk-free rate and volatility.sigma (
float) – Spot, strike, maturity, risk-free rate and volatility.q (
float) – Continuous dividend yield.option_type (
str)digital (
bool) – Price a cash-or-nothing digital paying 1 instead of a vanilla.n_paths (
int) – Number of terminal prices; ignored if z is given.z (
ArrayLike|None) – Standard normals to use, e.g. common random numbers withpathwise_greeks().
- Return type:
- Returns:
MonteCarloGreeks
Examples
Black-Scholes delta of the Hull Example 15.6 call is 0.7791:
>>> g = likelihood_ratio_greeks(42.0, 40.0, 0.5, 0.10, 0.20, n_paths=400_000, rng=0) >>> abs(g.delta - 0.7791) < 4 * g.delta_se True
- financialkit.monte_carlo.mlmc_level(level, n_paths, S, K, T, r, sigma, q=0.0, option_type='call', M=2, rng=None)[source]#
Samples of the level-
lcorrection \(P_l - P_{l-1}\) for a European option.\(P_l\) is the discounted payoff of the Euler-Maruyama solution of risk-neutral GBM with
M**levelsteps; the coarse solution sums the fine Brownian increments in groups of M, so the two are coupled. Level 0 returns \(P_0\) itself (one Euler step).- Parameters:
level (
int) – Level \(l \ge 0\).n_paths (
int) – Number of coupled samples.S (
float) – Spot, strike, maturity, risk-free rate and volatility.K (
float) – Spot, strike, maturity, risk-free rate and volatility.T (
float) – Spot, strike, maturity, risk-free rate and volatility.r (
float) – Spot, strike, maturity, risk-free rate and volatility.sigma (
float) – Spot, strike, maturity, risk-free rate and volatility.q (
float) – Continuous dividend yield.option_type (
str)M (
int) – Refinement factor between levels.
- Return type:
- Returns:
ndarray, shape (n_paths,)
Examples
The corrections vary far less than the payoff itself:
>>> y0 = mlmc_level(0, 20_000, 100.0, 100.0, 1.0, 0.05, 0.2, rng=0) >>> y4 = mlmc_level(4, 20_000, 100.0, 100.0, 1.0, 0.05, 0.2, rng=0) >>> bool(y4.var() < 0.01 * y0.var()) True
- financialkit.monte_carlo.moment_match(z)[source]#
Rescale normal draws to sample mean exactly 0 and sample variance exactly 1, per column (Barraquand 1995).
Examples
>>> import numpy as np >>> w = moment_match(np.random.default_rng(0).standard_normal((500, 3))) >>> bool(np.allclose(w.mean(axis=0), 0.0) and np.allclose(w.std(axis=0), 1.0)) True
- financialkit.monte_carlo.moment_matching_estimate(f, n, dim=None, n_replications=10, rng=None, confidence=0.95)[source]#
Estimate \(E[f(Z)]\) from moment-matched normals (Barraquand 1995).
Each replication draws
nnormals, appliesmoment_match()and averages \(f\). Matching makes the samples dependent and the estimator slightly biased (\(O(1/n)\)), so the standard error comes from n_replications independent replications.
- financialkit.monte_carlo.monte_carlo_integrate(f, n, dim=1, rng=None, confidence=0.95)[source]#
Crude Monte Carlo estimate of \(\int_{[0,1)^d} f(u)\,du\) (Metropolis-Ulam 1949).
- Parameters:
- Return type:
- Returns:
MonteCarloResult
Examples
The area of the quarter disc is \(\pi / 4\):
>>> import numpy as np >>> res = monte_carlo_integrate(lambda u: 4.0 * (u[:, 0] ** 2 + u[:, 1] ** 2 < 1), 100_000, dim=2, rng=0) >>> res.covers(np.pi) True
- financialkit.monte_carlo.multilevel_monte_carlo(S, K, T, r, sigma, q=0.0, option_type='call', eps=0.05, M=2, L_min=2, L_max=10, N0=2000, rng=None, confidence=0.95)[source]#
European option by Giles’s adaptive multilevel Monte Carlo, to root-mean-square error eps.
The algorithm of Giles (2008), Sec. 5: start with levels
0 .. L_minand N0 samples each; set the optimal sample sizes\[N_l = \left\lceil 2\varepsilon^{-2} \sqrt{V_l / C_l} \sum_k \sqrt{V_k C_k} \right\rceil,\]which split \(\varepsilon^2\) evenly between variance and squared bias; top up any level that is short; and once none is, add a level unless the bias estimate \(\max(|\hat Y_L|, |\hat Y_{L-1}| / M) / (M - 1)\) (weak order 1) is below \(\varepsilon / \sqrt 2\).
- Parameters:
S (
float) – Spot, strike, maturity, risk-free rate and volatility.K (
float) – Spot, strike, maturity, risk-free rate and volatility.T (
float) – Spot, strike, maturity, risk-free rate and volatility.r (
float) – Spot, strike, maturity, risk-free rate and volatility.sigma (
float) – Spot, strike, maturity, risk-free rate and volatility.q (
float) – Continuous dividend yield.option_type (
str)eps (
float) – Target root-mean-square error.M (
int) – Refinement factor.L_min (
int) – Smallest and largest finest level.L_max (
int) – Smallest and largest finest level.N0 (
int) – Initial samples on a new level.confidence (
float) – Coverage of the confidence interval (sampling error only).
- Return type:
- Returns:
MLMCResult
- financialkit.monte_carlo.neyman_allocation(weights, stds, n, min_per_stratum=2)[source]#
Neyman’s optimal allocation \(n_h \propto p_h \sigma_h\) of n samples to strata.
It minimises \(\sum_h p_h^2 \sigma_h^2 / n_h\) subject to \(\sum_h n_h = n\). Counts are rounded by largest remainder and floored at min_per_stratum so every stratum keeps a variance estimate.
- Parameters:
- Return type:
- Returns:
ndarray of int, shape (H,) – Sums to n (if
n >= H * min_per_stratum).
Examples
>>> neyman_allocation([0.5, 0.5], [1.0, 3.0], 100) array([25, 75])
- financialkit.monte_carlo.pathwise_greeks(S, K, T, r, sigma, q=0.0, option_type='call', n_paths=100000, rng=None, z=None)[source]#
Price, delta and vega of a vanilla option by pathwise differentiation (Broadie-Glasserman 1996).
The payoff is Lipschitz in \(S\) and \(\sigma\), so the derivative and the expectation commute; its derivative jumps at the strike, so there is no pathwise gamma (returned as
nan). Digital payoffs are discontinuous and are refused.- Parameters:
S (
float) – Spot, strike, maturity, risk-free rate and volatility.K (
float) – Spot, strike, maturity, risk-free rate and volatility.T (
float) – Spot, strike, maturity, risk-free rate and volatility.r (
float) – Spot, strike, maturity, risk-free rate and volatility.sigma (
float) – Spot, strike, maturity, risk-free rate and volatility.q (
float) – Continuous dividend yield.option_type (
str)n_paths (
int) – Number of terminal prices; ignored if z is given.
- Return type:
- Returns:
MonteCarloGreeks
Examples
>>> import numpy as np >>> g = pathwise_greeks(42.0, 40.0, 0.5, 0.10, 0.20, n_paths=100_000, rng=0) >>> abs(g.delta - 0.7791) < 4 * g.delta_se, bool(np.isnan(g.gamma)) (True, True)
- financialkit.monte_carlo.plot_allocation(result, ax=None)[source]#
Samples per stratum next to the within-stratum standard deviations.
- Parameters:
result (
StratifiedResult)ax (matplotlib.axes.Axes, optional)
- Returns:
matplotlib.axes.Axes
- financialkit.monte_carlo.plot_confidence_intervals(results, exact=None, ax=None)[source]#
Confidence intervals of repeated estimates, missed ones in red, with the exact value.
- Parameters:
results (
Sequence[MonteCarloResult])exact (
float|None) – True value; intervals that miss it are highlighted.ax (matplotlib.axes.Axes, optional)
- Returns:
matplotlib.axes.Axes
- financialkit.monte_carlo.plot_error_convergence(n, errors, ax=None, label='', rates=(0.5,))[source]#
Errors (or standard errors) against sample size on log-log axes, with reference slopes.
- Parameters:
n (array_like) – Sample sizes.
errors (array_like) – Absolute errors or standard errors, one per sample size.
ax (matplotlib.axes.Axes, optional) – Axes to draw on; a new figure is created if omitted.
label (
str) – Legend label of the data.rates (
Sequence[float]) – Draw \(n^{-a}\) guides through the first point for eacha.
- Returns:
matplotlib.axes.Axes
- financialkit.monte_carlo.plot_estimator_comparison(results, baseline=None, exact=None, ax=None)[source]#
Estimates with confidence intervals side by side, annotated with variance reduction factors.
- Parameters:
results (
Sequence[MonteCarloResult]) – Labelled by theirmethod.baseline (
MonteCarloResult|None) – Reference forvariance_ratio(); defaults to the first result.exact (
float|None) – True value, drawn as a horizontal line.ax (matplotlib.axes.Axes, optional)
- Returns:
matplotlib.axes.Axes
- financialkit.monte_carlo.plot_mlmc(result, axes=None)[source]#
Giles’s diagnostics: level variances and means (log2) and samples per level.
- Parameters:
result (
MLMCResult)axes (sequence of two matplotlib.axes.Axes, optional)
- Returns:
ndarray of matplotlib.axes.Axes
- financialkit.monte_carlo.plot_point_set(points, ax=None, title='', grid=0, dims=(0, 1))[source]#
A two-dimensional projection of a point set in the unit square.
- financialkit.monte_carlo.randomized_qmc_estimate(f, n, dim, n_replications=16, sequence='sobol', rng=None, confidence=0.95)[source]#
Randomized quasi-Monte Carlo estimate of \(\int_{[0,1)^d} f(u)\,du\) (Owen 1995).
Each replication integrates over an independently scrambled Sobol or Halton point set. Every scrambled point is uniform, so each replication is unbiased and their spread is an honest standard error; for smooth integrands Owen-scrambled nets reach \(O(n^{-3/2})\) root-mean-square error.
- Parameters:
- Return type:
- Returns:
MonteCarloResult
Examples
>>> import numpy as np >>> res = randomized_qmc_estimate(lambda u: np.prod(2 * u, axis=1), 1024, dim=3, rng=0) >>> res.covers(1.0), res.std_error < 1e-3 (True, True)
- financialkit.monte_carlo.sample_mean(samples, confidence=0.95, method='crude', n_evaluations=None)[source]#
The sample mean of i.i.d. draws with its standard error \(s / \sqrt{n}\).
- Parameters:
samples (
ArrayLike) – Independent draws of the quantity whose expectation is estimated.confidence (
float) – Coverage of the confidence interval.method (
str) – Label stored in the result.n_evaluations (
int|None) – Function evaluations behind the draws, if not one each (an antithetic pair or a replication costs several).
- Return type:
- Returns:
MonteCarloResult
Examples
>>> res = sample_mean([1.0, 2.0, 3.0, 4.0]) >>> res.estimate, round(res.std_error, 4) (2.5, 0.6455)
- financialkit.monte_carlo.sobol(n, dim, scramble=True, rng=None)[source]#
npoints of thedim-dimensional Sobol sequence (Sobol 1967), viascipy.stats.qmc.Sobol.The first \(2^m\) points form a \((t, m, s)\)-net: every elementary box of volume \(2^{t-m}\) holds exactly \(2^t\) of them. Use a power of two for n to keep that balance (scipy warns otherwise).
- Parameters:
- Return type:
- Returns:
ndarray, shape (n, dim)
Examples
The first eight unscrambled points stratify each axis into eighths:
>>> import numpy as np >>> x = sobol(8, 2, scramble=False) >>> np.sort(np.floor(8 * x[:, 1])).tolist() [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]
- financialkit.monte_carlo.star_discrepancy(points)[source]#
Exact star discrepancy \(D_N^*\) of a point set in one or two dimensions.
\[D_N^* = \sup_{a \in [0,1]^d} \left| \frac{\#\{x_i \in [0, a)\}}{N} - \prod_j a_j \right|.\]In one dimension \(D_N^* = \frac{1}{2N} + \max_i |x_{(i)} - \frac{2i - 1}{2N}|\) (Niederreiter 1992, Thm. 2.6). In two the supremum is attained on the grid of point coordinates (and 1), which is searched with closed and open boxes in \(O(N^2)\); coordinates are assumed distinct.
Examples
>>> import numpy as np >>> N = 10 >>> round(star_discrepancy((2 * np.arange(1, N + 1) - 1) / (2 * N)), 12) 0.05
- financialkit.monte_carlo.stratified_estimate(f, n, n_strata=10, allocation='neyman', n_pilot=200, rng=None, confidence=0.95)[source]#
Stratified estimate of \(E[f(U)]\), \(U \sim U[0, 1)\), over equal-probability strata.
Stratum \(h\) is \([h/H, (h+1)/H)\). The estimator and its variance are
\[\hat\mu = \sum_h p_h \bar f_h, \qquad \operatorname{Var}\hat\mu = \sum_h p_h^2 \frac{s_h^2}{n_h},\]which never exceeds the crude variance under proportional allocation and is smallest under Neyman’s allocation (Neyman 1934). A one-dimensional uniform drives any Black-Scholes payoff through
norm.ppf, so this stratifies the terminal price.- Parameters:
f (
Callable[[ndarray],ArrayLike]) – Vectorised integrand of uniforms of shape(n,).n (
int) – Samples after the pilot.n_strata (
int) – Number of strata \(H\).allocation (
str) – Neyman’s allocation uses the pilot’s within-stratum standard deviations; proportional gives \(n_h = n / H\).n_pilot (
int) – Pilot samples for Neyman’s allocation, spread evenly over the strata; their cost is included inn_samples.confidence (
float) – Coverage of the confidence interval.
- Return type:
- Returns:
StratifiedResult
- financialkit.monte_carlo.terminal_prices(S, T, r, sigma, z, q=0.0)[source]#
Risk-neutral Black-Scholes prices \(S_T = S e^{(r - q - \sigma^2/2)T + \sigma\sqrt{T} z}\).
- Parameters:
S (
float) – Spot, maturity, risk-free rate and volatility.T (
float) – Spot, maturity, risk-free rate and volatility.r (
float) – Spot, maturity, risk-free rate and volatility.sigma (
float) – Spot, maturity, risk-free rate and volatility.z (
ArrayLike) – Standard normal draws, from any source.q (
float) – Continuous dividend yield.
- Return type:
- Returns:
ndarray, shape (n,)
- financialkit.monte_carlo.van_der_corput(n, base=2, start=1)[source]#
The radical-inverse (van der Corput) sequence in base.
The \(k\)-th point mirrors the base-b digits of \(k\) about the radix point, \(k = \sum_j a_j b^j \mapsto \sum_j a_j b^{-j-1}\), so each new point falls in the largest gap left so far.
- Parameters:
- Return type:
- Returns:
ndarray, shape (n,)
Examples
>>> van_der_corput(7).tolist() [0.5, 0.25, 0.75, 0.125, 0.625, 0.375, 0.875]