Short-Maturity Option Pricing with the Heston Model

Calibrating Heston and Pricing Short-Maturity Options

This notebook prices and calibrates short-maturity options on SM Energy Company, with emphasis on 15-day contracts as a very short-maturity product.

The analysis is organized into three parts:

1. Calibrate the Heston model using the Lewis (2001) Fourier pricing formula.
2. Repeat the calibration using the Carr-Madan (1999) Fourier pricing formula.
3. Use the calibrated Heston parameters to price a 20-day ATM Asian call option by Monte Carlo simulation.

Main assumptions used throughout:

• Current stock price: S0 = 232.90
• Annual risk-free rate: r = 0.015
• One year has 250 trading days
• 15-day maturity: T = 15 / 250
• 20-day maturity: T = 20 / 250
• There are no dividends in the market data, so put-call parity is used without dividends

import warnings

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from scipy.integrate import quad, IntegrationWarning
from scipy.optimize import minimize
from IPython.display import display

# Integration warnings can appear for poor trial parameters during optimization.
# The objective function below penalizes invalid prices, so we hide these warnings to keep the notebook readable.
warnings.filterwarnings("ignore", category=IntegrationWarning)

plt.style.use("seaborn-v0_8-whitegrid")
pd.set_option("display.float_format", lambda x: f"{x:,.6f}")

1. Load and Inspect the Option Data

First we load the Excel file and check that the data has the expected columns and maturities. This is a simple step, but it matters because the calibration results depend on reading the market quotes correctly.

S0 = 232.90
r = 0.015
trading_days_per_year = 250

T_15 = 15 / trading_days_per_year
T_20 = 20 / trading_days_per_year

excel_file = "MScFE 622_Stochastic Modeling_GWP1_Option data.xlsx"

option_data = pd.read_excel(excel_file)
option_data.columns = option_data.columns.str.strip()

print(f"Number of rows: {len(option_data)}")
print(f"Columns: {list(option_data.columns)}")
print(f"Available maturities: {sorted(option_data['Days to maturity'].unique())}")
print(f"Option types: {sorted(option_data['Type'].unique())}")

display(option_data.head(10))

Number of rows: 30
Columns: ['Days to maturity', 'Strike', 'Price', 'Type']
Available maturities: [15, 60, 120]
Option types: ['C', 'P']

	Days to maturity	Strike	Price	Type
0	15	227.500000	10.520000	C
1	15	230.000000	10.050000	C
2	15	232.500000	7.750000	C
3	15	235.000000	6.010000	C
4	15	237.500000	4.750000	C
5	60	227.500000	16.780000	C
6	60	230.000000	17.650000	C
7	60	232.500000	16.860000	C
8	60	235.000000	16.050000	C
9	60	237.500000	15.100000	C

# This summary checks how many call and put quotes are available at each maturity.
summary_by_maturity = (
    option_data
    .groupby(["Days to maturity", "Type"])
    .agg(number_of_quotes=("Price", "count"), min_strike=("Strike", "min"), max_strike=("Strike", "max"))
    .reset_index()
)

display(summary_by_maturity)

	Days to maturity	Type	number_of_quotes	min_strike	max_strike
0	15	C	5	227.500000	237.500000
1	15	P	5	227.500000	237.500000
2	60	C	5	227.500000	237.500000
3	60	P	5	227.500000	237.500000
4	120	C	5	227.500000	237.500000
5	120	P	5	227.500000	237.500000

2. Keep Only the 15-Day Options

For Step 1A and Step 1B we calibrate the Heston model only to the 15-day options. We use both calls and puts, as required.

options_15 = (
    option_data[option_data["Days to maturity"] == 15]
    .copy()
    .sort_values(["Type", "Strike"])
    .reset_index(drop=True)
)

# To make sure the option type labels are clean.
options_15["Type"] = options_15["Type"].str.upper().str.strip()

display(options_15)

	Days to maturity	Strike	Price	Type
0	15	227.500000	10.520000	C
1	15	230.000000	10.050000	C
2	15	232.500000	7.750000	C
3	15	235.000000	6.010000	C
4	15	237.500000	4.750000	C
5	15	227.500000	4.320000	P
6	15	230.000000	5.200000	P
7	15	232.500000	6.450000	P
8	15	235.000000	7.560000	P
9	15	237.500000	8.780000	P

Before calibration, we also check put-call parity for the market quotes. The Heston model prices one no-arbitrage surface, so model call and put prices must satisfy parity. If the market quotes do not exactly satisfy parity, the calibration cannot fit every call and put perfectly at the same time.

calls_15 = options_15[options_15["Type"] == "C"][["Strike", "Price"]].rename(columns={"Price": "Call Price"})
puts_15 = options_15[options_15["Type"] == "P"][["Strike", "Price"]].rename(columns={"Price": "Put Price"})

parity_check = calls_15.merge(puts_15, on="Strike")
parity_check["Market C - P"] = parity_check["Call Price"] - parity_check["Put Price"]
parity_check["Parity RHS"] = S0 - parity_check["Strike"] * np.exp(-r * T_15)
parity_check["Difference"] = parity_check["Market C - P"] - parity_check["Parity RHS"]

display(parity_check.round(4))

	Strike	Call Price	Put Price	Market C - P	Parity RHS	Difference
0	227.500000	10.520000	4.320000	6.200000	5.604700	0.595300
1	230.000000	10.050000	5.200000	4.850000	3.106900	1.743100
2	232.500000	7.750000	6.450000	1.300000	0.609200	0.690800
3	235.000000	6.010000	7.560000	-1.550000	-1.888600	0.338600
4	237.500000	4.750000	8.780000	-4.030000	-4.386300	0.356300

3. Step 1A - Heston Calibration with the Lewis Formula

The Heston model lets the variance move randomly over time. The five parameters we calibrate are:

• kappa: speed of mean reversion in variance
• theta: long-run variance
• sigma: volatility of variance
• rho: correlation between the stock shock and variance shock
• v0: initial variance

The characteristic function below is for the log return. The Lewis formula then uses this characteristic function to compute European call prices by numerical integration. Put prices are calculated from put-call parity.

def heston_characteristic_function(u, T, r, kappa, theta, sigma, rho, v0):
    """Heston characteristic function for the log return ln(S_T / S_0)."""
    i = 1j
    c1 = kappa * theta
    c2 = -np.sqrt((rho * sigma * u * i - kappa) ** 2 - sigma ** 2 * (-u * i - u ** 2))
    c3 = (kappa - rho * sigma * u * i + c2) / (kappa - rho * sigma * u * i - c2)

    h1 = r * u * i * T + (c1 / sigma ** 2) * (
        (kappa - rho * sigma * u * i + c2) * T
        - 2 * np.log((1 - c3 * np.exp(c2 * T)) / (1 - c3))
    )
    h2 = ((kappa - rho * sigma * u * i + c2) / sigma ** 2) * (
        (1 - np.exp(c2 * T)) / (1 - c3 * np.exp(c2 * T))
    )

    return np.exp(h1 + h2 * v0)


def lewis_call_price(S0, K, T, r, params, integration_limit=100):
    """European call price from the Lewis (2001) Heston integral."""
    kappa, theta, sigma, rho, v0 = params
    log_moneyness = np.log(S0 / K)

    def integrand(u):
        shifted_u = u - 0.5j
        numerator = np.exp(1j * u * log_moneyness) * heston_characteristic_function(
            shifted_u, T, r, kappa, theta, sigma, rho, v0
        )
        return np.real(numerator / (u ** 2 + 0.25))

    integral_value, _ = quad(
        integrand,
        0,
        integration_limit,
        epsabs=1e-5,
        epsrel=1e-5,
        limit=100,
    )

    call_price = S0 - np.exp(-r * T) * np.sqrt(S0 * K) * integral_value / np.pi
    return float(call_price)


def put_price_from_parity(call_price, S0, K, T, r):
    """European put price from put-call parity with no dividends."""
    return call_price - S0 + K * np.exp(-r * T)


def price_call_or_put(S0, K, T, r, option_type, params, call_pricer):
    """Price one option. Calls are priced directly; puts use put-call parity."""
    call_price = call_pricer(S0, K, T, r, params)

    if option_type == "C":
        return call_price
    if option_type == "P":
        return put_price_from_parity(call_price, S0, K, T, r)

    raise ValueError("Option type must be 'C' or 'P'.")

Calibration Setup

We use regular mean squared error (MSE): the average squared difference between the market prices and the model prices.

For the optimization, we use L-BFGS-B because it allows simple parameter bounds. The bounds are deliberately broad, but they keep variances and volatility of variance positive and keep correlation inside a realistic range. We do not force the Feller condition because short-maturity option data can be difficult to fit if that extra restriction is imposed.

To reduce dependence on one starting guess, we try three reasonable starting points and keep the calibration with the smallest MSE. This is still simple, but it makes the result less arbitrary.

parameter_names = ["kappa", "theta", "sigma", "rho", "v0"]

parameter_bounds = [
    (0.01, 10.00),   # kappa must be positive
    (0.001, 1.00),   # theta must be positive
    (0.01, 5.00),    # sigma must be positive
    (-0.95, 0.95),   # rho is a correlation, kept away from exactly +/- 1
    (0.001, 1.00),   # v0 must be positive
]

starting_guesses = [
    np.array([1.50, 0.04, 0.30, -0.50, 0.04]),
    np.array([2.00, 0.10, 0.50, -0.40, 0.10]),
    np.array([0.80, 0.08, 1.00, -0.70, 0.08]),
]

