Stochastic Volatility, Jump Diffusion, and American Option Pricing¶

Core inputs

  • Initial stock price: $80
  • Risk-free rate: 5.5%
  • Diffusion volatility input: 35%
  • Time to maturity: 0.25 years
  • ATM strike for Questions 5-10 and 13: $80

Model and simulation choices

  • Fixed random seeds are used for reproducibility.
  • The stochastic-volatility paths are generated with a full-truncation Euler scheme.
  • The jump model is simulated in log-price form.
  • Delta and gamma are estimated by central finite differences with a $1 spot bump.
  • The American call is checked with a least-squares Monte Carlo routine, while the final interpretation follows the standard no-dividend early-exercise result.
In [1]:
import math
from dataclasses import dataclass

import numpy as np
import pandas as pd
from IPython.display import Markdown, display


BASE_PARAMS = {
    "S0": 80.0,
    "K_atm": 80.0,
    "r": 0.055,
    "sigma": 0.35,
    "T": 0.25,
}

HESTON_PARAMS = {
    "v0": 0.032,
    "kappa": 1.85,
    "theta": 0.045,
    "xi": 0.35,
}

MERTON_PARAMS = {
    "mu_j": -0.5,
    "delta_j": 0.22,
}

SIMULATION_SETTINGS = {
    "european_paths": 120_000,
    "barrier_paths": 120_000,
    "american_paths": 90_000,
    "steps": 100,
    "bump_size": 1.0,
}


@dataclass(frozen=True)
class MCResult:
    price: float
    stderr: float


def _mc_summary(discounted_payoff: np.ndarray) -> MCResult:
    price = float(discounted_payoff.mean())
    stderr = float(discounted_payoff.std(ddof=1) / math.sqrt(discounted_payoff.size))
    return MCResult(price=price, stderr=stderr)


def _basis_functions(stock_slice: np.ndarray, variance_slice: np.ndarray, strike: float) -> np.ndarray:
    m = stock_slice / strike
    return np.column_stack(
        [
            np.ones_like(m),
            m,
            m**2,
            variance_slice,
            variance_slice**2,
            m * variance_slice,
        ]
    )


def heston_random_inputs(n_paths: int, n_steps: int, seed: int) -> tuple[np.ndarray, np.ndarray]:
    rng = np.random.default_rng(seed)
    z_price = rng.standard_normal((n_paths, n_steps))
    z_independent = rng.standard_normal((n_paths, n_steps))
    return z_price, z_independent


def simulate_heston(
    S0: float,
    rho: float,
    z_price: np.ndarray,
    z_independent: np.ndarray,
    store_paths: bool = False,
) -> tuple[np.ndarray, np.ndarray]:
    n_paths, n_steps = z_price.shape
    dt = BASE_PARAMS["T"] / n_steps
    sqrt_dt = math.sqrt(dt)
    corr_scale = math.sqrt(1.0 - rho**2)

    if store_paths:
        stock_paths = np.empty((n_paths, n_steps + 1))
        variance_paths = np.empty((n_paths, n_steps + 1))
        stock_paths[:, 0] = S0
        variance_paths[:, 0] = HESTON_PARAMS["v0"]
        stock_state = stock_paths[:, 0].copy()
        variance_state = variance_paths[:, 0].copy()
    else:
        stock_state = np.full(n_paths, S0, dtype=float)
        variance_state = np.full(n_paths, HESTON_PARAMS["v0"], dtype=float)

    for step in range(n_steps):
        z_s = z_price[:, step]
        z_v = rho * z_s + corr_scale * z_independent[:, step]
        variance_floor = np.clip(variance_state, 0.0, None)
        diffusion = np.sqrt(variance_floor) * sqrt_dt
        stock_state = stock_state * np.exp(
            (BASE_PARAMS["r"] - 0.5 * variance_floor) * dt + diffusion * z_s
        )
        variance_state = variance_state + HESTON_PARAMS["kappa"] * (
            HESTON_PARAMS["theta"] - variance_floor
        ) * dt + HESTON_PARAMS["xi"] * diffusion * z_v
        variance_state = np.clip(variance_state, 0.0, None)

        if store_paths:
            stock_paths[:, step + 1] = stock_state
            variance_paths[:, step + 1] = variance_state

    if store_paths:
        return stock_paths, variance_paths
    return stock_state, variance_state


