Quant Finance: Option Pricing beyond Black-Scholes via Quantum Mechanics (C++ Code & Research)

Implementing wavefunction-derived probability densities and effective volatility modeling under quantum force potentials in C++17.

Download The source code using the button at the end of this article.

The paper did not include any official implementation, so this code was written by me from scratch. While I have tested several parts of it, it may still produce imperfect results. Consider it a demonstration baseline that can be extended and improved further.

Photo by Markus Spiske on Unsplash

The paper proposes a quantum mechanics approach to option pricing by mapping the Black-Scholes standard normal density to the ground state of a harmonic oscillator. Additional potentials model market forces, leading to modified probability densities. Option prices are computed by integrating the new density and adjusting the effective volatility. Several force types are analyzed, and the resulting option prices are compared numerically.

Implementation Assumptions

• ℏω = 1, m = 1, ω = 1 for all wavefunction-derived densities (consistent with Appendix Sec. 5.1).

• Numerical integration domain: [-12, 12] with 20001 points (Simpson's rule) for all unbounded densities.

• For quantum well, integration domain is [-a, a] with proportionally adjusted resolution.

• The x² force model includes a compile-time flag USECORRECTEDX2 to optionally replace the paper's erroneous linear-β term with a correct second-order perturbation (derived externally from standard QM).

• All option pricing uses the Black-Scholes formula structure with σeff and modified Neff(d_eff±).

• No external dependencies beyond C++17 standard library.

• Output is CSV files for external plotting; no built-in plotting library.

Introduction: The Black-Scholes Model and Its Limitations

Option pricing is about predicting the future value of an asset. The Black-Scholes model is the classic formula that assumes asset prices follow a random walk with a bell-curve distribution. This section introduces that formula and explains why it works, using simple terms like 'expected payoff' and 'discounting'. We will also look at the C++ implementation that serves as the baseline for the quantum extensions later in the tutorial.

The Black-Scholes Formula

The Black-Scholes model gives the theoretical price of a European call or put option under the assumption that the underlying asset price follows a geometric Brownian motion with constant volatility. The call price c and put price p are given by the following equations.

The call option price is

\(c = S_0 N(d_+) - K e^{-rT} N(d_-)\)

where S_0 is the current asset price, K is the strike price, r is the risk-free interest rate, T is the time to maturity, and N(·) is the cumulative distribution function of the standard normal distribution. The put option price is

\(p = K e^{-rT} N(-d_-) - S_0 N(-d_+)\)

These formulas express the option value as the discounted expected payoff under the risk-neutral measure. The terms N(d_+) and N(d_-) can be interpreted as the risk-adjusted probabilities that the option finishes in-the-money.

The Quantiles d_+ and d_-

The arguments d_+ and d_- are computed from the model parameters as

\(d_{\pm} = \frac{\ln(S_0/K) + (r \pm \sigma^2/2)T}{\sigma\sqrt{T}}\)

where σ is the annualized volatility of the asset. These quantiles arise from the log-normal distribution of the asset price at maturity. The term ln(S_0/K) is the log-moneyness, rT is the drift, and σ²T/2 is the Itô correction. The denominator σ√T standardises the distance to the strike in units of standard deviation.

The Standard Normal CDF

The cumulative distribution function N(d) is defined as the integral of the standard normal probability density P_BS(x):

\(N(d_{\pm}) = \int_{-\infty}^{d_{\pm}} P_{\text{BS}}(x) \, dx, \quad P_{\text{BS}}(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2/2}\)

This density is the familiar bell curve. In the Black-Scholes world, the log-return is normally distributed, so N(d_+) and N(d_-) are the probabilities that the option expires in-the-money under the risk-neutral measure (adjusted for the drift).

Implementation in C++

The baseline Black-Scholes pricing is implemented in the header file src/black_scholes_base.h. The code provides functions for the probability density, the cumulative distribution, the computation of d_±, and the option pricing itself.

First, the standard normal PDF and CDF are implemented using the standard library:

inline double normalPDF(double x) {
    constexpr double inv_sqrt_2pi = 0.3989422804014327; // 1/sqrt(2*pi)
    return inv_sqrt_2pi * std::exp(-0.5 * x * x);
}

inline double normalCDF(double x) {
    return 0.5 * std::erfc(-x * 0.7071067811865476); // 0.7071... = 1/sqrt(2)
}

The CDF uses the complementary error function std::erfc for high precision. The constant 0.7071... is 1/√2. This implementation is monotonic and returns values in [0,1].

Next, the quantiles d_+ and d_- are computed by the function compute_dpm:

inline std::pair<double, double> compute_dpm(double S0, double K, double r,
                                              double sigma, double T) {
    if (sigma <= 0.0) {
        throw std::invalid_argument("compute_dpm: sigma must be positive");
    }
    if (T <= 0.0) {
        throw std::invalid_argument("compute_dpm: T must be positive");
    }
    if (S0 <= 0.0 || K <= 0.0) {
        throw std::invalid_argument("compute_dpm: S0 and K must be positive");
    }

    double sigma_sqrt_T = sigma * std::sqrt(T);
    double log_moneyness = std::log(S0 / K);
    double drift = r * T;
    double half_var = 0.5 * sigma * sigma * T;

    double d_plus  = (log_moneyness + drift + half_var) / sigma_sqrt_T;
    double d_minus = (log_moneyness + drift - half_var) / sigma_sqrt_T;

    return {d_plus, d_minus};
}

This function validates the inputs (positive volatility, time, and prices) and then directly translates the mathematical formula into code. The half_var term is σ²T/2.

Finally, the option prices are computed by black_scholes_base, which returns a BSResult struct containing the call and put prices, the quantiles, and the CDF values:

struct BSResult {
    double call;      // European call option price
    double put;       // European put option price
    double d_plus;    // d_+ quantile
    double d_minus;   // d_- quantile
    double N_d_plus;  // N(d_+)
    double N_d_minus; // N(d_-)
};

inline BSResult black_scholes_base(double S0, double K, double r, double T,
                                    double sigma) {
    // ... input validation ...
    auto [d_plus, d_minus] = compute_dpm(S0, K, r, sigma, T);
    double N_d_plus  = normalCDF(d_plus);
    double N_d_minus = normalCDF(d_minus);
    double discount = std::exp(-r * T);
    double discounted_strike = K * discount;

    double call = S0 * N_d_plus - discounted_strike * N_d_minus;
    double N_neg_d_plus  = 1.0 - N_d_plus;
    double N_neg_d_minus = 1.0 - N_d_minus;
    double put = discounted_strike * N_neg_d_minus - S0 * N_neg_d_plus;

    return {call, put, d_plus, d_minus, N_d_plus, N_d_minus};
}

The call price is computed exactly as in the equation, and the put price uses the identity N(-d) = 1 - N(d) to avoid redundant CDF evaluations.

Worked Example

Using the paper's synthetic parameters—S0 = 20, K = 20, r = 0.1 (10%), T = 1 year, σ = 0.25 (25%)—the baseline call price is approximately 3.022 and the put price is approximately 1.119. The code computes these values with a simple call:

auto [call, put, dp, dm, Np, Nm] = black_scholes_base(20.0, 20.0, 0.1, 1.0, 0.25);
// call ≈ 3.022, put ≈ 1.119

Verification: Put-Call Parity

A fundamental consistency check for any option pricing model is put-call parity: c + K e^{-rT} = p + S_0. This relationship holds regardless of the distribution of the underlying, as long as the pricing is arbitrage-free. The implementation includes a static verification function verify_black_scholes_base() that checks this parity (and other properties) to ensure the code is correct. For the sample parameters, the left-hand side and right-hand side agree to within 1e-12.

Limitations of the Black-Scholes Model

The Black-Scholes model is elegant and widely used, but it rests on assumptions that are often violated in real markets: constant volatility, continuous trading, no transaction costs, and normally distributed log-returns. In practice, asset returns exhibit fat tails, skewness, and volatility clustering. The quantum mechanics approach introduced in the paper does not relax the continuous-trading or frictionless-market assumptions; instead, it provides a new way to generate non-Gaussian probability densities by adding "market potentials" to the harmonic oscillator that underlies the standard normal distribution. The following sections will build on this baseline to show how those modifications are implemented in C++.

Summary

We have reviewed the standard Black-Scholes formulas for European call and put options, the definition of d_±, and the role of the cumulative normal distribution. The C++ implementation in black_scholes_base.h provides a clean, verified baseline that we will later extend with quantum-mechanical probability densities. The next section will establish the analogy between the Black-Scholes density and the ground state of a quantum harmonic oscillator.

Quantum Mechanics Analogy: Mapping the Standard Normal Density to a Harmonic Oscillator Ground State

Imagine a marble rolling inside a smooth bowl. It tends to stay near the bottom, and if you recorded its position many times, you would find a bell-shaped pattern of positions. The Black-Scholes model assumes that the log-returns of an asset follow a similar bell-shaped (normal) distribution. We can think of the market as having a hidden "potential" that shapes the probability of price moves, just as the bowl determines where the marble is likely to be. This section builds the bridge between the financial model and the quantum mechanical harmonic oscillator.

The Harmonic Oscillator Potential

In quantum mechanics, a particle of mass m moving in a potential V(x) is described by the stationary Schrödinger equation. The baseline potential chosen in the paper is the harmonic oscillator:

\(V(x) = \frac{1}{2} m \omega^2 x^2\)

Here V(x) is the potential energy, m is the mass, omega is the angular frequency, and x is the position. This potential is shaped like a parabola. The corresponding Schrödinger equation reads

\(-\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2} + \frac{1}{2} m \omega^2 x^2 \psi = E \psi\)

where hbar is the reduced Planck constant, psi(x) is the wavefunction, and E is the energy eigenvalue. The ground state (lowest energy) wavefunction of this system is known exactly:

\(\psi_g(x) = \frac{\alpha}{\sqrt{\pi}} e^{-\alpha^2 x^2/2}, \quad \alpha = \sqrt{\frac{m\omega}{\hbar}}\)

The parameter alpha is the inverse length scale. The square modulus of the wavefunction gives the probability density of finding the particle at position x:

\(P_{g,HO}(x) = |\psi_g(x)|^2 = \frac{\alpha^2}{\pi} e^{-\alpha^2 x^2}\)

This is a Gaussian (bell curve). Its width is controlled by alpha.

Correspondence with the Black-Scholes Density

The Black-Scholes model uses the standard normal probability density, which can be written as

\(P_{g,HO}(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2/2} = P_{\text{BS}}(x) \quad \text{with } \alpha = 1/\sqrt{2}\)

Notice that both densities are Gaussians. To make the quantum density match the Black-Scholes density we would like the exponential part to be -x^2/2, which requires alpha^2 = 1/2, i.e., alpha = 1/sqrt(2). However, with this choice the prefactor becomes (1/2)/pi = 1/(2*pi), while the standard normal prefactor is 1/sqrt(2*pi). The normalization constants do not agree. The paper's derivation therefore contains an inconsistency in the normalization of the wavefunction.

For the purpose of option pricing, we side-step this issue and directly adopt the standard normal probability density as the baseline. This is the density that the Black-Scholes model uses and that will be modified later by market forces. In other words, we treat P_BS(x) as our reference ground state probability density, and we do not derive it from a wavefunction that would require a different normalization constant.

In the code, this baseline density is implemented by the function normalPDF.

Code: The Baseline Probability Density

The header src/black_scholes_base.h provides the baseline normal PDF:

inline double normalPDF(double x) {
    constexpr double inv_sqrt_2pi = 0.3989422804014327; // 1/sqrt(2*pi)
    return inv_sqrt_2pi * std::exp(-0.5 * x * x);
}

This function returns P_BS(x) exactly. It uses a precomputed constant 1/sqrt(2*pi) for efficiency.

Worked Example

Calling normalPDF(0.0) returns approximately 0.39894228. This is the peak of the standard normal bell curve, and it corresponds to the probability density at the mean. If we were to construct the wavefunction psi_g with alpha = 1/sqrt(2), we would obtain a squared modulus with a different peak value; the adopted normalPDF avoids that discrepancy and matches the Black-Scholes model directly.

Why This Analogy Matters

By mapping the standard normal density to the ground state of a harmonic oscillator, we gain a physical intuition: the parabolic potential keeps the distribution centered and bell-shaped. Adding extra potentials (market forces) will distort this baseline, leading to new probability densities that can capture skewness, fat tails, or hard bounds. The next sections show how to implement those modifications.

Advanced Note

The harmonic oscillator ground state is a Gaussian because the potential is quadratic. The paper's Eq. (8) provides the wavefunction, but we only ever need the probability density P(x) = |psi(x)|^2. Because of the normalization inconsistency, the code never uses a wavefunction directly; it computes probability densities using the target functional forms, normalizing them numerically when necessary.

Market Forces as Potentials: How Hidden Forces Modify the Probability Density

Just as a magnet can pull a marble away from the centre of a bowl, hidden “market forces” can pull the price distribution away from the standard bell curve. In the quantum analogy, these forces are represented by adding extra terms to the bowl’s shape—the potential energy function V(x). By choosing different potentials, we can generate a rich family of probability densities that go beyond the log‑normal assumption of Black–Scholes.

From Potential to Force

In classical physics, a conservative force is the negative gradient of a potential. The paper adopts the same relationship for the market force F(x):

\(F(x) = -\frac{dV}{dx}\)

Here F(x) is the market force acting at the dimensionless coordinate x (which corresponds to the log‑return variable in the Black–Scholes framework), and V(x) is the potential energy. A positive force pushes the distribution to the right; a negative force pushes it to the left. The shape of V(x) determines how the force changes with x.

The Baseline Potential

The starting point is the harmonic oscillator potential that reproduces the standard normal probability density of Black–Scholes. Its canonical form is V(x) = ½ m ω² x² (Eq. eqhopotential). Any additional potential V_extra(x) is added to this baseline, so the total potential becomes V_total(x) = ½ m ω² x² + V_extra(x). The new ground state of the modified Schrödinger equation then yields a different probability density.

Five Types of Market Forces