bounds_table = pd.DataFrame(parameter_bounds, columns=["Lower Bound", "Upper Bound"], index=parameter_names)
starting_guess_table = pd.DataFrame(starting_guesses, columns=parameter_names)
starting_guess_table.index = ["Start 1", "Start 2", "Start 3"]

display(bounds_table)
display(starting_guess_table)

	Lower Bound	Upper Bound
kappa	0.010000	10.000000
theta	0.001000	1.000000
sigma	0.010000	5.000000
rho	-0.950000	0.950000
v0	0.001000	1.000000

	kappa	theta	sigma	rho	v0
Start 1	1.500000	0.040000	0.300000	-0.500000	0.040000
Start 2	2.000000	0.100000	0.500000	-0.400000	0.100000
Start 3	0.800000	0.080000	1.000000	-0.700000	0.080000

def model_prices_for_options(option_rows, T, params, call_pricer):
    """Return model prices in the same order as the option data table."""
    call_price_by_strike = {}

    # Calls and puts have the same strikes, so this avoids pricing the same call twice.
    for strike in sorted(option_rows["Strike"].unique()):
        call_price_by_strike[strike] = call_pricer(S0, strike, T, r, params)

    model_prices = []
    for row in option_rows.itertuples(index=False):
        call_price = call_price_by_strike[row.Strike]
        if row.Type == "C":
            model_prices.append(call_price)
        else:
            model_prices.append(put_price_from_parity(call_price, S0, row.Strike, T, r))

    return np.array(model_prices)


def heston_mse(params, option_rows, T, call_pricer):
    """Mean squared error between market prices and Heston model prices."""
    kappa, theta, sigma, rho, v0 = params

    # This guard keeps the optimizer away from invalid Heston parameters.
    if kappa <= 0 or theta <= 0 or sigma <= 0 or v0 <= 0 or abs(rho) >= 1:
        return 1e8

    try:
        model_prices = model_prices_for_options(option_rows, T, params, call_pricer)
    except Exception:
        return 1e8

    if not np.all(np.isfinite(model_prices)):
        return 1e8

    market_prices = option_rows["Price"].to_numpy()
    return float(np.mean((market_prices - model_prices) ** 2))


def run_heston_calibration(label, call_pricer):
    """Run bounded MSE calibration from a small set of starting points."""
    calibration_results = []
    best_result = None

    for start_number, start in enumerate(starting_guesses, start=1):
        result = minimize(
            heston_mse,
            start,
            args=(options_15, T_15, call_pricer),
            method="L-BFGS-B",
            bounds=parameter_bounds,
            options={"maxiter": 80, "ftol": 1e-8},
        )

        calibration_results.append({
            "Method": label,
            "Start": start_number,
            "MSE": result.fun,
            "Success": result.success,
            "Function Evaluations": result.nfev,
        })

        if best_result is None or result.fun < best_result.fun:
            best_result = result

    display(pd.DataFrame(calibration_results))
    return best_result


def make_parameter_table(params, mse_value, method):
    table = pd.DataFrame({
        "Parameter": parameter_names,
        "Value": params,
    })
    table.loc[len(table)] = ["MSE", mse_value]
    table.insert(0, "Method", method)
    return table


def make_comparison_table(option_rows, T, params, call_pricer):
    comparison = option_rows.copy()
    comparison = comparison.rename(columns={"Price": "Market Price"})
    comparison["Model Price"] = model_prices_for_options(option_rows, T, params, call_pricer)
    comparison["Pricing Error"] = comparison["Model Price"] - comparison["Market Price"]
    return comparison


def plot_market_vs_model(comparison_table, title):
    fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)

    for ax, option_type, name in zip(axes, ["C", "P"], ["Calls", "Puts"]):
        subset = comparison_table[comparison_table["Type"] == option_type]
        ax.plot(subset["Strike"], subset["Market Price"], marker="o", label="Market price")
        ax.plot(subset["Strike"], subset["Model Price"], marker="s", linestyle="--", label="Model price")
        ax.set_title(name)
        ax.set_xlabel("Strike")
        ax.set_ylabel("Option price")
        ax.legend()

    fig.suptitle(title)
    plt.tight_layout()
    plt.show()

lewis_result = run_heston_calibration("Lewis", lewis_call_price)
lewis_params = lewis_result.x
lewis_mse = lewis_result.fun

lewis_parameter_table = make_parameter_table(lewis_params, lewis_mse, "Lewis")
display(lewis_parameter_table)

	Method	Start	MSE	Success	Function Evaluations
0	Lewis	1	0.246154	True	186
1	Lewis	2	0.246161	True	246
2	Lewis	3	0.246145	True	162

	Method	Parameter	Value
0	Lewis	kappa	0.794732
1	Lewis	theta	0.001000
2	Lewis	sigma	1.348435
3	Lewis	rho	-0.869343
4	Lewis	v0	0.105330
5	Lewis	MSE	0.246145

lewis_comparison = make_comparison_table(options_15, T_15, lewis_params, lewis_call_price)
display(lewis_comparison.round(4))

plot_market_vs_model(lewis_comparison, "Lewis Calibration: 15-Day Market Prices vs Heston Model Prices")

	Days to maturity	Strike	Market Price	Type	Model Price	Pricing Error
0	15	227.500000	10.520000	C	10.515100	-0.004900
1	15	230.000000	10.050000	C	8.858500	-1.191500
2	15	232.500000	7.750000	C	7.319400	-0.430600
3	15	235.000000	6.010000	C	5.906400	-0.103600
4	15	237.500000	4.750000	C	4.628200	-0.121800
5	15	227.500000	4.320000	P	4.910400	0.590400
6	15	230.000000	5.200000	P	5.751600	0.551600
7	15	232.500000	6.450000	P	6.710200	0.260200
8	15	235.000000	7.560000	P	7.795000	0.235000
9	15	237.500000	8.780000	P	9.014500	0.234500

Lewis Calibration Interpretation

The Lewis calibration gives one Heston parameter set for both calls and puts. Since model put prices are calculated by put-call parity, the model is internally consistent, but it cannot exactly match market calls and puts if the observed quotes do not satisfy parity perfectly.

For this short 15-day maturity, the calibration mainly identifies the short-run variance level and the skew-related parameters. Long-run parameters such as kappa and theta are harder to identify from only very short maturity data, so we interpret them with more caution than v0, sigma, and rho.

The long-run variance parameter theta reached the lower bound used in the calibration. This does not invalidate the calibration, but it means theta should be interpreted cautiously. Since the calibration uses only 15-day options, the short-maturity data is more informative about current variance and volatility skew than about the long-run variance level.

4. Step 1B - Heston Calibration with the Carr-Madan Formula

Now we repeat the same calibration using the Carr-Madan approach. The model is still Heston; only the Fourier pricing formula for European calls changes.

Carr-Madan uses a damping parameter alpha to make the call price transform integrable. We use alpha = 1.5, which is a common choice for equity call pricing.

def carr_madan_call_price(S0, K, T, r, params, alpha=1.5, integration_limit=100):
    """European call price using the Carr-Madan (1999) Fourier formula."""
    kappa, theta, sigma, rho, v0 = params
    log_strike = np.log(K)

    def log_price_characteristic_function(u):
        # The Heston function above is for log return, so we multiply by exp(iu log(S0))
        # to get the characteristic function of log(S_T).
        return np.exp(1j * u * np.log(S0)) * heston_characteristic_function(
            u, T, r, kappa, theta, sigma, rho, v0
        )

    def integrand(u):
        shifted_u = u - (alpha + 1) * 1j
        numerator = np.exp(-r * T) * log_price_characteristic_function(shifted_u)
        denominator = alpha ** 2 + alpha - u ** 2 + 1j * (2 * alpha + 1) * u
        return np.real(np.exp(-1j * u * log_strike) * numerator / denominator)

    integral_value, _ = quad(
        integrand,
        0,
        integration_limit,
        epsabs=1e-5,
        epsrel=1e-5,
        limit=100,
    )

    call_price = np.exp(-alpha * log_strike) * integral_value / np.pi
    return float(call_price)

carr_madan_result = run_heston_calibration("Carr-Madan", carr_madan_call_price)
carr_madan_params = carr_madan_result.x
carr_madan_mse = carr_madan_result.fun

carr_madan_parameter_table = make_parameter_table(carr_madan_params, carr_madan_mse, "Carr-Madan")
display(carr_madan_parameter_table)

	Method	Start	MSE	Success	Function Evaluations
0	Carr-Madan	1	0.246199	True	186
1	Carr-Madan	2	0.246206	True	240
2	Carr-Madan	3	0.246188	True	174

	Method	Parameter	Value
0	Carr-Madan	kappa	0.794212
1	Carr-Madan	theta	0.001000
2	Carr-Madan	sigma	1.359562
3	Carr-Madan	rho	-0.861178
4	Carr-Madan	v0	0.105490
5	Carr-Madan	MSE	0.246188

carr_madan_comparison = make_comparison_table(options_15, T_15, carr_madan_params, carr_madan_call_price)
display(carr_madan_comparison.round(4))

plot_market_vs_model(carr_madan_comparison, "Carr-Madan Calibration: 15-Day Market Prices vs Heston Model Prices")

	Days to maturity	Strike	Market Price	Type	Model Price	Pricing Error
0	15	227.500000	10.520000	C	10.515300	-0.004700
1	15	230.000000	10.050000	C	8.858400	-1.191600
2	15	232.500000	7.750000	C	7.319200	-0.430800
3	15	235.000000	6.010000	C	5.906300	-0.103700
4	15	237.500000	4.750000	C	4.628300	-0.121700
5	15	227.500000	4.320000	P	4.910600	0.590600
6	15	230.000000	5.200000	P	5.751500	0.551500
7	15	232.500000	6.450000	P	6.710000	0.260000
8	15	235.000000	7.560000	P	7.794900	0.234900
9	15	237.500000	8.780000	P	9.014600	0.234600