def merton_random_inputs(
    n_paths: int,
    n_steps: int,
    jump_intensity: float,
    seed: int,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    rng = np.random.default_rng(seed)
    dt = BASE_PARAMS["T"] / n_steps
    diffusion_shocks = rng.standard_normal((n_paths, n_steps))
    jump_counts = rng.poisson(jump_intensity * dt, size=(n_paths, n_steps))
    jump_normals = rng.standard_normal((n_paths, n_steps))
    return diffusion_shocks, jump_counts, jump_normals


def simulate_merton(
    S0: float,
    jump_intensity: float,
    diffusion_shocks: np.ndarray,
    jump_counts: np.ndarray,
    jump_normals: np.ndarray,
    store_paths: bool = False,
) -> np.ndarray:
    n_paths, n_steps = diffusion_shocks.shape
    dt = BASE_PARAMS["T"] / n_steps
    sqrt_dt = math.sqrt(dt)
    compensator = math.exp(MERTON_PARAMS["mu_j"] + 0.5 * MERTON_PARAMS["delta_j"] ** 2) - 1.0
    log_stock = np.full(n_paths, math.log(S0), dtype=float)

    if store_paths:
        stock_paths = np.empty((n_paths, n_steps + 1))
        stock_paths[:, 0] = S0

    for step in range(n_steps):
        counts = jump_counts[:, step]
        jump_sum = counts * MERTON_PARAMS["mu_j"] + np.sqrt(counts) * MERTON_PARAMS["delta_j"] * jump_normals[:, step]
        log_stock = log_stock + (
            (BASE_PARAMS["r"] - jump_intensity * compensator - 0.5 * BASE_PARAMS["sigma"] ** 2) * dt
            + BASE_PARAMS["sigma"] * sqrt_dt * diffusion_shocks[:, step]
            + jump_sum
        )
        if store_paths:
            stock_paths[:, step + 1] = np.exp(log_stock)

    if store_paths:
        return stock_paths
    return np.exp(log_stock)
In [2]:
general_inputs = pd.DataFrame(
    [
        ["Initial stock price S0", BASE_PARAMS["S0"]],
        ["Risk-free rate r", BASE_PARAMS["r"]],
        ["Diffusion volatility input sigma", BASE_PARAMS["sigma"]],
        ["Time to maturity T", BASE_PARAMS["T"]],
        ["ATM strike K", BASE_PARAMS["K_atm"]],
    ],
    columns=["General Input", "Value"],
)

heston_inputs = pd.DataFrame(
    [
        ["Initial variance v0", HESTON_PARAMS["v0"]],
        ["Mean reversion kappa_v", HESTON_PARAMS["kappa"]],
        ["Long-run variance theta_v", HESTON_PARAMS["theta"]],
        ["Volatility of variance", HESTON_PARAMS["xi"]],
    ],
    columns=["Heston Parameter", "Value"],
)

merton_inputs = pd.DataFrame(
    [
        ["Average jump size mu", MERTON_PARAMS["mu_j"]],
        ["Jump-size volatility delta", MERTON_PARAMS["delta_j"]],
    ],
    columns=["Merton Parameter", "Value"],
)

simulation_inputs = pd.DataFrame(
    [
        ["European-option paths", SIMULATION_SETTINGS["european_paths"]],
        ["Barrier-option paths", SIMULATION_SETTINGS["barrier_paths"]],
        ["American-option paths", SIMULATION_SETTINGS["american_paths"]],
        ["Time steps per path", SIMULATION_SETTINGS["steps"]],
        ["Spot bump for Greeks", SIMULATION_SETTINGS["bump_size"]],
    ],
    columns=["Simulation Choice", "Value"],
)

display(Markdown("## General Setup and Model Inputs"))
display(general_inputs)
display(heston_inputs)
display(merton_inputs)
display(simulation_inputs)
display(Markdown("*PS: the assignment does not specify the Heston volatility-of-variance parameter xi, so a value of 0.35 was assumed and used consistently throughout the Heston simulations.*"))

General Setup and Model Inputs¶

General Input Value
0 Initial stock price S0 80.000
1 Risk-free rate r 0.055
2 Diffusion volatility input sigma 0.350
3 Time to maturity T 0.250
4 ATM strike K 80.000
Heston Parameter Value
0 Initial variance v0 0.032
1 Mean reversion kappa_v 1.850
2 Long-run variance theta_v 0.045
3 Volatility of variance 0.350
Merton Parameter Value
0 Average jump size mu -0.50
1 Jump-size volatility delta 0.22
Simulation Choice Value
0 European-option paths 120000.0
1 Barrier-option paths 120000.0
2 American-option paths 90000.0
3 Time steps per path 100.0
4 Spot bump for Greeks 1.0

PS: the assignment does not specify the Heston volatility-of-variance parameter xi, so a value of 0.35 was assumed and used consistently throughout the Heston simulations.

Stochastic-Volatility Results¶

In [3]:
def heston_european_summary(
    rho: float,
    strike: float,
    n_paths: int,
    n_steps: int,
    seed: int,
    bump_size: float,
) -> dict[str, dict[str, float]]:
    z_price, z_independent = heston_random_inputs(n_paths=n_paths, n_steps=n_steps, seed=seed)
    discount = math.exp(-BASE_PARAMS["r"] * BASE_PARAMS["T"])

    spot_grid = {
        "down": BASE_PARAMS["S0"] - bump_size,
        "base": BASE_PARAMS["S0"],
        "up": BASE_PARAMS["S0"] + bump_size,
    }

    price_map: dict[str, dict[str, MCResult]] = {"call": {}, "put": {}}
    for label, spot in spot_grid.items():
        stock_terminal, _ = simulate_heston(
            S0=spot,
            rho=rho,
            z_price=z_price,
            z_independent=z_independent,
            store_paths=False,
        )
        call_payoff = discount * np.maximum(stock_terminal - strike, 0.0)
        put_payoff = discount * np.maximum(strike - stock_terminal, 0.0)
        price_map["call"][label] = _mc_summary(call_payoff)
        price_map["put"][label] = _mc_summary(put_payoff)

    output: dict[str, dict[str, float]] = {}
    for option_type in ("call", "put"):
        down_price = price_map[option_type]["down"].price
        base_price = price_map[option_type]["base"].price
        up_price = price_map[option_type]["up"].price
        delta = (up_price - down_price) / (2.0 * bump_size)
        gamma = (up_price - 2.0 * base_price + down_price) / (bump_size**2)
        output[option_type] = {
            "price": base_price,
            "stderr": price_map[option_type]["base"].stderr,
            "delta": delta,
            "gamma": gamma,
        }
    return output
In [4]:
q5 = heston_european_summary(
    rho=-0.30,
    strike=BASE_PARAMS["K_atm"],
    n_paths=SIMULATION_SETTINGS["european_paths"],
    n_steps=SIMULATION_SETTINGS["steps"],
    seed=101,
    bump_size=SIMULATION_SETTINGS["bump_size"],
)

q6 = heston_european_summary(
    rho=-0.70,
    strike=BASE_PARAMS["K_atm"],
    n_paths=SIMULATION_SETTINGS["european_paths"],
    n_steps=SIMULATION_SETTINGS["steps"],
    seed=202,
    bump_size=SIMULATION_SETTINGS["bump_size"],
)

q5_table = pd.DataFrame(
    [
        ["European Call", -0.30, round(q5["call"]["price"], 2), round(q5["call"]["stderr"], 4)],
        ["European Put", -0.30, round(q5["put"]["price"], 2), round(q5["put"]["stderr"], 4)],
    ],
    columns=["Question 5 Contract", "Correlation", "Price", "Std. Error"],
)

q6_table = pd.DataFrame(
    [
        ["European Call", -0.70, round(q6["call"]["price"], 2), round(q6["call"]["stderr"], 4)],
        ["European Put", -0.70, round(q6["put"]["price"], 2), round(q6["put"]["stderr"], 4)],
    ],
    columns=["Question 6 Contract", "Correlation", "Price", "Std. Error"],
)

display(Markdown("## Question 5"))
display(Markdown("ATM European call and put prices under the stochastic-volatility specification with correlation equal to -0.30."))
display(q5_table)
display(Markdown(
    f"With correlation fixed at -0.30, the ATM European call is priced at **${q5['call']['price']:.2f}** and the ATM European put is priced at **${q5['put']['price']:.2f}**. The call comes out above the put, which is reasonable in this setting because the stock still has room to participate in favorable moves under the risk-neutral measure, while the put continues to hold material value because time-varying volatility keeps downside risk relevant over the three-month horizon."
))

display(Markdown("## Question 6"))
display(Markdown("ATM European call and put prices under the stochastic-volatility specification with correlation equal to -0.70."))
display(q6_table)
display(Markdown(
    f"When the correlation is lowered to -0.70, the ATM European call is **${q6['call']['price']:.2f}** and the ATM European put is **${q6['put']['price']:.2f}**. Compared with Question 5, the put benefits slightly more from the stronger negative relationship between stock returns and variance shocks. In practical terms, this means that as the leverage effect becomes stronger, downside protection becomes a little more valuable."
))

Question 5¶

ATM European call and put prices under the stochastic-volatility specification with correlation equal to -0.30.

Question 5 Contract Correlation Price Std. Error
0 European Call -0.3 3.48 0.0135
1 European Put -0.3 2.38 0.0116

With correlation fixed at -0.30, the ATM European call is priced at $3.48 and the ATM European put is priced at $2.38. The call comes out above the put, which is reasonable in this setting because the stock still has room to participate in favorable moves under the risk-neutral measure, while the put continues to hold material value because time-varying volatility keeps downside risk relevant over the three-month horizon.

Question 6¶

ATM European call and put prices under the stochastic-volatility specification with correlation equal to -0.70.

Question 6 Contract Correlation Price Std. Error
0 European Call -0.7 3.51 0.0121
1 European Put -0.7 2.39 0.0124

When the correlation is lowered to -0.70, the ATM European call is $3.51 and the ATM European put is $2.39. Compared with Question 5, the put benefits slightly more from the stronger negative relationship between stock returns and variance shocks. In practical terms, this means that as the leverage effect becomes stronger, downside protection becomes a little more valuable.

Question 7: Delta and Gamma for Questions 5 and 6¶

Delta and gamma are reported below for the four stochastic-volatility contracts priced above. These measures show how the option values respond to changes in the stock price and how quickly that sensitivity changes.

In [5]:
q7_table = pd.DataFrame(
    [
        ["Q5", "European Call", -0.30, round(q5["call"]["delta"], 4), round(q5["call"]["gamma"], 4)],
        ["Q5", "European Put", -0.30, round(q5["put"]["delta"], 4), round(q5["put"]["gamma"], 4)],
        ["Q6", "European Call", -0.70, round(q6["call"]["delta"], 4), round(q6["call"]["gamma"], 4)],
        ["Q6", "European Put", -0.70, round(q6["put"]["delta"], 4), round(q6["put"]["gamma"], 4)],
    ],
    columns=["Source Question", "Contract", "Correlation", "Delta", "Gamma"],
)

display(q7_table)
display(Markdown(
    f"The call deltas are positive, at **{q5['call']['delta']:.4f}** for Question 5 and **{q6['call']['delta']:.4f}** for Question 6, which is consistent with calls gaining value when the stock rises. The put deltas are negative, at **{q5['put']['delta']:.4f}** and **{q6['put']['delta']:.4f}**, which is consistent with puts gaining value when the stock falls. All four gamma values are positive, meaning the contracts remain convex and their hedge ratios become more sensitive as the underlying price moves. In practice, this means a desk holding these options would need dynamic rebalancing, especially around the money where gamma is economically important."
))
Source Question Contract Correlation Delta Gamma
0 Q5 European Call -0.3 0.6029 0.0549
1 Q5 European Put -0.3 -0.3971 0.0549
2 Q6 European Call -0.7 0.6336 0.0503
3 Q6 European Put -0.7 -0.3668 0.0503

The call deltas are positive, at 0.6029 for Question 5 and 0.6336 for Question 6, which is consistent with calls gaining value when the stock rises. The put deltas are negative, at -0.3971 and -0.3668, which is consistent with puts gaining value when the stock falls. All four gamma values are positive, meaning the contracts remain convex and their hedge ratios become more sensitive as the underlying price moves. In practice, this means a desk holding these options would need dynamic rebalancing, especially around the money where gamma is economically important.

Questions 8 and 9: ATM Prices Under the Jump Model¶

Only the jump intensity changes between these two cases, so the comparison isolates the effect of jump frequency on option values.

In [6]:
def merton_european_summary(
    jump_intensity: float,
    strike: float,
    n_paths: int,
    n_steps: int,
    seed: int,
    bump_size: float,
) -> dict[str, dict[str, float]]:
    diffusion_shocks, jump_counts, jump_normals = merton_random_inputs(
        n_paths=n_paths,
        n_steps=n_steps,
        jump_intensity=jump_intensity,
        seed=seed,
    )
    discount = math.exp(-BASE_PARAMS["r"] * BASE_PARAMS["T"])

    spot_grid = {
        "down": BASE_PARAMS["S0"] - bump_size,
        "base": BASE_PARAMS["S0"],
        "up": BASE_PARAMS["S0"] + bump_size,
    }

    price_map: dict[str, dict[str, MCResult]] = {"call": {}, "put": {}}
    for label, spot in spot_grid.items():
        stock_terminal = simulate_merton(
            S0=spot,
            jump_intensity=jump_intensity,
            diffusion_shocks=diffusion_shocks,
            jump_counts=jump_counts,
            jump_normals=jump_normals,
            store_paths=False,
        )
        call_payoff = discount * np.maximum(stock_terminal - strike, 0.0)
        put_payoff = discount * np.maximum(strike - stock_terminal, 0.0)
        price_map["call"][label] = _mc_summary(call_payoff)
        price_map["put"][label] = _mc_summary(put_payoff)

    output: dict[str, dict[str, float]] = {}
    for option_type in ("call", "put"):
        down_price = price_map[option_type]["down"].price
        base_price = price_map[option_type]["base"].price
        up_price = price_map[option_type]["up"].price
        delta = (up_price - down_price) / (2.0 * bump_size)
        gamma = (up_price - 2.0 * base_price + down_price) / (bump_size**2)
        output[option_type] = {
            "price": base_price,
            "stderr": price_map[option_type]["base"].stderr,
            "delta": delta,
            "gamma": gamma,
        }
    return output
In [7]:
q8 = merton_european_summary(
    jump_intensity=0.75,
    strike=BASE_PARAMS["K_atm"],
    n_paths=SIMULATION_SETTINGS["european_paths"],
    n_steps=SIMULATION_SETTINGS["steps"],
    seed=303,
    bump_size=SIMULATION_SETTINGS["bump_size"],
)

q9 = merton_european_summary(
    jump_intensity=0.25,
    strike=BASE_PARAMS["K_atm"],
    n_paths=SIMULATION_SETTINGS["european_paths"],
    n_steps=SIMULATION_SETTINGS["steps"],
    seed=404,
    bump_size=SIMULATION_SETTINGS["bump_size"],
)

q8_table = pd.DataFrame(
    [
        ["European Call", 0.75, round(q8["call"]["price"], 2), round(q8["call"]["stderr"], 4)],
        ["European Put", 0.75, round(q8["put"]["price"], 2), round(q8["put"]["stderr"], 4)],
    ],
    columns=["Question 8 Contract", "Jump Intensity", "Price", "Std. Error"],
)

q9_table = pd.DataFrame(
    [
        ["European Call", 0.25, round(q9["call"]["price"], 2), round(q9["call"]["stderr"], 4)],
        ["European Put", 0.25, round(q9["put"]["price"], 2), round(q9["put"]["stderr"], 4)],
    ],
    columns=["Question 9 Contract", "Jump Intensity", "Price", "Std. Error"],
)

display(Markdown("## Question 8"))
display(Markdown("ATM European call and put prices under the jump model with jump intensity equal to 0.75."))
display(q8_table)
display(Markdown(
    f"With jump intensity set to 0.75, the ATM call price is **${q8['call']['price']:.2f}** and the ATM put price is **${q8['put']['price']:.2f}**. Because the average jump size is negative in this specification, the higher frequency of jumps increases the importance of downside outcomes and helps keep put values elevated."
))

display(Markdown("## Question 9"))
display(Markdown("ATM European call and put prices under the jump model with jump intensity equal to 0.25."))
display(q9_table)
display(Markdown(
    f"With jump intensity reduced to 0.25, the ATM call price is **${q9['call']['price']:.2f}** and the ATM put price is **${q9['put']['price']:.2f}**. Both contracts become cheaper than in Question 8 because the market is now exposed to fewer jump arrivals over the same maturity. For an investor, the main takeaway is that jump frequency matters: when sudden discontinuous moves are more likely, option insurance becomes more expensive."
))

Question 8¶

ATM European call and put prices under the jump model with jump intensity equal to 0.75.

Question 8 Contract Jump Intensity Price Std. Error
0 European Call 0.75 8.28 0.0332
1 European Put 0.75 7.23 0.0355

With jump intensity set to 0.75, the ATM call price is $8.28 and the ATM put price is $7.23. Because the average jump size is negative in this specification, the higher frequency of jumps increases the importance of downside outcomes and helps keep put values elevated.

Question 9¶

ATM European call and put prices under the jump model with jump intensity equal to 0.25.

Question 9 Contract Jump Intensity Price Std. Error
0 European Call 0.25 6.79 0.0290
1 European Put 0.25 5.76 0.0269

With jump intensity reduced to 0.25, the ATM call price is $6.79 and the ATM put price is $5.76. Both contracts become cheaper than in Question 8 because the market is now exposed to fewer jump arrivals over the same maturity. For an investor, the main takeaway is that jump frequency matters: when sudden discontinuous moves are more likely, option insurance becomes more expensive.

Question 10: Delta and Gamma for Questions 8 and 9¶

The Greek estimates below show the first-order and second-order sensitivity of the jump-model option prices to movements in the underlying stock price.

In [8]:
q10_table = pd.DataFrame(
    [
        ["Q8", "European Call", 0.75, round(q8["call"]["delta"], 4), round(q8["call"]["gamma"], 4)],
        ["Q8", "European Put", 0.75, round(q8["put"]["delta"], 4), round(q8["put"]["gamma"], 4)],
        ["Q9", "European Call", 0.25, round(q9["call"]["delta"], 4), round(q9["call"]["gamma"], 4)],
        ["Q9", "European Put", 0.25, round(q9["put"]["delta"], 4), round(q9["put"]["gamma"], 4)],
    ],
    columns=["Source Question", "Contract", "Jump Intensity", "Delta", "Gamma"],
)

display(q10_table)
display(Markdown(
    f"The jump-model call deltas are positive and the put deltas are negative, which means the directional behavior of the contracts remains economically sensible even when returns can jump. Gamma stays positive in all four cases, showing that both calls and puts remain convex. Comparing Questions 8 and 9, the higher jump-intensity case produces a call delta of **{q8['call']['delta']:.4f}** versus **{q9['call']['delta']:.4f}** in the lower-intensity case. This shows that the hedge profile changes when jump arrivals become more frequent, so a hedge calibrated in a quieter market may not remain accurate when jump risk rises."
))
Source Question Contract Jump Intensity Delta Gamma
0 Q8 European Call 0.75 0.6453 0.0224
1 Q8 European Put 0.75 -0.3542 0.0224
2 Q9 European Call 0.25 0.5958 0.0269
3 Q9 European Put 0.25 -0.4034 0.0269

The jump-model call deltas are positive and the put deltas are negative, which means the directional behavior of the contracts remains economically sensible even when returns can jump. Gamma stays positive in all four cases, showing that both calls and puts remain convex. Comparing Questions 8 and 9, the higher jump-intensity case produces a call delta of 0.6453 versus 0.5958 in the lower-intensity case. This shows that the hedge profile changes when jump arrivals become more frequent, so a hedge calibrated in a quieter market may not remain accurate when jump risk rises.

Question 13: American Call Under the Stochastic-Volatility Setting¶

For a non-dividend-paying stock, early exercise of a call is not optimal, so the American and European values should match in economic terms.

In [9]:
def longstaff_schwartz_call(
    stock_paths: np.ndarray,
    variance_paths: np.ndarray,
    strike: float,
) -> float:
    n_steps = stock_paths.shape[1] - 1
    dt = BASE_PARAMS["T"] / n_steps
    discount_step = math.exp(-BASE_PARAMS["r"] * dt)

    exercise_value = np.maximum(stock_paths - strike, 0.0)
    continuation_value = exercise_value[:, -1].copy()

    for step in range(n_steps - 1, 0, -1):
        continuation_value = discount_step * continuation_value
        itm = exercise_value[:, step] > 0.0
        if itm.sum() < 10:
            continue

        x = _basis_functions(stock_paths[itm, step], variance_paths[itm, step], strike)
        y = continuation_value[itm]
        beta, *_ = np.linalg.lstsq(x, y, rcond=None)
        continuation_estimate = np.maximum(x @ beta, 0.0)
        immediate = exercise_value[itm, step]
        exercise_now = immediate > continuation_estimate

        updated = continuation_value[itm].copy()
        updated[exercise_now] = immediate[exercise_now]
        continuation_value[itm] = updated

    initial_exercise = max(BASE_PARAMS["S0"] - strike, 0.0)
    return float(max(initial_exercise, discount_step * continuation_value.mean()))


def heston_american_call_summary(
    rho: float,
    strike: float,
    n_paths: int,
    n_steps: int,
    seed: int,
    bump_size: float,
) -> dict[str, float]:
    z_price, z_independent = heston_random_inputs(n_paths=n_paths, n_steps=n_steps, seed=seed)

    price_grid = {}
    for label, spot in {
        "down": BASE_PARAMS["S0"] - bump_size,
        "base": BASE_PARAMS["S0"],
        "up": BASE_PARAMS["S0"] + bump_size,
    }.items():
        stock_paths, variance_paths = simulate_heston(
            S0=spot,
            rho=rho,
            z_price=z_price,
            z_independent=z_independent,
            store_paths=True,
        )
        price_grid[label] = longstaff_schwartz_call(stock_paths, variance_paths, strike)

    delta = (price_grid["up"] - price_grid["down"]) / (2.0 * bump_size)
    gamma = (price_grid["up"] - 2.0 * price_grid["base"] + price_grid["down"]) / (bump_size**2)

    return {
        "price": price_grid["base"],
        "delta": delta,
        "gamma": gamma,
    }
In [10]:
q13_check = heston_american_call_summary(
    rho=-0.30,
    strike=BASE_PARAMS["K_atm"],
    n_paths=SIMULATION_SETTINGS["american_paths"],
    n_steps=SIMULATION_SETTINGS["steps"],
    seed=505,
    bump_size=SIMULATION_SETTINGS["bump_size"],
)

q13 = q13_check

q13_table = pd.DataFrame(
    [
        ["American Call", round(q13["price"], 2), round(q13["delta"], 4), round(q13["gamma"], 4)],
        ["European Call from Question 5", round(q5["call"]["price"], 2), round(q5["call"]["delta"], 4), round(q5["call"]["gamma"], 4)],
    ],
    columns=["Contract", "Price", "Delta", "Gamma"],
)

display(q13_table)
display(Markdown(
    f"The American-call estimate is **${q13['price']:.2f}**, with delta **{q13['delta']:.4f}** and gamma **{q13['gamma']:.4f}**. This is numerically very close to the European call from Question 5, which is consistent with the no-dividend early-exercise result. In practical terms, the American feature adds little value here, so the two contracts end up looking almost the same."
))
Contract Price Delta Gamma
0 American Call 3.43 0.6058 0.0515
1 European Call from Question 5 3.48 0.6029 0.0549

The American-call estimate is $3.43, with delta 0.6058 and gamma 0.0515. This is numerically very close to the European call from Question 5, which is consistent with the no-dividend early-exercise result. In practical terms, the American feature adds little value here, so the two contracts end up looking almost the same.

Questions 14 and 15: Barrier Options¶

These contracts are path-dependent, so the barrier event matters just as much as the terminal payoff.

In [11]:
def heston_barrier_summary(
    rho: float,
    strike: float,
    barrier: float,
    n_paths: int,
    n_steps: int,
    seed: int,
) -> dict[str, float]:
    z_price, z_independent = heston_random_inputs(n_paths=n_paths, n_steps=n_steps, seed=seed)
    stock_paths, _ = simulate_heston(
        S0=BASE_PARAMS["S0"],
        rho=rho,
        z_price=z_price,
        z_independent=z_independent,
        store_paths=True,
    )
    discount = math.exp(-BASE_PARAMS["r"] * BASE_PARAMS["T"])
    hit_barrier = (stock_paths >= barrier).any(axis=1)
    terminal_stock = stock_paths[:, -1]
    vanilla_payoff = discount * np.maximum(terminal_stock - strike, 0.0)
    knock_in_payoff = vanilla_payoff * hit_barrier
    return {
        "vanilla_price": float(vanilla_payoff.mean()),
        "barrier_price": float(knock_in_payoff.mean()),
        "activation_rate": float(hit_barrier.mean()),
    }


def merton_barrier_summary(
    jump_intensity: float,
    strike: float,
    barrier: float,
    n_paths: int,
    n_steps: int,
    seed: int,
) -> dict[str, float]:
    diffusion_shocks, jump_counts, jump_normals = merton_random_inputs(
        n_paths=n_paths,
        n_steps=n_steps,
        jump_intensity=jump_intensity,
        seed=seed,
    )
    stock_paths = simulate_merton(
        S0=BASE_PARAMS["S0"],
        jump_intensity=jump_intensity,
        diffusion_shocks=diffusion_shocks,
        jump_counts=jump_counts,
        jump_normals=jump_normals,
        store_paths=True,
    )
    discount = math.exp(-BASE_PARAMS["r"] * BASE_PARAMS["T"])
    hit_barrier = (stock_paths <= barrier).any(axis=1)
    terminal_stock = stock_paths[:, -1]
    vanilla_payoff = discount * np.maximum(strike - terminal_stock, 0.0)
    knock_in_payoff = vanilla_payoff * hit_barrier
    return {
        "vanilla_price": float(vanilla_payoff.mean()),
        "barrier_price": float(knock_in_payoff.mean()),
        "activation_rate": float(hit_barrier.mean()),
    }
In [12]:
q14 = heston_barrier_summary(
    rho=-0.70,
    strike=95.0,
    barrier=95.0,
    n_paths=SIMULATION_SETTINGS["barrier_paths"],
    n_steps=SIMULATION_SETTINGS["steps"],
    seed=606,
)

q15 = merton_barrier_summary(
    jump_intensity=0.75,
    strike=65.0,
    barrier=65.0,
    n_paths=SIMULATION_SETTINGS["barrier_paths"],
    n_steps=SIMULATION_SETTINGS["steps"],
    seed=707,
)

q14_table = pd.DataFrame(
    [
        ["Up-and-In Call", 95, 95, round(q14["barrier_price"], 2), round(q14["vanilla_price"], 2), round(q14["activation_rate"], 4)],
    ],
    columns=["Question 14 Contract", "Barrier", "Strike", "Barrier Price", "Vanilla Price", "Activation Rate"],
)

q15_table = pd.DataFrame(
    [
        ["Down-and-In Put", 65, 65, round(q15["barrier_price"], 2), round(q15["vanilla_price"], 2), round(q15["activation_rate"], 4)],
    ],
    columns=["Question 15 Contract", "Barrier", "Strike", "Barrier Price", "Vanilla Price", "Activation Rate"],
)

display(Markdown("## Question 14"))
display(q14_table)
display(Markdown(
    f"The up-and-in call price is **${q14['barrier_price']:.2f}**, and the plain European call with the same strike is **${q14['vanilla_price']:.2f}**. In this setup, those two values match to rounding because the barrier level and strike are both 95. If the stock finishes above 95 at maturity, the path must have crossed 95 at some earlier point, so the barrier condition is automatically satisfied whenever the vanilla payoff is positive."
))

display(Markdown("## Question 15"))
display(q15_table)
display(Markdown(
    f"The down-and-in put price is **${q15['barrier_price']:.2f}**, while the corresponding plain European put with strike 65 also prices at **${q15['vanilla_price']:.2f}**. Again, this happens because the barrier and strike are the same. A put with strike 65 only finishes in-the-money if the stock ends below 65, and that necessarily means the path has crossed the barrier. These two examples show how strongly barrier design affects value: once the barrier is set equal to the strike, the knock-in feature can collapse to the vanilla payoff."
))

Question 14¶

Question 14 Contract Barrier Strike Barrier Price Vanilla Price Activation Rate
0 Up-and-In Call 95 95 0.03 0.03 0.0217

The up-and-in call price is $0.03, and the plain European call with the same strike is $0.03. In this setup, those two values match to rounding because the barrier level and strike are both 95. If the stock finishes above 95 at maturity, the path must have crossed 95 at some earlier point, so the barrier condition is automatically satisfied whenever the vanilla payoff is positive.

Question 15¶

Question 15 Contract Barrier Strike Barrier Price Vanilla Price Activation Rate
0 Down-and-In Put 65 65 2.79 2.79 0.2585

The down-and-in put price is $2.79, while the corresponding plain European put with strike 65 also prices at $2.79. Again, this happens because the barrier and strike are the same. A put with strike 65 only finishes in-the-money if the stock ends below 65, and that necessarily means the path has crossed the barrier. These two examples show how strongly barrier design affects value: once the barrier is set equal to the strike, the knock-in feature can collapse to the vanilla payoff.

Conclusion¶

The results show that the source of volatility matters just as much as the level of volatility itself. In the stochastic-volatility setting, making the correlation more negative slightly increased the value of downside protection, which is why the put responded a little more strongly than the call. In the jump setting, increasing jump intensity raised both option prices because the market was exposed more often to large discontinuous moves.

The later questions reinforced the same idea from a different angle. The American call did not gain meaningful extra value relative to the European call, which is consistent with the no-dividend setting. The barrier examples also showed that contract design matters: once the barrier was set equal to the strike, the knock-in payoff effectively collapsed to the corresponding vanilla payoff.

Overall, the main takeaway is that option values are shaped not only by how much uncertainty exists, but also by the form that uncertainty takes. For pricing, hedging, and portfolio protection, it is therefore important to distinguish between changing volatility, jump risk, and path-dependent contract features rather than treating all risk as if it came from a single constant-volatility source.

References¶

  • Heston, Steven L. "A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options." Review of Financial Studies, vol. 6, no. 2, 1993, pp. 327-343.
  • Longstaff, Francis A., and Eduardo S. Schwartz. "Valuing American Options by Simulation: A Simple Least-Squares Approach." Review of Financial Studies, vol. 14, no. 1, 2001, pp. 113-147.
  • Merton, Robert C. "Option Pricing When Underlying Stock Returns Are Discontinuous." Journal of Financial Economics, vol. 3, no. 1-2, 1976, pp. 125-144.