The paper studies five distinct modifications to the potential, each modelling a different kind of market influence. The table below summarises them; the following paragraphs explain each one qualitatively.

Constant Force

The simplest modification is a linear potential. The canonical form from the paper is:

\(V_k(x) = k x\)

This corresponds to a constant force F = -k. It tilts the harmonic oscillator bowl, shifting the equilibrium point to a new position. The probability density remains a Gaussian with the same variance, but its peak moves to -x_k (where x_k is proportional to k). In financial terms, a constant force represents a persistent buying or selling pressure that biases the expected return.

Linear Force

A quadratic addition to the potential models a force proportional to the displacement. The canonical form is:

\(V_\lambda(x) = \lambda x^2, \quad \lambda > 0\)

This force F = -2λ x acts like a spring: it pulls the distribution back toward the origin more strongly than the baseline harmonic oscillator. The result is a Gaussian with a narrower variance. A larger λ makes the bowl steeper, reducing the probability of large moves and therefore lowering option prices.

x² Force (Cubic Potential)

A cubic potential introduces an asymmetry. The canonical form is:

\(V_\beta(x) = \beta x^3\)

The force F = -3β x² is always negative (for β > 0), pushing the distribution to the left, but the strength grows quadratically with distance. This breaks the left‑right symmetry of the Gaussian, creating a skewed probability density. (The paper’s perturbative treatment of this case contains a known error, which we will discuss in detail later.)

x³ Force (Quartic Potential)

A quartic potential is symmetric but can change the shape dramatically. The canonical form is:

\(V_\gamma(x) = \gamma x^4\)

For small γ, the force F = -4γ x³ is weak near the origin and grows rapidly for large |x|. As γ increases, the total potential ½ m ω² x² + γ x⁴ can develop a double‑well shape—two minima separated by a barrier. The ground‑state probability density then becomes bimodal, with peaks on either side of zero. This can model a market that expects a large move in either direction (e.g., ahead of a binary event).

Quantum Well (Hard Price Bounds)

The infinite square well potential imposes absolute limits on the variable x. The canonical form is:

\(V(x) = \begin{cases} 0 & |x| < a \\ \infty & |x| \ge a \end{cases}\)

Inside the well the particle is free; at the walls the potential is infinite, so the probability density must vanish at x = ±a and is exactly zero outside. This models a market with hard price bounds—for example, a currency peg or a circuit breaker that prevents trading beyond a certain range. The width a controls how much the distribution is squeezed.

How the Potentials Modify the Probability Density

Each extra potential is added to the harmonic oscillator baseline. The new stationary Schrödinger equation is solved (exactly or via perturbation theory) to obtain the ground‑state wavefunction ψ(x). The probability density is then P(x) = |ψ(x)|². Because the potential shapes the wavefunction, different potentials produce different densities:

• A constant force shifts the Gaussian without changing its width.

• A linear force narrows or broadens the Gaussian.

• An x² force skews the distribution.

• An x³ force can split the distribution into two peaks.

• A quantum well truncates the distribution at hard boundaries.

In every case, the new density P(x) replaces the standard normal PDF in the option pricing formula. The effective volatility σ_eff is then computed from the quantum standard deviation of x under P(x), and the cumulative distribution is integrated numerically. The next sections will walk through each model in detail, showing the exact mathematical forms and the corresponding C++ code.

Worked example (qualitative): For a constant force F = -k, the potential is V(x) = k x. This tilts the harmonic oscillator bowl, shifting the most likely position to -x_k. The bell curve keeps its shape but moves left or right, changing the probability that the option finishes in‑the‑money.

Effective Volatility: Computing σ_eff from Quantum Mechanical Expectation Values

Volatility measures how much the price jumps around. In the standard Black–Scholes world, the volatility σ is a single number that controls the width of the bell‑shaped distribution of log‑returns. When we introduce a market force, the probability density P(x) changes shape—it may shift, narrow, or even develop multiple peaks. The paper captures the effect of this new shape on option prices through an effective volatility. This section explains what the effective volatility is, how it is computed from the quantum mechanical probability density, and how the C++ code implements the calculation.

The Definition of Effective Volatility

The paper defines the effective volatility as a rescaling of the original Black–Scholes volatility by the standard deviation of the quantum probability density. The canonical equation from the paper (Eq. 11) is

\(\sigma_{\text{eff}} = \sigma \sigma_{\text{QM}}, \quad \sigma_{\text{QM}} = \sqrt{\langle x^2 \rangle_{\text{QM}} - \langle x \rangle_{\text{QM}}^2}\)

In this formula σ is the original volatility (for example 0.25), σ_QM is the quantum mechanical standard deviation of the dimensionless variable x under the market‑force‑modified probability density P(x), and σ_eff is the effective volatility that will be used in the modified Black–Scholes formula. The notation ⟨·⟩_QM denotes the expectation value with respect to P(x): ⟨x⟩_QM is the mean of x, and ⟨x²⟩_QM is the second moment.

In plain terms σ_QM is simply the standard deviation of the new probability curve. If the market force squeezes the distribution, σ_QM becomes smaller than 1 and the effective volatility drops. If the force stretches the distribution, σ_QM becomes larger than 1 and the effective volatility rises. For the baseline normal density P_BS(x) = (1/√(2π)) e^{-x²/2} the standard deviation is exactly 1, so σ_eff = σ and we recover the original Black–Scholes model.

Computing Expectation Values Numerically

The paper does not provide analytical formulas for ⟨x⟩ and ⟨x²⟩ for every force model. Instead the implementation computes them by numerical integration. The algorithm follows the method card effective_volatility_computation:

1. Choose an integration domain [a, b] that captures virtually all probability mass. For Gaussian‑like densities [-12, 12] is sufficient.

2. Build a uniform grid of N points (default 20001, which is odd and suitable for Simpson’s rule).

Use the button or URL below to download the source code

Evaluate P(x), x·P(x), and x²·P(x) at each grid point.

3. Compute ⟨x⟩ and ⟨x²⟩ using Simpson integration.

4. Compute σ_QM = √(⟨x²⟩ – ⟨x⟩²).

5. Return σ_eff = σ · σ_QM.

The choice of domain and number of points is a practical decision. The interval [-12, 12] covers more than 99.9% of the probability mass for any Gaussian‑like density. For densities with a non‑zero mean, such as the constant force model, the domain must be shifted accordingly to ensure the peak is well inside the integration range. The code uses 20001 points to achieve a relative error below 1e-12 for smooth densities.

The C++ Implementation

The computation lives in the header‑only file src/effective_volatility.h. It provides three main pieces: a general Simpson integrator, the VolResult struct, and the compute_eff_vol function.

Simpson Integration

The workhorse is simpson_integrate, which implements Simpson’s 1/3 rule for a uniformly spaced grid:

inline double simpson_integrate(const std::vector<double>& x,
                                const std::vector<double>& f) {
    const size_t n = x.size();
    if (n != f.size()) {
        throw std::invalid_argument(
            "simpson_integrate: x and f must have the same size");
    }
    if (n < 3 || n % 2 == 0) {
        throw std::invalid_argument(
            "simpson_integrate: size must be odd and >= 3");
    }

    const double dx = x[1] - x[0];
    double sum = f[0] + f[n - 1];

    for (size_t i = 1; i < n - 1; ++i) {
        sum += (i % 2 == 1) ? 4.0 * f[i] : 2.0 * f[i];
    }

    return (dx / 3.0) * sum;
}

This function is used throughout the codebase for expectation values, normalisation constants, and cumulative distribution functions. It requires an odd number of points (at least 3) and a constant step size dx. The pattern 4,2,4,2,… is the classic Simpson weight sequence.

The Result Struct

The VolResult struct bundles the outputs:

struct VolResult {
    double sigma_eff;  // Effective volatility: σ * σ_QM
    double sigma_QM;   // Quantum mechanical standard deviation of x
    double mean_x;     // Expectation value <x>
    double var_x;      // Variance <x^2> - <x>^2
};

All four fields are returned so that callers can inspect the mean, variance, and both volatility values.

The Main Function: compute_eff_vol

The function compute_eff_vol implements the algorithm described above. Here is the core logic (the full file contains input validation and a corrected version that properly assigns all fields):

inline VolResult compute_eff_vol(double sigma,
                                 std::function<double(double)> P,
                                 double a = -12.0,
                                 double b = 12.0,
                                 int N = 20001) {
    // ... input validation ...

    const double dx = (b - a) / static_cast<double>(N - 1);
    std::vector<double> x(N);
    std::vector<double> P_vals(N);
    std::vector<double> xP_vals(N);
    std::vector<double> x2P_vals(N);

    for (int i = 0; i < N; ++i) {
        x[i] = a + static_cast<double>(i) * dx;
        P_vals[i] = P(x[i]);
        xP_vals[i] = x[i] * P_vals[i];
        x2P_vals[i] = x[i] * x[i] * P_vals[i];
    }

    const double mean_x = simpson_integrate(x, xP_vals);
    const double mean_x2 = simpson_integrate(x, x2P_vals);
    const double var_x = mean_x2 - mean_x * mean_x;

    if (var_x <= 0.0) {
        throw std::runtime_error(
            "compute_eff_vol: variance is non-positive; "
            "check probability density");
    }

    const double sigma_QM = std::sqrt(var_x);
    const double sigma_eff = sigma * sigma_QM;

    return VolResult{sigma_eff, sigma_QM, mean_x, var_x};
}

Notice how the function builds three parallel arrays: P_vals, xP_vals, and x2P_vals. Simpson integration is then applied to each to obtain the mean and the second moment. (The full file also computes the integral of P itself for verification, though that value is not returned in the snippet above.) The variance is computed as ⟨x²⟩ – ⟨x⟩², and the quantum standard deviation is its square root. The effective volatility is simply the product σ * σ_QM.

Why σ_QM = 1 for the Baseline Normal Density

For the standard normal probability density P_BS(x) = (1/√(2π)) e^{-x²/2}, the mean is 0 and the variance is 1. Therefore σ_QM = 1 and σ_eff = σ. The code verifies this property in the verify_effective_volatility function:

VolResult res = compute_eff_vol(0.25, normalPDF, -12.0, 12.0, 20001);
// Check σ_QM ≈ 1.0
if (std::abs(res.sigma_QM - 1.0) > tol) { /* fail */ }
// Check mean_x ≈ 0.0
if (std::abs(res.mean_x) > tol) { /* fail */ }
// Check var_x ≈ 1.0
if (std::abs(res.var_x - 1.0) > tol) { /* fail */ }
// Check σ_eff = σ * σ_QM
if (std::abs(res.sigma_eff - 0.25) > tol) { /* fail */ }

These checks ensure that the numerical integration is accurate and that the baseline model is correctly recovered before any market force is applied.

Worked Example

Let’s walk through a concrete call. Suppose we have the original volatility σ = 0.25 and we want the effective volatility for the standard normal density. We call

VolResult res = compute_eff_vol(0.25, normalPDF);

Internally, the function builds a grid of 20001 points from -12.0 to 12.0, evaluates the normal PDF at each point, and computes the integrals. The result is res.sigma_QM ≈ 1.0, res.sigma_eff ≈ 0.25, res.mean_x ≈ 0.0, and res.var_x ≈ 1.0.

Now consider a market force that narrows the distribution, such as the linear force with λ = 0.5. The probability density becomes a Gaussian with variance 1/λ_ω = 1/√(1+λ) ≈ 0.816. The quantum standard deviation is then σ_QM ≈ 0.904, and the effective volatility drops to σ_eff ≈ 0.226. This reduction in effective volatility will lower the option price because extreme price moves become less likely.

Summary

The effective volatility σ_eff is the bridge between the quantum probability density and the Black–Scholes pricing formula. By computing the standard deviation of the modified density and scaling the original volatility, we preserve the structure of the Black–Scholes equations while incorporating the non‑Gaussian features introduced by market forces. The C++ implementation uses Simpson’s rule to evaluate the necessary expectation values numerically, with careful choice of integration domain and resolution to ensure high accuracy. The built‑in verification checks confirm that the baseline normal density yields σ_QM = 1, giving confidence that the machinery works correctly before we apply it to the force models in the following sections.

Modified Option Pricing: Integrating the New Density and Adjusting d±

Now that we have a market‑force‑driven probability density P(x) and an effective volatility sigma_eff, we need to combine them to obtain option prices. The key idea is to preserve the structure of the Black‑Scholes formulas but replace the standard normal cumulative distribution function N(·) with a new CDF derived from P(x), and to use sigma_eff when computing the integration limits d±. This section walks through the mathematics and the C++ implementation that brings all the pieces together.

Modified Integration Limits

The standard d± arguments (Eq. 3 from the paper) are defined by

\(d_{\pm} = \frac{\ln(S_0/K) + (r \pm \sigma^2/2)T}{\sigma\sqrt{T}}\)

To incorporate the market force we simply replace the original volatility σ with the effective volatility σ_eff obtained from the quantum mechanical expectation values. This gives the modified limits d_eff±. No new equation is introduced; the replacement is applied directly in the code by calling the same compute_dpm helper with sigma_eff instead of sigma.

Modified Black‑Scholes Formulas

The call and put pricing formulas (Eqs. 1 and 2) keep their algebraic form but now reference the new cumulative distribution N_eff(·) and the modified limits d_eff±. First, recall the original Black‑Scholes formulas:

\(c = S_0 N(d_+) - K e^{-rT} N(d_-)\)

\(p = K e^{-rT} N(-d_-) - S_0 N(-d_+)\)

In the modified version, N(·) is replaced by N_eff(·) and d± by d_eff±. The put formula is derived from put‑call parity and uses the same N_eff values.

Building the Cumulative Distribution Grid

Evaluating N_eff(d) on the fly for every option pricing call would be slow and numerically expensive. Instead, the implementation pre‑computes a cumulative distribution grid on a uniform partition of a sufficiently wide domain (the paper uses [-12, 12] with 20001 points) using Simpson’s rule. The function build_cdf_grid performs this task.