Lewis vs Carr-Madan Comparison

The Lewis and Carr-Madan results should be very close because both methods price the same Heston model. Small differences come from numerical integration and the different form of each Fourier inversion.

For the Asian option in Step 1C, we choose the calibration method with the lower MSE. If the two MSE values are almost the same, this choice is not economically dramatic, but it gives us a clear rule instead of choosing by eye.

method_comparison = pd.DataFrame({
    "Method": ["Lewis", "Carr-Madan"],
    "MSE": [lewis_mse, carr_madan_mse],
})

display(method_comparison)

if lewis_mse <= carr_madan_mse:
    selected_method = "Lewis"
    selected_params = lewis_params
    selected_mse = lewis_mse
else:
    selected_method = "Carr-Madan"
    selected_params = carr_madan_params
    selected_mse = carr_madan_mse

selected_parameter_table = make_parameter_table(selected_params, selected_mse, selected_method)
print(f"Selected calibration for Monte Carlo: {selected_method}")
display(selected_parameter_table)

	Method	MSE
0	Lewis	0.246145
1	Carr-Madan	0.246188

Selected calibration for Monte Carlo: Lewis

	Method	Parameter	Value
0	Lewis	kappa	0.794732
1	Lewis	theta	0.001000
2	Lewis	sigma	1.348435
3	Lewis	rho	-0.869343
4	Lewis	v0	0.105330
5	Lewis	MSE	0.246145

kappa, theta, sigma, rho, v0 = selected_params

feller_left = 2 * kappa * theta
feller_right = sigma ** 2
feller_satisfied = feller_left >= feller_right

feller_table = pd.DataFrame({
    "Quantity": ["2 * kappa * theta", "sigma^2", "Feller condition satisfied?"],
    "Value": [feller_left, feller_right, feller_satisfied]
})

display(feller_table)

	Quantity	Value
0	2 * kappa * theta	0.001589
1	sigma^2	1.818276
2	Feller condition satisfied?	False

Feller Condition Diagnostic

The Feller condition is checked as a diagnostic but not imposed during calibration. We did not impose it because the assignment focuses on fitting short-maturity market option prices, and imposing the condition can make the fit worse when the available maturity is very short.

5. Step 1C - Price a 20-Day ATM Asian Call Option

The Asian option payoff depends on the average stock price over the 20 trading days. We include the current stock price S0 in the average, so the average uses 21 stock prices: today plus 20 simulated daily prices.

For the simulation we use the risk-neutral Heston dynamics. The stock drift is the risk-free rate r. The variance process is simulated with a simple full-truncation Euler scheme, which keeps variance non-negative during simulation.

def simulate_heston_paths(S0, r, params, n_steps, n_paths, trading_days_per_year=250, seed=622):
    """Simulate stock paths under the risk-neutral Heston model using daily steps."""
    kappa, theta, sigma, rho, v0 = params
    dt = 1 / trading_days_per_year
    rng = np.random.default_rng(seed)

    stock_paths = np.empty((n_paths, n_steps + 1))
    stock_paths[:, 0] = S0
    variance = np.full(n_paths, v0)

    for step in range(1, n_steps + 1):
        z_variance = rng.standard_normal(n_paths)
        z_independent = rng.standard_normal(n_paths)
        z_stock = rho * z_variance + np.sqrt(max(1 - rho ** 2, 0)) * z_independent

        variance_positive = np.maximum(variance, 0)

        # Stock is simulated in log form so prices remain positive.
        stock_paths[:, step] = stock_paths[:, step - 1] * np.exp(
            (r - 0.5 * variance_positive) * dt
            + np.sqrt(variance_positive * dt) * z_stock
        )

        # Full truncation keeps the simulated variance from becoming negative.
        variance = (
            variance
            + kappa * (theta - variance_positive) * dt
            + sigma * np.sqrt(variance_positive * dt) * z_variance
        )
        variance = np.maximum(variance, 0)

    return stock_paths


def asian_call_payoffs(stock_paths, strike):
    """Arithmetic-average Asian call payoff, including S0 in the average."""
    average_stock_price = stock_paths.mean(axis=1)
    return np.maximum(average_stock_price - strike, 0)


def price_asian_call_heston(S0, strike, r, params, n_steps, n_paths, seed=622):
    stock_paths = simulate_heston_paths(
        S0=S0,
        r=r,
        params=params,
        n_steps=n_steps,
        n_paths=n_paths,
        trading_days_per_year=trading_days_per_year,
        seed=seed,
    )

    payoffs = asian_call_payoffs(stock_paths, strike)
    T = n_steps / trading_days_per_year
    discount_factor = np.exp(-r * T)

    fair_price = discount_factor * payoffs.mean()
    standard_error = discount_factor * payoffs.std(ddof=1) / np.sqrt(n_paths)

    return {
        "Fair Price": fair_price,
        "Monte Carlo Std. Error": standard_error,
        "Number of Simulations": n_paths,
        "Number of Daily Steps": n_steps,
        "Seed": seed,
    }

asian_strike = S0
asian_steps = 20
number_of_simulations = 100_000
random_seed = 622
bank_fee_rate = 0.04

asian_result = price_asian_call_heston(
    S0=S0,
    strike=asian_strike,
    r=r,
    params=selected_params,
    n_steps=asian_steps,
    n_paths=number_of_simulations,
    seed=random_seed,
)

fair_price = asian_result["Fair Price"]
fee_amount = bank_fee_rate * fair_price
client_price = fair_price + fee_amount

asian_summary = pd.DataFrame({
    "Item": [
        "Selected calibration method",
        "Asian option maturity",
        "Strike",
        "Simulations",
        "Random seed",
        "Fair price",
        "Monte Carlo standard error",
        "4% bank fee",
        "Final client price",
    ],
    "Value": [
        selected_method,
        "20 trading days",
        asian_strike,
        number_of_simulations,
        random_seed,
        fair_price,
        asian_result["Monte Carlo Std. Error"],
        fee_amount,
        client_price,
    ],
})

display(asian_summary)

	Item	Value
0	Selected calibration method	Lewis
1	Asian option maturity	20 trading days
2	Strike	232.900000
3	Simulations	100000
4	Random seed	622
5	Fair price	4.762997
6	Monte Carlo standard error	0.018475
7	4% bank fee	0.190520
8	Final client price	4.953517

mc_se = asian_result["Monte Carlo Std. Error"]

lower_95 = fair_price - 1.96 * mc_se
upper_95 = fair_price + 1.96 * mc_se

mc_interval_table = pd.DataFrame({
    "Item": ["Fair price", "Monte Carlo standard error", "95% lower bound", "95% upper bound"],
    "Value": [fair_price, mc_se, lower_95, upper_95]
})

display(mc_interval_table)

	Item	Value
0	Fair price	4.762997
1	Monte Carlo standard error	0.018475
2	95% lower bound	4.726785
3	95% upper bound	4.799209

Client Explanation

The Asian call option pays off if the average SM Energy stock price over the next 20 trading days is above the strike price of 232.90. Because it uses an average price instead of only the final stock price, the option is less sensitive to one large move on the last day.

Using the calibrated Heston model, we estimate the fair value by simulating many possible daily stock paths under the risk-neutral measure and discounting the average payoff back to today. The bank then adds a 4% fee to the fair value, and that gives the final price paid by the client.

6. Step 1 Summary and Recommended Action

This final table collects the main Step 1 outputs in one place for the report.

step1_summary = pd.DataFrame({
    "Item": [
        "Lewis MSE",
        "Carr-Madan MSE",
        "Selected method",
        "Fair Asian call price",
        "4% bank fee",
        "Final client price",
    ],
    "Value": [
        f"{lewis_mse:.6f}",
        f"{carr_madan_mse:.6f}",
        selected_method,
        f"{fair_price:.6f}",
        f"{fee_amount:.6f}",
        f"{client_price:.6f}",
    ]
})

display(step1_summary)

	Item	Value
0	Lewis MSE	0.246145
1	Carr-Madan MSE	0.246188
2	Selected method	Lewis
3	Fair Asian call price	4.762997
4	4% bank fee	0.190520
5	Final client price	4.953517

Recommended Action

Based on the selected calibration, the estimated fair value of the 20-day at-the-money Asian call is shown in the summary table above. After applying the bank's 4% fee, the final client quote is also shown above. We recommend quoting approximately $4.95 per option unit, subject to the assumptions of no dividends, constant risk-free rate, and the calibrated short-maturity volatility dynamics.

Step 2 - Bates Model with Jumps

In Step 2 we extend the Heston model by adding jumps in the stock price. This is useful because market option prices can reflect the risk of sudden large moves, not just continuous volatility.

For this step we use the 60-day market options, since the target put option has a 70-day maturity and 60 days is the closest available maturity in the data.

7. Keep Only the 60-Day Options

We now filter the same Excel data to the 60-day maturity. As in Step 1, we use both call and put prices in the calibration.

T_60 = 60 / trading_days_per_year

options_60 = (
    option_data[option_data["Days to maturity"] == 60]
    .copy()
    .sort_values(["Type", "Strike"])
    .reset_index(drop=True)
)

options_60["Type"] = options_60["Type"].str.upper().str.strip()

display(options_60)

	Days to maturity	Strike	Price	Type
0	60	227.500000	16.780000	C
1	60	230.000000	17.650000	C
2	60	232.500000	16.860000	C
3	60	235.000000	16.050000	C
4	60	237.500000	15.100000	C
5	60	227.500000	11.030000	P
6	60	230.000000	12.150000	P
7	60	232.500000	13.370000	P
8	60	235.000000	14.750000	P
9	60	237.500000	15.620000	P

