Mathematics (wraquant.math)¶
Advanced mathematical tools for quantitative finance: Levy processes, network analysis, optimal stopping, Hawkes processes, numerical methods, spectral analysis, signal processing, information theory, and ergodicity economics.
Quick Example¶
from wraquant.math import levy, network, optimal_stopping
# Simulate a Variance Gamma process
vg_paths = levy.variance_gamma(S0=100, sigma=0.2, theta=-0.1,
nu=0.5, T=1.0, n_paths=1000)
# Minimum spanning tree of asset correlations
mst = network.correlation_mst(returns_df)
print(f"MST edges: {mst['edges']}")
# American option optimal exercise boundary
boundary = optimal_stopping.exercise_boundary(
S_range=(80, 120), K=100, T=1.0, r=0.05, sigma=0.2
)
See also
Pricing (wraquant.price) – Levy process pricing (FFT, COS method)
Volatility Modeling (wraquant.vol) – Hawkes processes for volatility clustering
API Reference¶
Advanced mathematical tools for quantitative finance.
Provides specialized mathematical methods that go beyond standard statistics, serving as building blocks for sophisticated quant models. Covers financial network analysis, Levy process simulation and fitting, and optimal stopping theory – areas where standard libraries fall short of the finance-specific requirements.
Key sub-modules:
Network (
network) – Financial network analysis built on graph theory:correlation_network(construct a graph from a correlation matrix),minimum_spanning_tree(Mantegna’s MST for portfolio clustering),centrality_measures(identify systemically important assets),community_detection(find sector-like groupings),systemic_risk_score(aggregate contagion risk), andgranger_network(causal connectivity via Granger tests).Levy (
levy) – Fat-tailed stochastic processes beyond Gaussian models:variance_gamma_simulate/variance_gamma_pdf(VG process – captures excess kurtosis and skewness),nig_simulate/nig_pdf(Normal Inverse Gaussian),cgmy_simulate(tempered stable),levy_stable_simulate, andfit_variance_gamma/fit_nigfor parameter estimation.Optimal stopping (
optimal_stopping) – Decision-theoretic tools:longstaff_schwartz(American option pricing via regression),binomial_american(binomial lattice),optimal_exit_threshold(when to exit a mean-reverting trade),sequential_probability_ratio(Wald’s SPRT for hypothesis testing),cusum_stopping(cumulative sum changepoint detection), andsecretary_problem_threshold.
Example
>>> from wraquant.math import correlation_network, variance_gamma_simulate
>>> G = correlation_network(returns_df, threshold=0.5)
>>> paths = variance_gamma_simulate(n_paths=1000, n_steps=252, sigma=0.2)
Use wraquant.math for specialized mathematical modeling. For standard
statistical tests and regression, use wraquant.stats. For stochastic
process simulation in a pricing context, see wraquant.price.stochastic.
For network visualization, see wraquant.viz.plotly_network_graph.
- correlation_network(returns_df, threshold=0.3)[source]¶
Build an adjacency matrix from pairwise correlations.
Edges with absolute correlation below threshold are set to zero.
- Parameters:
- Returns:
adjacency– (n, n) adjacency matrix (0/1).correlation– (n, n) full correlation matrix.asset_names– list of column names.- Return type:
Example
>>> import pandas as pd, numpy as np >>> np.random.seed(42) >>> returns = pd.DataFrame(np.random.randn(252, 4), ... columns=['AAPL', 'MSFT', 'GOOG', 'AMZN']) >>> net = correlation_network(returns, threshold=0.1) >>> net['adjacency'].shape (4, 4) >>> net['asset_names'] ['AAPL', 'MSFT', 'GOOG', 'AMZN']
See also
minimum_spanning_treeBuild an MST from the correlation matrix.
granger_networkBuild a directed network from Granger causality.
centrality_measuresRank nodes by importance within the network.
- minimum_spanning_tree(correlation_matrix)[source]¶
Compute the minimum spanning tree of a correlation matrix.
Uses the distance metric
d = sqrt(2 * (1 - rho))and Prim’s algorithm.- Parameters:
correlation_matrix (
TypeAliasType) – (n, n) symmetric correlation matrix.- Returns:
(n, n) adjacency matrix of the MST (symmetric, 0/1 entries).
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.network import minimum_spanning_tree >>> corr = np.array([[1.0, 0.8, 0.3], ... [0.8, 1.0, 0.5], ... [0.3, 0.5, 1.0]]) >>> mst = minimum_spanning_tree(corr) >>> mst.shape (3, 3) >>> int(mst.sum() / 2) # MST of 3 nodes has 2 edges 2
See also
correlation_networkBuild a thresholded correlation network.
centrality_measuresCompute centrality on the resulting MST.
- centrality_measures(adjacency_matrix, asset_names=None)[source]¶
Compute degree, betweenness, and eigenvector centrality.
- Parameters:
- Returns:
degree– degree centrality (normalised). Higher valuesindicate nodes connected to more of the network.
betweenness– betweenness centrality. Higher values indicatenodes that bridge different clusters.
eigenvector– eigenvector centrality. Higher values indicatenodes connected to other well-connected nodes.
asset_names– node labels.- Return type:
Example
>>> import numpy as np >>> from wraquant.math.network import centrality_measures >>> adj = np.array([[0, 1, 1], [1, 0, 0], [1, 0, 0]]) >>> result = centrality_measures(adj, asset_names=['A', 'B', 'C']) >>> result['degree'][0] # node A connected to both B and C 1.0
See also
correlation_networkBuild the adjacency matrix from returns.
systemic_risk_scoreSystemic importance via MES or CoVaR.
- community_detection(adjacency_matrix, n_communities=2)[source]¶
Detect communities via spectral clustering on the graph Laplacian.
- Parameters:
adjacency_matrix (
TypeAliasType) – (n, n) adjacency matrix.n_communities (
int, default:2) – Number of communities to detect (default 2).
- Returns:
Integer community labels of length n.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.network import community_detection >>> # Block-diagonal adjacency (two clear communities) >>> adj = np.array([[0, 1, 1, 0, 0], ... [1, 0, 1, 0, 0], ... [1, 1, 0, 0, 0], ... [0, 0, 0, 0, 1], ... [0, 0, 0, 1, 0]], dtype=float) >>> labels = community_detection(adj, n_communities=2) >>> labels.shape (5,)
See also
correlation_networkBuild a network to cluster.
wraquant.ml.clustering.correlation_clusteringML-based asset clustering.
- systemic_risk_score(returns_df, method='mes', quantile=0.05)[source]¶
Compute a systemic importance score for each asset.
- Parameters:
returns_df (
DataFrame) – Asset returns, columns are assets.method (
str, default:'mes') –Scoring method: -
'mes'– Marginal Expected Shortfall: expected loss of assetconditional on market being in its lower quantile tail.
'covar'– Delta CoVaR: change in system VaR conditional on asset being at its VaR vs median.'connectedness'– Diebold-Yilmaz–style total connectedness from variance decomposition (simplified).
quantile (
float, default:0.05) – Tail quantile for MES / CoVaR (default 0.05).
- Returns:
Systemic risk score for each asset (higher = more systemic).
- Return type:
- Raises:
ValueError – If method is not recognised.
Example
>>> import pandas as pd, numpy as np >>> np.random.seed(42) >>> returns = pd.DataFrame(np.random.randn(252, 3) * 0.01, ... columns=['SPY', 'QQQ', 'IWM']) >>> scores = systemic_risk_score(returns, method='mes', quantile=0.05) >>> len(scores) == 3 True
See also
contagion_simulationSimulate shock propagation through a network.
centrality_measuresNetwork centrality as an alternative importance measure.
wraquant.risk.metricsVaR and CVaR risk metrics.
- contagion_simulation(adjacency_matrix, shock_node, shock_magnitude=1.0, threshold=0.5, max_rounds=100)[source]¶
Simulate shock propagation through a financial network.
At each round, a node that receives total stress above threshold “defaults” and propagates the shock to its neighbours weighted by the adjacency matrix.
- Parameters:
adjacency_matrix (
TypeAliasType) – (n, n) weighted adjacency matrix.shock_node (
int) – Index of the initially shocked node.shock_magnitude (
float, default:1.0) – Size of the initial shock (default 1.0).threshold (
float, default:0.5) – Stress threshold for contagion (default 0.5).max_rounds (
int, default:100) – Maximum propagation rounds (default 100).
- Returns:
defaulted– boolean array of which nodes defaulted.stress– cumulative stress on each node.rounds– number of propagation rounds.cascade_size– total number of defaults. A cascade_sizeclose to n indicates high systemic fragility.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.network import contagion_simulation >>> adj = np.array([[0, 0.6, 0.3], ... [0.4, 0, 0.5], ... [0.2, 0.3, 0]], dtype=float) >>> result = contagion_simulation(adj, shock_node=0, shock_magnitude=1.0, ... threshold=0.5) >>> result['defaulted'][0] True >>> result['cascade_size'] >= 1 True
See also
systemic_risk_scoreIdentify which nodes are most systemically important.
correlation_networkBuild the adjacency matrix from return data.
- granger_network(returns_df, max_lag=5, significance=0.05)[source]¶
Build a directed network from pairwise Granger causality tests.
For each pair (i, j) test whether lagged values of j help predict i beyond i’s own lags, using an F-test.
- Parameters:
- Returns:
adjacency– (n, n) directed adjacency matrix; entry (i, j) = 1means j Granger-causes i.
pvalues– (n, n) matrix of p-values.asset_names– column names.- Return type:
Example
>>> import pandas as pd, numpy as np >>> np.random.seed(42) >>> x = np.random.randn(200) >>> y = np.concatenate([[0], x[:-1]]) + np.random.randn(200) * 0.1 >>> df = pd.DataFrame({'X': x, 'Y': y}) >>> result = granger_network(df, max_lag=3, significance=0.05) >>> result['adjacency'].shape (2, 2)
See also
correlation_networkUndirected network from correlations.
wraquant.math.information.transfer_entropyInformation-theoretic causality.
- variance_gamma_pdf(x, sigma, nu, theta)[source]¶
Evaluate the Variance Gamma probability density function.
\[f(x) = \frac{2\,e^{\theta\,x / \sigma^2}} {\nu^{1/\nu}\,\sqrt{2\pi}\,\sigma\,\Gamma(1/\nu)} \left(\frac{x^2}{2\sigma^2/\nu + \theta^2}\right)^{1/(2\nu) - 1/4} K_{1/\nu - 1/2}\!\left(\frac{1}{\sigma^2} \sqrt{x^2(2\sigma^2/\nu + \theta^2)}\right)\]- Parameters:
- Returns:
PDF values.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.levy import variance_gamma_pdf >>> x = np.linspace(-0.05, 0.05, 5) >>> pdf = variance_gamma_pdf(x, sigma=0.01, nu=0.5, theta=-0.001) >>> pdf.shape (5,) >>> pdf[2] > pdf[0] # peak near center True
See also
variance_gamma_simulateSimulate paths from a VG process.
fit_variance_gammaCalibrate VG parameters to return data.
nig_pdfNormal Inverse Gaussian PDF (alternative fat-tailed model).
- variance_gamma_simulate(sigma, nu, theta, n_steps, dt=0.003968253968253968, seed=None)[source]¶
Simulate a Variance Gamma process via time-changed Brownian motion.
X(t) = theta * G(t) + sigma * W(G(t))
where G is a Gamma subordinator with shape = dt/nu and scale = nu.
- Parameters:
- Returns:
Cumulative VG process values of length n_steps + 1 (starts at 0).
- Return type:
Example
>>> from wraquant.math.levy import variance_gamma_simulate >>> path = variance_gamma_simulate(sigma=0.01, nu=0.5, theta=-0.001, ... n_steps=252, seed=42) >>> path.shape (253,) >>> path[0] 0.0
See also
variance_gamma_pdfEvaluate the VG density.
nig_simulateSimulate a Normal Inverse Gaussian process.
cgmy_simulateSimulate a CGMY process (generalises VG).
- nig_pdf(x, alpha, beta, mu, delta)[source]¶
Evaluate the Normal Inverse Gaussian probability density function.
\[f(x) = \frac{\alpha\,\delta}{\pi} \exp\!\left(\delta\sqrt{\alpha^2 - \beta^2} + \beta(x - \mu)\right) \frac{K_1\!\left(\alpha\sqrt{\delta^2 + (x-\mu)^2}\right)} {\sqrt{\delta^2 + (x - \mu)^2}}\]- Parameters:
- Returns:
PDF values.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.levy import nig_pdf >>> x = np.linspace(-0.05, 0.05, 5) >>> pdf = nig_pdf(x, alpha=50.0, beta=-5.0, mu=0.0, delta=0.01) >>> pdf.shape (5,) >>> (pdf >= 0).all() True
See also
nig_simulateSimulate paths from an NIG process.
fit_nigCalibrate NIG parameters to return data.
variance_gamma_pdfVG density (alternative fat-tailed model).
- nig_simulate(alpha, beta, mu, delta, n_steps, dt=0.003968253968253968, seed=None)[source]¶
Simulate a Normal Inverse Gaussian process.
NIG is obtained as a normal variance-mean mixture with an Inverse Gaussian subordinator.
- Parameters:
- Returns:
Cumulative NIG process of length n_steps + 1 (starts at 0).
- Return type:
Example
>>> from wraquant.math.levy import nig_simulate >>> path = nig_simulate(alpha=50.0, beta=-5.0, mu=0.0, delta=0.01, ... n_steps=252, seed=42) >>> path.shape (253,) >>> path[0] 0.0
See also
nig_pdfEvaluate the NIG density.
variance_gamma_simulateSimulate a VG process.
- cgmy_simulate(C, G, M, Y, n_steps, dt=0.003968253968253968, seed=None)[source]¶
Simulate a CGMY process via series representation (shot-noise).
For Y < 0 the process has finite activity and can be approximated by compound Poisson. For 0 < Y < 2 the process has infinite activity and is approximated by truncating the small-jump component at epsilon and adding a Brownian correction.
- Parameters:
C (
float) – Overall activity level (> 0).G (
float) – Rate of exponential decay for negative jumps (> 0).M (
float) – Rate of exponential decay for positive jumps (> 0).Y (
float) – Fine structure parameter (Y < 2).n_steps (
int) – Number of time steps.dt (
float, default:0.003968253968253968) – Time increment (default 1/252).
- Returns:
Cumulative CGMY process of length n_steps + 1 (starts at 0).
- Return type:
Example
>>> from wraquant.math.levy import cgmy_simulate >>> path = cgmy_simulate(C=1.0, G=5.0, M=10.0, Y=0.5, ... n_steps=252, seed=42) >>> path.shape (253,) >>> path[0] 0.0
Notes
The CGMY process generalises the Variance Gamma process (VG corresponds to Y = 0). Higher Y values produce more small jumps; Y < 0 yields finite activity (compound Poisson).
See also
variance_gamma_simulateSimulate VG (special case Y = 0).
levy_stable_simulateSimulate a stable Levy process.
- fit_variance_gamma(returns)[source]¶
Fit a Variance Gamma distribution to return data via MLE.
- Parameters:
returns (
TypeAliasType) – Observed returns.- Returns:
sigma– fitted volatility.nu– fitted variance rate.theta– fitted drift.log_likelihood– maximised log-likelihood.- Return type:
Example
>>> import numpy as np >>> from wraquant.math.levy import fit_variance_gamma, variance_gamma_simulate >>> # Generate synthetic VG returns and re-fit >>> path = variance_gamma_simulate(sigma=0.01, nu=0.5, theta=-0.001, ... n_steps=2000, seed=42) >>> returns = np.diff(path) >>> result = fit_variance_gamma(returns) >>> result['sigma'] > 0 True
See also
variance_gamma_pdfEvaluate the fitted density.
fit_nigFit an NIG distribution instead.
- fit_nig(returns)[source]¶
Fit a Normal Inverse Gaussian distribution to return data via MLE.
- Parameters:
returns (
TypeAliasType) – Observed returns.- Returns:
alpha– fitted tail parameter.beta– fitted skewness parameter.mu– fitted location.delta– fitted scale.log_likelihood– maximised log-likelihood.- Return type:
Example
>>> import numpy as np >>> from wraquant.math.levy import fit_nig >>> rng = np.random.default_rng(42) >>> returns = rng.standard_t(df=5, size=1000) * 0.01 >>> result = fit_nig(returns) >>> result['alpha'] > 0 and result['delta'] > 0 True
See also
nig_pdfEvaluate the fitted density.
fit_variance_gammaFit a VG distribution instead.
- levy_stable_simulate(alpha, beta, n_steps, seed=None)[source]¶
Simulate a stable Lévy process using the Chambers-Mallows-Stuck method.
- Parameters:
- Returns:
Cumulative process values of length n_steps + 1 (starts at 0).
- Return type:
Example
>>> from wraquant.math.levy import levy_stable_simulate >>> path = levy_stable_simulate(alpha=1.7, beta=0.0, n_steps=500, seed=42) >>> path.shape (501,) >>> path[0] 0.0
Notes
alpha = 2 recovers Brownian motion; alpha = 1 with beta = 0 yields a Cauchy process. Smaller alpha produces heavier tails.
See also
variance_gamma_simulateVG process (finite variance, unlike stable).
cgmy_simulateCGMY process (tempered stable, finite moments).
- characteristic_function_vg(u, sigma, nu, theta)[source]¶
Characteristic function of the Variance Gamma distribution at time 1.
\[\phi(u) = \left(1 - i\,\theta\,\nu\,u + \tfrac{1}{2}\sigma^2\nu\,u^2\right)^{-1/\nu}\]- Parameters:
- Returns:
Complex-valued characteristic function values.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.levy import characteristic_function_vg >>> phi = characteristic_function_vg(np.array([0.0, 1.0]), sigma=0.01, ... nu=0.5, theta=-0.001) >>> abs(phi[0]) # phi(0) = 1 for any characteristic function 1.0
See also
variance_gamma_pdfDensity evaluation (Fourier-inversion counterpart).
- longstaff_schwartz(paths, strike, rf_rate, dt, basis_functions='laguerre', option_type='put')[source]¶
Longstaff-Schwartz algorithm for American option pricing.
- Parameters:
paths (
TypeAliasType) – Simulated price paths of shape(n_paths, n_steps + 1). Column 0 is the initial price.strike (
float) – Strike price.rf_rate (
float) – Risk-free rate (annualised, continuously compounded).dt (
float) – Time step (in years).basis_functions (
str, default:'laguerre') – Basis for the continuation-value regression (default'laguerre').option_type (
str, default:'put') – Option type (default'put').
- Returns:
Estimated American option price.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.optimal_stopping import longstaff_schwartz >>> rng = np.random.default_rng(42) >>> S0, K, r, vol, T, n_steps = 100, 100, 0.05, 0.2, 1.0, 50 >>> dt = T / n_steps >>> # Simulate GBM paths >>> z = rng.standard_normal((10000, n_steps)) >>> paths = np.zeros((10000, n_steps + 1)) >>> paths[:, 0] = S0 >>> for i in range(n_steps): ... paths[:, i+1] = paths[:, i] * np.exp((r - 0.5*vol**2)*dt + vol*np.sqrt(dt)*z[:, i]) >>> price = longstaff_schwartz(paths, strike=K, rf_rate=r, dt=dt) >>> price > 0 True
Notes
Reference: Longstaff, F. A. & Schwartz, E. S. (2001). “Valuing American Options by Simulation.” Review of Financial Studies, 14(1), 113-147.
See also
binomial_americanBinomial tree pricing (exact, but slower for many steps).
wraquant.priceOptions pricing module with Black-Scholes and SDEs.
- binomial_american(spot, strike, rf_rate, vol, T, n_steps, option_type='put')[source]¶
Price an American option using the CRR binomial tree.
- Parameters:
spot (
float) – Current underlying price.strike (
float) – Strike price.rf_rate (
float) – Risk-free rate (annualised).vol (
float) – Volatility (annualised).T (
float) – Time to expiration (years).n_steps (
int) – Number of tree steps.option_type (
str, default:'put') – Option type (default'put').
- Returns:
American option price.
- Return type:
Example
>>> from wraquant.math.optimal_stopping import binomial_american >>> price = binomial_american(spot=100, strike=100, rf_rate=0.05, ... vol=0.2, T=1.0, n_steps=200) >>> price > 0 True >>> # American put is worth at least European put >>> price >= 4.0 # approximate lower bound True
See also
longstaff_schwartzMonte Carlo pricing for path-dependent payoffs.
- optimal_exit_threshold(mu, sigma, transaction_cost)[source]¶
Compute the optimal exit threshold for an OU mean-reverting process.
For a process \(dX = -\mu\,X\,dt + \sigma\,dW\), the optimal exit threshold balances expected profit against the cost of trading.
Uses the analytical approximation:
\[x^* \approx \sigma \sqrt{\frac{2}{\mu}} \ln\!\left(\frac{\sigma}{\sqrt{2\mu}\,c}\right)\]when the threshold is large relative to transaction_cost.
- Parameters:
- Returns:
entry_threshold– optimal entry (in units of process value).exit_threshold– optimal exit threshold. Higher valuesmean waiting for a larger dislocation before exiting.
expected_profit– approximate expected profit per trade.- Return type:
Example
>>> from wraquant.math.optimal_stopping import optimal_exit_threshold >>> result = optimal_exit_threshold(mu=5.0, sigma=0.5, ... transaction_cost=0.001) >>> result['exit_threshold'] > 0 True >>> result['expected_profit'] > 0 True
See also
wraquant.stats.cointegrationTest for cointegration before pairs trading.
wraquant.execution.optimalOptimal execution with Almgren-Chriss.
- sequential_probability_ratio(observations, h0_dist, h1_dist, alpha=0.05, beta=0.05)[source]¶
Sequential Probability Ratio Test (SPRT).
Sequentially tests H0 vs H1 with controlled error probabilities.
- Parameters:
observations (
TypeAliasType) – Sequentially observed data.h0_dist (
tuple[str,dict[str,float]]) –(distribution_name, params)for H0. Supported:('normal', {'mu': ..., 'sigma': ...}).h1_dist (
tuple[str,dict[str,float]]) –(distribution_name, params)for H1.alpha (
float, default:0.05) – Type I error probability (default 0.05).beta (
float, default:0.05) – Type II error probability (default 0.05).
- Returns:
decision–'reject_h0','accept_h0', or'inconclusive'.stopping_time– index at which the decision was reached. Fewerobservations needed indicates stronger signal.
log_ratio– final log-likelihood ratio.- Return type:
Example
>>> import numpy as np >>> from wraquant.math.optimal_stopping import sequential_probability_ratio >>> rng = np.random.default_rng(42) >>> # Generate data from H1 (mu=0.5) and test against H0 (mu=0) >>> obs = rng.normal(0.5, 1.0, size=100) >>> result = sequential_probability_ratio( ... obs, ... h0_dist=('normal', {'mu': 0.0, 'sigma': 1.0}), ... h1_dist=('normal', {'mu': 0.5, 'sigma': 1.0}), ... ) >>> result['decision'] 'reject_h0'
See also
cusum_stoppingCUSUM-based change detection (non-parametric alternative).
- cusum_stopping(observations, target_mean, threshold)[source]¶
CUSUM stopping rule for detecting a mean shift.
Monitors cumulative sum of deviations from target_mean and triggers a stop when the CUSUM statistic exceeds threshold.
- Parameters:
- Returns:
detected– whether a change was detected.stopping_time– index of detection (or length of data if none).cusum_pos– final positive CUSUM statistic.cusum_neg– final negative CUSUM statistic.cusum_values– array of max(cusum_pos, |cusum_neg|) at each step.- Return type:
Example
>>> import numpy as np >>> from wraquant.math.optimal_stopping import cusum_stopping >>> # Returns with a mean shift at index 50 >>> obs = np.concatenate([np.random.randn(50) * 0.01, ... np.random.randn(50) * 0.01 + 0.05]) >>> result = cusum_stopping(obs, target_mean=0.0, threshold=0.5) >>> result['detected'] True >>> result['stopping_time'] > 50 # detected after the shift True
See also
sequential_probability_ratioParametric sequential test.
wraquant.ts.changepointsChange-point detection for time series.
- secretary_problem_threshold(n_candidates)[source]¶
Compute the optimal stopping threshold for the secretary problem.
The classical 1/e rule: reject the first r - 1 candidates, then hire the first one who is better than all previously seen.
- Parameters:
n_candidates (
int) – Total number of candidates.- Returns:
threshold– number of candidates to unconditionally reject.optimal_fraction– threshold / n_candidates. Converges to1/e (approximately 0.368) as n_candidates grows.
success_probability– probability of selecting the best candidate.Converges to 1/e for large n_candidates.
- Return type:
Example
>>> from wraquant.math.optimal_stopping import secretary_problem_threshold >>> result = secretary_problem_threshold(100) >>> 30 <= result['threshold'] <= 40 # close to 1/e * 100 ~ 37 True >>> result['success_probability'] > 0.3 True
See also
sequential_probability_ratioSequential testing framework.
Levy Processes¶
Lévy processes for fat-tailed asset models.
Provides simulation, density evaluation, and calibration for: - Variance Gamma (VG) - Normal Inverse Gaussian (NIG) - CGMY - Stable Lévy processes
All implementations are pure numpy/scipy.
- variance_gamma_pdf(x, sigma, nu, theta)[source]¶
Evaluate the Variance Gamma probability density function.
\[f(x) = \frac{2\,e^{\theta\,x / \sigma^2}} {\nu^{1/\nu}\,\sqrt{2\pi}\,\sigma\,\Gamma(1/\nu)} \left(\frac{x^2}{2\sigma^2/\nu + \theta^2}\right)^{1/(2\nu) - 1/4} K_{1/\nu - 1/2}\!\left(\frac{1}{\sigma^2} \sqrt{x^2(2\sigma^2/\nu + \theta^2)}\right)\]- Parameters:
- Returns:
PDF values.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.levy import variance_gamma_pdf >>> x = np.linspace(-0.05, 0.05, 5) >>> pdf = variance_gamma_pdf(x, sigma=0.01, nu=0.5, theta=-0.001) >>> pdf.shape (5,) >>> pdf[2] > pdf[0] # peak near center True
See also
variance_gamma_simulateSimulate paths from a VG process.
fit_variance_gammaCalibrate VG parameters to return data.
nig_pdfNormal Inverse Gaussian PDF (alternative fat-tailed model).
- variance_gamma_simulate(sigma, nu, theta, n_steps, dt=0.003968253968253968, seed=None)[source]¶
Simulate a Variance Gamma process via time-changed Brownian motion.
X(t) = theta * G(t) + sigma * W(G(t))
where G is a Gamma subordinator with shape = dt/nu and scale = nu.
- Parameters:
- Returns:
Cumulative VG process values of length n_steps + 1 (starts at 0).
- Return type:
Example
>>> from wraquant.math.levy import variance_gamma_simulate >>> path = variance_gamma_simulate(sigma=0.01, nu=0.5, theta=-0.001, ... n_steps=252, seed=42) >>> path.shape (253,) >>> path[0] 0.0
See also
variance_gamma_pdfEvaluate the VG density.
nig_simulateSimulate a Normal Inverse Gaussian process.
cgmy_simulateSimulate a CGMY process (generalises VG).
- nig_pdf(x, alpha, beta, mu, delta)[source]¶
Evaluate the Normal Inverse Gaussian probability density function.
\[f(x) = \frac{\alpha\,\delta}{\pi} \exp\!\left(\delta\sqrt{\alpha^2 - \beta^2} + \beta(x - \mu)\right) \frac{K_1\!\left(\alpha\sqrt{\delta^2 + (x-\mu)^2}\right)} {\sqrt{\delta^2 + (x - \mu)^2}}\]- Parameters:
- Returns:
PDF values.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.levy import nig_pdf >>> x = np.linspace(-0.05, 0.05, 5) >>> pdf = nig_pdf(x, alpha=50.0, beta=-5.0, mu=0.0, delta=0.01) >>> pdf.shape (5,) >>> (pdf >= 0).all() True
See also
nig_simulateSimulate paths from an NIG process.
fit_nigCalibrate NIG parameters to return data.
variance_gamma_pdfVG density (alternative fat-tailed model).
- nig_simulate(alpha, beta, mu, delta, n_steps, dt=0.003968253968253968, seed=None)[source]¶
Simulate a Normal Inverse Gaussian process.
NIG is obtained as a normal variance-mean mixture with an Inverse Gaussian subordinator.
- Parameters:
- Returns:
Cumulative NIG process of length n_steps + 1 (starts at 0).
- Return type:
Example
>>> from wraquant.math.levy import nig_simulate >>> path = nig_simulate(alpha=50.0, beta=-5.0, mu=0.0, delta=0.01, ... n_steps=252, seed=42) >>> path.shape (253,) >>> path[0] 0.0
See also
nig_pdfEvaluate the NIG density.
variance_gamma_simulateSimulate a VG process.
- cgmy_simulate(C, G, M, Y, n_steps, dt=0.003968253968253968, seed=None)[source]¶
Simulate a CGMY process via series representation (shot-noise).
For Y < 0 the process has finite activity and can be approximated by compound Poisson. For 0 < Y < 2 the process has infinite activity and is approximated by truncating the small-jump component at epsilon and adding a Brownian correction.
- Parameters:
C (
float) – Overall activity level (> 0).G (
float) – Rate of exponential decay for negative jumps (> 0).M (
float) – Rate of exponential decay for positive jumps (> 0).Y (
float) – Fine structure parameter (Y < 2).n_steps (
int) – Number of time steps.dt (
float, default:0.003968253968253968) – Time increment (default 1/252).
- Returns:
Cumulative CGMY process of length n_steps + 1 (starts at 0).
- Return type:
Example
>>> from wraquant.math.levy import cgmy_simulate >>> path = cgmy_simulate(C=1.0, G=5.0, M=10.0, Y=0.5, ... n_steps=252, seed=42) >>> path.shape (253,) >>> path[0] 0.0
Notes
The CGMY process generalises the Variance Gamma process (VG corresponds to Y = 0). Higher Y values produce more small jumps; Y < 0 yields finite activity (compound Poisson).
See also
variance_gamma_simulateSimulate VG (special case Y = 0).
levy_stable_simulateSimulate a stable Levy process.
- fit_variance_gamma(returns)[source]¶
Fit a Variance Gamma distribution to return data via MLE.
- Parameters:
returns (
TypeAliasType) – Observed returns.- Returns:
sigma– fitted volatility.nu– fitted variance rate.theta– fitted drift.log_likelihood– maximised log-likelihood.- Return type:
Example
>>> import numpy as np >>> from wraquant.math.levy import fit_variance_gamma, variance_gamma_simulate >>> # Generate synthetic VG returns and re-fit >>> path = variance_gamma_simulate(sigma=0.01, nu=0.5, theta=-0.001, ... n_steps=2000, seed=42) >>> returns = np.diff(path) >>> result = fit_variance_gamma(returns) >>> result['sigma'] > 0 True
See also
variance_gamma_pdfEvaluate the fitted density.
fit_nigFit an NIG distribution instead.
- fit_nig(returns)[source]¶
Fit a Normal Inverse Gaussian distribution to return data via MLE.
- Parameters:
returns (
TypeAliasType) – Observed returns.- Returns:
alpha– fitted tail parameter.beta– fitted skewness parameter.mu– fitted location.delta– fitted scale.log_likelihood– maximised log-likelihood.- Return type:
Example
>>> import numpy as np >>> from wraquant.math.levy import fit_nig >>> rng = np.random.default_rng(42) >>> returns = rng.standard_t(df=5, size=1000) * 0.01 >>> result = fit_nig(returns) >>> result['alpha'] > 0 and result['delta'] > 0 True
See also
nig_pdfEvaluate the fitted density.
fit_variance_gammaFit a VG distribution instead.
- levy_stable_simulate(alpha, beta, n_steps, seed=None)[source]¶
Simulate a stable Lévy process using the Chambers-Mallows-Stuck method.
- Parameters:
- Returns:
Cumulative process values of length n_steps + 1 (starts at 0).
- Return type:
Example
>>> from wraquant.math.levy import levy_stable_simulate >>> path = levy_stable_simulate(alpha=1.7, beta=0.0, n_steps=500, seed=42) >>> path.shape (501,) >>> path[0] 0.0
Notes
alpha = 2 recovers Brownian motion; alpha = 1 with beta = 0 yields a Cauchy process. Smaller alpha produces heavier tails.
See also
variance_gamma_simulateVG process (finite variance, unlike stable).
cgmy_simulateCGMY process (tempered stable, finite moments).
- characteristic_function_vg(u, sigma, nu, theta)[source]¶
Characteristic function of the Variance Gamma distribution at time 1.
\[\phi(u) = \left(1 - i\,\theta\,\nu\,u + \tfrac{1}{2}\sigma^2\nu\,u^2\right)^{-1/\nu}\]- Parameters:
- Returns:
Complex-valued characteristic function values.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.levy import characteristic_function_vg >>> phi = characteristic_function_vg(np.array([0.0, 1.0]), sigma=0.01, ... nu=0.5, theta=-0.001) >>> abs(phi[0]) # phi(0) = 1 for any characteristic function 1.0
See also
variance_gamma_pdfDensity evaluation (Fourier-inversion counterpart).
Network Analysis¶
Financial network analysis.
Build and analyse networks from correlation structures, Granger causality, and other dependence measures. Includes systemic risk scoring and contagion simulation.
- correlation_network(returns_df, threshold=0.3)[source]¶
Build an adjacency matrix from pairwise correlations.
Edges with absolute correlation below threshold are set to zero.
- Parameters:
- Returns:
adjacency– (n, n) adjacency matrix (0/1).correlation– (n, n) full correlation matrix.asset_names– list of column names.- Return type:
Example
>>> import pandas as pd, numpy as np >>> np.random.seed(42) >>> returns = pd.DataFrame(np.random.randn(252, 4), ... columns=['AAPL', 'MSFT', 'GOOG', 'AMZN']) >>> net = correlation_network(returns, threshold=0.1) >>> net['adjacency'].shape (4, 4) >>> net['asset_names'] ['AAPL', 'MSFT', 'GOOG', 'AMZN']
See also
minimum_spanning_treeBuild an MST from the correlation matrix.
granger_networkBuild a directed network from Granger causality.
centrality_measuresRank nodes by importance within the network.
- minimum_spanning_tree(correlation_matrix)[source]¶
Compute the minimum spanning tree of a correlation matrix.
Uses the distance metric
d = sqrt(2 * (1 - rho))and Prim’s algorithm.- Parameters:
correlation_matrix (
TypeAliasType) – (n, n) symmetric correlation matrix.- Returns:
(n, n) adjacency matrix of the MST (symmetric, 0/1 entries).
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.network import minimum_spanning_tree >>> corr = np.array([[1.0, 0.8, 0.3], ... [0.8, 1.0, 0.5], ... [0.3, 0.5, 1.0]]) >>> mst = minimum_spanning_tree(corr) >>> mst.shape (3, 3) >>> int(mst.sum() / 2) # MST of 3 nodes has 2 edges 2
See also
correlation_networkBuild a thresholded correlation network.
centrality_measuresCompute centrality on the resulting MST.
- centrality_measures(adjacency_matrix, asset_names=None)[source]¶
Compute degree, betweenness, and eigenvector centrality.
- Parameters:
- Returns:
degree– degree centrality (normalised). Higher valuesindicate nodes connected to more of the network.
betweenness– betweenness centrality. Higher values indicatenodes that bridge different clusters.
eigenvector– eigenvector centrality. Higher values indicatenodes connected to other well-connected nodes.
asset_names– node labels.- Return type:
Example
>>> import numpy as np >>> from wraquant.math.network import centrality_measures >>> adj = np.array([[0, 1, 1], [1, 0, 0], [1, 0, 0]]) >>> result = centrality_measures(adj, asset_names=['A', 'B', 'C']) >>> result['degree'][0] # node A connected to both B and C 1.0
See also
correlation_networkBuild the adjacency matrix from returns.
systemic_risk_scoreSystemic importance via MES or CoVaR.
- community_detection(adjacency_matrix, n_communities=2)[source]¶
Detect communities via spectral clustering on the graph Laplacian.
- Parameters:
adjacency_matrix (
TypeAliasType) – (n, n) adjacency matrix.n_communities (
int, default:2) – Number of communities to detect (default 2).
- Returns:
Integer community labels of length n.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.network import community_detection >>> # Block-diagonal adjacency (two clear communities) >>> adj = np.array([[0, 1, 1, 0, 0], ... [1, 0, 1, 0, 0], ... [1, 1, 0, 0, 0], ... [0, 0, 0, 0, 1], ... [0, 0, 0, 1, 0]], dtype=float) >>> labels = community_detection(adj, n_communities=2) >>> labels.shape (5,)
See also
correlation_networkBuild a network to cluster.
wraquant.ml.clustering.correlation_clusteringML-based asset clustering.
- systemic_risk_score(returns_df, method='mes', quantile=0.05)[source]¶
Compute a systemic importance score for each asset.
- Parameters:
returns_df (
DataFrame) – Asset returns, columns are assets.method (
str, default:'mes') –Scoring method: -
'mes'– Marginal Expected Shortfall: expected loss of assetconditional on market being in its lower quantile tail.
'covar'– Delta CoVaR: change in system VaR conditional on asset being at its VaR vs median.'connectedness'– Diebold-Yilmaz–style total connectedness from variance decomposition (simplified).
quantile (
float, default:0.05) – Tail quantile for MES / CoVaR (default 0.05).
- Returns:
Systemic risk score for each asset (higher = more systemic).
- Return type:
- Raises:
ValueError – If method is not recognised.
Example
>>> import pandas as pd, numpy as np >>> np.random.seed(42) >>> returns = pd.DataFrame(np.random.randn(252, 3) * 0.01, ... columns=['SPY', 'QQQ', 'IWM']) >>> scores = systemic_risk_score(returns, method='mes', quantile=0.05) >>> len(scores) == 3 True
See also
contagion_simulationSimulate shock propagation through a network.
centrality_measuresNetwork centrality as an alternative importance measure.
wraquant.risk.metricsVaR and CVaR risk metrics.
- contagion_simulation(adjacency_matrix, shock_node, shock_magnitude=1.0, threshold=0.5, max_rounds=100)[source]¶
Simulate shock propagation through a financial network.
At each round, a node that receives total stress above threshold “defaults” and propagates the shock to its neighbours weighted by the adjacency matrix.
- Parameters:
adjacency_matrix (
TypeAliasType) – (n, n) weighted adjacency matrix.shock_node (
int) – Index of the initially shocked node.shock_magnitude (
float, default:1.0) – Size of the initial shock (default 1.0).threshold (
float, default:0.5) – Stress threshold for contagion (default 0.5).max_rounds (
int, default:100) – Maximum propagation rounds (default 100).
- Returns:
defaulted– boolean array of which nodes defaulted.stress– cumulative stress on each node.rounds– number of propagation rounds.cascade_size– total number of defaults. A cascade_sizeclose to n indicates high systemic fragility.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.network import contagion_simulation >>> adj = np.array([[0, 0.6, 0.3], ... [0.4, 0, 0.5], ... [0.2, 0.3, 0]], dtype=float) >>> result = contagion_simulation(adj, shock_node=0, shock_magnitude=1.0, ... threshold=0.5) >>> result['defaulted'][0] True >>> result['cascade_size'] >= 1 True
See also
systemic_risk_scoreIdentify which nodes are most systemically important.
correlation_networkBuild the adjacency matrix from return data.
- granger_network(returns_df, max_lag=5, significance=0.05)[source]¶
Build a directed network from pairwise Granger causality tests.
For each pair (i, j) test whether lagged values of j help predict i beyond i’s own lags, using an F-test.
- Parameters:
- Returns:
adjacency– (n, n) directed adjacency matrix; entry (i, j) = 1means j Granger-causes i.
pvalues– (n, n) matrix of p-values.asset_names– column names.- Return type:
Example
>>> import pandas as pd, numpy as np >>> np.random.seed(42) >>> x = np.random.randn(200) >>> y = np.concatenate([[0], x[:-1]]) + np.random.randn(200) * 0.1 >>> df = pd.DataFrame({'X': x, 'Y': y}) >>> result = granger_network(df, max_lag=3, significance=0.05) >>> result['adjacency'].shape (2, 2)
See also
correlation_networkUndirected network from correlations.
wraquant.math.information.transfer_entropyInformation-theoretic causality.
Optimal Stopping¶
Optimal stopping theory for finance.
Longstaff-Schwartz American option pricing, binomial American options, optimal exit thresholds, sequential testing, and change-point detection.
- longstaff_schwartz(paths, strike, rf_rate, dt, basis_functions='laguerre', option_type='put')[source]¶
Longstaff-Schwartz algorithm for American option pricing.
- Parameters:
paths (
TypeAliasType) – Simulated price paths of shape(n_paths, n_steps + 1). Column 0 is the initial price.strike (
float) – Strike price.rf_rate (
float) – Risk-free rate (annualised, continuously compounded).dt (
float) – Time step (in years).basis_functions (
str, default:'laguerre') – Basis for the continuation-value regression (default'laguerre').option_type (
str, default:'put') – Option type (default'put').
- Returns:
Estimated American option price.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.optimal_stopping import longstaff_schwartz >>> rng = np.random.default_rng(42) >>> S0, K, r, vol, T, n_steps = 100, 100, 0.05, 0.2, 1.0, 50 >>> dt = T / n_steps >>> # Simulate GBM paths >>> z = rng.standard_normal((10000, n_steps)) >>> paths = np.zeros((10000, n_steps + 1)) >>> paths[:, 0] = S0 >>> for i in range(n_steps): ... paths[:, i+1] = paths[:, i] * np.exp((r - 0.5*vol**2)*dt + vol*np.sqrt(dt)*z[:, i]) >>> price = longstaff_schwartz(paths, strike=K, rf_rate=r, dt=dt) >>> price > 0 True
Notes
Reference: Longstaff, F. A. & Schwartz, E. S. (2001). “Valuing American Options by Simulation.” Review of Financial Studies, 14(1), 113-147.
See also
binomial_americanBinomial tree pricing (exact, but slower for many steps).
wraquant.priceOptions pricing module with Black-Scholes and SDEs.
- binomial_american(spot, strike, rf_rate, vol, T, n_steps, option_type='put')[source]¶
Price an American option using the CRR binomial tree.
- Parameters:
spot (
float) – Current underlying price.strike (
float) – Strike price.rf_rate (
float) – Risk-free rate (annualised).vol (
float) – Volatility (annualised).T (
float) – Time to expiration (years).n_steps (
int) – Number of tree steps.option_type (
str, default:'put') – Option type (default'put').
- Returns:
American option price.
- Return type:
Example
>>> from wraquant.math.optimal_stopping import binomial_american >>> price = binomial_american(spot=100, strike=100, rf_rate=0.05, ... vol=0.2, T=1.0, n_steps=200) >>> price > 0 True >>> # American put is worth at least European put >>> price >= 4.0 # approximate lower bound True
See also
longstaff_schwartzMonte Carlo pricing for path-dependent payoffs.
- optimal_exit_threshold(mu, sigma, transaction_cost)[source]¶
Compute the optimal exit threshold for an OU mean-reverting process.
For a process \(dX = -\mu\,X\,dt + \sigma\,dW\), the optimal exit threshold balances expected profit against the cost of trading.
Uses the analytical approximation:
\[x^* \approx \sigma \sqrt{\frac{2}{\mu}} \ln\!\left(\frac{\sigma}{\sqrt{2\mu}\,c}\right)\]when the threshold is large relative to transaction_cost.
- Parameters:
- Returns:
entry_threshold– optimal entry (in units of process value).exit_threshold– optimal exit threshold. Higher valuesmean waiting for a larger dislocation before exiting.
expected_profit– approximate expected profit per trade.- Return type:
Example
>>> from wraquant.math.optimal_stopping import optimal_exit_threshold >>> result = optimal_exit_threshold(mu=5.0, sigma=0.5, ... transaction_cost=0.001) >>> result['exit_threshold'] > 0 True >>> result['expected_profit'] > 0 True
See also
wraquant.stats.cointegrationTest for cointegration before pairs trading.
wraquant.execution.optimalOptimal execution with Almgren-Chriss.
- sequential_probability_ratio(observations, h0_dist, h1_dist, alpha=0.05, beta=0.05)[source]¶
Sequential Probability Ratio Test (SPRT).
Sequentially tests H0 vs H1 with controlled error probabilities.
- Parameters:
observations (
TypeAliasType) – Sequentially observed data.h0_dist (
tuple[str,dict[str,float]]) –(distribution_name, params)for H0. Supported:('normal', {'mu': ..., 'sigma': ...}).h1_dist (
tuple[str,dict[str,float]]) –(distribution_name, params)for H1.alpha (
float, default:0.05) – Type I error probability (default 0.05).beta (
float, default:0.05) – Type II error probability (default 0.05).
- Returns:
decision–'reject_h0','accept_h0', or'inconclusive'.stopping_time– index at which the decision was reached. Fewerobservations needed indicates stronger signal.
log_ratio– final log-likelihood ratio.- Return type:
Example
>>> import numpy as np >>> from wraquant.math.optimal_stopping import sequential_probability_ratio >>> rng = np.random.default_rng(42) >>> # Generate data from H1 (mu=0.5) and test against H0 (mu=0) >>> obs = rng.normal(0.5, 1.0, size=100) >>> result = sequential_probability_ratio( ... obs, ... h0_dist=('normal', {'mu': 0.0, 'sigma': 1.0}), ... h1_dist=('normal', {'mu': 0.5, 'sigma': 1.0}), ... ) >>> result['decision'] 'reject_h0'
See also
cusum_stoppingCUSUM-based change detection (non-parametric alternative).
- cusum_stopping(observations, target_mean, threshold)[source]¶
CUSUM stopping rule for detecting a mean shift.
Monitors cumulative sum of deviations from target_mean and triggers a stop when the CUSUM statistic exceeds threshold.
- Parameters:
- Returns:
detected– whether a change was detected.stopping_time– index of detection (or length of data if none).cusum_pos– final positive CUSUM statistic.cusum_neg– final negative CUSUM statistic.cusum_values– array of max(cusum_pos, |cusum_neg|) at each step.- Return type:
Example
>>> import numpy as np >>> from wraquant.math.optimal_stopping import cusum_stopping >>> # Returns with a mean shift at index 50 >>> obs = np.concatenate([np.random.randn(50) * 0.01, ... np.random.randn(50) * 0.01 + 0.05]) >>> result = cusum_stopping(obs, target_mean=0.0, threshold=0.5) >>> result['detected'] True >>> result['stopping_time'] > 50 # detected after the shift True
See also
sequential_probability_ratioParametric sequential test.
wraquant.ts.changepointsChange-point detection for time series.
- secretary_problem_threshold(n_candidates)[source]¶
Compute the optimal stopping threshold for the secretary problem.
The classical 1/e rule: reject the first r - 1 candidates, then hire the first one who is better than all previously seen.
- Parameters:
n_candidates (
int) – Total number of candidates.- Returns:
threshold– number of candidates to unconditionally reject.optimal_fraction– threshold / n_candidates. Converges to1/e (approximately 0.368) as n_candidates grows.
success_probability– probability of selecting the best candidate.Converges to 1/e for large n_candidates.
- Return type:
Example
>>> from wraquant.math.optimal_stopping import secretary_problem_threshold >>> result = secretary_problem_threshold(100) >>> 30 <= result['threshold'] <= 40 # close to 1/e * 100 ~ 37 True >>> result['success_probability'] > 0.3 True
See also
sequential_probability_ratioSequential testing framework.
Hawkes Processes¶
Hawkes process for modelling event clustering.
Useful for order-flow analysis, jump modelling, and trade-arrival processes.
- hawkes_intensity(times, mu, alpha, beta)[source]¶
Compute the conditional intensity of a univariate Hawkes process.
\[\begin{split}\\lambda(t_i) = \\mu + \\alpha \\sum_{t_j < t_i} \\beta\\, e^{-\\beta (t_i - t_j)}\end{split}\]The intensity is evaluated at each event time in times.
- Parameters:
- Returns:
Intensity evaluated at each event time. Values above mu indicate self-excitation from recent events.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.hawkes import hawkes_intensity >>> times = np.array([0.1, 0.2, 0.25, 0.5, 1.0]) >>> intensity = hawkes_intensity(times, mu=1.0, alpha=0.5, beta=2.0) >>> intensity[0] # first event: just background 1.0 >>> intensity[2] > intensity[0] # cluster of events excites intensity True
See also
simulate_hawkesGenerate event times from a Hawkes process.
fit_hawkesCalibrate parameters from observed event times.
- simulate_hawkes(mu, alpha, beta, T, seed=None)[source]¶
Simulate event times from a univariate Hawkes process.
Uses the Ogata thinning algorithm.
- Parameters:
- Returns:
Sorted array of event times in
[0, T].- Return type:
- Raises:
ValueError – If the branching ratio
alpha / beta >= 1(non-stationary).
Example
>>> from wraquant.math.hawkes import simulate_hawkes >>> events = simulate_hawkes(mu=1.0, alpha=0.5, beta=2.0, T=100.0, ... seed=42) >>> len(events) > 0 True >>> events[0] >= 0 and events[-1] <= 100.0 True
See also
hawkes_intensityCompute intensity at observed event times.
fit_hawkesFit parameters from event data.
wraquant.microstructure.toxicityToxicity models for order flow.
- fit_hawkes(times, T=None)[source]¶
Fit a univariate Hawkes process via maximum likelihood.
- Parameters:
- Returns:
mu– fitted background intensity.alpha– fitted excitation magnitude.beta– fitted decay rate.log_likelihood– maximised log-likelihood value.branching_ratio–alpha / beta. Values close to 1indicate strong self-excitation; close to 0 means events are nearly independent.
- Return type:
Example
>>> from wraquant.math.hawkes import simulate_hawkes, fit_hawkes >>> events = simulate_hawkes(mu=1.0, alpha=0.5, beta=2.0, T=1000.0, ... seed=42) >>> result = fit_hawkes(events, T=1000.0) >>> 0 < result['branching_ratio'] < 1 True
See also
hawkes_intensityEvaluate intensity with fitted parameters.
hawkes_branching_ratioQuick branching ratio computation.
- hawkes_branching_ratio(alpha, beta)[source]¶
Compute the branching ratio of a Hawkes process.
The branching ratio
alpha / betadetermines stationarity: the process is stationary if and only if the ratio is strictly less than 1.- Parameters:
- Returns:
Branching ratio
alpha / beta. Must be < 1 for stationarity.- Return type:
Example
>>> from wraquant.math.hawkes import hawkes_branching_ratio >>> hawkes_branching_ratio(alpha=0.5, beta=2.0) 0.25
See also
fit_hawkesEstimate parameters (includes branching ratio in output).
Numerical Methods¶
General-purpose numerical methods useful for quantitative finance.
- finite_difference_gradient(fn, x, dx=1e-07)[source]¶
Numerical gradient via central finite differences.
- Parameters:
- Returns:
Gradient vector of the same shape as x.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.numerical import finite_difference_gradient >>> # Gradient of f(x) = x[0]^2 + x[1]^2 at (3, 4) should be (6, 8) >>> grad = finite_difference_gradient(lambda x: x[0]**2 + x[1]**2, ... np.array([3.0, 4.0])) >>> np.allclose(grad, [6.0, 8.0], atol=1e-5) True
See also
finite_difference_hessianCompute the second-derivative matrix.
wraquant.math.information.fisher_informationFisher information via Hessian.
- finite_difference_hessian(fn, x, dx=1e-05)[source]¶
Numerical Hessian matrix via central finite differences.
- Parameters:
- Returns:
Hessian matrix of shape
(len(x), len(x)).- Return type:
Example
>>> import numpy as np >>> from wraquant.math.numerical import finite_difference_hessian >>> # Hessian of f(x) = x[0]^2 + 3*x[1]^2 is diag(2, 6) >>> H = finite_difference_hessian(lambda x: x[0]**2 + 3*x[1]**2, ... np.array([1.0, 1.0])) >>> np.allclose(np.diag(H), [2.0, 6.0], atol=1e-3) True
See also
finite_difference_gradientFirst-order derivatives.
- newton_raphson(fn, x0, tol=1e-08, max_iter=100, dfn=None)[source]¶
Find a root of fn using the Newton-Raphson method.
- Parameters:
fn (
Callable[[float],float]) – Functionf(x) -> floatwhose root is sought.x0 (
float) – Initial guess.tol (
float, default:1e-08) – Convergence tolerance on|f(x)|(default 1e-8).max_iter (
int, default:100) – Maximum number of iterations (default 100).dfn (
Optional[Callable[[float],float]], default:None) – Derivativef'(x). IfNone, a central-difference approximation is used.
- Returns:
Approximate root.
- Return type:
- Raises:
RuntimeError – If the method does not converge within max_iter iterations.
Example
>>> from wraquant.math.numerical import newton_raphson >>> # Find implied volatility: solve BS_price(vol) - market_price = 0 >>> root = newton_raphson(lambda x: x**2 - 2, x0=1.5) >>> abs(root - 2**0.5) < 1e-8 True
See also
bisectionGuaranteed convergence but slower (no derivative needed).
- bisection(fn, a, b, tol=1e-08, max_iter=100)[source]¶
Find a root of fn on
[a, b]using the bisection method.- Parameters:
fn (
Callable[[float],float]) – Continuous functionf(x) -> float.a (
float) – Left endpoint of the bracket.b (
float) – Right endpoint of the bracket.tol (
float, default:1e-08) – Convergence tolerance on the bracket width (default 1e-8).max_iter (
int, default:100) – Maximum number of iterations (default 100).
- Returns:
Approximate root.
- Return type:
- Raises:
ValueError – If
f(a)andf(b)have the same sign.RuntimeError – If the method does not converge within max_iter iterations.
Example
>>> from wraquant.math.numerical import bisection >>> root = bisection(lambda x: x**3 - 1, a=0.0, b=2.0) >>> abs(root - 1.0) < 1e-8 True
See also
newton_raphsonFaster convergence when a derivative is available.
- trapezoidal_integration(fn, a, b, n=1000)[source]¶
Numerical integration using the trapezoidal rule.
- Parameters:
- Returns:
Approximate value of the definite integral.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.numerical import trapezoidal_integration >>> # Integrate sin(x) from 0 to pi (exact answer = 2) >>> result = trapezoidal_integration(np.sin, 0, np.pi, n=1000) >>> abs(result - 2.0) < 0.001 True
See also
monte_carlo_integrationStochastic integration for high dimensions.
- monte_carlo_integration(fn, bounds, n_samples=100000, seed=None)[source]¶
Monte Carlo integration over a hyper-rectangular domain.
- Parameters:
fn (
Callable[...,float]) – Function to integrate. Should accept a 1-D array of lengthlen(bounds)(or individual floats for 1-D).bounds (
list[tuple[float,float]]) – Integration bounds for each dimension, e.g.[(0, 1), (0, 2)].n_samples (
int, default:100000) – Number of random samples (default 100 000).seed (
int|None, default:None) – Random seed for reproducibility.
- Returns:
estimate– estimated value of the integral.std_error– standard error of the estimate. Decreases as1 / sqrt(n_samples).- Return type:
Example
>>> import numpy as np >>> from wraquant.math.numerical import monte_carlo_integration >>> # Integrate x^2 from 0 to 1 (exact answer = 1/3) >>> result = monte_carlo_integration(lambda x: x[0]**2, ... bounds=[(0, 1)], seed=42) >>> abs(result['estimate'] - 1/3) < 0.01 True
See also
trapezoidal_integrationDeterministic integration for 1-D problems.
Spectral Methods¶
Spectral analysis / FFT tools for financial data.
- fft_decompose(data, n_components=None)[source]¶
Perform FFT decomposition of a 1-D signal.
- Parameters:
- Returns:
frequencies– normalised frequencies (cycles / sample).amplitudes– amplitude of each frequency component.phases– phase angle (radians) of each component.power– power spectrum (amplitude squared).- Return type:
Example
>>> import numpy as np >>> from wraquant.math.spectral import fft_decompose >>> t = np.arange(256) >>> signal = np.sin(2 * np.pi * t / 20) # period=20 bars >>> result = fft_decompose(signal, n_components=3) >>> result['frequencies'].shape[0] 3
See also
dominant_frequenciesExtract the strongest cyclical components.
spectral_densityPower spectral density estimation.
bandpass_filterIsolate a frequency band from the signal.
- dominant_frequencies(data, n_top=5)[source]¶
Identify the dominant cyclical frequencies in data.
- Parameters:
data (
TypeAliasType) – Input time series.n_top (
int, default:5) – Number of top frequencies to return (default 5).
- Returns:
frequency– normalised frequencies of the top components.period– corresponding period in bars (1 / frequency).Periods near common trading horizons (5, 21, 63 bars) may indicate weekly, monthly, or quarterly cycles.
amplitude– amplitude of each component.- Return type:
Example
>>> import numpy as np >>> from wraquant.math.spectral import dominant_frequencies >>> t = np.arange(500) >>> signal = np.sin(2 * np.pi * t / 21) + 0.5 * np.sin(2 * np.pi * t / 63) >>> result = dominant_frequencies(signal, n_top=2) >>> # The two dominant periods should be near 21 and 63 >>> sorted(result['period'])[:2] [21.0, 63.0]
See also
fft_decomposeFull FFT decomposition.
spectral_entropyMeasure spectral complexity.
- spectral_density(data, method='periodogram')[source]¶
Estimate the power spectral density (PSD) of data.
- Parameters:
data (
TypeAliasType) – Input time series.method (
str, default:'periodogram') – Estimation method (default'periodogram').
- Returns:
frequencies– frequency axis.psd– power spectral density values. Higher PSD atlow frequencies indicates trending behaviour; flat PSD indicates white noise.
- Return type:
- Raises:
ValueError – If method is not recognised.
Example
>>> import numpy as np >>> from wraquant.math.spectral import spectral_density >>> data = np.random.randn(500) >>> result = spectral_density(data, method='welch') >>> len(result['frequencies']) == len(result['psd']) True
See also
fft_decomposeRaw FFT decomposition.
spectral_entropySingle-number summary of spectral shape.
- bandpass_filter(data, low_freq, high_freq, sampling_rate=1.0)[source]¶
Apply a bandpass filter to isolate a frequency band.
Uses a 5th-order Butterworth filter applied forward-backward (zero phase shift) via
scipy.signal.sosfiltfilt().- Parameters:
- Returns:
Filtered signal.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.spectral import bandpass_filter >>> t = np.arange(500) >>> # Signal with a 21-bar cycle plus high-frequency noise >>> signal = np.sin(2 * np.pi * t / 21) + 0.5 * np.random.randn(500) >>> filtered = bandpass_filter(signal, low_freq=0.04, high_freq=0.06) >>> filtered.shape (500,)
See also
wraquant.math.signals.savitzky_golayPolynomial smoothing filter.
dominant_frequenciesIdentify which frequencies to isolate.
- spectral_entropy(data)[source]¶
Compute the spectral entropy of data.
Spectral entropy measures the “flatness” of the power spectrum. A perfectly periodic signal has low spectral entropy; white noise has high spectral entropy (close to
log2(N/2)).- Parameters:
data (
TypeAliasType) – Input time series.- Returns:
Normalised spectral entropy in [0, 1]. Values near 0 indicate a strongly periodic signal; values near 1 indicate noise-like behaviour with no dominant frequency.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.spectral import spectral_entropy >>> # Pure sine wave has low spectral entropy >>> t = np.arange(256) >>> se_sine = spectral_entropy(np.sin(2 * np.pi * t / 20)) >>> se_noise = spectral_entropy(np.random.randn(256)) >>> se_sine < se_noise True
See also
wraquant.math.information.entropyShannon entropy of a data distribution.
spectral_densityFull PSD estimation.
Signal Processing¶
Signal processing filters for financial time series.
- savitzky_golay(data, window=11, polyorder=3)[source]¶
Apply a Savitzky-Golay smoothing filter.
- Parameters:
- Returns:
Smoothed signal.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.signals import savitzky_golay >>> noisy = np.sin(np.linspace(0, 4*np.pi, 200)) + np.random.randn(200)*0.3 >>> smooth = savitzky_golay(noisy, window=15, polyorder=3) >>> smooth.shape (200,) >>> np.std(smooth) < np.std(noisy) # smoothing reduces variance True
See also
kalman_smoothKalman-filter-based smoothing (adaptive).
exponential_smoothSimple exponential smoothing.
- kalman_smooth(data, process_var=1e-05, measurement_var=0.1)[source]¶
Simple 1-D Kalman smoother (forward pass).
Implements a constant-level model with Gaussian noise.
- Parameters:
- Returns:
Kalman-smoothed estimate (same length as data).
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.signals import kalman_smooth >>> true_signal = np.cumsum(np.random.randn(200) * 0.01) >>> noisy = true_signal + np.random.randn(200) * 0.05 >>> smooth = kalman_smooth(noisy, process_var=1e-5, measurement_var=0.01) >>> smooth.shape (200,)
Notes
Larger process_var makes the filter more responsive to changes; larger measurement_var makes it smoother. For trending markets, increase process_var; for noisy mean-reverting signals, increase measurement_var.
See also
savitzky_golayPolynomial smoothing (no tuning parameters beyond window).
wraquant.regimesKalman filter with state-space models.
- wavelet_denoise(data, wavelet='db4', level=None, threshold='soft')[source]¶
Wavelet denoising of a time series.
Requires the optional
pywaveletspackage (install viapip install wraquant[timeseries]).- Parameters:
- Returns:
Denoised signal.
- Return type:
- Raises:
ImportError – If
pywtis not installed.
Example
>>> import numpy as np >>> from wraquant.math.signals import wavelet_denoise >>> noisy = np.sin(np.linspace(0, 4*np.pi, 256)) + np.random.randn(256)*0.5 >>> clean = wavelet_denoise(noisy, wavelet='db4') >>> np.std(clean) < np.std(noisy) True
See also
savitzky_golayPolynomial smoothing (no optional dependencies).
median_filterRobust to outlier spikes.
- median_filter(data, kernel_size=5)[source]¶
Median filter for spike removal.
- Parameters:
data (
TypeAliasType) – Input time series.kernel_size (
int, default:5) – Size of the median-filter kernel (default 5).
- Returns:
Filtered signal.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.signals import median_filter >>> data = np.array([1.0, 1.1, 50.0, 1.05, 0.95]) # spike at index 2 >>> filtered = median_filter(data, kernel_size=3) >>> filtered[2] < 10 # spike removed True
See also
savitzky_golayPolynomial smoothing for continuous signals.
- exponential_smooth(data, alpha=0.3)[source]¶
Simple exponential smoothing.
\[\begin{split}s_t = \\alpha \\, x_t + (1 - \\alpha) \\, s_{t-1}\end{split}\]- Parameters:
data (
TypeAliasType) – Input time series.alpha (
float, default:0.3) – Smoothing factor in(0, 1](default 0.3).
- Returns:
Smoothed signal.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.signals import exponential_smooth >>> data = np.array([1.0, 2.0, 3.0, 2.5, 2.0]) >>> smooth = exponential_smooth(data, alpha=0.5) >>> smooth[0] 1.0 >>> smooth[-1] < data[-2] # smoothed toward the end True
See also
kalman_smoothAdaptive smoothing with noise model.
savitzky_golayLocal polynomial smoothing.
- hodrick_prescott(data, lamb=1600.0)[source]¶
Hodrick-Prescott filter decomposition.
Decomposes data into a trend and a cyclical component.
- Parameters:
data (
TypeAliasType) – Input time series.lamb (
float, default:1600.0) – Smoothing parameter (default 1600, standard for quarterly data). Use 6.25 for annual data, 129600 for monthly data.
- Returns:
(trend, cycle)wheretrend + cycle == data.- Return type:
Example
>>> import numpy as np >>> from wraquant.math.signals import hodrick_prescott >>> data = np.cumsum(np.random.randn(100) * 0.01) # random walk >>> trend, cycle = hodrick_prescott(data, lamb=1600) >>> trend.shape == cycle.shape == (100,) True >>> np.allclose(trend + cycle, data) True
Notes
Common lambda values: 6.25 (annual), 1600 (quarterly), 129600 (monthly). The HP filter is controversial for financial data due to end-point bias and spurious cyclicality.
See also
wraquant.ts.decompositionSeasonal decomposition of time series.
wraquant.math.spectral.bandpass_filterFrequency-domain detrending.
Information Theory¶
Information-theoretic measures for financial analysis.
- fisher_information(log_likelihood_fn, params, dx=1e-05)[source]¶
Numerical Fisher information matrix via second derivatives.
Computes the negative of the Hessian of log_likelihood_fn evaluated at params using central finite differences.
- Parameters:
- Returns:
Fisher information matrix of shape
(len(params), len(params)). Larger diagonal entries indicate parameters that are more precisely estimable from the data.- Return type:
Example
>>> import numpy as np >>> from wraquant.math.information import fisher_information >>> # Log-likelihood of Normal(mu, sigma=1): -0.5 * sum((x - mu)^2) >>> data = np.array([1.0, 2.0, 3.0]) >>> ll_fn = lambda p: -0.5 * np.sum((data - p[0])**2) >>> fim = fisher_information(ll_fn, np.array([2.0])) >>> fim.shape (1, 1) >>> fim[0, 0] > 0 # positive definite True
See also
wraquant.math.numerical.finite_difference_hessianGeneral Hessian computation.
- mutual_information(x, y, bins=20)[source]¶
Mutual information between two series (discretised).
\[I(X; Y) = H(X) + H(Y) - H(X, Y)\]- Parameters:
x (
TypeAliasType) – First data series.y (
TypeAliasType) – Second data series.bins (
int, default:20) – Number of histogram bins per dimension (default 20).
- Returns:
Mutual information in nats (>= 0). Zero indicates independence; higher values indicate stronger dependence (linear or nonlinear).
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.information import mutual_information >>> rng = np.random.default_rng(42) >>> x = rng.standard_normal(1000) >>> y = x + rng.standard_normal(1000) * 0.1 # highly dependent >>> z = rng.standard_normal(1000) # independent >>> mutual_information(x, y) > mutual_information(x, z) True
See also
transfer_entropyDirected (causal) information flow.
conditional_entropyH(X | Y) = H(X) - I(X; Y).
- transfer_entropy(source, target, lag=1, bins=10)[source]¶
Transfer entropy from source to target.
Measures the directional information flow from source to target beyond what target’s own past explains.
\[\begin{split}TE_{X \\to Y} = H(Y_t | Y_{t-k}) - H(Y_t | Y_{t-k}, X_{t-k})\end{split}\]- Parameters:
- Returns:
Transfer entropy in nats (>= 0). Higher values indicate stronger directional information flow from source to target.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.information import transfer_entropy >>> rng = np.random.default_rng(42) >>> x = rng.standard_normal(500) >>> y = np.concatenate([[0], x[:-1]]) + rng.standard_normal(500) * 0.1 >>> te_x_to_y = transfer_entropy(x, y, lag=1) >>> te_y_to_x = transfer_entropy(y, x, lag=1) >>> te_x_to_y > te_y_to_x # x drives y, not vice versa True
Notes
Reference: Schreiber, T. (2000). “Measuring Information Transfer.” Physical Review Letters, 85(2), 461-464.
See also
mutual_informationSymmetric (undirected) dependence measure.
wraquant.math.network.granger_networkLinear Granger causality.
- entropy(data, bins=20, method='histogram')[source]¶
Shannon entropy of a data series.
- Parameters:
- Returns:
Shannon entropy in nats. Higher values indicate more dispersed (uncertain) distributions; lower values indicate concentrated distributions.
- Return type:
- Raises:
ValueError – If method is not recognised.
Example
>>> import numpy as np >>> from wraquant.math.information import entropy >>> uniform = np.random.uniform(size=1000) >>> peaked = np.random.normal(0, 0.01, size=1000) >>> entropy(uniform) > entropy(peaked) True
See also
conditional_entropyEntropy of X given Y.
mutual_informationShared information between two variables.
wraquant.math.spectral.spectral_entropyEntropy of the power spectrum.
- kl_divergence(p, q, bins=20)[source]¶
KL divergence D_KL(P || Q) estimated from samples.
- Parameters:
p (
TypeAliasType) – Samples from distribution P.q (
TypeAliasType) – Samples from distribution Q.bins (
int, default:20) – Number of histogram bins (default 20).
- Returns:
KL divergence in nats (>= 0). Zero when P and Q are identical; larger values indicate greater distributional difference. Note: KL divergence is asymmetric – D_KL(P||Q) != D_KL(Q||P).
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.information import kl_divergence >>> rng = np.random.default_rng(42) >>> p = rng.normal(0, 1, size=5000) >>> q = rng.normal(0.5, 1, size=5000) >>> kl_divergence(p, q) > 0 True
See also
entropyShannon entropy of a single distribution.
mutual_informationSymmetric dependence measure.
- conditional_entropy(x, y, bins=20)[source]¶
Conditional entropy H(X | Y).
\[H(X | Y) = H(X, Y) - H(Y)\]- Parameters:
x (
TypeAliasType) – First data series.y (
TypeAliasType) – Second data series (the conditioning variable).bins (
int, default:20) – Number of histogram bins per dimension (default 20).
- Returns:
Conditional entropy in nats. Lower values mean Y is more informative about X. When H(X|Y) = 0, X is fully determined by Y.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.information import conditional_entropy >>> rng = np.random.default_rng(42) >>> x = rng.standard_normal(1000) >>> y = x + rng.standard_normal(1000) * 0.1 # y almost determines x >>> h_x_given_y = conditional_entropy(x, y) >>> h_x_given_y < 1.0 # low because y is informative about x True
See also
entropyUnconditional Shannon entropy H(X).
mutual_informationI(X;Y) = H(X) - H(X|Y).
Ergodicity¶
Ergodicity economics (Ole Peters framework).
Tools for distinguishing ensemble averages from time averages, computing Kelly-optimal fractions, and measuring ergodicity.
- ensemble_average(returns)[source]¶
Arithmetic mean of returns (the ensemble average).
- Parameters:
returns (
TypeAliasType) – Simple (arithmetic) returns, e.g.[0.05, -0.02, 0.03].- Returns:
Arithmetic mean of the returns.
- Return type:
Example
>>> from wraquant.math.ergodicity import ensemble_average >>> ensemble_average([0.05, -0.02, 0.03]) 0.02
See also
time_averageGeometric mean (the time-average growth rate).
ergodicity_gapDifference between ensemble and time averages.
- time_average(returns)[source]¶
Time-average growth rate (geometric mean return).
Computes the annualised-per-period geometric growth rate:
\[\begin{split}g = \\left(\\prod_{i=1}^{N}(1 + r_i)\\right)^{1/N} - 1\end{split}\]- Parameters:
returns (
TypeAliasType) – Simple (arithmetic) returns.- Returns:
Geometric mean return per period. Always less than or equal to the arithmetic mean (ensemble average) due to Jensen’s inequality. This is the growth rate actually experienced by a single investor over time.
- Return type:
Example
>>> from wraquant.math.ergodicity import time_average >>> ta = time_average([0.10, -0.10]) >>> ta < 0.0 # net loss despite zero arithmetic mean True
See also
ensemble_averageArithmetic mean (the ensemble expectation).
ergodicity_gapQuantifies the difference between the two averages.
- ergodicity_gap(returns)[source]¶
Difference between the ensemble average and the time average.
A positive gap means the ensemble (arithmetic) average overstates the realised long-run growth.
- Parameters:
returns (
TypeAliasType) – Simple returns.- Returns:
ensemble_average(returns) - time_average(returns). A larger gap indicates more volatility drag and greater non-ergodicity.- Return type:
Example
>>> from wraquant.math.ergodicity import ergodicity_gap >>> # High volatility creates a large gap >>> gap = ergodicity_gap([0.50, -0.50, 0.50, -0.50]) >>> gap > 0 True
See also
ergodicity_ratioRatio version of the same concept.
kelly_fractionOptimal leverage accounting for the gap.
- kelly_fraction(returns)[source]¶
Optimal Kelly criterion fraction for a simple binary-style bet.
For a series of returns this computes the leverage that maximises the expected log-growth rate via a simple numerical optimisation over a grid.
- Parameters:
returns (
TypeAliasType) – Simple returns for each period.- Returns:
Optimal Kelly fraction (leverage). A value of 1.0 means full investment; values > 1.0 indicate levered positions. Values near 0 suggest the edge is too small relative to risk.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.ergodicity import kelly_fraction >>> rng = np.random.default_rng(42) >>> # Strategy with positive expected return >>> returns = rng.normal(0.001, 0.02, size=1000) >>> f = kelly_fraction(returns) >>> f > 0 # positive edge -> positive Kelly fraction True
Notes
Reference: Kelly, J. L. (1956). “A New Interpretation of Information Rate.” Bell System Technical Journal, 35(4), 917-926.
See also
growth_optimal_leverageKelly with a risk-free rate.
ergodicity_gapUnderstand why Kelly < 1 / edge.
- growth_optimal_leverage(returns, risk_free=0.0)[source]¶
Leverage that maximises the time-average growth rate.
Maximises
E[log(1 + risk_free + f * (r - risk_free))]over f.- Parameters:
returns (
TypeAliasType) – Simple returns per period.risk_free (
float, default:0.0) – Risk-free rate per period (default 0.0).
- Returns:
Growth-optimal leverage.
- Return type:
Example
>>> import numpy as np >>> from wraquant.math.ergodicity import growth_optimal_leverage >>> rng = np.random.default_rng(42) >>> returns = rng.normal(0.001, 0.02, size=1000) >>> lev = growth_optimal_leverage(returns, risk_free=0.0001) >>> lev > 0 True
See also
kelly_fractionSimpler version without risk-free rate.
- ergodicity_ratio(returns)[source]¶
Ratio of the time-average to the ensemble-average growth rate.
A ratio of 1.0 indicates an ergodic process; ratios below 1.0 indicate that time averaging yields lower growth than the ensemble expectation.
- Parameters:
returns (
TypeAliasType) – Simple returns.- Returns:
time_average(returns) / ensemble_average(returns). Returns1.0if the ensemble average is zero.- Return type:
Example
>>> import numpy as np >>> from wraquant.math.ergodicity import ergodicity_ratio >>> # Low-volatility returns are nearly ergodic >>> low_vol = np.random.randn(1000) * 0.001 + 0.0005 >>> ratio_lv = ergodicity_ratio(low_vol) >>> 0.9 < ratio_lv <= 1.0 True
See also
ergodicity_gapAdditive version of the same concept.