struct ModifiedBSResult {
    double call;         // Modified call option price
    double put;          // Modified put option price
    double d_eff_plus;   // d_eff+ using effective volatility
    double d_eff_minus;  // d_eff- using effective volatility
    double N_eff_plus;   // N_eff(d_eff+)
    double N_eff_minus;  // N_eff(d_eff-)
};

// Build a cumulative distribution grid from a probability density P(x).
// Returns a pair (x_grid, cdf_grid) where cdf_grid[i] = ∫_{a}^{x_grid[i]} P(t) dt.
std::pair<std::vector<double>, std::vector<double>>
build_cdf_grid(std::function<double(double)> P, double a, double b, int N);

Inside build_cdf_grid, the density P(x) is evaluated at every grid point, and a cumulative Simpson integration is applied. The first element of cdf_grid is always 0.0, and for a properly normalised P, the last element should be very close to 1.0.

Interpolating the CDF at Arbitrary Points

Once the grid is built, N_eff(d) can be evaluated for any d by linear interpolation between the two nearest grid points. The helper interpolate_cdf handles that, including clamping to 0.0 for values below the domain and to cdf_grid.back() for values above it.

// Interpolate the CDF at an arbitrary point d.
double interpolate_cdf(double d,
                       const std::vector<double>& x_grid,
                       const std::vector<double>& cdf_grid);

Binary search locates the interval containing d, and a linear blend yields the interpolated value. This keeps the option pricing fast even when the grid is very fine.

Bringing It All Together: modifiedoptionpricing

The main function modified_option_pricing accepts the financial parameters (S0, K, r, T), the effective volatility sigma_eff, the density function P (kept only for documentation, not called inside the function), and the pre‑computed grids. It then:

1. Computes d_eff± by calling the standard compute_dpm with sigma_eff.

2. Evaluates N_eff(d_eff±) via interpolate_cdf.

3. Applies the modified call and put formulas.

ModifiedBSResult modified_option_pricing(
    double S0,
    double K,
    double r,
    double T,
    double sigma_eff,
    std::function<double(double)> /*P*/,
    const std::vector<double>& x_grid,
    const std::vector<double>& cdf_grid)
{
    // Compute d_eff± using sigma_eff
    auto [d_eff_plus, d_eff_minus] = compute_dpm(S0, K, r, sigma_eff, T);

    // Evaluate modified CDF
    double N_eff_plus  = interpolate_cdf(d_eff_plus,  x_grid, cdf_grid);
    double N_eff_minus = interpolate_cdf(d_eff_minus, x_grid, cdf_grid);

    double discount = std::exp(-r * T);
    double call = S0 * N_eff_plus - K * discount * N_eff_minus;
    double put  = K * discount * (1.0 - N_eff_minus) - S0 * (1.0 - N_eff_plus);

    return {call, put, d_eff_plus, d_eff_minus, N_eff_plus, N_eff_minus};
}

Verification: Recovering the Baseline

A crucial sanity check is that when we supply the standard normal PDF and use sigma_eff = sigma, the modified pricing must reproduce the original Black‑Scholes prices. The implementation includes a verification function verify_modified_option_pricing that performs this test alongside monotonicity and put‑call parity checks. Using the paper’s numeric example (S0 = 20, K = 20, r = 0.1, T = 1, σ = 0.25) the modified call and put prices match the baseline values within a tolerance of 1e-10.

Putting It All Together

To compute an option price under a given market force model, you would:

1. Choose a force strength (e.g., λ for the linear force) and obtain the corresponding probability density P(x) (e.g., linear_force_PDF(x, λ)).

2. Compute the effective volatility by calling compute_eff_vol(sigma, P, -12.0, 12.0, 20001), which returns sigma_eff.

3. Build the CDF grid with build_cdf_grid(P, -12.0, 12.0, 20001).

4. Call modified_option_pricing(S0, K, r, T, sigma_eff, P, x_grid, cdf_grid) to obtain the modified call and put prices.

When P is the standard normal PDF, the result collapses identically to the original Black‑Scholes price, confirming the internal consistency of the approach.

Model 1: Constant Market Force (Shifted Gaussian)

A constant market force is the simplest departure from the Black–Scholes baseline. Imagine a steady wind that pushes the marble to one side of the bowl. The bell‑shaped probability density keeps its shape but shifts left or right. In financial terms, this corresponds to a net buying or selling pressure that changes the most likely outcome without altering the size of typical fluctuations. This section builds the constant force model from the paper, explains why the effective volatility remains unchanged, and walks through the C++ implementation that reproduces Figure 1.

The Potential and the Force

The paper introduces a linear potential V_k(x) = k x to model a constant force, where k is the force strength constant. A positive k gives a force F(x) = -dV/dx = -k, which pushes the distribution to the left (negative direction). A negative k pushes it to the right. The potential is added to the harmonic oscillator baseline, tilting the bowl.

The Shifted Gaussian Density

Solving the Schrödinger equation with the linear potential yields a ground state wavefunction whose squared modulus is a shifted Gaussian. The paper’s Eq. (15) is garbled in the original text; we reconstruct the probability density as a standard normal distribution shifted by -x_k:

\(P_k(x) = \frac{1}{\sqrt{2\pi}} e^{-(x + x_k)^2/2}\)

Here P_k(x) is the probability density at the dimensionless coordinate x, and x_k is the peak shift. The shift is related to the force strength by x_k = 2k/(\hbar\omega). Following the Appendix (Sec. 5.1) we set the energy scale ℏω = 1, so the shift simplifies to x_k = 2k. The density is a Gaussian with mean -x_k and variance 1.

Effective Volatility

The paper defines the effective volatility as the product of the original volatility and the quantum standard deviation:

\(\sigma_{\text{eff}} = \sigma \sigma_{\text{QM}}, \quad \sigma_{\text{QM}} = \sqrt{\langle x^2 \rangle_{\text{QM}} - \langle x \rangle_{\text{QM}}^2}\)

For the constant force model, the variance of the shifted Gaussian is unchanged, so σ_QM = 1 and therefore σ_eff = σ. This is a crucial property: a constant force shifts the distribution but does not change its width, so the effective volatility is the same as the Black–Scholes volatility. The only effect on option prices comes from the altered moneyness through the modified integration limits d_eff±.

Implementation in C++

The constant force model is implemented in the header src/constant_force_model.h. Two functions capture the mathematics.

Computing the shift – The function constant_force_xk returns x_k = 2k:

inline double constant_force_xk(double k) {
    // Eq. (14): x_k = 2k / (hbar * omega), with hbar*omega = 1
    return 2.0 * k;
}

Evaluating the density – The function constant_force_PDF implements the reconstructed Eq. (15). It computes the shifted argument x + x_k and evaluates the standard normal PDF:

inline double constant_force_PDF(double x, double k) {
    const double xk = constant_force_xk(k);
    const double shifted = x + xk;  // mean is -x_k, so (x - (-x_k)) = x + x_k
    constexpr double inv_sqrt_2pi = 0.3989422804014327;  // 1 / sqrt(2*pi)
    return inv_sqrt_2pi * std::exp(-0.5 * shifted * shifted);
}

Notice that the density is always positive and integrates to 1 over the real line. For k = 0 the shift x_k is zero, and the function returns exactly the standard normal PDF used in the baseline Black–Scholes model.

How the Shift Affects Option Prices

Because σ_eff = σ, the modified integration limits d_eff± are computed with the original volatility. However, the cumulative distribution function N_eff(d_eff±) is obtained by integrating the shifted density P_k(x) from -∞ to d_eff±. The shift changes the probability mass to the left or right of the strike, altering the option’s moneyness.

Consider a call option with the paper’s synthetic parameters: S0 = 20, K = 20, r = 0.1, T = 1, σ = 0.25. For k = 0.5 we have x_k = 1.0. The density is centred at -1.0, meaning the distribution is shifted to the left. The probability of finishing in‑the‑money decreases, so the call price falls below the baseline value of approximately 3.022. Conversely, a negative k shifts the density to the right, increasing the call price. The put price moves in the opposite direction.

Verification

The header also provides a verification function verify_constant_force_model that checks three invariants:

1. For k = 0, P_k(x) equals the standard normal PDF at several test points.

2. The integral of P_k over the domain [-12, 12] is 1.0 within a tight tolerance for several values of k.

3. The quantum standard deviation σ_QM computed numerically is 1.0 regardless of the shift.

These checks ensure that the implementation faithfully reproduces the paper’s shifted Gaussian and that the effective volatility remains unchanged. The verification is designed to be called from the main program after all headers are included.

Worked Example

To see the model in action, suppose we set k = 0.5. Then x_k = 1.0. The density P_k(x) is a Gaussian centred at -1.0. Using the modified pricing pipeline (building a CDF grid from constant_force_PDF, computing σ_eff = σ, and evaluating modified_option_pricing), the call price drops from the baseline 3.022 to a lower value, while the put price increases. The parameter sweep function sweep_constant_force in src/paper_figures.h automates this process over a range of k values, producing the data for Figure 1 of the paper.

In summary, the constant force model is the simplest illustration of how a market potential can alter option prices. It shifts the probability density without changing its shape, and the implementation is a direct translation of the reconstructed Eq. (15). The next section will examine the linear force, which modifies the variance and therefore the effective volatility.

Model 2: Linear Market Force (Modified Variance Gaussian)

A linear force is like making the bowl steeper or shallower. A steeper bowl (λ > 0) squeezes the bell curve, making extreme moves less likely. This lowers the option price because big payoffs are rarer. In this section we implement the linear market force model from the paper, which modifies the variance of the probability density while keeping its Gaussian shape. We will see how the effective frequency parameter λ_ω emerges from the quadratic potential, how the analytical standard deviation σ_QM avoids numerical integration, and how the C++ code in linear_force_model.h brings everything together.

The Quadratic Potential and the Modified Frequency

The paper adds a quadratic potential to the harmonic oscillator baseline: V_λ(x) = λ x² with λ > 0. The corresponding market force is F(x) = -dV/dx = -2λ x, a restoring force proportional to the displacement from equilibrium. Adding this term changes the effective frequency of the system. With the conventions ℏω = 1, m = 1, and ω = 1 (Appendix Sec. 5.1), the new frequency becomes ω' = ω √(1 + λ/ω) = √(1 + λ). The paper defines the effective frequency parameter as λ_ω = √(1 + λ/ω), which simplifies to λ_ω = √(1 + λ) when ω = 1. This parameter controls the width of the ground state probability density.

The Modified Gaussian Density

The ground state wavefunction for the modified harmonic oscillator yields a probability density that is still a Gaussian, but with a variance that depends on λ. The canonical equation from the paper (Eq. 18) is

\(P_\lambda(x) = \frac{\lambda_\omega}{\sqrt{2\pi}} e^{-\lambda_\omega x^2/2}\)

Here P_λ(x) is the probability density at the dimensionless coordinate x, and λ_ω = √(1 + λ). The density is always positive and integrates to 1 over the real line. When λ = 0, λ_ω = 1 and P_λ(x) reduces to the standard normal density P_BS(x). As λ increases, the variance 1/λ_ω decreases, making the distribution narrower.

Effective Volatility for the Linear Force

Because the density is a zero‑mean Gaussian, the quantum mechanical standard deviation σ_QM can be computed analytically. From the definition in Eq. (11),

\(\sigma_{\text{eff}} = \sigma \sigma_{\text{QM}}, \quad \sigma_{\text{QM}} = \sqrt{\langle x^2 \rangle_{\text{QM}} - \langle x \rangle_{\text{QM}}^2}\)

For P_λ(x) the mean is zero and the second moment is 1/λ_ω, so σ_QM = 1/λ_ω. The effective volatility becomes σ_eff = σ / λ_ω = σ / √(1 + λ). This analytical result avoids numerical integration for the linear force model. A larger λ (stronger restoring force) reduces σ_eff, which in turn lowers option prices because the distribution of log‑returns is less spread out.

C++ Implementation

The header‑only file src/linear_force_model.h provides three core functions and a verification routine. Let’s examine each.

Computing `λ_ω` and `σ_QM`

The functions linear_force_lambda_omega and linear_force_sigma_QM encapsulate the analytical formulas. They validate that λ > 0 and return the exact values.

inline double linear_force_lambda_omega(double lambda) {
    if (lambda <= 0.0) {
        throw std::domain_error("linear_force_lambda_omega: lambda must be > 0, got "
                                + std::to_string(lambda));
    }
    return std::sqrt(1.0 + lambda);
}

inline double linear_force_sigma_QM(double lambda) {
    if (lambda <= 0.0) {
        throw std::domain_error("linear_force_sigma_QM: lambda must be > 0, got "
                                + std::to_string(lambda));
    }
    return 1.0 / std::sqrt(1.0 + lambda);
}

Notice that linear_force_sigma_QM simply returns 1.0 / λ_ω, matching the analytical derivation. Both functions throw a std::domain_error if λ is not strictly positive, as required by the paper’s Eq. (16).

Evaluating the Probability Density

The function linear_force_PDF implements Eq. (18) directly.

inline double linear_force_PDF(double x, double lambda) {
    if (lambda <= 0.0) {
        throw std::domain_error("linear_force_PDF: lambda must be > 0, got "
                                + std::to_string(lambda));
    }
    const double lambda_omega = std::sqrt(1.0 + lambda);
    const double norm = lambda_omega / std::sqrt(2.0 * M_PI);
    return norm * std::exp(-0.5 * lambda_omega * x * x);
}

The normalisation constant norm = λ_ω / √(2π) ensures the integral over the real line equals 1. The exponential uses -0.5 * λ_ω * x², which is the standard Gaussian exponent with variance 1/λ_ω. For large |x| the exponential underflows to zero, which is correct.

Verification

The function verify_linear_force_model runs several static checks:

• For a very small λ (e.g. 1e-10), P_λ(x) must match the standard normal PDF at several test points.

• The integral of P_λ over [-12, 12] must be approximately 1.0 for a range of λ values.

• The numerical σ_QM (computed via Simpson integration of x² P_λ(x)) must agree with the analytical linear_force_sigma_QM.