Before calibrating, we do the same quick market sanity check as before. This helps explain why the model may not match every quote exactly: the model prices a smooth no-arbitrage curve, while the observed market quotes can be a little noisy.

calls_60 = options_60[options_60["Type"] == "C"][["Strike", "Price"]].rename(columns={"Price": "Call Price"})
puts_60 = options_60[options_60["Type"] == "P"][["Strike", "Price"]].rename(columns={"Price": "Put Price"})

parity_check_60 = calls_60.merge(puts_60, on="Strike")
parity_check_60["Market C - P"] = parity_check_60["Call Price"] - parity_check_60["Put Price"]
parity_check_60["Parity RHS"] = S0 - parity_check_60["Strike"] * np.exp(-r * T_60)
parity_check_60["Difference"] = parity_check_60["Market C - P"] - parity_check_60["Parity RHS"]

display(parity_check_60.round(4))

	Strike	Call Price	Put Price	Market C - P	Parity RHS	Difference
0	227.500000	16.780000	11.030000	5.750000	6.217500	-0.467500
1	230.000000	17.650000	12.150000	5.500000	3.726500	1.773500
2	232.500000	16.860000	13.370000	3.490000	1.235500	2.254500
3	235.000000	16.050000	14.750000	1.300000	-1.255500	2.555500
4	237.500000	15.100000	15.620000	-0.520000	-3.746500	3.226500

8. Bates Characteristic Function and Pricing Functions

The Bates model combines the Heston stochastic variance process with lognormal jumps in the stock price. We use eight parameters:

• kappa, theta, sigma, rho, v0: the same Heston variance parameters as Step 1
• lambda_j: annual jump intensity
• mu_j: average log jump size
• delta_j: volatility of the log jump size

The jump compensator keeps the stock price process risk-neutral, so the expected growth rate is still tied to the risk-free rate.

def bates_characteristic_function(u, T, r, kappa, theta, sigma, rho, v0, lambda_j, mu_j, delta_j):
    """Bates characteristic function for the log return ln(S_T / S_0)."""
    jump_compensator = np.exp(mu_j + 0.5 * delta_j ** 2) - 1
    adjusted_r = r - lambda_j * jump_compensator

    heston_part = heston_characteristic_function(
        u, T, adjusted_r, kappa, theta, sigma, rho, v0
    )

    jump_part = np.exp(
        lambda_j * T * (
            np.exp(1j * u * mu_j - 0.5 * delta_j ** 2 * u ** 2) - 1
        )
    )

    return heston_part * jump_part


def bates_lewis_call_price(S0, K, T, r, params, integration_limit=100):
    """European call price under Bates using the Lewis integral."""
    kappa, theta, sigma, rho, v0, lambda_j, mu_j, delta_j = params
    log_moneyness = np.log(S0 / K)

    def integrand(u):
        shifted_u = u - 0.5j
        numerator = np.exp(1j * u * log_moneyness) * bates_characteristic_function(
            shifted_u, T, r, kappa, theta, sigma, rho, v0, lambda_j, mu_j, delta_j
        )
        return np.real(numerator / (u ** 2 + 0.25))

    integral_value, _ = quad(
        integrand,
        0,
        integration_limit,
        epsabs=1e-5,
        epsrel=1e-5,
        limit=100,
    )

    call_price = S0 - np.exp(-r * T) * np.sqrt(S0 * K) * integral_value / np.pi
    return float(call_price)


def bates_carr_madan_call_price(S0, K, T, r, params, alpha=1.5, integration_limit=100):
    """European call price under Bates using the Carr-Madan formula."""
    kappa, theta, sigma, rho, v0, lambda_j, mu_j, delta_j = params
    log_strike = np.log(K)

    def log_price_characteristic_function(u):
        return np.exp(1j * u * np.log(S0)) * bates_characteristic_function(
            u, T, r, kappa, theta, sigma, rho, v0, lambda_j, mu_j, delta_j
        )

    def integrand(u):
        shifted_u = u - (alpha + 1) * 1j
        numerator = np.exp(-r * T) * log_price_characteristic_function(shifted_u)
        denominator = alpha ** 2 + alpha - u ** 2 + 1j * (2 * alpha + 1) * u
        return np.real(np.exp(-1j * u * log_strike) * numerator / denominator)

    integral_value, _ = quad(
        integrand,
        0,
        integration_limit,
        epsabs=1e-5,
        epsrel=1e-5,
        limit=100,
    )

    call_price = np.exp(-alpha * log_strike) * integral_value / np.pi
    return float(call_price)


def bates_price_by_type(S0, K, T, r, option_type, params, call_pricer):
    """Price a Bates call directly and a Bates put through put-call parity."""
    call_price = call_pricer(S0, K, T, r, params)

    if option_type == "C":
        return call_price
    if option_type == "P":
        return put_price_from_parity(call_price, S0, K, T, r)

    raise ValueError("Option type must be 'C' or 'P'.")

Bates Calibration Setup

We keep the same calibration idea as Step 1: minimize the MSE between market prices and model prices. The bounds are broad enough to let the optimizer move, but they avoid impossible values such as negative variance or correlation outside [-1, 1].

We do not force the Feller condition here either. The main goal is to fit the 60-day market option prices, and adding the jump terms already makes the calibration more flexible.

bates_parameter_names = [
    "kappa",
    "theta",
    "sigma",
    "rho",
    "v0",
    "lambda_j",
    "mu_j",
    "delta_j",
]

bates_parameter_bounds = [
    (0.01, 10.00),    # kappa must be positive
    (0.001, 1.00),    # theta must be positive
    (0.01, 5.00),     # sigma must be positive
    (-0.95, 0.95),    # rho is a correlation
    (0.001, 1.00),    # v0 must be positive
    (0.00, 3.00),     # lambda_j is the annual jump intensity
    (-0.50, 0.50),    # mu_j is the average log jump size
    (0.001, 1.00),    # delta_j must be positive
]

bates_starting_guesses = [
    np.array([1.20, 0.06, 0.60, -0.50, 0.06, 0.20, -0.05, 0.20]),
    np.array([2.00, 0.08, 0.80, -0.40, 0.08, 0.50, -0.08, 0.25]),
    np.array([0.80, 0.10, 1.00, -0.70, 0.10, 0.10, 0.00, 0.15]),
]

bates_bounds_table = pd.DataFrame(
    bates_parameter_bounds,
    columns=["Lower Bound", "Upper Bound"],
    index=bates_parameter_names,
)

bates_starting_guess_table = pd.DataFrame(bates_starting_guesses, columns=bates_parameter_names)
bates_starting_guess_table.index = ["Start 1", "Start 2", "Start 3"]

display(bates_bounds_table)
display(bates_starting_guess_table)

	Lower Bound	Upper Bound
kappa	0.010000	10.000000
theta	0.001000	1.000000
sigma	0.010000	5.000000
rho	-0.950000	0.950000
v0	0.001000	1.000000
lambda_j	0.000000	3.000000
mu_j	-0.500000	0.500000
delta_j	0.001000	1.000000

	kappa	theta	sigma	rho	v0	lambda_j	mu_j	delta_j
Start 1	1.200000	0.060000	0.600000	-0.500000	0.060000	0.200000	-0.050000	0.200000
Start 2	2.000000	0.080000	0.800000	-0.400000	0.080000	0.500000	-0.080000	0.250000
Start 3	0.800000	0.100000	1.000000	-0.700000	0.100000	0.100000	0.000000	0.150000

def model_prices_for_bates_options(option_rows, T, params, call_pricer):
    """Return Bates model prices in the same order as the option data table."""
    call_price_by_strike = {}

    for strike in sorted(option_rows["Strike"].unique()):
        call_price_by_strike[strike] = call_pricer(S0, strike, T, r, params)

    model_prices = []
    for row in option_rows.itertuples(index=False):
        call_price = call_price_by_strike[row.Strike]
        if row.Type == "C":
            model_prices.append(call_price)
        else:
            model_prices.append(put_price_from_parity(call_price, S0, row.Strike, T, r))

    return np.array(model_prices)


def bates_mse(params, option_rows, T, call_pricer):
    """Mean squared error between market prices and Bates model prices."""
    kappa, theta, sigma, rho, v0, lambda_j, mu_j, delta_j = params

    if (
        kappa <= 0
        or theta <= 0
        or sigma <= 0
        or v0 <= 0
        or abs(rho) >= 1
        or lambda_j < 0
        or delta_j <= 0
    ):
        return 1e8

    try:
        model_prices = model_prices_for_bates_options(option_rows, T, params, call_pricer)
    except Exception:
        return 1e8

    if not np.all(np.isfinite(model_prices)) or np.any(model_prices < 0):
        return 1e8

    market_prices = option_rows["Price"].to_numpy()
    return float(np.mean((market_prices - model_prices) ** 2))


def run_bates_calibration(label, call_pricer):
    """Run Bates MSE calibration and keep the best successful start when available."""
    results = []
    rows = []

    for start_number, start in enumerate(bates_starting_guesses, start=1):
        result = minimize(
            bates_mse,
            start,
            args=(options_60, T_60, call_pricer),
            method="L-BFGS-B",
            bounds=bates_parameter_bounds,
            options={"maxiter": 80, "ftol": 1e-8},
        )

        results.append((start_number, result))
        rows.append({
            "Method": label,
            "Start": start_number,
            "MSE": result.fun,
            "Success": result.success,
            "Function Evaluations": result.nfev,
        })

    successful_results = [(start, result) for start, result in results if result.success]
    if successful_results:
        selected_start, selected_result = min(successful_results, key=lambda item: item[1].fun)
        selection_note = "best successful start"
    else:
        selected_start, selected_result = min(results, key=lambda item: item[1].fun)
        selection_note = "lowest MSE start because no run fully converged"

    calibration_table = pd.DataFrame(rows)
    display(calibration_table)
    print(f"Using Start {selected_start} for {label}: {selection_note}.")

    return selected_result, calibration_table


def make_bates_parameter_table(params, mse_value, method):
    table = pd.DataFrame({
        "Parameter": bates_parameter_names,
        "Value": params,
    })
    table.loc[len(table)] = ["MSE", mse_value]
    table.insert(0, "Method", method)
    return table


def make_bates_comparison_table(option_rows, T, params, call_pricer):
    comparison = option_rows.copy()
    comparison = comparison.rename(columns={"Price": "Market Price"})
    comparison["Model Price"] = model_prices_for_bates_options(option_rows, T, params, call_pricer)
    comparison["Pricing Error"] = comparison["Model Price"] - comparison["Market Price"]
    return comparison


def plot_pricing_error_by_strike(comparison_table, title):
    fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)

    for ax, option_type, name in zip(axes, ["C", "P"], ["Calls", "Puts"]):
        subset = comparison_table[comparison_table["Type"] == option_type]
        ax.axhline(0, color="black", linewidth=1)
        ax.plot(subset["Strike"], subset["Pricing Error"], marker="o", color="tab:red")
        ax.set_title(name)
        ax.set_xlabel("Strike")
        ax.set_ylabel("Model price - market price")

    fig.suptitle(title)
    plt.tight_layout()
    plt.show()

9. Step 2A - Bates Calibration with the Lewis Formula

We first calibrate the Bates model to the 60-day calls and puts using the Lewis pricing formula. Puts are still priced through put-call parity, so the model remains internally consistent.

bates_lewis_result, bates_lewis_runs = run_bates_calibration(
    "Bates-Lewis",
    bates_lewis_call_price,
)

bates_lewis_params = bates_lewis_result.x
bates_lewis_mse = bates_lewis_result.fun

bates_lewis_parameter_table = make_bates_parameter_table(
    bates_lewis_params,
    bates_lewis_mse,
    "Bates-Lewis",
)

display(bates_lewis_parameter_table)

	Method	Start	MSE	Success	Function Evaluations
0	Bates-Lewis	1	1.341693	False	1134
1	Bates-Lewis	2	1.367897	True	630
2	Bates-Lewis	3	1.596102	True	243

Using Start 2 for Bates-Lewis: best successful start.

	Method	Parameter	Value
0	Bates-Lewis	kappa	2.034435
1	Bates-Lewis	theta	0.183842
2	Bates-Lewis	sigma	0.446345
3	Bates-Lewis	rho	0.288915
4	Bates-Lewis	v0	0.001010
5	Bates-Lewis	lambda_j	0.711281
6	Bates-Lewis	mu_j	0.313091
7	Bates-Lewis	delta_j	0.070289
8	Bates-Lewis	MSE	1.367897

bates_lewis_comparison = make_bates_comparison_table(
    options_60,
    T_60,
    bates_lewis_params,
    bates_lewis_call_price,
)

display(bates_lewis_comparison.round(4))

plot_market_vs_model(
    bates_lewis_comparison,
    "Bates-Lewis Calibration: 60-Day Market Prices vs Model Prices",
)

plot_pricing_error_by_strike(
    bates_lewis_comparison,
    "Bates-Lewis Pricing Error by Strike",
)

	Days to maturity	Strike	Market Price	Type	Model Price	Pricing Error
0	60	227.500000	16.780000	C	17.416500	0.636500
1	60	230.000000	17.650000	C	16.410600	-1.239400
2	60	232.500000	16.860000	C	15.483600	-1.376400
3	60	235.000000	16.050000	C	14.628300	-1.421700
4	60	237.500000	15.100000	C	13.837700	-1.262300
5	60	227.500000	11.030000	P	11.198900	0.168900
6	60	230.000000	12.150000	P	12.684100	0.534100
7	60	232.500000	13.370000	P	14.248100	0.878100
8	60	235.000000	14.750000	P	15.883800	1.133800
9	60	237.500000	15.620000	P	17.584300	1.964300

10. Step 2B - Bates Calibration with the Carr-Madan Formula

Next we repeat the same calibration using the Carr-Madan call pricing formula. The data, parameter bounds, starting guesses, and MSE objective are kept the same so the comparison is fair.

bates_carr_madan_result, bates_carr_madan_runs = run_bates_calibration(
    "Bates-Carr-Madan",
    bates_carr_madan_call_price,
)

bates_carr_madan_params = bates_carr_madan_result.x
bates_carr_madan_mse = bates_carr_madan_result.fun

bates_carr_madan_parameter_table = make_bates_parameter_table(
    bates_carr_madan_params,
    bates_carr_madan_mse,
    "Bates-Carr-Madan",
)

display(bates_carr_madan_parameter_table)

	Method	Start	MSE	Success	Function Evaluations
0	Bates-Carr-Madan	1	1.326551	False	1170
1	Bates-Carr-Madan	2	1.369157	True	882
2	Bates-Carr-Madan	3	1.596138	True	243

Using Start 2 for Bates-Carr-Madan: best successful start.

	Method	Parameter	Value
0	Bates-Carr-Madan	kappa	2.015882
1	Bates-Carr-Madan	theta	0.232440
2	Bates-Carr-Madan	sigma	0.381825
3	Bates-Carr-Madan	rho	0.616397
4	Bates-Carr-Madan	v0	0.002850
5	Bates-Carr-Madan	lambda_j	0.379818
6	Bates-Carr-Madan	mu_j	0.415602
7	Bates-Carr-Madan	delta_j	0.029798
8	Bates-Carr-Madan	MSE	1.369157

bates_carr_madan_comparison = make_bates_comparison_table(
    options_60,
    T_60,
    bates_carr_madan_params,
    bates_carr_madan_call_price,
)

display(bates_carr_madan_comparison.round(4))

plot_market_vs_model(
    bates_carr_madan_comparison,
    "Bates-Carr-Madan Calibration: 60-Day Market Prices vs Model Prices",
)

plot_pricing_error_by_strike(
    bates_carr_madan_comparison,
    "Bates-Carr-Madan Pricing Error by Strike",
)

	Days to maturity	Strike	Market Price	Type	Model Price	Pricing Error
0	60	227.500000	16.780000	C	17.458900	0.678900
1	60	230.000000	17.650000	C	16.430500	-1.219500
2	60	232.500000	16.860000	C	15.482200	-1.377800
3	60	235.000000	16.050000	C	14.608500	-1.441500
4	60	237.500000	15.100000	C	13.803600	-1.296400
5	60	227.500000	11.030000	P	11.241300	0.211300
6	60	230.000000	12.150000	P	12.704000	0.554000
7	60	232.500000	13.370000	P	14.246700	0.876700
8	60	235.000000	14.750000	P	15.864000	1.114000
9	60	237.500000	15.620000	P	17.550200	1.930200

11. Compare the Two Bates Calibrations

Both pricing formulas are valuing the same Bates model, so their fitted prices should be reasonably close. Small differences can still happen because the Fourier formulas are written differently and the numerical optimizer may settle in slightly different places.

The Bates-Lewis MSE is 1.367897, while the Bates-Carr-Madan MSE is 1.369157. The difference is small, so both methods provide similar calibration quality. Bates-Lewis was selected because it has the slightly lower MSE.

For the final 70-day put price, we select the calibration with the lower MSE among the successful calibration runs. If the MSE values are very close, this rule keeps the decision simple and transparent.

bates_method_comparison = pd.DataFrame({
    "Method": ["Bates-Lewis", "Bates-Carr-Madan"],
    "MSE": [bates_lewis_mse, bates_carr_madan_mse],
})

display(bates_method_comparison)

if bates_lewis_mse <= bates_carr_madan_mse:
    selected_bates_method = "Bates-Lewis"
    selected_bates_params = bates_lewis_params
    selected_bates_mse = bates_lewis_mse
    selected_bates_call_pricer = bates_lewis_call_price
else:
    selected_bates_method = "Bates-Carr-Madan"
    selected_bates_params = bates_carr_madan_params
    selected_bates_mse = bates_carr_madan_mse
    selected_bates_call_pricer = bates_carr_madan_call_price

print(f"Selected Bates calibration for Step 2C: {selected_bates_method}")
display(make_bates_parameter_table(selected_bates_params, selected_bates_mse, selected_bates_method))

	Method	MSE
0	Bates-Lewis	1.367897
1	Bates-Carr-Madan	1.369157

Selected Bates calibration for Step 2C: Bates-Lewis

	Method	Parameter	Value
0	Bates-Lewis	kappa	2.034435
1	Bates-Lewis	theta	0.183842
2	Bates-Lewis	sigma	0.446345
3	Bates-Lewis	rho	0.288915
4	Bates-Lewis	v0	0.001010
5	Bates-Lewis	lambda_j	0.711281
6	Bates-Lewis	mu_j	0.313091
7	Bates-Lewis	delta_j	0.070289
8	Bates-Lewis	MSE	1.367897

def bates_bound_diagnostic_table(params, bounds, names):
    rows = []
    for name, value, (lower, upper) in zip(names, params, bounds):
        rows.append({
            "Parameter": name,
            "Value": value,
            "At Lower Bound?": np.isclose(value, lower, atol=1e-5),
            "At Upper Bound?": np.isclose(value, upper, atol=1e-5),
        })
    return pd.DataFrame(rows)

selected_bates_bound_check = bates_bound_diagnostic_table(
    selected_bates_params,
    bates_parameter_bounds,
    bates_parameter_names,
)

display(selected_bates_bound_check)

	Parameter	Value	At Lower Bound?	At Upper Bound?
0	kappa	2.034435	False	False
1	theta	0.183842	False	False
2	sigma	0.446345	False	False
3	rho	0.288915	False	False
4	v0	0.001010	True	False
5	lambda_j	0.711281	False	False
6	mu_j	0.313091	False	False
7	delta_j	0.070289	False	False

Bates Calibration Interpretation

The model is trying to fit both calls and puts with one consistent pricing curve. The 60-day market quotes are not perfectly smooth, so a non-zero MSE is expected. What matters most here is whether the calibrated prices track the overall level and shape of the market prices.

Some optimization starts can reach a lower MSE without fully converging. For the reported parameters, we use the best successful calibration run for each pricing method so the final comparison is more stable.

The bound diagnostic above is included so we do not hide any calibration pressure. In the selected Bates calibration, v0 is very close to the lower bound, so the initial variance estimate should be interpreted cautiously. This does not invalidate the price, but it suggests the 60-day data is putting more weight on the jump and medium-term variance parameters than on a high starting variance level.

The selected Bates calibration produces a positive correlation parameter and a positive average jump size. These values should be interpreted as calibration parameters rather than literal forecasts. With a small 60-day option dataset, the optimizer may use the jump component and the variance component together to match the observed option prices.

12. Step 2C - Price the 70-Day 95% Moneyness Put

The target derivative is a 70-day put option with 95% moneyness. Since moneyness is 95%, the strike is 0.95 * S0. This put pays max(K - S_T, 0) at maturity.

put_maturity_days = 70
T_70 = put_maturity_days / trading_days_per_year
put_moneyness = 0.95
put_strike = put_moneyness * S0

selected_bates_call_price_70 = selected_bates_call_pricer(
    S0,
    put_strike,
    T_70,
    r,
    selected_bates_params,
)

estimated_bates_put_price_70 = put_price_from_parity(
    selected_bates_call_price_70,
    S0,
    put_strike,
    T_70,
    r,
)

bates_put_summary = pd.DataFrame({
    "Item": [
        "Current stock price",
        "Strike",
        "Maturity in days",
        "Maturity in years",
        "Selected calibration method",
        "Estimated put price",
    ],
    "Value": [
        f"{S0:.6f}",
        f"{put_strike:.6f}",
        put_maturity_days,
        f"{T_70:.6f}",
        selected_bates_method,
        f"{estimated_bates_put_price_70:.6f}",
    ]
})

display(bates_put_summary)

	Item	Value
0	Current stock price	232.900000
1	Strike	221.255000
2	Maturity in days	70
3	Maturity in years	0.280000
4	Selected calibration method	Bates-Lewis
5	Estimated put price	9.664595

Fee Note

For Step 2C, we report the model-implied fair value of the 70-day put option. The assignment explicitly states the 4% bank fee for the Step 1 OTC Asian option. Since Step 2 describes the new product as a put option rather than an OTC Asian instrument and does not explicitly restate the fee, the value reported here is the fair premium before any additional commercial markup. If the same 4% bank fee were applied, the client quote would be fair price x 1.04.

put_fee_rate = 0.04
put_fee_amount = put_fee_rate * estimated_bates_put_price_70
put_client_price_with_fee = estimated_bates_put_price_70 + put_fee_amount

step2_fee_sensitivity = pd.DataFrame({
    "Item": [
        "Fair put price",
        "Optional 4% fee",
        "Put price with optional 4% fee",
    ],
    "Value": [
        f"{estimated_bates_put_price_70:.6f}",
        f"{put_fee_amount:.6f}",
        f"{put_client_price_with_fee:.6f}",
    ],
})

display(step2_fee_sensitivity)

	Item	Value
0	Fair put price	9.664595
1	Optional 4% fee	0.386584
2	Put price with optional 4% fee	10.051179

Technical Interpretation

The 70-day put is out of the money because the strike is below the current stock price. It still has value because there is a chance the stock price falls below the strike before maturity.

Using 60-day options to price a 70-day option is a reasonable short extrapolation because the maturities are close, but it is still an assumption. The final price should therefore be read as a model-based estimate under the calibrated 60-day market conditions, not as a guaranteed trading price.

print(
    f"Based on the selected 60-day calibration, the estimated value of the 70-day "
    f"95% moneyness put is {estimated_bates_put_price_70:.6f}. "
    f"We recommend quoting approximately ${estimated_bates_put_price_70:.2f} per option unit, "
    f"subject to the assumptions of no dividends, a constant risk-free rate, and the fitted medium-maturity option dynamics."
)

Based on the selected 60-day calibration, the estimated value of the 70-day 95% moneyness put is 9.664595. We recommend quoting approximately $9.66 per option unit, subject to the assumptions of no dividends, a constant risk-free rate, and the fitted medium-maturity option dynamics.

Client Explanation

From a client perspective, this option gives protection if SM Energy falls below 221.26 over about 70 trading days. The strike is below today's stock price, so the option is cheaper than protection that starts exactly at today's price, but it still has value because a meaningful drop is possible before maturity.

The estimated fair price is $9.66. If the bank applies the same 4% fee convention used in Step 1, the client quote becomes approximately $10.05. The quote should be reviewed if market prices, interest rates, or the stock price move materially before the trade is executed.

13. Step 2 Summary

This compact table collects the main Step 2 outputs for the report.

step2_summary = pd.DataFrame({
    "Item": [
        "Bates-Lewis MSE",
        "Bates-Carr-Madan MSE",
        "Selected method",
        "Current stock price",
        "Put strike",
        "Put maturity",
        "Estimated put price",
        "Optional put price with 4% fee",
    ],
    "Value": [
        f"{bates_lewis_mse:.6f}",
        f"{bates_carr_madan_mse:.6f}",
        selected_bates_method,
        f"{S0:.6f}",
        f"{put_strike:.6f}",
        f"{put_maturity_days} trading days",
        f"{estimated_bates_put_price_70:.6f}",
        f"{put_client_price_with_fee:.6f}",
    ],
})

display(step2_summary)

	Item	Value
0	Bates-Lewis MSE	1.367897
1	Bates-Carr-Madan MSE	1.369157
2	Selected method	Bates-Lewis
3	Current stock price	232.900000
4	Put strike	221.255000
5	Put maturity	70 trading days
6	Estimated put price	9.664595
7	Optional put price with 4% fee	10.051179

Step 3 - Euribor Term Structure and Interest Rate Simulation

In Step 3 we model future interest rates because changes in rates can affect future derivative prices. The bank operates mostly in a European setting, so we use the Euribor rates given in the assignment.

The workflow is:

1. Build and interpolate the Euribor term structure.
2. Fit a CIR interest rate model to the interpolated weekly rates.
3. Simulate one year of future rates and compare the expected future 12-month Euribor with today's 12-month Euribor.

14. Step 3A - Euribor Term Structure

We start by putting the given Euribor quotes into a small table. The rates are entered in decimal form, but we also show percentages because they are easier to read.

from scipy.interpolate import CubicSpline

euribor_term_structure = pd.DataFrame({
    "Maturity": ["1 week", "1 month", "3 months", "6 months", "12 months"],
    "Maturity (Years)": [1 / 52, 1 / 12, 3 / 12, 6 / 12, 1.0],
    "Euribor Rate": [0.00648, 0.00679, 0.01173, 0.01809, 0.02556],
})

euribor_term_structure["Euribor Rate (%)"] = 100 * euribor_term_structure["Euribor Rate"]

display(euribor_term_structure)

	Maturity	Maturity (Years)	Euribor Rate	Euribor Rate (%)
0	1 week	0.019231	0.006480	0.648000
1	1 month	0.083333	0.006790	0.679000
2	3 months	0.250000	0.011730	1.173000
3	6 months	0.500000	0.018090	1.809000
4	12 months	1.000000	0.025560	2.556000

fig, ax = plt.subplots(figsize=(8, 4))

ax.plot(
    euribor_term_structure["Maturity (Years)"],
    100 * euribor_term_structure["Euribor Rate"],
    marker="o",
    label="Observed Euribor rates",
)

ax.set_title("Original Euribor Term Structure")
ax.set_xlabel("Maturity in years")
ax.set_ylabel("Annual rate (%)")
ax.legend()
plt.tight_layout()
plt.show()

Weekly Interpolation

The market quotes are only available at five maturities. To fit a smoother curve, we use cubic spline interpolation to estimate weekly rates from week 1 to week 52. This does not create new market data; it just gives a smooth curve between the observed points.

weeks = np.arange(1, 53)
weekly_maturities = weeks / 52

euribor_spline = CubicSpline(
    euribor_term_structure["Maturity (Years)"],
    euribor_term_structure["Euribor Rate"],
    bc_type="natural",
)

weekly_euribor_rates = euribor_spline(weekly_maturities)

weekly_euribor_curve = pd.DataFrame({
    "Week": weeks,
    "Maturity (Years)": weekly_maturities,
    "Interpolated Euribor Rate": weekly_euribor_rates,
    "Interpolated Euribor Rate (%)": 100 * weekly_euribor_rates,
})

display(weekly_euribor_curve.head(10))
display(weekly_euribor_curve.tail(5))

	Week	Maturity (Years)	Interpolated Euribor Rate	Interpolated Euribor Rate (%)
0	1	0.019231	0.006480	0.648000
1	2	0.038462	0.006506	0.650640
2	3	0.057692	0.006572	0.657231
3	4	0.076923	0.006717	0.671728
4	5	0.096154	0.006978	0.697792
5	6	0.115385	0.007357	0.735714
6	7	0.134615	0.007836	0.783612
7	8	0.153846	0.008396	0.839566
8	9	0.173077	0.009017	0.901657
9	10	0.192308	0.009680	0.967965

	Week	Maturity (Years)	Interpolated Euribor Rate	Interpolated Euribor Rate (%)
47	48	0.923077	0.024576	2.457603
48	49	0.942308	0.024823	2.482334
49	50	0.961538	0.025070	2.506952
50	51	0.980769	0.025315	2.531495
51	52	1.000000	0.025560	2.556000

fig, ax = plt.subplots(figsize=(9, 4))

ax.plot(
    weekly_euribor_curve["Maturity (Years)"],
    100 * weekly_euribor_curve["Interpolated Euribor Rate"],
    label="Interpolated weekly curve",
)
ax.scatter(
    euribor_term_structure["Maturity (Years)"],
    100 * euribor_term_structure["Euribor Rate"],
    color="tab:red",
    zorder=3,
    label="Observed Euribor rates",
)

ax.set_title("Observed and Interpolated Euribor Term Structure")
ax.set_xlabel("Maturity in years")
ax.set_ylabel("Annual rate (%)")
ax.legend()
plt.tight_layout()
plt.show()

15. CIR Term Structure Calibration

The CIR model keeps interest rates non-negative through a square-root process. For calibration, we compare the interpolated weekly Euribor rates with the zero rates implied by the CIR zero-coupon bond formula.

We treat the Euribor rates as annualized zero rates for this fitting exercise. This is a simplifying assumption, but it keeps the analysis transparent and consistent with the assignment goal.

def cir_zero_coupon_price(maturity, kappa, theta, sigma, r0):
    """CIR zero-coupon bond price P(0,T)."""
    maturity = np.asarray(maturity)
    gamma = np.sqrt(kappa ** 2 + 2 * sigma ** 2)
    exp_gamma_t = np.exp(gamma * maturity)

    denominator = (gamma + kappa) * (exp_gamma_t - 1) + 2 * gamma
    B = 2 * (exp_gamma_t - 1) / denominator
    A = (
        2 * gamma * np.exp((kappa + gamma) * maturity / 2) / denominator
    ) ** (2 * kappa * theta / sigma ** 2)

    return A * np.exp(-B * r0)


def cir_implied_zero_rate(maturity, kappa, theta, sigma, r0):
    """Convert CIR zero-coupon prices into continuously compounded zero rates."""
    prices = cir_zero_coupon_price(maturity, kappa, theta, sigma, r0)
    prices = np.maximum(prices, 1e-12)
    return -np.log(prices) / maturity


def cir_mse(params, market_maturities, market_rates):
    """MSE between interpolated market rates and CIR fitted rates."""
    kappa, theta, sigma, r0 = params

    if kappa <= 0 or theta <= 0 or sigma <= 0 or r0 < 0:
        return 1e8

    try:
        fitted_rates = cir_implied_zero_rate(market_maturities, kappa, theta, sigma, r0)
    except Exception:
        return 1e8

    if not np.all(np.isfinite(fitted_rates)):
        return 1e8

    return float(np.mean((market_rates - fitted_rates) ** 2))

CIR Calibration Setup

The four calibrated parameters are:

• kappa: speed of mean reversion
• theta: long-run average rate level
• sigma: rate volatility
• r0: starting short rate

The bounds keep the parameters positive and keep the starting rate in a reasonable range for the Euribor curve. We use a few starting guesses and keep the calibration with the lowest MSE. This is still simple, but it avoids relying too much on one starting point.

cir_parameter_names = ["kappa", "theta", "sigma", "r0"]

cir_parameter_bounds = [
    (0.001, 5.00),    # kappa must be positive
    (0.0001, 0.10),   # theta must be positive
    (0.001, 0.50),    # sigma must be positive
    (0.0001, 0.05),   # r0 must be positive and realistic for this curve
]

cir_starting_guesses = [
    np.array([0.50, 0.030, 0.030, 0.0065]),
    np.array([1.00, 0.040, 0.080, 0.0065]),
    np.array([2.00, 0.050, 0.150, 0.0050]),
    np.array([4.00, 0.040, 0.180, 0.0040]),
]

cir_bounds_table = pd.DataFrame(
    cir_parameter_bounds,
    columns=["Lower Bound", "Upper Bound"],
    index=cir_parameter_names,
)

cir_starting_guess_table = pd.DataFrame(cir_starting_guesses, columns=cir_parameter_names)
cir_starting_guess_table.index = ["Start 1", "Start 2", "Start 3", "Start 4"]

display(cir_bounds_table)
display(cir_starting_guess_table)

	Lower Bound	Upper Bound
kappa	0.001000	5.000000
theta	0.000100	0.100000
sigma	0.001000	0.500000
r0	0.000100	0.050000

	kappa	theta	sigma	r0
Start 1	0.500000	0.030000	0.030000	0.006500
Start 2	1.000000	0.040000	0.080000	0.006500
Start 3	2.000000	0.050000	0.150000	0.005000
Start 4	4.000000	0.040000	0.180000	0.004000

cir_calibration_runs = []
best_cir_result = None

for start_number, start in enumerate(cir_starting_guesses, start=1):
    result = minimize(
        cir_mse,
        start,
        args=(weekly_maturities, weekly_euribor_rates),
        method="L-BFGS-B",
        bounds=cir_parameter_bounds,
        options={"maxiter": 500, "ftol": 1e-14},
    )

    cir_calibration_runs.append({
        "Start": start_number,
        "MSE": result.fun,
        "Success": result.success,
        "Function Evaluations": result.nfev,
    })

    if best_cir_result is None or result.fun < best_cir_result.fun:
        best_cir_result = result

cir_calibration_runs_table = pd.DataFrame(cir_calibration_runs)
display(cir_calibration_runs_table)

cir_params = best_cir_result.x
cir_mse_value = best_cir_result.fun

cir_parameter_table = pd.DataFrame({
    "Parameter": cir_parameter_names + ["MSE"],
    "Value": list(cir_params) + [cir_mse_value],
})

display(cir_parameter_table)

	Start	MSE	Success	Function Evaluations
0	1	0.000001	True	50
1	2	0.000000	True	35
2	3	0.000000	True	30
3	4	0.000001	True	30

	Parameter	Value
0	kappa	1.999904
1	theta	0.041966
2	sigma	0.150004
3	r0	0.003946
4	MSE	0.000000

cir_fitted_rates = cir_implied_zero_rate(weekly_maturities, *cir_params)

cir_comparison = pd.DataFrame({
    "Week": weeks,
    "Maturity (Years)": weekly_maturities,
    "Interpolated Market Rate": weekly_euribor_rates,
    "CIR Fitted Rate": cir_fitted_rates,
    "Fitting Error": cir_fitted_rates - weekly_euribor_rates,
})

display(cir_comparison.round(6))

	Week	Maturity (Years)	Interpolated Market Rate	CIR Fitted Rate	Fitting Error
0	1	0.019231	0.006480	0.004667	-0.001813
1	2	0.038462	0.006506	0.005371	-0.001135
2	3	0.057692	0.006572	0.006057	-0.000515
3	4	0.076923	0.006717	0.006726	0.000008
4	5	0.096154	0.006978	0.007378	0.000400
5	6	0.115385	0.007357	0.008013	0.000656
6	7	0.134615	0.007836	0.008633	0.000797
7	8	0.153846	0.008396	0.009238	0.000842
8	9	0.173077	0.009017	0.009827	0.000811
9	10	0.192308	0.009680	0.010402	0.000723
10	11	0.211538	0.010366	0.010963	0.000598
11	12	0.230769	0.011056	0.011511	0.000455
12	13	0.250000	0.011730	0.012044	0.000314
13	14	0.269231	0.012373	0.012565	0.000192
14	15	0.288462	0.012984	0.013074	0.000090
15	16	0.307692	0.013564	0.013570	0.000006
16	17	0.326923	0.014115	0.014054	-0.000061
17	18	0.346154	0.014639	0.014527	-0.000112
18	19	0.365385	0.015138	0.014988	-0.000150
19	20	0.384615	0.015614	0.015438	-0.000176
20	21	0.403846	0.016068	0.015878	-0.000191
21	22	0.423077	0.016504	0.016307	-0.000197
22	23	0.442308	0.016922	0.016726	-0.000195
23	24	0.461538	0.017324	0.017136	-0.000188
24	25	0.480769	0.017713	0.017535	-0.000177
25	26	0.500000	0.018090	0.017926	-0.000164
26	27	0.519231	0.018457	0.018307	-0.000150
27	28	0.538462	0.018815	0.018680	-0.000135
28	29	0.557692	0.019164	0.019044	-0.000120
29	30	0.576923	0.019504	0.019400	-0.000104
30	31	0.596154	0.019835	0.019747	-0.000088
31	32	0.615385	0.020159	0.020087	-0.000072
32	33	0.634615	0.020476	0.020419	-0.000057
33	34	0.653846	0.020785	0.020744	-0.000041
34	35	0.673077	0.021088	0.021061	-0.000027
35	36	0.692308	0.021384	0.021371	-0.000013
36	37	0.711538	0.021674	0.021674	0.000001
37	38	0.730769	0.021958	0.021971	0.000013
38	39	0.750000	0.022238	0.022261	0.000023
39	40	0.769231	0.022512	0.022545	0.000033
40	41	0.788462	0.022782	0.022822	0.000040
41	42	0.807692	0.023047	0.023093	0.000046
42	43	0.826923	0.023309	0.023359	0.000049
43	44	0.846154	0.023568	0.023618	0.000050
44	45	0.865385	0.023824	0.023873	0.000049
45	46	0.884615	0.024077	0.024121	0.000045
46	47	0.903846	0.024327	0.024365	0.000037
47	48	0.923077	0.024576	0.024603	0.000027
48	49	0.942308	0.024823	0.024836	0.000013
49	50	0.961538	0.025070	0.025064	-0.000005
50	51	0.980769	0.025315	0.025288	-0.000027
51	52	1.000000	0.025560	0.025506	-0.000054

