Stochastic Computing & Probabilistic Derivatives

Stochastic computing extends numerical differentiation to probabilistic settings, enabling uncertainty quantification and risk analysis. By applying stochastic link functions and calculus corrections, derivatives can reflect the underlying probability distributions of the data, making this essential for financial modeling, risk management, and scientific applications with inherent uncertainty.

🎲 Core Concepts

Stochastic Link Functions transform deterministic derivatives to account for the probabilistic nature of the underlying data. Instead of computing derivatives of f(x), we compute derivatives that reflect the expected behavior when data follows a specific probability distribution.

Itô’s Lemma provides correction terms for stochastic calculus when the underlying process follows a stochastic differential equation. For a function f(X_t) where X_t follows a diffusion process:

\[df = f'(X_t)dX_t + \frac{1}{2}f''(X_t)(dX_t)^2\]

Stratonovich Integral offers an alternative interpretation of stochastic integration that preserves the classical chain rule, useful when the noise is “physical” rather than mathematical.

Financial Applications: Essential for option pricing, risk management, and volatility modeling where asset prices follow geometric Brownian motion or other stochastic processes.

💰 Example 1: Stock Price Derivatives (Geometric Brownian Motion)

Classic Application: Computing option Greeks with stochastic corrections.

import numpy as np
import matplotlib.pyplot as plt
from pydelt.interpolation import SplineInterpolator

# Simulate geometric Brownian motion stock price
# dS_t = μS_t dt + σS_t dW_t
np.random.seed(42)
T = 1.0          # 1 year
N = 252          # Daily observations
dt = T / N
mu = 0.05        # Expected return (5%)
sigma = 0.2      # Volatility (20%)
S0 = 100         # Initial stock price

# Generate price path
t = np.linspace(0, T, N+1)
W = np.random.randn(N+1).cumsum() * np.sqrt(dt)  # Brownian motion
S = S0 * np.exp((mu - 0.5*sigma**2)*t + sigma*W)  # GBM solution

# Fit interpolator
spline = SplineInterpolator(smoothing=0.01)
spline.fit(t, S)

# Compare regular vs stochastic derivatives
regular_deriv_func = spline.differentiate(order=1)
regular_derivatives = regular_deriv_func(t)

# Set log-normal stochastic link (appropriate for stock prices)
spline.set_stochastic_link('lognormal', sigma=sigma, method='ito')
stochastic_deriv_func = spline.differentiate(order=1)
stochastic_derivatives = stochastic_deriv_func(t)

# Also try Stratonovich method
spline.set_stochastic_link('lognormal', sigma=sigma, method='stratonovich')
stratonovich_deriv_func = spline.differentiate(order=1)
stratonovich_derivatives = stratonovich_deriv_func(t)

# Analysis
print("Stock Price Derivative Analysis:")
print(f"Regular derivative mean: {np.mean(regular_derivatives):.2f}")
print(f"Itô stochastic derivative mean: {np.mean(stochastic_derivatives):.2f}")
print(f"Stratonovich derivative mean: {np.mean(stratonovich_derivatives):.2f}")

# Theoretical expectation: E[dS/dt] = μS for GBM
theoretical_mean = mu * np.mean(S)
print(f"Theoretical mean (μS): {theoretical_mean:.2f}")

# Compute differences
ito_correction = np.mean(stochastic_derivatives - regular_derivatives)
stratonovich_correction = np.mean(stratonovich_derivatives - regular_derivatives)

print(f"\nStochastic Corrections:")
print(f"Itô correction: {ito_correction:.2f}")
print(f"Stratonovich correction: {stratonovich_correction:.2f}")

# Option Greeks approximation
# Delta (price sensitivity) ≈ derivative w.r.t. underlying
current_price = S[-1]
current_delta_regular = regular_derivatives[-1] / current_price
current_delta_stochastic = stochastic_derivatives[-1] / current_price

print(f"\nOption Greeks Approximation:")
print(f"Regular Delta: {current_delta_regular:.4f}")
print(f"Stochastic Delta: {current_delta_stochastic:.4f}")

🏦 Example 2: Interest Rate Modeling

Application: Modeling interest rate derivatives with mean reversion.

# Vasicek interest rate model simulation
# dr_t = κ(θ - r_t)dt + σ dW_t
# κ: mean reversion speed, θ: long-term mean, σ: volatility

def vasicek_simulation(r0, kappa, theta, sigma, T, N):
    """Simulate Vasicek interest rate model"""
    dt = T / N
    t = np.linspace(0, T, N+1)
    r = np.zeros(N+1)
    r[0] = r0

    for i in range(N):
        dW = np.random.randn() * np.sqrt(dt)
        r[i+1] = r[i] + kappa*(theta - r[i])*dt + sigma*dW

    return t, r

# Model parameters
r0 = 0.03      # Initial rate (3%)
kappa = 0.5    # Mean reversion speed
theta = 0.04   # Long-term mean (4%)
sigma = 0.01   # Volatility (1%)
T = 5.0        # 5 years
N = 1000       # Time steps

# Simulate interest rate path
np.random.seed(123)
t, r = vasicek_simulation(r0, kappa, theta, sigma, T, N)

# Fit with different interpolators
spline_rates = SplineInterpolator(smoothing=0.001)
spline_rates.fit(t, r)

# Normal distribution appropriate for interest rates (can be negative)
spline_rates.set_stochastic_link('normal', sigma=sigma, method='ito')

# Compute rate derivatives (duration-like measures)
rate_deriv_func = spline_rates.differentiate(order=1)
rate_derivatives = rate_deriv_func(t)

# Second derivatives (convexity-like measures)
rate_second_deriv_func = spline_rates.differentiate(order=2)
rate_second_derivatives = rate_second_deriv_func(t)

# Theoretical drift: E[dr/dt] = κ(θ - r)
theoretical_drift = kappa * (theta - r)

print("Interest Rate Analysis:")
print(f"Mean rate: {np.mean(r):.4f}")
print(f"Rate volatility: {np.std(r):.4f}")
print(f"Mean derivative: {np.mean(rate_derivatives):.6f}")
print(f"Theoretical mean drift: {np.mean(theoretical_drift):.6f}")

# Duration and convexity approximations
current_rate = r[-1]
duration_approx = -rate_derivatives[-1] / current_rate
convexity_approx = rate_second_derivatives[-1] / current_rate

print(f"\nBond Risk Measures (approximations):")
print(f"Modified duration: {duration_approx:.4f}")
print(f"Convexity: {convexity_approx:.4f}")

🧬 Example 3: Population Dynamics with Uncertainty

Application: Biological population modeling with environmental stochasticity.

# Stochastic logistic growth with environmental noise
# dN/dt = rN(1 - N/K) + σN ξ(t)
# where ξ(t) is white noise

def stochastic_logistic(N0, r, K, sigma, T, N_steps):
    """Simulate stochastic logistic growth"""
    dt = T / N_steps
    t = np.linspace(0, T, N_steps+1)
    N = np.zeros(N_steps+1)
    N[0] = N0

    for i in range(N_steps):
        # Deterministic growth
        growth = r * N[i] * (1 - N[i]/K) * dt
        # Stochastic perturbation
        noise = sigma * N[i] * np.random.randn() * np.sqrt(dt)
        N[i+1] = max(0, N[i] + growth + noise)  # Prevent negative population

    return t, N

# Population parameters
N0 = 10        # Initial population
r = 0.5        # Growth rate
K = 1000       # Carrying capacity
sigma = 0.1    # Environmental noise strength
T = 20         # Time horizon
N_steps = 500  # Time steps

# Simulate population
np.random.seed(456)
t, N = stochastic_logistic(N0, r, K, sigma, T, N_steps)

# Fit interpolator
population_interp = SplineInterpolator(smoothing=0.1)
population_interp.fit(t, N)

# Use gamma distribution (positive values, flexible shape)
# Gamma parameters chosen to match population characteristics
alpha = 4.0  # Shape parameter
beta = alpha / np.mean(N)  # Rate parameter (scale = 1/beta)

population_interp.set_stochastic_link('gamma', alpha=alpha, beta=beta, method='ito')

# Compute growth rates
growth_rate_func = population_interp.differentiate(order=1)
stochastic_growth_rates = growth_rate_func(t)

# Compare with deterministic logistic growth rate
deterministic_growth_rates = r * N * (1 - N/K)

# Per capita growth rates
per_capita_stochastic = stochastic_growth_rates / N
per_capita_deterministic = deterministic_growth_rates / N

print("Population Dynamics Analysis:")
print(f"Final population: {N[-1]:.0f}")
print(f"Carrying capacity: {K}")
print(f"Mean growth rate: {np.mean(stochastic_growth_rates):.2f}")
print(f"Growth rate volatility: {np.std(stochastic_growth_rates):.2f}")

# Find maximum growth rate period
max_growth_idx = np.argmax(stochastic_growth_rates)
max_growth_time = t[max_growth_idx]
max_growth_pop = N[max_growth_idx]

print(f"\nMaximum growth rate:")
print(f"Time: {max_growth_time:.1f}")
print(f"Population: {max_growth_pop:.0f}")
print(f"Growth rate: {stochastic_growth_rates[max_growth_idx]:.2f}")

# Environmental impact analysis
stochastic_correction = np.mean(stochastic_growth_rates - deterministic_growth_rates)
print(f"\nEnvironmental stochasticity effect:")
print(f"Mean correction: {stochastic_correction:.3f}")

⚙️ Advanced Stochastic Features

Custom Link Functions

Create custom probability distributions:

from pydelt.stochastic import StochasticLinkFunction

class WeibullLink(StochasticLinkFunction):
    def __init__(self, shape, scale):
        self.shape = shape
        self.scale = scale

    def transform(self, x):
        return self.scale * (-np.log(1 - x))**(1/self.shape)

    def inverse_transform(self, y):
        return 1 - np.exp(-(y/self.scale)**self.shape)

    def derivative_transform(self, x, f_prime, f_double_prime=None, method='ito'):
        # Implement Weibull-specific corrections
        pass

# Use custom link
weibull_link = WeibullLink(shape=2.0, scale=1.0)
interpolator.set_stochastic_link(weibull_link, method='ito')

Method Comparison

Compare Itô vs Stratonovich corrections:

# Fit same data with both methods
interp_ito = SplineInterpolator(smoothing=0.1)
interp_ito.fit(data_x, data_y)
interp_ito.set_stochastic_link('lognormal', sigma=0.2, method='ito')

interp_strat = SplineInterpolator(smoothing=0.1)
interp_strat.fit(data_x, data_y)
interp_strat.set_stochastic_link('lognormal', sigma=0.2, method='stratonovich')

# Compare derivatives
ito_deriv = interp_ito.differentiate(order=1)(eval_points)
strat_deriv = interp_strat.differentiate(order=1)(eval_points)

difference = np.mean(np.abs(ito_deriv - strat_deriv))
print(f"Itô vs Stratonovich difference: {difference:.4f}")

Parameter Sensitivity Analysis

Analyze sensitivity to distribution parameters:

# Test different volatility levels
sigmas = [0.1, 0.2, 0.3, 0.4]
derivative_means = []

for sigma in sigmas:
    interp = SplineInterpolator(smoothing=0.1)
    interp.fit(data_x, data_y)
    interp.set_stochastic_link('lognormal', sigma=sigma, method='ito')

    deriv_func = interp.differentiate(order=1)
    derivatives = deriv_func(eval_points)
    derivative_means.append(np.mean(derivatives))

print("Sensitivity to volatility parameter:")
for sigma, mean_deriv in zip(sigmas, derivative_means):
    print(f"σ = {sigma:.1f}: Mean derivative = {mean_deriv:.3f}")

🎓 Best Practices

Distribution Selection: 1. Stock Prices: Log-normal (positive, multiplicative noise) 2. Interest Rates: Normal (can be negative, additive noise) 3. Waiting Times: Exponential or Gamma (positive, memoryless or aging) 4. Proportions: Beta (bounded between 0 and 1) 5. Count Data: Poisson (discrete, non-negative)

Method Selection: - Itô Calculus: Mathematical finance, theoretical models - Stratonovich Calculus: Physical systems, engineering applications - When in doubt: Try both and compare results

Parameter Estimation: Use maximum likelihood or method of moments to estimate distribution parameters from data:

from scipy import stats

# Estimate log-normal parameters
sigma_est = np.std(np.log(price_data))
print(f"Estimated volatility: {sigma_est:.3f}")

# Use estimated parameter
interpolator.set_stochastic_link('lognormal', sigma=sigma_est, method='ito')

Validation: Always compare with known analytical solutions when available:

# For geometric Brownian motion: E[dS/dt] = μS
theoretical_drift = mu * stock_prices
numerical_drift = stochastic_derivatives

error = np.sqrt(np.mean((numerical_drift - theoretical_drift)**2))
print(f"Drift estimation error: {error:.4f}")

⚠️ Limitations & Considerations

Numerical Stability: - Second derivatives required for Itô corrections may amplify noise - Use appropriate smoothing parameters - Consider robust interpolation methods for noisy data

Model Assumptions: - Stochastic link functions assume specific probability distributions - Validate distribution assumptions with data - Consider model uncertainty in critical applications

Computational Cost: - Stochastic corrections require additional derivative computations - Itô method needs second derivatives (more expensive) - Consider computational budget for real-time applications

Interpretation: - Stochastic derivatives reflect expected behavior under distributional assumptions - Results depend on chosen probability distribution and parameters - Always validate against domain knowledge and empirical data

🔬 Research Applications

Financial Engineering: - Option pricing with stochastic volatility - Risk-neutral measure transformations - Credit risk modeling with jump processes

Scientific Computing: - Uncertainty quantification in differential equations - Stochastic partial differential equations - Climate modeling with random forcing

Machine Learning: - Bayesian neural networks with stochastic derivatives - Uncertainty-aware optimization - Robust control with distributional assumptions

🔗 Integration with Other Features

Stochastic computing combines powerfully with other pydelt features:

Multivariate Stochastic Derivatives: .. code-block:: python

# Apply stochastic links to multivariate functions from pydelt.multivariate import MultivariateDerivatives

mv = MultivariateDerivatives(SplineInterpolator, smoothing=0.1) mv.fit(input_data, output_data)

# Set stochastic link for underlying interpolator mv.interpolator_class.set_stochastic_link(‘lognormal’, sigma=0.2)

# Gradients now include stochastic corrections stochastic_gradient = mv.gradient()

Neural Networks with Stochastic Links: .. code-block:: python

# Combine automatic differentiation with stochastic corrections nn_stochastic = NeuralNetworkInterpolator(hidden_layers=[128, 64]) nn_stochastic.fit(data_x, data_y) nn_stochastic.set_stochastic_link(‘gamma’, alpha=2.0, beta=1.0)

# Exact derivatives with probabilistic corrections exact_stochastic_deriv = nn_stochastic.differentiate(order=1)

This completes the progressive learning path from basic interpolation to advanced stochastic computing, providing a comprehensive framework for numerical differentiation with uncertainty quantification.