• The identity σ_QM = 1/λ_ω must hold exactly.

These checks confirm that the implementation is consistent with the paper’s equations and that the density is properly normalised.

Worked Example

Let’s use the paper’s synthetic parameters: S0 = 20, K = 20, r = 0.1, T = 1, σ = 0.25. For λ = 0.5 we compute:

• λ_ω = √(1 + 0.5) ≈ 1.2247

• σ_QM = 1 / 1.2247 ≈ 0.8165

• σ_eff = 0.25 × 0.8165 ≈ 0.2041

The baseline Black‑Scholes call price (with σ = 0.25) is approximately 3.022. With the reduced effective volatility, the call price drops to about 2.4. This makes sense: a narrower distribution means the asset is less likely to finish far in‑the‑money, so the option is worth less.

In code, you would obtain the density at any point with linear_force_PDF(x, 0.5). To reproduce Figure 2 from the paper, the sweep_linear_force function (in paper_figures.h) loops over a range of λ values, builds the CDF grid from linear_force_PDF, computes σ_eff analytically, and calls modified_option_pricing. The resulting call and put prices are stored in a SweepResult struct for plotting.

Key Takeaways

• The linear force model corresponds to a quadratic potential V(x) = λ x² that changes the harmonic oscillator frequency.

• The probability density remains Gaussian, but its variance becomes 1/λ_ω with λ_ω = √(1 + λ).

• The quantum standard deviation σ_QM and the effective volatility σ_eff can be computed analytically, avoiding numerical integration.

• As λ increases, the distribution narrows, reducing option prices.

• The implementation in linear_force_model.h is self‑contained, validates inputs, and includes built‑in verification checks that ensure correctness against the paper’s equations.

Model 3: x² Market Force (Cubic Potential, Perturbative) – with Error Discussion

An x² force is like a bowl that is steeper on one side than the other. The paper attempts to approximate the new bell‑shaped probability density using a mathematical technique called perturbation theory, but the approximation contains a mistake. In this section we will implement the model exactly as the paper presents it, explain why the linear term in β is incorrect, and then show how to switch to a corrected second‑order expansion using a compile‑time flag. The C++ code lives in src/x2_force_model.h and provides both the paper's formula and the optional corrected version.

The Cubic Potential and the x² Force

The paper introduces a cubic potential V_β(x) = β x³ to model a force that grows quadratically with the displacement x. Here β is the force strength parameter. The corresponding market force is F(x) = -dV/dx = -3β x². A positive β creates a force that pushes the distribution to the left for positive x and to the right for negative x, skewing the probability density. Because the potential is an odd function, it breaks the left‑right symmetry of the harmonic oscillator.

Perturbation Theory and the Paper's Formula

To obtain the ground state wavefunction under the cubic potential, the paper uses non‑degenerate perturbation theory. The general first‑order correction to a wavefunction is given by