fig, ax = plt.subplots(figsize=(9, 4))

ax.plot(
    cir_comparison["Maturity (Years)"],
    100 * cir_comparison["Interpolated Market Rate"],
    label="Interpolated market rate",
)
ax.plot(
    cir_comparison["Maturity (Years)"],
    100 * cir_comparison["CIR Fitted Rate"],
    linestyle="--",
    label="CIR fitted rate",
)

ax.set_title("CIR Fit to Interpolated Weekly Euribor Rates")
ax.set_xlabel("Maturity in years")
ax.set_ylabel("Annual rate (%)")
ax.legend()
plt.tight_layout()
plt.show()

CIR Calibration Interpretation

The CIR curve follows the upward-sloping Euribor term structure closely, and the MSE is small because the fitted curve is close to the interpolated weekly rates. The calibrated r0 is below the current 12-month Euribor, which is consistent with the original term structure being upward sloping.

None of the selected CIR parameters is at its stated calibration bound, so the fit does not appear to be driven by a boundary solution. The volatility parameter sigma is useful for simulation, but it should still be interpreted cautiously because we are fitting only a term structure of rates, not interest-rate option prices. In practice, volatility would usually be supported by cap, floor, or swaption data.

16. Step 3B - Simulate 12-Month Euribor Rates for One Year

Using the calibrated CIR parameters, we simulate daily rates for one trading year. We use 250 daily steps and 100,000 Monte Carlo paths, matching the trading-day convention used earlier in the assignment.

The simulation uses full truncation, meaning negative simulated values are pushed back to zero. This keeps the square-root process numerically stable and consistent with the non-negative rate structure of the CIR model.

def simulate_cir_rate_paths(r0, kappa, theta, sigma, n_steps, n_paths, dt, seed=622):
    """Simulate CIR rate paths with full truncation to keep rates non-negative."""
    rng = np.random.default_rng(seed)
    rate_paths = np.empty((n_paths, n_steps + 1))
    rate_paths[:, 0] = r0

    for step in range(1, n_steps + 1):
        z = rng.standard_normal(n_paths)
        previous_rate = rate_paths[:, step - 1]
        positive_rate = np.maximum(previous_rate, 0)

        rate_paths[:, step] = (
            previous_rate
            + kappa * (theta - positive_rate) * dt
            + sigma * np.sqrt(positive_rate * dt) * z
        )
        rate_paths[:, step] = np.maximum(rate_paths[:, step], 0)

    return rate_paths

cir_simulation_steps = trading_days_per_year
cir_number_of_simulations = 100_000
cir_random_seed = 622
cir_dt = 1 / trading_days_per_year

kappa_cir, theta_cir, sigma_cir, r0_cir = cir_params

cir_rate_paths = simulate_cir_rate_paths(
    r0=r0_cir,
    kappa=kappa_cir,
    theta=theta_cir,
    sigma=sigma_cir,
    n_steps=cir_simulation_steps,
    n_paths=cir_number_of_simulations,
    dt=cir_dt,
    seed=cir_random_seed,
)

terminal_cir_rates = cir_rate_paths[:, -1]

print(f"Simulated paths: {cir_number_of_simulations:,}")
print(f"Daily steps: {cir_simulation_steps}")

Simulated paths: 100,000
Daily steps: 250

fig, ax = plt.subplots(figsize=(9, 4))

sample_path_count = 20
time_grid_days = np.arange(cir_simulation_steps + 1)

for path_number in range(sample_path_count):
    ax.plot(time_grid_days, 100 * cir_rate_paths[path_number], linewidth=1, alpha=0.8)

ax.set_title("Sample Simulated Euribor Rate Paths")
ax.set_xlabel("Trading day")
ax.set_ylabel("Rate (%)")
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(8, 4))

ax.hist(100 * terminal_cir_rates, bins=50, color="tab:blue", alpha=0.75, edgecolor="white")
ax.axvline(100 * terminal_cir_rates.mean(), color="tab:red", linestyle="--", label="Mean")
ax.set_title("Simulated 12-Month Euribor in One Year")
ax.set_xlabel("Terminal rate (%)")
ax.set_ylabel("Frequency")
ax.legend()
plt.tight_layout()
plt.show()

expected_terminal_rate = terminal_cir_rates.mean()
median_terminal_rate = np.median(terminal_cir_rates)
std_terminal_rate = terminal_cir_rates.std(ddof=1)
lower_95_rate, upper_95_rate = np.percentile(terminal_cir_rates, [2.5, 97.5])
min_terminal_rate = terminal_cir_rates.min()
max_terminal_rate = terminal_cir_rates.max()
current_12m_euribor = euribor_term_structure.loc[
    euribor_term_structure["Maturity"] == "12 months",
    "Euribor Rate",
].iloc[0]

simulation_summary = pd.DataFrame({
    "Item": [
        "Current 12-month Euribor",
        "Expected 12-month Euribor in 1 year",
        "Median terminal rate",
        "Standard deviation",
        "95% lower percentile",
        "95% upper percentile",
        "Minimum terminal rate",
        "Maximum terminal rate",
    ],
    "Rate": [
        current_12m_euribor,
        expected_terminal_rate,
        median_terminal_rate,
        std_terminal_rate,
        lower_95_rate,
        upper_95_rate,
        min_terminal_rate,
        max_terminal_rate,
    ],
})

simulation_summary["Rate (%)"] = 100 * simulation_summary["Rate"]

display(simulation_summary)

	Item	Rate	Rate (%)
0	Current 12-month Euribor	0.025560	2.556000
1	Expected 12-month Euribor in 1 year	0.036852	3.685171
2	Median terminal rate	0.035255	3.525459
3	Standard deviation	0.013507	1.350734
4	95% lower percentile	0.015333	1.533329
5	95% upper percentile	0.067600	6.760039
6	Minimum terminal rate	0.002346	0.234625
7	Maximum terminal rate	0.135022	13.502184

rate_difference = expected_terminal_rate - current_12m_euribor

if rate_difference > 0:
    rate_direction = "higher than"
elif rate_difference < 0:
    rate_direction = "lower than"
else:
    rate_direction = "about the same as"

print(
    f"The expected 12-month Euribor in one year is {100 * expected_terminal_rate:.4f}%, "
    f"which is {rate_direction} the current 12-month Euribor of {100 * current_12m_euribor:.4f}%."
)

The expected 12-month Euribor in one year is 3.6852%, which is higher than the current 12-month Euribor of 2.5560%.

Simulation Interpretation

The simulated expected 12-month Euribor is higher than the current 12-month Euribor. This follows from the upward-sloping Euribor curve and the calibrated mean-reversion level. The percentile range is also important because future rates are uncertain, not just a single number.

For derivative pricing, higher future interest rates can change discounting and forward pricing assumptions. The exact effect depends on the product payoff, but the main point is that future client quotes may need to be updated if the rate environment moves materially.

Non-Technical Explanation

The current Euribor curve suggests that rates may move higher over the next year. In the simulation, the average 12-month Euribor in one year is above today's 12-month Euribor rate of 2.556%.

For the bank, this means future product prices may not be the same as today's prices. If rates rise, discount factors and forward pricing inputs change, so future client quotes should be refreshed using the latest rate curve rather than relying on today's assumptions.

17. Step 3 Summary and Recommendation

This final table collects the main Step 3 results for the report.

step3_summary = pd.DataFrame({
    "Item": [
        "Calibrated kappa",
        "Calibrated theta",
        "Calibrated sigma",
        "Calibrated r0",
        "CIR calibration MSE",
        "Current 12-month Euribor",
        "Expected 12-month Euribor in 1 year",
        "95% lower bound",
        "95% upper bound",
        "Number of simulations",
        "Random seed",
    ],
    "Value": [
        f"{kappa_cir:.6f}",
        f"{theta_cir:.6f}",
        f"{sigma_cir:.6f}",
        f"{r0_cir:.6f}",
        f"{cir_mse_value:.10f}",
        f"{current_12m_euribor:.6f}",
        f"{expected_terminal_rate:.6f}",
        f"{lower_95_rate:.6f}",
        f"{upper_95_rate:.6f}",
        f"{cir_number_of_simulations:,}",
        cir_random_seed,
    ],
})

display(step3_summary)

	Item	Value
0	Calibrated kappa	1.999904
1	Calibrated theta	0.041966
2	Calibrated sigma	0.150004
3	Calibrated r0	0.003946
4	CIR calibration MSE	0.0000001739
5	Current 12-month Euribor	0.025560
6	Expected 12-month Euribor in 1 year	0.036852
7	95% lower bound	0.015333
8	95% upper bound	0.067600
9	Number of simulations	100,000
10	Random seed	622

Recommendation

Based on the calibrated rate curve and simulation, we should expect the rate input used for future pricing to be higher than today's 12-month Euribor on average. We recommend updating derivative quotes with the latest Euribor curve at the time of pricing and checking rate sensitivity, especially for products with longer maturities or payoff structures that are sensitive to discounting.