PYTHON ALGO

[fnoor1]

import numpy as np from scipy.stats import norm from scipy.optimize import least_squares import pyswarms as ps

Black-Scholes Model

def black_scholes(S, K, T, r, sigma): d1 = (np.log(S / K) + (r + 0.5 * sigma * 2) T) / (sigma np.sqrt(T)) d2 = d1 - sigma np.sqrt(T) call_price = S norm.cdf(d1) - K np.exp(-r T) norm.cdf(d2) return call_price

Binomial Tree Model

def binomial_tree(S, K, T, r, sigma, steps): dt = T / steps u = np.exp(sigma np.sqrt(dt)) d = 1 / u p = (np.exp(r dt) - d) / (u - d)

# Calculate the option price at each node
option_prices = np.zeros((steps + 1, steps + 1))
for i in range(steps + 1):
    option_prices[i, -1] = max(S * u ** (2 * i - steps) - K, 0)

# Work backward through the tree
for j in range(steps - 1, -1, -1):
    for i in range(j + 1):
        option_prices[i, j] = np.exp(-r * dt) * (p * option_prices[i, j + 1] + (1 - p) * option_prices[i + 1, j + 1])

return option_prices[0, 0]

Heston Model

def heston_call_price(S, K, T, r, sigma, kappa, theta, rho, v0, num_simulations=100000):

Implement the Heston model for option pricing using Monte Carlo simulations

def heston_call_price(S, K, T, r, sigma, kappa, theta, rho, v0, num_simulations=100000, num_steps=100): np.random.seed(42) dt = T / num_steps prices = np.zeros(num_simulations)

for sim in range(num_simulations):
    S_t = S
    v_t = v0

    for step in range(num_steps):
        dW_S = np.random.normal(0, np.sqrt(dt))
        dW_v = rho * dW_S + np.sqrt(1 - rho ** 2) * np.random.normal(0, np.sqrt(dt))

        S_t *= np.exp((r - 0.5 * v_t) * dt + np.sqrt(v_t) * dW_S)
        v_t = max(v_t + kappa * (theta - v_t) * dt + sigma * np.sqrt(v_t) * dW_v, 0)

    prices[sim] = S_t

call_prices = np.maximum(prices - K, 0)
call_price = np.exp(-r * T) * np.mean(call_prices)

return call_price

Heston Model Calibration

def heston_squared_errors(params, S, K, T, r, market_prices): kappa, theta, rho, v0 = params errors = np.zeros_like(market_prices) for i, K_i in enumerate(K): estimated_price = heston_call_price(S, K_i, T, r, sigma, kappa, theta, rho, v0) errors[i] = (estimated_price - market_prices[i]) ** 2 return errors

Custom Model Calibration

def custom_model_weights(weights, S, K, T, r, sigma, steps, kappa, theta, rho, v0, num_simulations=100000): weight_bs, weight_bt, weight_heston = weights price_bs = black_scholes(S, K, T, r, sigma) price_bt = binomial_tree(S, K, T, r, sigma, steps) price_heston = heston_call_price(S, K, T, r, sigma, kappa, theta, rho, v0, num_simulations) return weight_bs price_bs + weight_bt price_bt + weight_heston * price_heston

Objective function for Particle Swarm Optimization (PSO)

def pso_objective_function(params, S, K, T, r, sigma, steps, market_prices, num_simulations=100000): weights, heston_params = np.split(params, [3]) kappa, theta, rho, v0 = heston_params errors = np.zeros_like(market_prices) for i, K_i in enumerate(K): estimated_price = custom_model_weights(weights, S, K_i, T, r, sigma, steps, kappa, theta, rho, v0, num_simulations) errors[i] = (estimated_price - market_prices[i]) ** 2 return np.sum(errors)

def fars_custom_algo(S, K, T, r, sigma, steps, observed_strike_prices, observed_market_prices):

Particle Swarm Optimization (PSO) for calibration

bounds = ([0, 0, 0, 0, 0, -1, 0], [1, 1, 1, np.inf, np.inf, 1, np.inf])
pso_options = {"c1": 0.5, "c2": 0.3, "w": 0.9}
optimizer = ps.single.GlobalBestPSO(n_particles=50, dimensions=7, options=pso_options, bounds=bounds)
cost, pos = optimizer.optimize(pso_objective_function, 100, S=S, K=observed_strike_prices, T=T, r=r, sigma=sigma, steps=steps, market_prices=observed_market_prices)

# Extract the calibrated weights and Heston parameters
weights, heston_params = np.split(pos, [3])
kappa, theta, rho, v0 = heston_params

# Calculate the custom call price using calibrated weights and Heston parameters
custom_call_price = custom_model_weights(weights, S, K, T, r, sigma, steps, kappa, theta, rho, v0)
return custom_call_price

if name == "main":

Replace these variables with your own data

observed_market_prices = np.array([8.0, 6.0, 4.0])  # Example observed option prices
observed_strike_prices = np.array([95.0, 100.0, 105.0])  # Example observed strike prices
S = 100
K = 100
T = 1
r = 0.05
sigma = 0.2
steps = 100

# Run the custom algo
custom_call_price = fars_custom_algo(S, K, T, r, sigma, steps, observed_strike_prices, observed_market_prices)
print(f"Fars Custom Algo call price: {custom_call_price:.2f}")

[fnoor1]

import numpy as np from scipy.stats import norm from sobol_seq import i4_sobol_generate

def heston_call_price_quasi_monte_carlo(S, K, T, r, sigma, kappa, theta, rho, v0, num_simulations=100000, steps=100): dt = T / steps S_t = np.zeros(num_simulations) v_t = np.zeros(num_simulations)

# Generate Sobol sequences for the two Brownian motions
sobol_points = i4_sobol_generate(2 * steps, num_simulations)

for i in range(num_simulations):
    S_prev = S
    v_prev = v0
    for t in range(steps):
        # Convert Sobol points to standard normal using the inverse CDF
        W_S = norm.ppf(sobol_points[i, 2*t])
        W_v = norm.ppf(sobol_points[i, 2*t + 1])

        # Scale the normal variates to match the desired covariance structure
        W_S *= np.sqrt(dt)
        W_v = W_v * np.sqrt(dt) * rho + W_S * np.sqrt(1 - rho**2)

        S_prev = S_prev * np.exp((r - 0.5 * v_prev) * dt + np.sqrt(v_prev) * W_S)
        v_prev = v_prev + kappa * (theta - v_prev) * dt + sigma * np.sqrt(v_prev) * W_v
        v_prev = max(v_prev, 0)  # Apply truncation or reflection to prevent negative variance
    S_t[i] = S_prev
    v_t[i] = v_prev

# Calculate the call option prices
call_payoffs = np.maximum(S_t - K, 0)
call_price = np.exp(-r * T) * np.mean(call_payoffs)
return call_price

Continue with the rest of the custom model implementation from before

...