\(\psi_n = \psi_n^{(0)} + \sum_{m \neq n} \frac{H'_{nm}}{E_n^{(0)} - E_m^{(0)}} \psi_m^{(0)} + \ldots\)

where ψ_n^{(0)} are the unperturbed eigenstates, E_n^{(0)} the unperturbed energies, and H'_{nm} the matrix elements of the perturbation V_β(x). The paper then squares the approximate wavefunction to obtain the probability density. The claimed result (Eq. 22) is

\(P_\beta(x) = \frac{C^2}{\sqrt{2\pi}} e^{-x^2/2} \left[1 - \frac{\beta}{\hbar\omega}\frac{1}{3}x^3\right]^2\)

In this expression C² is a normalisation constant that must be computed numerically because the perturbation expansion does not preserve the unit integral. The factor ℏω is the energy quantum; following Appendix Sec. 5.1 we set ℏω = 1 throughout the implementation, so the term simplifies to (β/3) x³.

The Error in the Paper's Result

There is a fundamental problem with the formula above. The unperturbed ground state of the harmonic oscillator is an even function, while the cubic perturbation β x³ is odd. The first‑order correction to the wavefunction involves matrix elements H'_{nm} that couple the ground state to odd‑parity excited states. However, the first‑order correction to the probability density (the squared modulus of a real wavefunction) contains a cross term 2 ψ^{(0)} ψ^{(1)}. Because ψ^{(0)} is even and ψ^{(1)} is odd, their product is odd, and the integral of the cross term over the whole real line vanishes. Therefore the first‑order correction to the normalised probability density is zero. The leading correction is of order β², not linear in β. The paper's Eq. (22) erroneously retains a linear β term, which is spurious.

Despite this error, we provide the paper's formula for reproducibility. A compile‑time flag USE_CORRECTED_X2 can be defined to replace it with a correct second‑order expansion derived from standard quantum mechanics.

Implementation of the Paper's Formula

The header src/x2_force_model.h contains three essential functions for the paper's model. First, the unnormalised raw density is computed by x2_force_PDF_raw:

inline double x2_force_PDF_raw(double x, double beta) {
    const double x3 = x * x * x;
    const double correction = 1.0 - (beta / 3.0) * x3;
    const double gaussian = std::exp(-0.5 * x * x);
    return gaussian * correction * correction;
}

This directly implements the bracket [1 - (β/3) x³]² multiplied by the Gaussian factor exp(-x²/2). Notice that the factor 1/√(2π) is omitted because the normalisation constant C² will absorb it.

Next, the normalisation constant C² is obtained by numerical integration over the domain [-12, 12] with 20001 points (Simpson's rule):

inline double compute_x2_normalization(double beta, double a = -12.0, double b = 12.0, int N = 20001) {
    if (N < 3 || N % 2 == 0) {
        throw std::invalid_argument("compute_x2_normalization: N must be odd and >= 3, got " + std::to_string(N));
    }
    if (a >= b) {
        throw std::invalid_argument("compute_x2_normalization: a must be less than b");
    }

    const double dx = (b - a) / (N - 1);
    std::vector<double> x(N), f(N);
    for (int i = 0; i < N; ++i) {
        x[i] = a + i * dx;
        f[i] = x2_force_PDF_raw(x[i], beta);
    }

    double integral = x2_force_forward::simpson_integrate(x, f);

    if (integral <= 0.0) {
        throw std::runtime_error("compute_x2_normalization: integral of P_raw is non-positive ("
                                 + std::to_string(integral) + "). The perturbation expansion "
                                 "may have broken down for beta = " + std::to_string(beta));
    }

    return 1.0 / integral;
}

The function validates its inputs, evaluates the raw density on a uniform grid, integrates with Simpson's rule, and throws an exception if the integral is non‑positive—which can happen for large |β| when the truncated expansion breaks down and the raw density becomes negative over most of the domain.

Finally, the normalised density is a simple scaling:

inline double x2_force_PDF(double x, double beta, double C2) {
    return C2 * x2_force_PDF_raw(x, beta);
}

The Corrected Second‑Order Perturbation (Optional)

When USE_CORRECTED_X2 is defined, the header provides an alternative set of functions that use the correct O(β²) expansion. The corrected raw density is

inline double x2_force_PDF_corrected_raw(double x, double beta) {
    const double x2 = x * x;
    const double x4 = x2 * x2;
    const double x6 = x4 * x2;
    const double Q = (4.0 * x6 - 36.0 * x4 + 81.0 * x2 - 15.0) / 72.0;
    const double beta2 = beta * beta;
    const double gaussian = std::exp(-0.5 * x2);
    return gaussian * (1.0 + beta2 * Q);
}

This formula is not from the paper; it is derived from standard quantum mechanics. The polynomial Q(x) ensures that the density remains symmetric and non‑negative for sufficiently small β. The corresponding normalisation function compute_x2_normalization_corrected and the normalised accessor x2_force_PDF_corrected follow the same pattern as the paper's version.

Worked Example

Let us set β = 0.1. First we compute the normalisation constant:

double C2 = compute_x2_normalization(0.1);

This integrates the raw density over [-12, 12] and returns C² ≈ 1.0 (the deviation from unity is small because the perturbation is weak). Then we can evaluate the normalised density at any point:

double p = x2_force_PDF(0.0, 0.1, C2);  // near 0.3989

To price an option, we would build a CDF grid from this density, compute σ_eff via compute_eff_vol, and call modified_option_pricing. The call price decreases slightly compared to the Black‑Scholes baseline because the distribution becomes skewed, reducing the probability of large positive returns.

Verification

The file includes a verify_x2_force_model() function that runs several static checks:

• For β = 0, the density must equal the standard normal PDF pointwise.

• The normalised density must integrate to 1.0 within a tight tolerance for several small β values.

• A warning is emitted if the density becomes negative anywhere in the domain, indicating that the perturbation expansion has broken down.

If USE_CORRECTED_X2 is defined, analogous checks are performed for the corrected density. These verifications help ensure that the implementation faithfully reproduces the paper's formulas (with the known caveat) and that the optional corrected version behaves correctly.

Summary

The x² force model illustrates both the power and the pitfalls of perturbation theory. The paper's formula, while algebraically simple, contains a linear β term that should not be present. Our implementation provides the paper's version for exact reproducibility and an optional corrected second‑order expansion for those who need a physically consistent density. In the next section we will examine the x³ force, where the first‑order correction is non‑zero and leads to a richer, non‑monotonic behaviour.

Model 4: x³ Market Force (Quartic Potential, Perturbative)

An x³ force is like a bowl that becomes flatter at the bottom and then rises steeply. For strong forces, the marble can sit in two places, creating a double‑humped probability distribution. This can make options cheaper or more expensive depending on where the strike price falls relative to the two peaks. In this section we implement the x³ market force model from the paper using first‑order perturbation theory, walk through the C++ code in src/x3_force_model.h, and explain why the option price does not change monotonically with the force strength γ.

The Quartic Potential

The paper introduces a quartic potential to model an x³ force. The potential is V_γ(x) = γ * x^4 with γ > 0. The corresponding market force is F(x) = -dV/dx = -4 * γ * x^3. Because the potential is symmetric (an even function of x), the force is an odd function: it pulls the distribution toward the origin for positive γ, making the bottom of the effective bowl flatter than the pure harmonic oscillator. For large enough γ the total potential ½ m ω² x² + γ x⁴ develops a double‑well shape, and the probability density can become bimodal—the marble can sit on either side of the origin.

Perturbation Theory for the Wavefunction

The paper uses non‑degenerate perturbation theory (Eq. 33 in the Appendix) to obtain an approximate ground state wavefunction. The unperturbed system is the harmonic oscillator, and the perturbation is H' = γ x⁴. The first‑order correction to the wavefunction is given by the canonical perturbation expansion:

\(\psi_n = \psi_n^{(0)} + \sum_{m \neq n} \frac{H'_{nm}}{E_n^{(0)} - E_m^{(0)}} \psi_m^{(0)} + \ldots\)

Here ψ_n^(0) are the unperturbed harmonic oscillator eigenstates, E_n^(0) are the unperturbed energies, and H'_{nm} is the matrix element of the perturbation between states n and m. For the ground state (n = 0), the quartic potential couples it to the m = 2 and m = 4 states. The paper gives the resulting ground state wavefunction (Eq. 25) as

\(\psi_{g,\gamma}(x) = \frac{C}{(2\pi)^{1/4}} e^{-x^2/4} \left[1 - \frac{\gamma}{4\sqrt{2}\hbar\omega}(x^4 - 9)\right]\)

With the convention ℏω = 1 (Appendix Sec. 5.1), the factor simplifies to γ / (4√2).

The Probability Density

The probability density is the squared modulus of the wavefunction. Because the ground state is real, this is simply P_γ(x) = |ψ_{g,γ}(x)|². The paper writes it as (Eq. 26)

\(P_\gamma(x) = \frac{C^2}{\sqrt{2\pi}} e^{-x^2/2} \left[1 - \frac{\gamma}{4\sqrt{2}\hbar\omega}(x^4 - 9)\right]^2\)

This is a Gaussian multiplied by a quartic correction. The constant 9 emerges from the matrix elements of the perturbation; its origin is derived in the paper's Appendix. For γ = 0 the expression reduces exactly to the standard normal PDF P_BS(x).

With ℏω = 1 the density simplifies to the raw unnormalized form used in the code:

\(P_{\text{raw}}(x) = e^{-x^2/2} \left[1 - \frac{\gamma}{4\sqrt{2}}(x^4 - 9)\right]^2\)

Important caveat: The perturbation expansion is valid only for small γ. For large γ the correction can overcome the Gaussian envelope, making P_γ(x) negative for some x. The C++ code detects this and emits a warning, because a negative probability density is unphysical and signals the breakdown of the first‑order approximation.

Normalisation

Perturbation theory does not preserve the normalization of the wavefunction. The constant C² must therefore be computed numerically so that ∫ P_γ(x) dx = 1 over the chosen integration domain. The paper does not give an analytical value for C²; we compute it by Simpson integration.

Implementation in src/x3forcemodel.h

The header‑only file src/x3_force_model.h provides three public functions and a verification routine. Let us examine each one.

The Raw Unnormalized Density

The function x3_force_PDF_raw implements the bracket part of Eq. 26 without the normalization constant:

inline double x3_force_PDF_raw(double x, double gamma) {
    const double factor = gamma / (4.0 * std::sqrt(2.0));  // γ/(4√2)
    const double x2 = x * x;
    const double x4 = x2 * x2;
    const double bracket = 1.0 - factor * (x4 - 9.0);
    return std::exp(-0.5 * x2) * bracket * bracket;
}

Notice how ℏω is set to 1, consistent with Appendix Sec. 5.1, so the denominator 4√2 ℏω becomes simply 4√2. The function computes x², then x⁴, forms the bracket [1 - γ/(4√2) (x⁴ - 9)], squares it, and multiplies by the Gaussian envelope exp(-x²/2). The squaring guarantees the raw density is non‑negative, which is a mathematical identity from the squared modulus. However, the normalized density C² * P_raw can still become negative if C² is negative—which we prevent by validating the integral.

Numerical Normalization

The function compute_x3_normalization computes C²:

inline double compute_x3_normalization(double gamma, double a, double b, int N) {
    // Input validation ...
    const double dx = (b - a) / (N - 1);
    std::vector<double> x(N), f(N);
    for (int i = 0; i < N; ++i) {
        x[i] = a + i * dx;
        f[i] = x3_force_PDF_raw(x[i], gamma);
    }
    const double integral = simpson_integrate(x, f);
    if (integral <= 0.0) {
        throw std::runtime_error("Integral of P_raw is non-positive ...");
    }
    return 1.0 / integral;
}

It creates a uniform grid on [a, b] (default [-12, 12] with 20001 points), evaluates the raw density at every grid point, and integrates with Simpson's rule. If the integral is non‑positive, a runtime_error is thrown because the perturbation expansion has broken down. Otherwise it returns C² = 1 / integral. For γ = 0, this returns 1/√(2π), exactly the normal distribution's normalization.

The Normalized Density

inline double x3_force_PDF(double x, double gamma, double C2) {
    return C2 * x3_force_PDF_raw(x, gamma);
}

This multiplies the raw density by the pre‑computed C². All subsequent calculations—effective volatility, CDF grid construction, and option pricing—call this function.

How the Option Price Behaves

The paper's Figure 4 shows non‑monotonic behavior: as γ increases from zero, the call option price first decreases, reaches a minimum, and then increases again. This can be understood physically. For small γ, the quartic potential flattens the bottom of the well, which broadens the probability density slightly but also changes its shape. The effective standard deviation σ_QM may decrease, reducing σ_eff and lowering the option price (similar to the linear force case). For larger γ, the potential develops two minima, and the probability density becomes bimodal with peaks at non‑zero x. If the strike price happens to lie near one of the peaks, the probability of finishing in‑the‑money can actually increase, raising the option price. The non‑monotonicity is a signature of the double‑well physics.

Worked Example

Let γ = 0.2. First compute the normalization:

double C2 = compute_x3_normalization(0.2, -12.0, 12.0, 20001);

For γ = 0.2, C2 is approximately 0.3989 (close to 1/√(2π) ≈ 0.398942). To evaluate the density at, say, x = 1.5:

double density = x3_force_PDF(1.5, 0.2, C2);

The raw bracket at x = 1.5 is 1 - 0.2/(4√2) (5.0625 - 9) = 1 + 0.03536 × 3.9375 ≈ 1.139. Squaring gives ≈ 1.298. Multiplying by exp(-1.125) ≈ 0.3247 yields a raw density of ≈ 0.421. After multiplying by C2 ≈ 0.3989, the normalized density is roughly 0.168.

To compute the full option price, you would follow the same pipeline described in earlier sections: build the CDF grid from x3_force_PDF, compute σ_eff via compute_eff_vol (which numerically integrates x * P(x) and x² * P(x)), obtain d_eff± using σ_eff, interpolate the CDF at d_eff±, and finally apply the Black‑Scholes formula with the modified N_eff values.

When the Model Breaks Down

The verification function verify_x3_force_model (included in the same header) checks the γ = 0 recovery, the integral normalization, and whether the density becomes negative. For γ larger than about 1.0, the correction term γ/(4√2) (x⁴ - 9) can exceed 1 for large |x|, making the bracket negative. Although squaring makes the raw density positive, the normalized density can still dip below zero at extreme x because numerical errors in the normalization constant magnify tiny negative raw values. The verification emits a warning in this case:

WARNING: x3_force_model - For gamma=2.0, P_gamma becomes negative (min=-1.23e-06 at x=...). 
Perturbation expansion may be invalid for this gamma.

This is the signal that you have pushed the perturbative model beyond its regime of validity. For quantitative work, you would restrict γ to values where the density is everywhere non‑negative.

Summary

The x³ market force model introduces a quartic potential γ x⁴ that flattens the harmonic oscillator well and eventually creates a double‑well structure. The probability density is a Gaussian multiplied by a quartic correction, obtained from first‑order perturbation theory. The C++ implementation computes the raw density x3_force_PDF_raw, normalizes it numerically with compute_x3_normalization, and produces the final density x3_force_PDF. The resulting option price shows non‑monotonic behaviour—a minimum followed by an increase—which reflects the onset of bimodality in the underlying distribution. The model is a striking example of how quantum mechanical potentials can generate option‑price patterns that are impossible in the standard Black‑Scholes world.

Model 5: Quantum Well (Hard Price Bounds)

Imagine the price of an asset is trapped between two walls. It cannot go beyond a certain range, no matter what. This is like a currency that is pegged to a narrow band, or a stock subject to circuit breakers that halt trading when the price moves too far. In such a market, the probability of extreme returns is exactly zero, and the shape of the distribution inside the allowed band is not a simple bell curve. The paper models this situation with an infinite square well potential, the simplest quantum mechanical system that enforces hard boundaries. This section builds the quantum well model, derives its probability density and effective volatility, and walks through the C++ implementation in src/quantum_well_model.h.

The infinite square well potential

The potential that creates hard price bounds is defined piecewise in Eq. (27) of the paper:

\(V(x) = \begin{cases} 0 & |x| < a \\ \infty & |x| \ge a \end{cases}\)

Here a is the half‑width of the well, a positive number that controls how wide the allowed price range is. V(x) is the potential energy. Inside the well the potential is zero; outside it is infinite, meaning the price cannot exist there. The parameter a is qualitatively proportional to the price range ΔS; the paper does not give an exact mapping, but larger a means a wider allowed band.

Ground state wavefunction and probability density

Solving the stationary Schrödinger equation inside the well with the boundary condition that the wavefunction vanishes at x = ±a gives the ground state wavefunction. The infinite walls force the wavefunction to zero at the boundaries, which selects a standing‑wave solution—a sine function with exactly one half‑wavelength fitting inside the well. The result is Eq. (28):

\(\psi_g(x) = \frac{1}{\sqrt{a}} \sin\left[\frac{\pi}{2a}(x + a)\right], \quad x \in [-a, a]\)

ψ_g(x) is the ground state wavefunction. The factor 1/√a ensures the wavefunction is normalized. The argument π(x + a) / (2a) shifts and scales the sine so that it is zero at x = -a, reaches its maximum at x = 0, and returns to zero at x = a. Outside the well the wavefunction is identically zero.

The probability density is the squared modulus of the wavefunction, which yields Eq. (29):

\(P_{\text{QW}}(x) = \begin{cases} \frac{1}{a} \sin^2\left[\frac{\pi}{2a}(x + a)\right] & x \in [-a, a] \\ 0 & \text{otherwise} \end{cases}\)

P_QW(x) is the probability density for the quantum well model. This density is zero outside [-a, a], eliminating any chance of extreme returns. Inside the well it has a sin² shape that peaks at x = 0 with value 1/a and goes smoothly to zero at the boundaries. The density is symmetric and integrates to 1 over the interval.

Effective volatility

The effective volatility is defined by Eq. (11):

\(\sigma_{\text{eff}} = \sigma \sigma_{\text{QM}}, \quad \sigma_{\text{QM}} = \sqrt{\langle x^2 \rangle_{\text{QM}} - \langle x \rangle_{\text{QM}}^2}\)

σ_eff is the effective volatility that replaces the original σ in the Black‑Scholes formula. σ_QM is the quantum standard deviation computed from the probability density P_QW(x). The expectation values ⟨x⟩_QM and ⟨x²⟩_QM are integrals over the density: ⟨x⟩_QM = ∫ x P_QW(x) dx and ⟨x²⟩_QM = ∫ x² P_QW(x) dx.

For the quantum well, the density is symmetric about zero, so the mean ⟨x⟩_QM is zero. The second moment can be computed exactly by integrating x² times the density over [-a, a]:

\(\langle x^2 \rangle_{\text{QM}} = \int_{-a}^{a} x^2 \cdot \frac{1}{a} \sin^2\!\left[\frac{\pi(x+a)}{2a}\right] dx = \frac{a^2}{3} - \frac{2a^2}{\pi^2}\)

Consequently, the quantum standard deviation is σ_QM = √(a²/3 - 2a²/π²). This analytical formula is exact and avoids numerical integration for the effective volatility. The variance is guaranteed to be non‑negative because 1/3 ≈ 0.3333 and 2/π² ≈ 0.2026, so the difference is positive for any a > 0.

For small a, σ_QM is small, so the effective volatility is low and option prices are depressed. As a grows, σ_QM increases, and the option price rises. However, the option price saturates for large a rather than converging to the Black‑Scholes price. This is because the well always imposes hard bounds: the density is exactly zero outside [-a, a], so the tails of the distribution are always truncated, no matter how wide the well becomes. The cumulative distribution N_eff(d_eff±) therefore never matches the standard normal CDF, and the option price remains below the unbounded Black‑Scholes value.

Modified option pricing with the quantum well

To price an option under the quantum well model, the effective volatility σ_eff is used in place of σ when computing the integration limits d_eff± via the standard formula (Eq. 3, reconstructed):

\(d_{\pm} = \frac{\ln(S_0/K) + (r \pm \sigma^2/2)T}{\sigma\sqrt{T}}\)

With σ replaced by σ_eff, the modified limits become:

\(d_{\text{eff}\pm} = \frac{\ln(S_0/K) + (r \pm \sigma_{\text{eff}}^2/2)T}{\sigma_{\text{eff}}\sqrt{T}}\)

The modified cumulative distribution N_eff(d_eff±) is then obtained by numerically integrating P_QW(x) from -∞ to d_eff±. Because P_QW(x) is zero outside [-a, a], the integration domain is effectively truncated to [-a, min(d_eff±, a)]. The option prices are computed using the standard Black‑Scholes formula structure (Eqs. 1‑2) with σ_eff and N_eff replacing σ and N.

C++ implementation

The header src/quantum_well_model.h provides two core functions: one for the probability density and one for the analytical σ_QM. Both are straightforward translations of the equations above.

inline double quantum_well_sigma_QM(double a) {
    if (a <= 0.0) {
        throw std::invalid_argument("quantum_well_sigma_QM: a must be positive, got " + std::to_string(a));
    }
    const double a2 = a * a;
    const double var = a2 / 3.0 - 2.0 * a2 / (M_PI * M_PI);
    return std::sqrt(var);
}

The function quantum_well_sigma_QM computes the exact standard deviation. Input validation ensures a > 0, throwing std::invalid_argument otherwise. The variance is guaranteed to be non‑negative because 1/3 ≈ 0.3333 and 2/π² ≈ 0.2026, so the difference is positive for any a > 0.

inline double quantum_well_PDF(double x, double a) {
    if (a <= 0.0) {
        throw std::invalid_argument("quantum_well_PDF: a must be positive, got " + std::to_string(a));
    }
    if (x <= -a || x >= a) {
        return 0.0;
    }
    const double argument = M_PI * (x + a) / (2.0 * a);
    const double sin_val = std::sin(argument);
    return (1.0 / a) * sin_val * sin_val;
}

The PDF function first validates a > 0, then checks whether x lies outside the well; if so, it returns zero immediately. Inside the well, it evaluates the sin² expression exactly as given in Eq. (29).

Verification

The header file includes a comprehensive verification function verify_quantum_well_model() that runs several static checks. Here is an excerpt showing the normalization check and the analytical‑vs‑numerical σ_QM check:

// Check 1: Integral of P_QW over [-a, a] ≈ 1.0
for (double a : a_values) {
    const int N = 20001;
    std::vector<double> x(N);
    std::vector<double> f(N);
    const double dx = 2.0 * a / (N - 1);
    for (int i = 0; i < N; ++i) {
        x[i] = -a + i * dx;
        f[i] = quantum_well_PDF(x[i], a);
    }
    const double integral = simpson_integrate(x, f);
    const bool passed = std::abs(integral - 1.0) < tolerance;
    // ... pass/fail reporting ...
}

// Check 2: Analytical σ_QM vs numerical integration
for (double a : a_values) {
    const double sigma_analytical = quantum_well_sigma_QM(a);
    // Numerical computation of σ_QM via Simpson integration of x*f and x²*f
    // ...
    const bool passed = std::abs(sigma_analytical - sigma_numerical) < tolerance;
    // ... pass/fail reporting ...
}

The full verification function performs seven checks:

• Normalization: The integral of P_QW over [-a, a] equals 1.0 for several test values of a.

• Analytical vs. numerical σ_QM: The analytical formula matches a numerical integration of ⟨x²⟩ and ⟨x⟩.

• Zero outside the well: P_QW(x) returns exactly 0 for |x| ≥ a.

• Non‑negativity: The density is never negative inside the well.

• Small‑a limit: As a → 0, σ_QM → 0.

• Peak value: P_QW(0, a) = 1/a.

• Boundary condition: P_QW(±a, a) = 0.

These checks ensure that the implementation faithfully reproduces the paper’s model and that the probability density is a valid distribution. When all checks pass, the model is ready to be used in the parameter sweeps that reproduce Figure 5 of the paper.

Worked example

Take the paper’s standard parameters: S0 = 20, K = 20, r = 0.1, T = 1, σ = 0.25. For a well half‑width a = 2.0:

1. Quantum standard deviation: σ_QM = quantum_well_sigma_QM(2.0) ≈ 0.737.

2. Effective volatility: σ_eff = 0.25 × 0.737 ≈ 0.184.

3. Modified d_eff±: Using σ_eff in the d± formula:

• d_eff+ = [ln(20/20) + (0.1 + 0.184²/2) × 1] / (0.184 × √1) ≈ 0.635

• d_eff- = d_eff+ - 0.184 ≈ 0.451

4. Modified CDF: Integrate P_QW(x) from -∞ to d_eff±. Since d_eff± are within [-2, 2], the integration uses the sin² density. The resulting N_eff(d_eff+) and N_eff(d_eff-) are smaller than their standard normal counterparts because the distribution is narrower.

5. Option price: The call price computed via the modified Black‑Scholes formula is lower than the baseline 3.022.

The code call quantum_well_PDF(x, 2.0) returns the density at any x, and quantum_well_sigma_QM(2.0) returns approximately 0.737.

Reproducing Figure 5

Figure 5 of the paper plots the option price against the well half‑width a. The parameter sweep function sweep_quantum_well in src/paper_figures.h varies a over a range (e.g., 0.1 to 5.0) while keeping S0, K, r, T, and σ fixed at the standard values. For each a, it builds a CDF grid from quantum_well_PDF, computes σ_eff via quantum_well_sigma_QM, and calls modified_option_pricing. The resulting call price increases with a and saturates, never reaching the Black‑Scholes baseline.

Limitations

The quantum well model is a stylized representation of hard price bounds. The abrupt cutoff at x = ±a is unrealistic for most markets, where price limits are rarely absolute. The parameter a is only qualitatively related to the price range ΔS; the paper provides no calibration procedure. Because the density is always zero outside the well, the model never converges to the Black‑Scholes price, even for arbitrarily large a. Despite these limitations, the model illustrates how a simple potential can capture the effect of price boundaries on option valuation, completing the set of five market force models presented in the paper.

Code Walkthrough: File Structure, Key Functions, and Data Flow

Now that we have studied every market force model and the machinery that turns a probability density into an option price, it is time to look at how all the pieces fit together inside the C++ codebase. The implementation is organised as a collection of header‑only libraries, making it easy to include and reuse individual components. This section presents the file dependency hierarchy, walks through the purpose of each module, and traces the exact path a single force parameter takes from input to a pair of call and put prices. We will also see how the sweep functions that reproduce the paper’s figures are built and how the final results are written to CSV files for external plotting.

File Dependency Hierarchy

The project follows a strict, layered build order. Every file depends only on modules that are listed above it, and no circular dependencies exist. The recommended inclusion order is:

1. `src/black_scholes_base.h` – the baseline Black‑Scholes formulas and the standard normal CDF/PDF. No dependencies.

2. `src/effective_volatility.h` – the Simpson integrator and the compute_eff_vol function. The primary functions in this header (simpson_integrate and compute_eff_vol) have no dependency on normalPDF. The verification function verify_effective_volatility defines a local lambda that matches the standard normal PDF for its own checks, so the header does not need to include black_scholes_base.h.

3. `src/modified_option_pricing.h` – the modified Black‑Scholes pricing that uses a precomputed CDF grid and effective volatility. It includes both black_scholes_base.h (for compute_dpm and normalPDF) and effective_volatility.h (for the VolResult struct). Note that build_cdf_grid implements its own cumulative integration internally and does not call simpson_integrate from effective_volatility.h; simpson_integrate is used elsewhere, such as in compute_eff_vol and the normalisation functions of the perturbative models.

4. Market force model headers – each file (constant_force_model.h, linear_force_model.h, x2_force_model.h, x3_force_model.h, quantum_well_model.h) defines one probability density function and its helpers. They do not include the core headers directly. Instead they use forward declarations (e.g., extern double normalPDF(double x); or namespace wrappers) for the symbols they need during verification. The symbols are resolved at link time when all headers are included together in a translation unit such as main.cpp.

5. `src/paper_figures.h` – the sweep functions that reproduce Figures 1‑5. This header depends on all five model headers plus the three core modules. It uses forward declarations to avoid including everything at once, but the user must include the required headers before this file.

6. `src/main.cpp` – the entry point that calls all verification routines and then generates the CSV data for the figures. It includes every header listed above.

If you are compiling the project from scratch, you only need to compile main.cpp; all other code is pulled in through the headers.

Core Modules

black_scholes_base.h – The Baseline

This file provides the four building blocks that the rest of the code depends on:

• normalPDF(double x) – the standard normal density P_BS(x).

• normalCDF(double x) – the standard normal cumulative distribution N(x) implemented via std::erfc.

• compute_dpm(S0, K, r, sigma, T) – returns d_+ and d_- from the reconstructed Eq. (3).

• black_scholes_base(…) – returns a BSResult struct with call, put, and the intermediate d and N(d) values.

A small excerpt shows the interface:

struct BSResult {
    double call, put, d_plus, d_minus, N_d_plus, N_d_minus;
};

BSResult black_scholes_base(double S0, double K, double r, double T, double sigma);

Every subsequent module uses normalPDF and compute_dpm; the CDF normalCDF is used only for verification and as a baseline reference.

effective_volatility.h – The Integrator and σ_QM

This header defines the VolResult struct and the two workhorses of numerical integration:

struct VolResult {
    double sigma_eff, sigma_QM, mean_x, var_x;
};

double simpson_integrate(const std::vector<double>& x, const std::vector<double>& f);
VolResult compute_eff_vol(double sigma, std::function<double(double)> P,
                          double a = -12.0, double b = 12.0, int N = 20001);

simpson_integrate applies Simpson’s 1/3 rule on a uniform grid. It is used by compute_eff_vol and by the normalisation functions of the perturbative models. The function compute_eff_vol builds a grid, evaluates P(x), x·P(x), and x²·P(x), and returns the effective volatility σ_eff = σ · σ_QM.

modified_option_pricing.h – The Modified Black‑Scholes Engine

This header is the glue that combines the new probability density with the Black‑Scholes formula structure. It provides:

struct ModifiedBSResult {
    double call, put, d_eff_plus, d_eff_minus, N_eff_plus, N_eff_minus;
};

std::pair<std::vector<double>, std::vector<double>>
build_cdf_grid(std::function<double(double)> P, double a, double b, int N);

double interpolate_cdf(double d, const std::vector<double>& x_grid,
                       const std::vector<double>& cdf_grid);

ModifiedBSResult modified_option_pricing(double S0, double K, double r, double T,
                                         double sigma_eff,
                                         std::function<double(double)> P,
                                         const std::vector<double>& x_grid,
                                         const std::vector<double>& cdf_grid);

• build_cdf_grid precomputes the CDF N_eff(x) = ∫_{a}^{x} P(t) dt on a uniform grid. It uses a mixed integration scheme: Simpson’s rule for completed panels (even indices) and the trapezoidal rule for the first interval and for half‑panels at odd indices. The first element is 0.0, and the last should be approximately 1.0 for a properly normalised density.

• interpolate_cdf performs linear interpolation between the two nearest grid points to evaluate the CDF at arbitrary d_eff± values. When d falls outside the grid, it clamps the result to 0.0 for d ≤ x_grid[0] and to cdf_grid.back() (which is approximately 1.0 for a normalised density) for d ≥ x_grid[n-1].

• modified_option_pricing calls compute_dpm with sigma_eff to obtain d_eff±, then uses interpolate_cdf to get N_eff(d_eff±), and finally plugs everything into the call and put formulas. The P parameter is not used inside the function body; it is retained in the signature for interface consistency and documentation.

When P is the standard normal PDF and sigma_eff equals the original σ, the function reproduces the baseline Black‑Scholes prices to within 1e-10.

Market Force Model Headers

Each of the five model files follows a uniform pattern:

• One or two helper functions that compute analytical parameters (e.g., constant_force_xk, linear_force_sigma_QM, quantum_well_sigma_QM).

• A function *_PDF(x, param) that returns the value of the probability density at x.

• For the perturbative models (x2_force_model.h and x3_force_model.h), a separate compute_*_normalization function and a *_PDF that takes the precomputed normalisation constant C2.

• A verify_* function that checks normalisation, zero‑force recovery, and (for the x^2 and x^3 models) warns when the density becomes negative.

For example, the constant force model exposes:

double constant_force_xk(double k);
double constant_force_PDF(double x, double k);

The linear force model provides analytical σ_QM:

double linear_force_lambda_omega(double lambda);
double linear_force_sigma_QM(double lambda);
double linear_force_PDF(double x, double lambda);

The perturbative models require a normalisation step:

// x² force
double compute_x2_normalization(double beta, double a = -12.0, double b = 12.0, int N = 20001);
double x2_force_PDF(double x, double beta, double C2);

// x³ force
double compute_x3_normalization(double gamma, double a = -12.0, double b = 12.0, int N = 20001);
double x3_force_PDF(double x, double gamma, double C2);

The quantum well model has an analytical σ_QM and a piecewise PDF:

double quantum_well_sigma_QM(double a);
double quantum_well_PDF(double x, double a);

All these functions are inline and defined in the headers, so no separate compilation is required.

Sweep Functions and Figure Generation

The file src/paper_figures.h contains five sweep functions, one for each figure in the paper. They all share the same signature pattern:

SweepResult sweep_constant_force(double S0, double K, double r, double T, double sigma,
                                 double k_min, double k_max, int n_steps);
SweepResult sweep_linear_force(double S0, double K, double r, double T, double sigma,
                                double lambda_min, double lambda_max, int n_steps);
// ... and similarly for x2, x3, and quantum_well

A SweepResult is simply a struct with three vectors:

struct SweepResult {
    std::vector<double> param_values;
    std::vector<double> call_prices;
    std::vector<double> put_prices;
};

Each sweep function iterates over the requested parameter range, performs the following steps for every parameter value, and records the resulting call and put prices.

Let us trace the sweep_constant_force function as a concrete example. The code, slightly simplified, looks like this:

SweepResult sweep_constant_force(double S0, double K, double r, double T, double sigma,
                                 double k_min, double k_max, int n_steps)
{
    SweepResult result;
    const double a = -12.0, b = 12.0;
    const int N = 20001;
    const double sigma_eff = sigma;   // variance unchanged for constant force

    for (int i = 0; i < n_steps; ++i) {
        double k = k_min + (k_max - k_min) * i / (n_steps - 1);

        auto P = [k](double x) { return constant_force_PDF(x, k); };
        auto [x_grid, cdf_grid] = build_cdf_grid(P, a, b, N);
        ModifiedBSResult mbs = modified_option_pricing(S0, K, r, T, sigma_eff, P, x_grid, cdf_grid);

        result.param_values.push_back(k);
        result.call_prices.push_back(mbs.call);
        result.put_prices.push_back(mbs.put);
    }
    return result;
}

• Step 1 – Create the density: A lambda captures the current k and calls constant_force_PDF(x, k). This lambda is a valid std::function<double(double)>.

• Step 2 – Build the CDF grid: build_cdf_grid evaluates the lambda on N points from -12 to 12 and accumulates the integral using the mixed Simpson/trapezoidal scheme described above. The result is a pair of grids: x_grid and cdf_grid.

• Step 3 – Compute effective volatility: For the constant force model we know analytically that σ_QM = 1, so sigma_eff = sigma. The code does not call compute_eff_vol here; it simply uses sigma directly. The same is true for sweep_linear_force and sweep_quantum_well, which use analytical σ_QM formulas (linear_force_sigma_QM and quantum_well_sigma_QM). Only the perturbative sweeps (sweep_x2_force and sweep_x3_force) call compute_eff_vol because their variance changes in a way that is not captured by a simple analytical expression.

• Step 4 – Price the option: modified_option_pricing takes the original market parameters, sigma_eff, the density lambda (retained in the signature for interface consistency but not used inside the function), and the precomputed grids. It computes d_eff±, interpolates the CDF, and returns the modified call and put prices.

• Step 5 – Store the results: The parameter value and the two prices are appended to the SweepResult vectors.

The other sweep functions follow the same logic, with the only differences being the density function and the method for obtaining sigma_eff.

The Entry Point: main.cpp

The file main.cpp orchestrates the entire workflow. It defines two helper functions—run_verification and generate_all_figures—and a main that calls them in sequence. The programme uses the paper’s sample parameters (S0=20, K=20, r=0.1, T=1, σ=0.25) throughout.

int main() {
    try {
        run_verification();
        generate_all_figures();
        return 0;
    } catch (const std::exception& e) {
        std::cerr << "Error: " << e.what() << "\n";
        return 1;
    }
}

run_verification calls every verify_* function from the corresponding headers. If any check fails, it sets a flag and, after all checks have run, throws std::runtime_error("Verification failure") to stop the programme before any figure data is written. This ensures that users are immediately alerted to inconsistencies.

generate_all_figures invokes the five sweep functions with parameter ranges chosen to reproduce the paper’s Figures 1‑5. For example, the constant force sweep is called as:

SweepResult res = sweep_constant_force(20.0, 20.0, 0.1, 1.0, 0.25, -1.0, 1.0, 50);
write_sweep_csv("fig1_constant_force.csv", res, "k");

Each sweep result is written to a CSV file by a helper function write_sweep_csv. The output format is straightforward:

k,call_price,put_price
-1.0000000000,3.022...,1.119...
-0.9591836735,3.020...,1.117...
...

The first column is the force parameter, the second is the modified call price, and the third is the modified put price. The values are printed with std::setprecision(10), which sets the number of significant digits (not decimal places). For values around 3.0 this gives about nine decimal places; for values around 0.1 it gives about eleven. These CSV files can be loaded into any plotting tool (Python’s matplotlib, Excel, gnuplot) to generate the figures.

Data Flow for a Single Sweep

To summarise the data flow, let us follow a single parameter, say k = 0.5, through the constant force sweep:

1. constant_force_PDF is called with x values from -12 to 12 and k = 0.5.

2. build_cdf_grid accumulates the integral of P_k(x) and returns x_grid and cdf_grid.

3. modified_option_pricing receives sigma_eff = 0.25, the density lambda, and the grids.

4. Inside modified_option_pricing, compute_dpm(20, 20, 0.1, 0.25, 1.0) returns d_eff+ ≈ 0.525 and d_eff- ≈ 0.275.

5. interpolate_cdf evaluates the CDF at d_eff+ and d_eff- using the precomputed grid.

6. The call price is computed as 20 * N_eff(0.525) - 20 * exp(-0.1) * N_eff(0.275).

7. The result is appended to SweepResult.call_prices.

This identical pattern is used for every force model, every parameter, and every point in the figures. The separation of concerns—density definition, CDF construction, effective volatility computation, and Black‑Scholes evaluation—makes the code easy to extend. Adding a new force model requires only writing a new *_PDF function and a corresponding sweep function, without touching the core pricing engine.

How to Plot the Results

After running the programme, you will have five CSV files. A minimal Python snippet to plot Figure 1 looks like:

import pandas as pd
import matplotlib.pyplot as plt
data = pd.read_csv('fig1_constant_force.csv')
plt.plot(data['k'], data['call_price'], label='Call')
plt.plot(data['k'], data['put_price'], label='Put')
plt.xlabel('k')
plt.ylabel('Option Price')
plt.legend()
plt.show()

No built‑in plotting library is used in the C++ code; the separation of computation and visualisation keeps the implementation lightweight and portable.

With this walkthrough, you now have a complete picture of how the code is structured and how the mathematical models from the paper are translated into a working C++ implementation. The next section will verify that every component behaves correctly and that the paper’s figures are faithfully reproduced.

Reproducing the Paper's Figures: Parameter Sweeps and Plotting

The previous sections built the core machinery that turns a probability density and an effective volatility into a modified option price. We now use that machinery to recreate the five qualitative graphs shown in the paper. Each graph plots the call option price against the strength parameter of one market force, while all standard Black‑Scholes inputs are held fixed.

The paper itself provides only illustrative sketches and does not give the exact parameter ranges. The sweep functions in src/paper_figures.h make the ranges explicit so you can reproduce the described behaviour and experiment with your own intervals. This section explains the fixed financial parameters, introduces the SweepResult data structure and the numerical integration domain, walks through each sweep function, and shows how to export the resulting data to CSV files for plotting with external tools.

Fixed Financial Parameters

All five sweeps share the same underlying asset and option parameters. They are chosen so the baseline Black‑Scholes call price is approximately 3.022 and the put price approximately 1.119:

• Current asset price S0 = 20

• Strike price K = 20 (at‑the‑money)

• Risk‑free rate r = 0.1 (10 % per annum)

• Time to maturity T = 1 year

• Volatility σ = 0.25 (25 % per annum)

A crucial consistency requirement is that when the force parameter is zero (or its equivalent neutral value), the modified option price must recover the baseline Black‑Scholes price. Every sweep function is designed to honour this property, and the built‑in verification routine verify_paper_figures formally checks it.

The SweepResult Structure

All sweep functions return a SweepResult struct, defined in src/paper_figures.h. It holds three std::vector<double> members of equal length:

• param_values – the force parameter at each step (k, λ, β, γ, or a),

• call_prices – the corresponding modified call option price,

• put_prices – the corresponding modified put option price.

The uniform structure lets you loop over a sweep and write rows to a CSV file without knowing which force model produced the data.

struct SweepResult {
    std::vector<double> param_values;
    std::vector<double> call_prices;
    std::vector<double> put_prices;
};

Every sweep function accepts the five fixed Black‑Scholes parameters (S0, K, r, T, σ), the minimum and maximum of the force parameter, and an integer n_steps (≥ 2) that controls the number of equally spaced points. The first entry always equals param_min and the last equals param_max.

Numerical Integration Domain

For unbounded densities the sweep functions rely on the same numerical integration domain used throughout the codebase. The constants are gathered in the SweepConfig namespace inside src/paper_figures.h:

• SweepConfig::INTEGRAL_A = -12.0 – lower bound

• SweepConfig::INTEGRAL_B = 12.0 – upper bound

• SweepConfig::INTEGRAL_N = 20001 – number of Simpson grid points (must be odd)

These values are passed to build_cdf_grid and compute_eff_vol for every sweep model except the quantum well, which uses its own bounded support [-a, a] with the same number of points. The domain [-12, 12] is wide enough to capture the tails of the standard normal distribution and all the perturbed densities used in the paper.

Constant Market Force (Figure 1)

The constant force model uses a shifted Gaussian whose variance is unchanged; therefore the standard deviation of the quantum density is one, σ_eff equals the original σ, and no numerical integration for the effective volatility is needed. The sweep function sweep_constant_force iterates over k, builds a probability density through constant_force_PDF(x, k), constructs a CDF grid with build_cdf_grid, and calls modified_option_pricing.

A call that reproduces Figure 1 is:

SweepResult fig1 = sweep_constant_force(
    20.0, 20.0, 0.1, 1.0, 0.25,   // S0, K, r, T, sigma
    -2.0, 2.0, 50                   // k_min, k_max, n_steps
);

The loop body inside the sweep function is shown below. Notice that sigma_eff is assigned the original sigma before the loop because the constant force only shifts the mean, leaving the width intact.

const double sigma_eff = sigma;   // constant force does not change variance
for (int i = 0; i < n_steps; ++i) {
    double k = k_min + (k_max - k_min) * static_cast<double>(i) / (n_steps - 1);

    auto P = [k](double x) { return constant_force_PDF(x, k); };
    auto [x_grid, cdf_grid] = build_cdf_grid(P, a, b, N);

    ModifiedBSResult mbs = modified_option_pricing(
        S0, K, r, T, sigma_eff, P, x_grid, cdf_grid);

    result.param_values.push_back(k);
    result.call_prices.push_back(mbs.call);
    result.put_prices.push_back(mbs.put);
}

The call price decreases monotonically as k increases. A positive k gives a force F = −k that pushes the probability mass to the left, away from the strike, thereby reducing the probability that the option finishes in‑the‑money. For negative k the distribution shifts to the right, increasing the call price. When k = 0 the density is exactly the standard normal and the sweep recovers the baseline Black‑Scholes price.

Linear Market Force (Figure 2)

The linear force modifies the variance of the Gaussian. The quantum standard deviation is available analytically as σ_QM = 1 / sqrt(1 + λ), provided by the helper linear_force_sigma_QM. The sweep function sweep_linear_force therefore bypasses numerical integration for effective volatility and computes sigma_eff = σ · σ_QM directly.

A call that reproduces Figure 2 might be:

SweepResult fig2 = sweep_linear_force(
    20.0, 20.0, 0.1, 1.0, 0.25,   // S0, K, r, T, sigma
    0.1, 2.0, 50                    // lambda_min, lambda_max, n_steps
);

Larger λ makes the bowl steeper, narrowing the distribution and reducing σ_eff. Because the probability of large price movements is lower, the at‑the‑money call option becomes less valuable. The call price therefore decreases monotonically with λ. At λ = 0 the density reduces to the standard normal, and the sweep returns the baseline Black‑Scholes price.

x² Market Force (Figure 3)

The x² force model uses the paper's perturbative density. As discussed earlier, that density contains a spurious linear term in β. The sweep function sweep_x2_force nevertheless implements the paper's formula so you can reproduce the original figure. For each β it recomputes the normalisation constant C² through compute_x2_normalization and obtains the effective volatility numerically from the normalised density via compute_eff_vol.

A typical call is:

SweepResult fig3 = sweep_x2_force(
    20.0, 20.0, 0.1, 1.0, 0.25,   // S0, K, r, T, sigma
    -0.5, 0.5, 50                   // beta_min, beta_max, n_steps
);

The call price decreases as β increases. This curve should be interpreted cautiously because of the perturbation error in the paper. The provided code does not contain a corrected second‑order density; if you wish to replace the paper's formula with the correct expansion you need to modify the source files directly (see the discussion of the perturbation error in the earlier section on the x² force model). At β = 0 the density is exactly the standard normal and the sweep recovers the baseline price.

x³ Market Force (Figure 4)

The x³ force model uses first‑order perturbation theory for a quartic potential. The sweep function sweep_x3_force is called with a range of γ values. The paper reports a non‑monotonic behaviour: the option price first decreases and then increases. This is reproduced with a call such as:

SweepResult fig4 = sweep_x3_force(
    20.0, 20.0, 0.1, 1.0, 0.25,   // S0, K, r, T, sigma
    0.0, 0.5, 50                    // gamma_min, gamma_max, n_steps
);

For small γ the quartic term steepens the potential slightly, narrowing the central peak and reducing the call price. As γ grows larger the potential develops a double‑well shape, which can create a bimodal density. The increased probability mass in the tails eventually raises the probability of finishing in‑the‑money, causing the call price to rise again. The density can become negative for large γ, signalling the breakdown of the perturbation expansion. The sweep function still produces numbers, but you should restrict γ to a range where the density stays positive. When γ = 0 the density is the standard normal and the sweep returns the baseline Black‑Scholes price.

Quantum Well (Figure 5)

The quantum well model confines the probability mass to the interval [-a, a]. The effective volatility is obtained analytically from quantum_well_sigma_QM(a), which returns √(a²/3 − 2a²/π²). The sweep function sweep_quantum_well builds the CDF grid only over the support [-a, a] to save computation.

To reproduce Figure 5, call:

SweepResult fig5 = sweep_quantum_well(
    20.0, 20.0, 0.1, 1.0, 0.25,   // S0, K, r, T, sigma
    0.5, 5.0, 50                    // a_min, a_max, n_steps
);

The call price increases with a and saturates for large a. A narrow well (small a) traps the probability mass near zero, making large movements impossible and drastically reducing the call option value. As the well widens, the distribution can spread further, and the option price rises. However, even for a very wide well the price never reaches the baseline Black‑Scholes value because the hard boundaries always eliminate the extreme returns that a log‑normal distribution would allow. The saturation level is therefore below the Black‑Scholes price.

Exporting to CSV and Plotting

The sweep functions return the data in memory. To plot the figures you can write a short loop that outputs rows of param_values[i], call_prices[i], put_prices[i] to a comma‑separated file. The src/main.cpp entry point follows this pattern—it calls all five sweep functions and writes the results to separate CSV files. A minimal helper you could write for your own project might look like this:

#include <fstream>
#include "paper_figures.h"

void write_sweep_to_csv(const SweepResult& sweep, const std::string& filename) {
    std::ofstream out(filename);
    out << "param,call,put\n";
    for (size_t i = 0; i < sweep.param_values.size(); ++i) {
        out << sweep.param_values[i] << ","
            << sweep.call_prices[i]  << ","
            << sweep.put_prices[i]   << "\n";
    }
}

This helper is not part of the supplied code files; it is meant as an illustration of how to convert a SweepResult into CSV format. Once the CSV files exist, any plotting tool (Python with matplotlib, gnuplot, Excel, etc.) can create the graphs. The horizontal axis is the force parameter and the vertical axis is the call price (or put price). Drawing a horizontal line at the baseline Black‑Scholes price makes the effect of the market force visually obvious.

Built‑in Verification

The header src/paper_figures.h also provides a function verify_paper_figures that runs static integrity checks on the sweep machinery. It confirms that for zero force parameters the modified prices match the baseline and that call prices stay within the theoretical bounds [0, S0]. Call it before generating the figures to catch configuration errors early:

if (!verify_paper_figures()) {
    std::cerr << "Paper figure verification failed.\n";
    return 1;
}

With the sweep functions described above, you have a complete toolset for exploring how different market forces alter standard Black‑Scholes prices. The next section covers the full verification strategy in more detail, and the final section summarises the key insights and limitations of the quantum approach.

Verification: Static Checks, Invariants, and Paper-Fidelity Review

Before we trust the option prices produced by our quantum market force models, we need to be confident that the code is correct. The implementation includes a comprehensive set of built‑in verification functions that test mathematical invariants, check consistency with the paper's equations, and ensure that edge cases behave as expected. This section explains the verification strategy, walks through the key invariants that every model must satisfy, and shows how the C++ code in src/main.cpp orchestrates all the checks before any figure data is generated.

Important: The verification described here is static and semantic only. No code execution is performed in this tutorial. The checks are designed to be compiled and run as part of a normal build process, but the tutorial itself does not execute them. The pass/fail messages shown are illustrative of what the code would produce when compiled and executed.

The Verification Strategy

The verification approach has three layers. First, the code is designed to compile cleanly with strict warnings enabled (-Wall -Wextra -Werror), catching potential bugs at the compiler level. Second, every header file contains a dedicated verify_* function that runs a battery of mathematical checks and prints pass/fail messages to the console. Third, the entry point main() calls all verification functions before generating any output; if any check fails, the program aborts with an error message rather than producing potentially incorrect results.

The checks are static and semantic—they test the mathematical properties of the code but do not require execution against external market data. The paper does not provide validation against real option prices, so the verification focuses on internal consistency: do the probability densities integrate to one? Does the modified pricing recover the baseline Black‑Scholes when the force is turned off? Does put‑call parity hold for every model?

The Orchestration in main.cpp

The run_verification function in src/main.cpp calls every model's verification routine and collects the results. Here is the relevant excerpt:

void run_verification() {
    std::cout << "=== Running verification checks ===\n";
    bool all_passed = true;

    auto check = [&](const std::string& name, bool result) {
        std::cout << "  " << name << ": " << (result ? "PASSED" : "FAILED") << "\n";
        if (!result) all_passed = false;
    };

    check("Black-Scholes baseline", verify_black_scholes_base());
    check("Effective volatility", verify_effective_volatility());
    check("Modified option pricing", verify_modified_option_pricing());
    check("Constant force model", verify_constant_force_model());
    check("Linear force model", verify_linear_force_model());
    check("x^2 force model", verify_x2_force_model());
    check("x^3 force model", verify_x3_force_model());
    check("Quantum well model", verify_quantum_well_model());
    check("Paper figures (sweep sanity)", verify_paper_figures());

    if (!all_passed) {
        std::cerr << "\nSome verification checks FAILED. Aborting figure generation.\n";
        throw std::runtime_error("Verification failure");
    }
    std::cout << "All verification checks passed.\n\n";
}

Notice the order: the baseline Black‑Scholes module is verified first because all other modules depend on it. The effective volatility and modified pricing modules come next, followed by each force model, and finally the sweep functions that reproduce the paper's figures. If any single check fails, the program throws an exception and stops. This guarantees that the CSV files written later are produced by code that has passed every internal consistency test.

Invariant 1: Probability Density Normalization

Every probability density P(x) used in the computation must integrate to one over the real line. This is a fundamental requirement: if the density is not normalized, the cumulative distribution N_eff(d) will not be a valid CDF, and the option prices will be wrong.

Each model's verification function builds a uniform grid over the integration domain (typically [-12, 12] with 20001 points) and uses Simpson's rule to compute the integral of P(x). The result is compared to 1.0 with a tolerance of 1e-10. For example, the constant force model verification in src/constant_force_model.h contains a loop that sweeps over several force strengths k and verifies that the shifted Gaussian integrates to one regardless of the shift. The same pattern appears in every model's verification function. For the perturbation‑based models (x² and x³ forces), the normalization constant C² is computed numerically, and the check confirms that the normalized density indeed integrates to one.

Invariant 2: Put‑Call Parity

Put‑call parity is a model‑independent relationship that must hold for any European option pricing formula that uses a valid cumulative distribution function. The identity is

\(c + K e^{-rT} = p + S_0\)

where c is the call price, p is the put price, K is the strike, r is the risk‑free rate, T is the time to maturity, and S_0 is the current asset price. This parity relies only on the definition of the CDF, not on the specific shape of the probability density. Therefore, it must hold for the baseline Black‑Scholes model and for every modified model we implement. Note that this identity is a standard financial invariant and is not explicitly numbered as an equation in the paper; it follows directly from the definitions of the call and put prices in Eqs. (1) and (2).

The baseline verification in src/black_scholes_base.h checks this explicitly:

double parity_lhs = res.call + K * std::exp(-r * T);
double parity_rhs = res.put + S0;
bool pass = std::abs(parity_lhs - parity_rhs) < eps;

The modified option pricing verification in src/modified_option_pricing.h repeats the same check for the modified prices, ensuring that the new CDF N_eff does not break the parity. If a market force model produced a density that was not a valid probability distribution, the parity check would catch it.

Invariant 3: Baseline Recovery When Force Parameters Are Zero

When the force strength parameter is zero, every model must reduce exactly to the standard Black‑Scholes baseline. This is the most important paper‑fidelity check: the quantum approach is an extension, not a replacement, and it must nest the original model.

Each model's verification function tests this at multiple points. For the constant force model, k = 0 must give P_k(x) = normalPDF(x). The verification in src/constant_force_model.h checks this with a loop over several test points:

const double k = 0.0;
const double test_points[] = {-3.0, -1.0, 0.0, 1.0, 3.0};
for (double x : test_points) {
    double pk = constant_force_PDF(x, k);
    double np = normalPDF(x);
    if (std::abs(pk - np) > tol) {
        all_passed = false;
    }
}

For the linear force model, the limit λ → 0 is tested with a very small positive λ (since the function throws for λ ≤ 0). The x² and x³ force models check that β = 0 and γ = 0 recover the standard normal PDF after normalization. The quantum well model does not have a zero‑force limit in the same sense, but its verification checks that the analytical σ_QM formula matches numerical integration.

The sweep verification in src/paper_figures.h takes this further: it runs a full pricing sweep with zero force and compares every call and put price to the baseline Black‑Scholes result, ensuring that the entire pipeline—density, effective volatility, CDF grid, and modified pricing—collapses to the original model.

Invariant 4: Effective Volatility Consistency

For the baseline standard normal PDF, the quantum standard deviation σ_QM must be exactly 1.0, so the effective volatility σ_eff must equal the input σ. The verification in src/effective_volatility.h checks this:

VolResult res = compute_eff_vol(0.25, normalPDF, -12.0, 12.0, 20001);
if (std::abs(res.sigma_QM - 1.0) > tol) { /* fail */ }
if (std::abs(res.sigma_eff - 0.25) > tol) { /* fail */ }

For the constant force model, the variance is unchanged by the shift, so σ_QM must remain 1.0 for any k. The verification confirms this numerically. For the linear force model, the analytical formula σ_QM = 1 / sqrt(1 + λ) is tested against numerical integration, and the identity σ_QM = 1 / λ_ω is checked exactly. For the quantum well, the analytical formula σ_QM = sqrt(a²/3 - 2a²/π²) is verified against Simpson integration over the well domain.

Invariant 5: CDF Monotonicity and Range

A valid cumulative distribution function must be monotonically non‑decreasing and must return values in [0, 1]. The modified option pricing verification in src/modified_option_pricing.h checks both properties for the CDF grid built from normalPDF:

bool monotonic = true;
for (size_t i = 1; i < cdf_grid.size(); ++i) {
    if (cdf_grid[i] < cdf_grid[i - 1] - 1e-15) {
        monotonic = false;
        break;
    }
}

It also verifies that the CDF grid ends at approximately 1.0 and that the interpolated values N_eff(d_eff±) lie in [0, 1]. These checks ensure that the numerical integration and interpolation routines are working correctly before they are used with the more exotic force‑model densities.

Invariant 6: Simpson Integration Accuracy

The Simpson integrator is the workhorse of the entire codebase. It is used for normalization, expectation values, and CDF construction. The effective volatility verification tests it against a known integral: the area under x² on [0, 1] must equal 1/3.

const int N = 101;
const double a = 0.0, b = 1.0;
// ... build x and f vectors with f[i] = x[i] * x[i] ...
const double integral = simpson_integrate(x, f);
if (std::abs(integral - 1.0 / 3.0) > tol) { /* fail */ }

This test is independent of any probability density and confirms that the Simpson rule implementation is correct to machine precision. The paper relies on numerical integration throughout (for expectation values in Eq. (11) and for the modified CDF), so the accuracy of this integrator is critical to the fidelity of the entire implementation.

Model‑Specific Warnings

The verification functions for the perturbation‑based models (x² and x³ forces) include an additional check that is not a hard failure but a warning. Because the perturbation expansions are truncated, the resulting probability density can become negative for large force strengths or large |x|. The verification scans the entire integration grid and reports the minimum value:

if (min_val < 0.0) {
    std::cout << "WARNING: x3_force_model - For gamma=" << gamma
              << ", P_gamma becomes negative (min=" << min_val
              << " at x=" << min_x
              << "). Perturbation expansion may be invalid for this gamma."
              << std::endl;
}

This warning alerts the user that the perturbation approximation has broken down and that the resulting option prices should not be trusted. It does not cause the verification to fail, because the paper's formula is known to be approximate, but it provides important guidance for choosing parameter ranges.

What the Verification Does Not Check

The verification suite is thorough but has deliberate boundaries. It does not check that the option prices match any market data, because the paper provides no such validation. It does not execute the code and compare output to expected numerical values from the paper, because the paper's figures are qualitative sketches without precise data points. It does not verify that the corrected x² force perturbation (enabled by USE_CORRECTED_X2) is mathematically correct, because that formula is derived externally and not part of the paper.

The verification is a paper‑fidelity review: it ensures that the implementation faithfully reproduces the equations and methods described in the paper, that the mathematical invariants are preserved, and that the code is internally consistent. When all checks pass, you can be confident that the CSV files produced by main() correctly implement the quantum option pricing framework as the paper intended.

Conclusion: Insights, Limitations, and Extensions

We have travelled from the familiar Black–Scholes formula to a family of option pricing models inspired by quantum mechanics. The journey began with a simple observation: the standard normal probability density that underpins Black–Scholes is exactly the squared ground‑state wavefunction of a quantum harmonic oscillator. From that starting point, the paper builds a framework in which hidden “market forces” are represented by additional potential energy terms. Each potential distorts the harmonic oscillator bowl, producing a new probability density that can capture skewness, fat tails, or hard price bounds—features that the log‑normal model cannot express.

What the Framework Achieves

The core contribution of the paper is a phenomenological recipe for generating non‑Gaussian distributions within a Black–Scholes‑like pricing structure. The recipe has three steps:

1. Choose a potential V(x) that encodes the market force you want to model.

2. Solve (or approximate) the stationary Schrödinger equation to obtain the ground‑state probability density P(x).

3. Compute the effective volatility σ_eff = σ · σ_QM from the standard deviation of P(x), then price options by integrating P(x) up to the modified limits d_eff±.

This procedure yields a family of models that are all anchored to the same baseline: when the extra potential vanishes, every model collapses back to the standard Black–Scholes result. The five forces studied in the paper—constant, linear, x², x³, and the quantum well—demonstrate how different potential shapes produce qualitatively different option price behaviours. A constant force shifts the distribution without changing its width. A linear force squeezes or stretches the variance. The x³ force can create a double‑humped density, leading to non‑monotonic price curves. The quantum well imposes hard boundaries, eliminating the possibility of extreme returns altogether.

The C++ implementation we have built mirrors this structure exactly. Each force model lives in its own header file, depends only on the core pricing and integration modules, and exposes a probability density function that plugs directly into the modified pricing pipeline. The sweep functions in src/paper_figures.h automate the parameter studies that reproduce the paper's qualitative figures, and the verification functions in every module check that mathematical invariants—normalisation, put‑call parity, zero‑force recovery—hold to high precision.

Limitations and Caveats

Despite its elegance, the quantum mechanics approach has important limitations that anyone using the code should keep in mind.

It is a phenomenological model, not a first‑principles derivation. The paper does not derive the potentials from market micro‑structure, agent behaviour, or any economic theory. The potentials are chosen because they are mathematically tractable and produce interesting probability densities. This means the force parameters k, λ, β, γ, and a do not have direct financial interpretations that can be read off a balance sheet or a ticker tape. They are dials that change the shape of the distribution, and their values would need to be calibrated to market data before the model could be used for real pricing.

There is no dynamics. The paper works entirely with stationary states of the Schrödinger equation. It does not propose a stochastic process for the underlying asset price, nor does it derive a time‑dependent evolution that could be used for path‑dependent options or hedging calculations. The mapping between the dimensionless quantum coordinate x and the financial log‑return is implicit; the paper never writes down a stochastic differential equation that corresponds to a given potential. The effective volatility σ_eff is a heuristic that preserves the algebraic form of the Black‑Scholes formula, but it is not derived from a dynamic hedging argument.

The perturbation treatment for the `x²` force is incorrect. As we discussed in detail in the section on Model 3, the paper's first‑order perturbation result for the cubic potential contains a linear term in β that should vanish by symmetry. The correct leading correction is of order β². The implementation provides the paper's formula for reproducibility, but it also includes a compile‑time flag USE_CORRECTED_X2 that switches to a proper second‑order expansion. Anyone using the x² force model for quantitative work should enable the corrected version or derive their own.

The constant force density was reconstructed. The original Eq. (15) in the paper is garbled. We reconstructed it as a shifted Gaussian with variance 1, which is the physically correct result for a linear potential added to the harmonic oscillator. This reconstruction is consistent with the paper's qualitative description and with standard quantum mechanics, but it is not a verbatim reproduction of the published text.

No market validation is provided. The paper's numerical examples use a single set of synthetic parameters (S0=20, K=20, r=10%, T=1 year, σ=25%). There is no comparison with real option prices, no calibration to implied volatility surfaces, and no statistical test of the model's predictive power. The figures are qualitative illustrations of how option prices change with the force parameters, not evidence that the model describes actual markets.

Numerical integration introduces small errors. The implementation uses Simpson's rule on a truncated domain [-12, 12] with 20001 points. For Gaussian‑like densities this captures more than 99.9% of the probability mass, but for densities with heavy tails or for very large force parameters, the truncation may miss some probability. The quantum well model avoids this issue by integrating only over its compact support [-a, a]. The normalisation constants for the perturbative densities are computed numerically, so they inherit the integration error.

Extensions and Future Directions

The codebase is designed to be a starting point for further exploration. Here are several directions that a motivated reader could pursue.

Calibration to market data. The most natural next step is to treat the force parameters as free parameters and calibrate them to observed option prices. For example, one could minimise the sum of squared differences between model prices and market prices across a range of strikes and maturities. The constant force model (which shifts the mean) could capture skewness in the implied volatility smile, while the linear force model (which changes the variance) could capture the overall level of implied volatility. The x³ force model, with its double‑well potential, might fit volatility smiles that have a pronounced curvature.

Exotic and path‑dependent options. Because the paper only provides stationary probability densities, the current framework cannot price barrier options, Asian options, or any contract whose payoff depends on the path of the underlying asset. Extending the model to a full stochastic process would require solving the time‑dependent Schrödinger equation or finding a classical stochastic process whose stationary distribution matches the quantum probability density. This is a non‑trivial problem, but it would greatly expand the scope of the approach.

Time‑dependent potentials. The paper assumes that the market force potential is static. A natural generalisation is to allow the potential to change over time, perhaps to model changing market conditions or approaching news events. In quantum mechanics, time‑dependent potentials lead to transitions between energy levels; in the financial analogy, this might correspond to regime shifts or stochastic volatility.

Other potentials. The five forces studied in the paper are only a small sample of the possible potentials one could add to the harmonic oscillator. A periodic potential (like a cosine) could model seasonality or support and resistance levels. A potential with a narrow dip and wide flat regions could model a market that is usually quiet but occasionally experiences large moves. The perturbation theory framework used for the x² and x³ forces can be applied to any small perturbation, and for larger perturbations one could use numerical methods to solve the Schrödinger equation directly.

Higher‑order perturbation theory. The x² force model highlights the danger of truncating a perturbation expansion too early. Implementing higher‑order corrections—either analytically with computer algebra or numerically by diagonalising the Hamiltonian in a truncated basis—would improve the accuracy of the perturbative densities and extend their range of validity to larger force parameters.

Alternative numerical methods. The current implementation relies on Simpson integration on a uniform grid. For densities with sharp features or heavy tails, adaptive quadrature or Gauss‑Hermite integration could provide better accuracy with fewer function evaluations. The modular design of the code makes it straightforward to swap in a different integrator.

Final Thoughts

The quantum mechanics approach to option pricing is a creative and thought‑provoking framework. It does not replace the Black‑Scholes model, nor does it claim to be a complete theory of financial markets. Instead, it offers a new language for thinking about non‑Gaussian distributions and a systematic way to generate them. The analogy between a market force and a potential energy term is intuitive: just as a physicist can predict how a particle's position distribution changes when you reshape the potential, a quant can explore how option prices change when you reshape the underlying return distribution.

The C++ implementation accompanying this tutorial makes the framework concrete. Every equation from the paper is translated into working code, every model is verified against mathematical invariants, and every figure can be reproduced with a single function call. Whether you use the code as a learning tool, a starting point for your own research, or a source of ideas for a production pricing system, we hope it deepens your understanding of both option pricing and the unexpected connections between finance and quantum physics.

Use the button below to download code: