Introduction

This project investigate two primary drivers of volatility beyond constant volatility models: stochastic volatility (via the Heston Model) and jump risk (via the Merton Model). Specifically, we (1) price ATM European calls and puts under the Heston Model using two different correlations ($\rho = –0.30$ and $\rho = –0.70$) and compute their deltas and gammas, (2) price ATM European calls and puts under the Merton jump‐diffusion model for two jump intensities ($\lambda = 0.75$ and $\lambda = 0.25$) and compute their deltas and gammas, (3) verify put–call parity for all four European‐option cases, (4) generate and compare Heston and Merton based call prices across seven equally spaced strikes (moneyness = {0.85, 0.90, 0.95, 1.00, 1.05, 1.10, 1.15}), (5) price an American call under the Heston Model via Longstaff–Schwartz and compare its price and Greeks to its European counterpart, and (6) price two barrier options (a Heston based up and in call and a Merton based down and in put) and compare each to its plain vanilla equivalent. By trading off stochastic variance and jump components, we isolate how model features mean reversion, leverage effect, jump amplitude, jump frequency alter option payoffs, sensitivities, and arbitrage relationships.

Team member A

Question 5

To solve Question 5, we price ATM European call and put options using the Heston model with Monte Carlo simulation, correlation $\rho = -0.30$, and the given parameters. The Heston model dynamics are given by two stochastic differential equations: $dS_t = rS_t\,dt + \sqrt{v_t}\,S_t\,dW_t^S \quad\text{and}\quad dv_t = \kappa(\theta - v_t)\,dt + \sigma_v\sqrt{v_t}\,dW_t^v,$ where the Brownian motions satisfy $dW_t^S \cdot dW_t^v = \rho\,dt$ (Heston 330). We implement a Monte Carlo simulation with 100 time steps (representing daily intervals over three months) and 500,000 paths to ensure high accuracy. The simulation uses the full‐truncation scheme to keep variance nonnegative by applying $\max(v_t,0)$ during updates (Andersen and Brotherton‐Ratcliffe 78), and it employs a log‐transformation to simulate $\ln(S_t)$, preventing negative stock prices (Kahl and Jäckel 115). Correlated Brownian motions are generated via Cholesky decomposition. The model parameters are: $S_0 = 80$, $K = 80$ (ATM), $r = 0.055$, $T = 0.25$, $v_0 = 0.032$, $\kappa = 1.85$, $\theta = 0.045$, $\sigma_v = 0.35$, and $\rho = -0.30$. The European call and put option prices are obtained from the discounted expected payoffs, using $e^{-rT}\max(S_T - K,0)$ and $e^{-rT}\max(K - S_T,0)$, respectively (Glasserman 345; Heston 332), that's option Pricing: $\text{Call Price} = e^{-rT} \mathbb{E}[\max(S_T - K, 0)]$ and $\text{Put Price} = e^{-rT} \mathbb{E}[\max(K - S_T, 0)]$, while for discretization, we used full-truncation scheme for variance $v_t = \max(v_t, 0)$ (Andersen and Brotherton-Ratcliffe 78) and log-transformation for $S_t$ (Kahl and Jäckel 115)

import numpy as np
np.random.seed(42)

# Parameters
S0 = 80.0
K = 80.0
r = 0.055
T = 0.25
v0 = 0.032
kappa = 1.85
theta = 0.045
sigma_v = 0.35
rho = -0.30
n_simulations = 500000
n_steps = 100
dt = T / n_steps

# Correlated Brownian motions
np.random.seed(42)
dW = np.random.normal(size=(n_simulations, n_steps, 2))
dW_v = dW[:, :, 0]
dW_S = rho * dW[:, :, 0] + np.sqrt(1 - rho**2) * dW[:, :, 1]

v = np.full(n_simulations, v0)
lnS = np.full(n_simulations, np.log(S0))

# Euler-Maruyama simulation
for t in range(n_steps):
    v_positive = np.maximum(v, 0)
    sqrt_v = np.sqrt(v_positive)

    # Update variance
    v += kappa * (theta - v_positive) * dt + sigma_v * sqrt_v * np.sqrt(dt) * dW_v[:, t]

    # Update log stock price
    lnS += (r - 0.5 * v_positive) * dt + sqrt_v * np.sqrt(dt) * dW_S[:, t]

# Terminal stock price
S_T = np.exp(lnS)

# Payoffs
call_payoff = np.maximum(S_T - K, 0)
put_payoff = np.maximum(K - S_T, 0)

# Discounted prices
call_price = np.exp(-r * T) * np.mean(call_payoff)
put_price = np.exp(-r * T) * np.mean(put_payoff)

call_price_rounded = round(call_price, 2)
put_price_rounded = round(put_price, 2)

print(f"ATM European Call Price: ${call_price_rounded:.2f}")
print(f"ATM European Put Price: ${put_price_rounded:.2f}")

ATM European Call Price: $3.46
ATM European Put Price: $2.39

The simulated Heston Model pricing with $\rho = -0.30$ yields an ATM European call price of $3.46 and put price of $2.39, reflecting the impact of negative correlation between stock returns and volatility. This inverse relationship amplifies downside movements through higher volatility when prices fall, boosting the put option's value beyond what would occur under positive correlation. Simultaneously, the call option is suppressed as price rallies correspond with lower volatility, limiting upside potential. The call's premium over the put ($1.07) aligns with the positive drift from the 5.5% risk-free rate, while maintaining approximate put-call parity with the theoretical difference of $1.09, confirming model consistency.

Question 6

For this question, all parameters, formulas and simulations are identical to Question 5, except that the correlation ($\rho$) is –0.70 and the number of time steps is 100.

import numpy as np

np.random.seed(42)
# Parameters (updated ρ = -0.70)
S0 = 80.0
K = 80.0
r = 0.055
T = 0.25
v0 = 0.032
kappa = 1.85
theta = 0.045
sigma_v = 0.35
rho = -0.70
n_simulations = 500000
n_steps = 100
dt = T / n_steps

# Correlated Brownian motions
np.random.seed(42)
dW = np.random.normal(size=(n_simulations, n_steps, 2))
dW_v = dW[:, :, 0]
dW_S = rho * dW[:, :, 0] + np.sqrt(1 - rho**2) * dW[:, :, 1]

# Initialize arrays
v = np.full(n_simulations, v0)
lnS = np.full(n_simulations, np.log(S0))

for t in range(n_steps):
    v_positive = np.maximum(v, 0)
    sqrt_v = np.sqrt(v_positive)

    # Update variance (full truncation)
    v += kappa * (theta - v_positive) * dt + sigma_v * sqrt_v * np.sqrt(dt) * dW_v[:, t]

    # Update log stock price
    lnS += (r - 0.5 * v_positive) * dt + sqrt_v * np.sqrt(dt) * dW_S[:, t]

# Terminal prices
S_T = np.exp(lnS)

# Payoffs and pricing
call_price = np.exp(-r*T) * np.mean(np.maximum(S_T - K, 0))
put_price = np.exp(-r*T) * np.mean(np.maximum(K - S_T, 0))


print(f"ATM European Call (ρ=-0.70): ${round(call_price, 2):.2f}")
print(f"ATM European Put (ρ=-0.70): ${round(put_price, 2):.2f}")

ATM European Call (ρ=-0.70): $3.49
ATM European Put (ρ=-0.70): $2.41

The simulated Heston Model pricing with $\rho = -0.70$ yields an ATM European call price of $3.49 and put price of $2.41, showing a muted response to the stronger negative correlation compared to the $\rho = -0.30$ case ($3.46 call, $2.39 put). This limited sensitivity suggests that while the increased inverse relationship between stock returns and volatility should theoretically amplify the leverage effect, the specific parameterization (particularly the volatility mean-reversion at $κ = 1.85$ and initial variance $v_0 = 3.2\%$ near long-term $\theta = 4.5\%$) dampens extreme volatility reactions. The call-put spread remains nearly identical at $1.08 versus $1.07 previously, indicating that the positive drift component dominates option valuations more than correlation effects in this configuration. The put-call parity holds well (3.49 - 2.41 = 1.08 vs. theoretical 1.09), confirming model consistency despite the parameter change.

Question 7

To calculate delta and gamma for the European call and put options from Questions 5 ($\rho = -0.30$) and 6 ($\rho = -0.70$) using the central difference method. We used a stock price perturbation of $dS = 0.01\times S0 = 0.80.$, stock price: $S_0 = \$80.00$, strike price (K): $80.00 (ATM), risk-free rate (r): 5.5% (0.055) and time to maturity (T): 0.25 years

Delta and Gamma formulas

• Delta ($\Delta$) $$\Delta \approx \frac{V(S_0 + dS) - V(S_0 - dS)}{2 \cdot dS}$$

• Gamma ($\Gamma$) $$\Gamma \approx \frac{V(S_0 + dS) - 2V(S_0) + V(S_0 - dS)}{dS^2}$$

import numpy as np
import pandas as pd

def heston_price(S0, K, T, r, v0, kappa_v, theta_v, sigma_v, rho, opt_type='call', num_sims=100000, num_steps=100):
    """
    Calculates European option price using Heston model via Monte Carlo simulation
    """
    dt = T / num_steps

    # Generate correlated Brownian motions
    dW_S = np.random.normal(0, np.sqrt(dt), (num_sims, num_steps))
    dZ = np.random.normal(0, np.sqrt(dt), (num_sims, num_steps))
    dW_v = rho * dW_S + np.sqrt(1 - rho**2) * dZ

    S = np.zeros((num_sims, num_steps + 1))
    v = np.zeros((num_sims, num_steps + 1))
    S[:, 0] = S0
    v[:, 0] = v0

    # Euler-Maruyama discretization
    for t in range(num_steps):
        v_current = np.maximum(v[:, t], 0)
        v_next = v[:, t] + kappa_v * (theta_v - v_current) * dt + sigma_v * np.sqrt(v_current) * dW_v[:, t]
        v[:, t + 1] = np.maximum(v_next, 0)
        S[:, t + 1] = S[:, t] * (1 + r * dt + np.sqrt(v_current) * dW_S[:, t])

    # Calculate payoff
    S_T = S[:, -1]
    if opt_type == 'call':
        payoff = np.maximum(S_T - K, 0)
    else:  # put
        payoff = np.maximum(K - S_T, 0)

    price = np.exp(-r * T) * np.mean(payoff)
    return round(price, 2)

def calculate_greeks():
    S0 = 80.0
    K = 80.0
    T = 0.25
    r = 0.055
    v0 = 0.032
    kappa_v = 1.85
    theta_v = 0.045
    sigma_v = 0.35
    dS = 0.80  # 1% of S0

    results = {}

    # For each correlation case
    for rho, case in zip([-0.30, -0.70], ['Q5', 'Q6']):
        np.random.seed(42 if case == 'Q5' else 123)

        # Base prices
        base_call = heston_price(S0, K, T, r, v0, kappa_v, theta_v, sigma_v, rho, 'call')
        base_put = heston_price(S0, K, T, r, v0, kappa_v, theta_v, sigma_v, rho, 'put')

        # Perturbed prices
        np.random.seed(42 if case == 'Q5' else 123)  # same seed path
        call_plus = heston_price(S0 + dS, K, T, r, v0, kappa_v, theta_v, sigma_v, rho, 'call')
        call_minus = heston_price(S0 - dS, K, T, r, v0, kappa_v, theta_v, sigma_v, rho, 'call')

        np.random.seed(42 if case == 'Q5' else 123)  # same seed path
        put_plus = heston_price(S0 + dS, K, T, r, v0, kappa_v, theta_v, sigma_v, rho, 'put')
        put_minus = heston_price(S0 - dS, K, T, r, v0, kappa_v, theta_v, sigma_v, rho, 'put')

        #  deltas
        call_delta = (call_plus - call_minus) / (2 * dS)
        put_delta = (put_plus - put_minus) / (2 * dS)

        #  gammas
        call_gamma = (call_plus - 2 * base_call + call_minus) / (dS ** 2)
        put_gamma = (put_plus - 2 * base_put + put_minus) / (dS ** 2)

        results[case] = {
            'call': {'price': base_call, 'delta': call_delta, 'gamma': call_gamma},
            'put': {'price': base_put, 'delta': put_delta, 'gamma': put_gamma}
        }

    return results


greek_results = calculate_greeks()

def results_to_dataframe(results):
    data = []
    for question, opt_types in results.items():
        for opt_type, metrics in opt_types.items():
            row = {
                'Question': question,
                'Option Type': opt_type.capitalize(),
                'Price': metrics['price'],
                'Delta': metrics['delta'],
                'Gamma': metrics['gamma']
            }
            data.append(row)
    df = pd.DataFrame(data)
    return df


df_greeks = results_to_dataframe(greek_results)
print("Table 1: Delta and gamma for the European call and put options")
print(df_greeks)

Table 1: Delta and gamma for the European call and put options
  Question Option Type  Price    Delta     Gamma
0       Q5        Call   3.45  0.60625  0.046875
1       Q5         Put   2.38 -0.41250  0.031250
2       Q6        Call   3.52  0.65000  0.031250
3       Q6         Put   2.39 -0.35625  0.078125

For $\rho = -0.30$, the call option with a price of $3.45 shows a delta of 0.60625, indicating it gains value as the stock price rises, while its gamma of 0.046875 reflects moderate sensitivity of delta to stock price changes. The put option, priced at $2.38, has a delta of -0.41250, confirming its inverse relationship with the stock price, and a gamma of 0.031250, suggesting less convexity than the call (Table 1).

For $\rho = -0.70$, the call option's price increases slightly to $3.52, with a higher delta of 0.65000, showing greater responsiveness to stock movements, while its gamma decreases to 0.031250, indicating reduced convexity. The put option's price remains stable at $2.39, but its delta becomes less negative (-0.35625), implying reduced sensitivity to downside moves, while its gamma more than doubles to 0.078125, highlighting increased convexity near the strike. These results demonstrate that more negative correlation ($\rho = -0.70$) leads to higher call deltas and put gammas, reflecting the asymmetric impact of correlation on option sensitivities calls become more responsive to upside moves while puts develop stronger convexity near expiration, consistent with the Heston model's volatility clustering behavior where negative correlation amplifies downside volatility spikes. The gamma values particularly show how correlation changes affect the curvature of option prices, with puts under $\rho = -0.70$ exhibiting nearly $2.5 \times$ the gamma of those under $\rho = -0.30$, a critical insight for volatility traders managing gamma exposure.

Team member B

Question 8

The Merton jump-diffusion model describes the stock price as: $\frac{dS_t}{S_{t-}} = (\mu - \lambda k)dt + \sigma dW_t + dJ_t$, where $\mu$ is the drift rate (risk-free rate $r = 5.5\%$), $\sigma = 35\%$ is the diffusion volatility, $J_t$ is a compound Poisson process (Merton 128) with intensity $\lambda = 0.75$.

Jump sizes are lognormal: $\ln(1 + J) \sim \mathcal{N}(\mu, \delta^2)\Longrightarrow \ln(1 + J) \sim \mathcal{N}(\mu = -0.5, \delta = 0.22$, and the terminal stock price is $S_T = S_0 \exp\left(\left(r - \frac{\sigma^2}{2} - \lambda \kappa\right)T + \sigma \sqrt{T}Z + \sum_{i=1}^{N_T} (\mu + \delta Z_i)\right)$, where $\kappa = e^{\mu + \delta^2/2} - 1$, drift adjustment: $r - \frac{\sigma^2}{2} - \lambda \kappa$ and terminal variance $\text{Var} = \sigma^2 T + \lambda T(\mu^2 + \delta^2)$ (Merton 130). For pricing we used same payoff expectations as Heston.

import numpy as np
import pandas as pd

np.random.seed(42)
def merton_jump_price(S0, K, T, r, sigma, lam, mu_j, delta_j, opt_type, num_sims=100000):
    """
    Prices European options using the Merton jump-diffusion model via Monte Carlo simulation.
    """
    # Calculate jump compensator κ
    kappa = np.exp(mu_j + 0.5 * delta_j**2) - 1

    # Adjust drift for risk-neutral measure
    drift_adj = r - 0.5 * sigma**2 - lam * kappa
    total_drift = drift_adj * T
    sigma_sqrtT = sigma * np.sqrt(T)

    np.random.seed(42)
    poisson_jumps = np.random.poisson(lam * T, num_sims)
    Z1 = np.random.normal(0, 1, num_sims)  # Diffusion shock
    Z2 = np.random.normal(0, 1, num_sims)  # Jump size shock

    jump_part = np.zeros(num_sims)
    non_zero_mask = poisson_jumps > 0
    jump_part[non_zero_mask] = (
        mu_j * poisson_jumps[non_zero_mask] +
        delta_j * np.sqrt(poisson_jumps[non_zero_mask]) * Z2[non_zero_mask]
    )

    #  total log return
    log_return = total_drift + sigma_sqrtT * Z1 + jump_part

    #  terminal stock prices
    S_T = S0 * np.exp(log_return)

    #  payoffs
    if opt_type == 'call':
        payoff = np.maximum(S_T - K, 0)
    else:  # put
        payoff = np.maximum(K - S_T, 0)

    # Discount payoffs and average
    price = np.exp(-r * T) * np.mean(payoff)
    return round(price, 2)

# Parameters for Question 8
params = {
    'S0': 80.0,       # Spot price
    'K': 80.0,        # ATM strike
    'T': 0.25,        # 3 months
    'r': 0.055,       # Risk-free rate
    'sigma': 0.35,    # Volatility
    'lam': 0.75,      # Jump intensity
    'mu_j': -0.5,     # Mean jump size
    'delta_j': 0.22   # Jump size volatility
}

# Price call and put options
call_price = merton_jump_price(opt_type='call', **params)
put_price = merton_jump_price(opt_type='put', **params)

results = pd.DataFrame({
    'Question': ['8', '8'],
    'Option Type': ['Call', 'Put'],
    'Price': [call_price, put_price]
})

print("Table 2: Merton Model with λ = 0.75")
print(results.to_string(index=False))

Table 2: Merton Model with λ = 0.75
Question Option Type  Price
       8        Call   8.30
       8         Put   7.19

From the Table 2 we see that the Merton model results ($\lambda=0.75$) show an $8.30 call and $7.19 put price, revealing how jump risk distorts option valuations. The near-parity pricing (just $1.11 difference) signals strong crash pricing from negative jumps ($\mu=-0.5$), inflating the put value while the call maintains worth through diffusion volatility. This creates a volatility smile both options trade rich to Black-Scholes, with puts disproportionately expensive due to their direct exposure to downward jumps. The high jump intensity ($\lambda=0.75$) amplifies these effects, making both options price in substantial tail risk.

Question 9

We'll price ATM European call and put options using the Merton jump-diffusion model with reduced jump intensity ($\lambda = 0.25$). This represents a market environment with less frequent jumps compared to Question 8. Note that we used same formulas as Question 8, except $\lambda = 0.25$.

import numpy as np
import pandas as pd

np.random.seed(42)
def merton_jump_price(S0, K, T, r, sigma, lam, mu_j, delta_j, opt_type, num_sims=100000):
    """Prices European options under Merton jump-diffusion model"""
    # Calculate jump compensator κ
    kappa = np.exp(mu_j + 0.5 * delta_j**2) - 1

    # Adjust drift for risk-neutral measure
    drift_adj = r - 0.5 * sigma**2 - lam * kappa
    total_drift = drift_adj * T
    sigma_sqrtT = sigma * np.sqrt(T)

    np.random.seed(42)
    poisson_jumps = np.random.poisson(lam * T, num_sims)
    Z1 = np.random.normal(0, 1, num_sims)  # Diffusion shock
    Z2 = np.random.normal(0, 1, num_sims)  # Jump size shock

    #  jump component
    jump_part = np.zeros(num_sims)
    non_zero_mask = poisson_jumps > 0
    jump_part[non_zero_mask] = (
        mu_j * poisson_jumps[non_zero_mask] +
        delta_j * np.sqrt(poisson_jumps[non_zero_mask]) * Z2[non_zero_mask]
    )

    # total log return
    log_return = total_drift + sigma_sqrtT * Z1 + jump_part

    # terminal stock prices
    S_T = S0 * np.exp(log_return)

    # payoffs
    if opt_type == 'call':
        payoff = np.maximum(S_T - K, 0)
    else:  # put
        payoff = np.maximum(K - S_T, 0)

    # Discount payoffs and  average
    price = np.exp(-r * T) * np.mean(payoff)
    return round(price, 2)

# Parameters for Question 9
params_q9 = {
    'S0': 80.0,       # Spot price
    'K': 80.0,        # ATM strike
    'T': 0.25,        # 3 months
    'r': 0.055,       # Risk-free rate
    'sigma': 0.35,    # Volatility
    'lam': 0.25,      # Reduced jump intensity
    'mu_j': -0.5,     # Mean jump size
    'delta_j': 0.22   # Jump size volatility
}

# Price call and put options
call_price_q9 = merton_jump_price(opt_type='call', **params_q9)
put_price_q9 = merton_jump_price(opt_type='put', **params_q9)

results_q9 = pd.DataFrame({
    'Question': ['9', '9'],
    'Option Type': ['Call', 'Put'],
    'Price': [call_price_q9, put_price_q9]
})

print("Table 3: Merton Model with λ = 0.25")
print(results_q9.to_string(index=False))

Table 3: Merton Model with λ = 0.25
Question Option Type  Price
       9        Call   6.85
       9         Put   5.68

From Table 3, we see that the Merton model results with $\lambda=0.25$ show a call price of $6.85 and put price of $5.68, revealing how reduced jump frequency moderates the extreme pricing seen in Question 8 (Q8). The narrower $1.17 call-put spread (versus Q8's $1.11) still reflects persistent crash pricing from negative jumps ($\mu=-0.5$), but the overall lower premiums demonstrate decreased sensitivity to tail events. The call maintains 78% of its Q8 value due to ongoing diffusion volatility ($\sigma=35\%$), while the put retains just 79% of its Q8 price as reduced jump intensity ($\lambda=0.25$) diminishes its crash protection premium. This pricing pattern suggests a flatter volatility smile compared to Q8, with options behaving more like Black-Scholes but still pricing in meaningful jump risk, particularly for puts which remain elevated relative to non-jump models. The results highlight how jump intensity directly controls the weight assigned to extreme moves in option valuation.

Question 10

We'll calculate delta and gamma for options priced in Questions 8 ($\lambda=0.75$) and 9 ($\lambda=0.25$) using central difference approximation. The implementation ensures consistent random paths across perturbations to minimize Monte Carlo noise. Note that we used identical central difference formulas as Question 7, applied to Merton model prices.

import numpy as np
import pandas as pd

np.random.seed(42)
def merton_jump_price(S0, K, T, r, sigma, lam, mu_j, delta_j, opt_type, num_sims=100000, seed=42):
    """Prices European options under Merton jump-diffusion model with seed control"""
    np.random.seed(seed)
    kappa = np.exp(mu_j + 0.5 * delta_j**2) - 1
    drift_adj = r - 0.5 * sigma**2 - lam * kappa
    total_drift = drift_adj * T
    sigma_sqrtT = sigma * np.sqrt(T)

    poisson_jumps = np.random.poisson(lam * T, num_sims)
    Z1 = np.random.normal(0, 1, num_sims)
    Z2 = np.random.normal(0, 1, num_sims)

    jump_part = np.zeros(num_sims)
    non_zero_mask = poisson_jumps > 0
    jump_part[non_zero_mask] = (
        mu_j * poisson_jumps[non_zero_mask] +
        delta_j * np.sqrt(poisson_jumps[non_zero_mask]) * Z2[non_zero_mask]
    )

    log_return = total_drift + sigma_sqrtT * Z1 + jump_part
    S_T = S0 * np.exp(log_return)

    if opt_type == 'call':
        payoff = np.maximum(S_T - K, 0)
    else:
        payoff = np.maximum(K - S_T, 0)

    return np.exp(-r * T) * np.mean(payoff)

def calculate_merton_greeks(lam, call_price_base, put_price_base):
    """Calculates delta and gamma for Merton model options"""
    dS = 0.80  # 1% of $80 stock price
    base_params = {
        'K': 80.0, 'T': 0.25, 'r': 0.055, 'sigma': 0.35,
        'mu_j': -0.5, 'delta_j': 0.22, 'num_sims': 100000
    }

    # Call option Greeks
    call_plus = merton_jump_price(S0=80.80, lam=lam, opt_type='call', seed=42, **base_params)
    call_minus = merton_jump_price(S0=79.20, lam=lam, opt_type='call', seed=42, **base_params)
    call_delta = (call_plus - call_minus) / (2 * dS)
    call_gamma = (call_plus - 2 * call_price_base + call_minus) / (dS ** 2)

    # Put option Greeks
    put_plus = merton_jump_price(S0=80.80, lam=lam, opt_type='put', seed=42, **base_params)
    put_minus = merton_jump_price(S0=79.20, lam=lam, opt_type='put', seed=42, **base_params)
    put_delta = (put_plus - put_minus) / (2 * dS)
    put_gamma = (put_plus - 2 * put_price_base + put_minus) / (dS ** 2)

    return {
        'call': {'delta': call_delta, 'gamma': call_gamma},
        'put': {'delta': put_delta, 'gamma': put_gamma}
    }

# Base prices from Q8 and Q9
q8_prices = {'call': 6.85, 'put': 5.68}  # λ=0.75
q9_prices = {'call': 5.78, 'put': 3.82}  # λ=0.25

# Greeks
q8_greeks = calculate_merton_greeks(lam=0.75, call_price_base=6.85, put_price_base=5.68)
q9_greeks = calculate_merton_greeks(lam=0.25, call_price_base=5.78, put_price_base=3.82)

results = pd.DataFrame({
    'Question': ['8', '8', '9', '9'],
    'Option Type': ['Call', 'Put', 'Call', 'Put'],
    'λ': [0.75, 0.75, 0.25, 0.25],
    'Price': [6.85, 5.68, 5.78, 3.82],
    'Delta': [
        q8_greeks['call']['delta'],
        q8_greeks['put']['delta'],
        q9_greeks['call']['delta'],
        q9_greeks['put']['delta']
    ],
    'Gamma': [
        q8_greeks['call']['gamma'],
        q8_greeks['put']['gamma'],
        q9_greeks['call']['gamma'],
        q9_greeks['put']['gamma']
    ]
})

results['Delta'] = results['Delta'].round(3)
results['Gamma'] = results['Gamma'].round(3)

print("Table 4: Delta and Gamma for Merton Options")
print(results.to_string(index=False))

Table 4: Delta and Gamma for Merton Options
Question Option Type     λ  Price  Delta  Gamma
       8        Call  0.75   6.85  0.650  4.538
       8         Put  0.75   5.68 -0.350  4.725
       9        Call  0.25   5.78  0.600  3.374
       9         Put  0.25   3.82 -0.401  5.834

From Table 4, we see that calls remain positive (about 0.65 at $\lambda = 0.75$ and 0.60 at $\lambda = 0.25$), reflecting bullish exposure that strengthens with more frequent jumps, while puts stay negative (around –0.35 at $\lambda = 0.75$ and –0.40 at $\lambda = 0.25$), indicating bearish exposure that becomes slightly more sensitive when jumps are rarer. At low jump intensity ($\lambda = 0.25$), put gamma ($\approx 5.83$) exceeds call gamma ($\approx 3.37$), meaning out-of-the-money puts develop greater curvature when jumps are infrequent; at higher jump intensity ($\lambda = 0.75$), call and put gammas converge ($\approx$ 4.54 vs. $\approx 4.73$), reflecting more symmetric tail-risk once jumps are common. Across both $\lambda’s$, gamma is well above zero, confirming that incorporating jump risk ($\mu = –0.5$) adds extra convexity compared to Black–Scholes especially for puts so the Merton model produces a steeper implied volatility smile for out-of-the-money puts as markets price in downside crash risk.

Team member C

Question 11

To verify put-call parity for European options priced under the Heston (Questions 5-6) and Merton (Questions 8-9) models, we use the fundamental put-call parity equation: $C - P = S_0 - K·e^{-rT}$ (Hull 250), where C is call price, P is put price, $S_0$ is current stock price (80), K is strike price (80, ATM), r is risk-free rate (5.5% = 0.055) and T is time to maturity (0.25 years). For validation, we compared $C - P (LHS)$ to $S_0 - K e^{-rT}$ (RHS) with tolerance $0.10, where RHS is the right-hand side and LHS is the left-hand side (LHS) from the model prices.

import numpy as np
import pandas as pd

np.random.seed(42)
# Put-call parity parameters
S0 = 80.0
K = 80.0
r = 0.055
T = 0.25
rhs = S0 - K * np.exp(-r * T)  # Exact theoretical value: 1.09544

results = {
    'Q5 (Heston ρ=-0.30)': (3.46, 2.39),
    'Q6 (Heston ρ=-0.70)': (3.49, 2.41),
    'Q8 (Merton λ=0.75)': (8.30, 7.19),
    'Q9 (Merton λ=0.25)': (6.85, 5.68)
}

# Validate parity
parity_data = []
for case, (call, put) in results.items():
    lhs = call - put
    difference = lhs - rhs
    satisfies = "Yes" if abs(difference) <= 0.10 else "No"
    parity_data.append({
        'Question': case,
        'Call Price': call,
        'Put Price': put,
        'C - P (LHS)': round(lhs, 3),
        'S₀ - K·e^{-rT} (RHS)': round(rhs, 3),
        'Difference': round(difference, 3),
        'Satisfies Parity?': satisfies
    })

df_parity = pd.DataFrame(parity_data)
print("Table 5: Put-Call Parity Validation")
print(df_parity.to_string(index=False))

Table 5: Put-Call Parity Validation
            Question  Call Price  Put Price  C - P (LHS)  S₀ - K·e^{-rT} (RHS)  Difference Satisfies Parity?
 Q5 (Heston ρ=-0.30)        3.46       2.39         1.07                 1.092      -0.022               Yes
 Q6 (Heston ρ=-0.70)        3.49       2.41         1.08                 1.092      -0.012               Yes
  Q8 (Merton λ=0.75)        8.30       7.19         1.11                 1.092       0.018               Yes
  Q9 (Merton λ=0.25)        6.85       5.68         1.17                 1.092       0.078               Yes

The results (Table 5) confirm that all models satisfy put-call parity within the acceptable tolerance, with the Heston model showing tighter alignment (differences of -$0.02247 for $\rho=-0.30$ and -$0.01247 for $\rho=-0.70$) due to lower simulation noise, while the Merton model exhibits slightly larger deviations ($0.01753 for $\lambda=0.75$ and $0.07753 for $\lambda=0.25$) attributable to the increased variance introduced by jump processes, particularly at lower jump intensities where the impact of infrequent but significant jumps amplifies pricing discrepancies. That's, all four cases satisfy put‐call parity within a few cents, confirming arbitrage‐free pricing. The Heston model’s small negative differences (–$0.02247 and –$0.01247) reflect slight call undervaluation due to negative correlation limiting upside volatility. The Merton model’s positive deviations (+$0.01753 and +$0.07753) stem from jump‐induced put inflation, especially when jumps are less frequent at λ=0.25. Overall, parity holds tightly, validating the Monte Carlo implementation despite increased variance from jumps.

Question 12

We'll price European call options for 7 strikes using Heston Model ($\rho = -0.30$ from Q5 and Merton Model ($\lambda = 0.75$ from Q8). Strikes are derived from moneyness values ($S_0/K=K_i = \frac{S_0}{\text{moneyness}_i}, \quad \text{moneyness} \in \{0.85, 0.90, \dots, 1.15\}$). Same Heston/Merton Monte Carlo as Questions 5/8, reused for all strikes.

We'll run one simulation per model (500,000 paths, 100 steps) and price all strikes using the same terminal stock prices for efficiency. Full-truncation scheme (Heston) and lognormal jumps (Merton) are implemented as before.

import numpy as np
import pandas as pd

# Parameters
S0 = 80.0
r = 0.055
T = 0.25
n_simulations = 500000
n_steps = 100
dt = T / n_steps
moneyness_values = [0.85, 0.90, 0.95, 1.00, 1.05, 1.10, 1.15]
strikes = [round(S0 / m, 4) for m in moneyness_values]  # Exact strikes

# =====================
# Heston Model (ρ=-0.30)
# =====================
np.random.seed(42)
v0, kappa, theta, sigma_v, rho = 0.032, 1.85, 0.045, 0.35, -0.30

# Correlated Brownian motions
dW = np.random.normal(size=(n_simulations, n_steps, 2))
dW_S = rho * dW[:, :, 0] + np.sqrt(1 - rho**2) * dW[:, :, 1]

v = np.full(n_simulations, v0)
lnS = np.log(S0) * np.ones(n_simulations)
for t in range(n_steps):
    v_positive = np.maximum(v, 0)
    sqrt_v = np.sqrt(v_positive)
    v += kappa * (theta - v_positive) * dt + sigma_v * sqrt_v * np.sqrt(dt) * dW[:, t, 0]
    lnS += (r - 0.5 * v_positive) * dt + sqrt_v * np.sqrt(dt) * dW_S[:, t]
S_T_heston = np.exp(lnS)

# =====================
# Merton Model
# =====================
np.random.seed(42)
lam, mu_j, delta_j = 0.75, -0.5, 0.22
kappa = np.exp(mu_j + 0.5 * delta_j**2) - 1
drift_adj = r - 0.5 * sigma_v**2 - lam * kappa

# Simulate jumps and diffusion
poisson_jumps = np.random.poisson(lam * T, n_simulations)
Z_diff = np.random.normal(0, 1, n_simulations)
Z_jump = np.random.normal(0, 1, n_simulations)
jump_part = mu_j * poisson_jumps + delta_j * np.sqrt(poisson_jumps) * Z_jump
log_return = drift_adj * T + sigma_v * np.sqrt(T) * Z_diff + jump_part
S_T_merton = S0 * np.exp(log_return)

# =====================
# Price options for all strikes
# =====================
def price_calls(S_T, strikes):
    prices = []
    for K in strikes:
        payoff = np.maximum(S_T - K, 0)
        price = np.exp(-r * T) * np.mean(payoff)
        prices.append(round(price, 2))
    return prices

heston_prices = price_calls(S_T_heston, strikes)
merton_prices = price_calls(S_T_merton, strikes)

results = pd.DataFrame({
    'Moneyness (S0/K)': moneyness_values,
    'Strike': strikes,
    'Heston Call Price': heston_prices,
    'Merton Call Price': merton_prices
})

print("Table 6: European Call Prices for Different Strikes")
print(results.to_string(index=False))

Table 6: European Call Prices for Different Strikes
 Moneyness (S0/K)   Strike  Heston Call Price  Merton Call Price
             0.85  94.1176               0.14               2.80
             0.90  88.8889               0.55               4.34
             0.95  84.2105               1.59               6.19
             1.00  80.0000               3.46               8.29
             1.05  76.1905               5.96              10.52
             1.10  72.7273               8.76              12.81
             1.15  69.5652              11.59              15.09

The results (Table 6) show that Merton call prices are uniformly higher than Heston across all strikes ($2.80 vs $0.14 at 94.12 OTM; $15.09 vs $11.59 at 69.57 ITM), reflecting how jump risk amplifies tail exposure despite negative mean jump size ($\mu=-0.5$). The Heston model's negative correlation ($\rho=-0.30$) suppresses OTM call valuations (only $0.14 at 94.12 strike) by linking price rallies to volatility damping, while its stochastic volatility generates moderate convexity—prices rise smoothly from $0.14 (94.12) to $11.59 (69.57). Conversely, in Merton model, OTM calls retain substantial value ($2.80 at 94.12) due to upward jump potential, while ITM calls command significant premiums ($15.09 at 69.57) as jumps override strike barriers.

All team members

Question 13

To solve Question 13, we price an American call option using the Heston model with correlation $\rho = -0.30$ (same parameters as Question 5) and compute its delta and gamma. Unlike European options, American options allow early exercise, requiring a different pricing approach. We use the Longstaff-Schwartz Monte Carlo (LSM) algorithm for valuation. Note that for longstaff-Schwartz Algorithm we used the formula: $V_t = \max\left(\text{exercise}_t, \mathbb{E}[e^{-r \Delta t} V_{t+1} | \mathcal{F}_t]\right)$, where continuation value is estimated via regression on $S_t$ and $v_t$ (Longstaff and Schwartz 113). For Greeks, same central differences as Question 7.

import numpy as np
import pandas as pd

np.random.seed(42)
# ====================
# Heston Model Parameters
# ====================
S0 = 80.0
K = 80.0
r = 0.055
T = 0.25
v0 = 0.032
kappa = 1.85
theta = 0.045
sigma_v = 0.35
rho = -0.30
n_simulations = 500000
n_steps = 100
dt = T / n_steps
dS = 0.80

# ====================
# Heston Paths
# ====================
def simulate_heston_paths(S0, seed=None):
    np.random.seed(seed)
    dW = np.random.normal(size=(n_simulations, n_steps, 2))
    dW_v = dW[:, :, 0]
    dW_S = rho * dW[:, :, 0] + np.sqrt(1 - rho**2) * dW[:, :, 1]

    v = np.full(n_simulations, v0)
    lnS = np.full(n_simulations, np.log(S0))

    S_paths = np.zeros((n_simulations, n_steps + 1))
    v_paths = np.zeros((n_simulations, n_steps + 1))
    S_paths[:, 0] = S0
    v_paths[:, 0] = v0

    # Euler‐Maruyama discretization
    for t in range(n_steps):
        v_positive = np.maximum(v, 0)
        sqrt_v = np.sqrt(v_positive)

        # Variance update
        v += kappa * (theta - v_positive) * dt + sigma_v * sqrt_v * np.sqrt(dt) * dW_v[:, t]

        # Log‐price update
        lnS += (r - 0.5 * v_positive) * dt + sqrt_v * np.sqrt(dt) * dW_S[:, t]

        # Store S and v at next time
        S_paths[:, t + 1] = np.exp(lnS)
        v_paths[:, t + 1] = v

    return S_paths, v_paths

# ====================
# Longstaff‐Schwartz for American Call
# ====================
def lsm_american_call(S_paths, v_paths):
    df = np.exp(-r * dt)
    payoff = np.maximum(S_paths[:, -1] - K, 0)
    cash_flow = payoff.copy()

    for t in range(n_steps - 1, 0, -1):
        in_money = S_paths[:, t] > K
        if np.any(in_money):
            X = S_paths[in_money, t]
            Y = cash_flow[in_money] * df

            v_t = v_paths[in_money, t]
            X_basis = np.vstack([
                np.ones_like(X),
                X,
                X**2,
                v_t,
                X * v_t
            ]).T

            coeffs = np.linalg.lstsq(X_basis, Y, rcond=None)[0]
            continuation = X_basis @ coeffs

            immediate_ex = X - K

            exercise_mask = immediate_ex > continuation
            exercise_idx_global = np.where(in_money)[0][exercise_mask]

            cash_flow[exercise_idx_global] = immediate_ex[exercise_mask]
            cash_flow[exercise_idx_global] *= np.exp(r * (t * dt))

    return np.mean(cash_flow * np.exp(-r * T))

# ====================
# Price American Call
# ====================
S_paths_base, v_paths_base = simulate_heston_paths(S0, seed=42)
am_call_price = lsm_american_call(S_paths_base, v_paths_base)
print(f"American Call Price: ${am_call_price:.2f}")

# ====================
# Delta and Gamma Calculation
# ====================
# S0 + dS
S_paths_plus, _ = simulate_heston_paths(S0 + dS, seed=42)
call_plus = lsm_american_call(S_paths_plus, v_paths_base)

# S0 - dS
S_paths_minus, _ = simulate_heston_paths(S0 - dS, seed=42)
call_minus = lsm_american_call(S_paths_minus, v_paths_base)

# Greeks
delta = (call_plus - call_minus) / (2 * dS)
gamma = (call_plus - 2 * am_call_price + call_minus) / (dS ** 2)

print(f"Delta: {delta:.3f}")
print(f"Gamma: {gamma:.3f}")

# ====================
# Comparison with European Call
# ====================
euro_price = 3.46
euro_delta = 0.606
euro_gamma = 0.047

euro = {
    'Price': euro_price,
    'Delta': euro_delta,
    'Gamma': euro_gamma
}

american = {
    'Price': round(am_call_price, 3),
    'Delta': round(delta, 3),
    'Gamma': round(gamma, 3)
}

df = pd.DataFrame({
    'European': euro,
    'American': american
})

df = df.round(3)

print("Table 7: Questions 5 and 7 for the case of an American call option")
print(df)

American Call Price: $3.48

The American call price from Question 13 is $3.48 (Table 7), compared with the European call price of $3.46, its delta of 0.60 is slightly below the European delta of 0.606, while its gamma of 0.05 exceeds the European gamma of 0.047, indicating that early‐exercise flexibility in the American option adds a small premium to the price and increases convexity (gamma) even as sensitivity to spot (delta) remains very similar. This, is because when there are no dividends, any early‐exercise premium is essentially zero, so this gap reflects Monte Carlo/Longstaff–Schwartz noise. The delta for the American call (0.603) is almost identical to the European value (0.606) , since early exercise has negligible impact at‐the‐money for a non‐dividend stock. The gamma for the American call (0.053) is slightly higher than the European gamma (0.047) , as optionality can increase convexity near the strike, though short maturity and no dividends keep this effect very small.

Question 14

To solve this question we used the following formulas: Payoff: $\text{Payoff} = \begin{cases} \max(S_T - K, 0) & \text{if } \max_{t \leq T} S_t \geq H \\ 0 & \text{otherwise} \end{cases}$
(Hull 530), and pricing as discounted expected payoff with barrier condition.

import numpy as np

# Parameters from Question 6 (Heston ρ=-0.70)
S0 = 80.0
K_barrier = 95.0
H = 95.0  # Barrier
r = 0.055
T = 0.25
v0 = 0.032
kappa = 1.85
theta = 0.045
sigma_v = 0.35
rho = -0.70
n_simulations = 500000
n_steps = 100
dt = T / n_steps

np.random.seed(42)

# Correlated Brownian motions
dW = np.random.normal(size=(n_simulations, n_steps, 2))
dW_v = dW[:, :, 0]
dW_S = rho * dW[:, :, 0] + np.sqrt(1 - rho**2) * dW[:, :, 1]

v = np.full(n_simulations, v0)
lnS = np.full(n_simulations, np.log(S0))
S_paths = np.zeros((n_simulations, n_steps + 1))
S_paths[:, 0] = S0

# Euler-Maruyama simulation with full truncation
for t in range(n_steps):
    v_positive = np.maximum(v, 0)
    sqrt_v = np.sqrt(v_positive)
    v += kappa * (theta - v_positive) * dt + sigma_v * sqrt_v * np.sqrt(dt) * dW_v[:, t]
    lnS += (r - 0.5 * v_positive) * dt + sqrt_v * np.sqrt(dt) * dW_S[:, t]
    S_paths[:, t + 1] = np.exp(lnS)

# Terminal prices
S_T = S_paths[:, -1]

# barrier breach (up-and-in condition)
max_S = np.max(S_paths, axis=1)
barrier_breached = max_S >= H

# Payoffs
cui_payoff = np.where(barrier_breached, np.maximum(S_T - K_barrier, 0), 0)
# European call (no barrier)
euro_payoff = np.maximum(S_T - K_barrier, 0)

# Discounted prices
cui_price = np.exp(-r * T) * np.mean(cui_payoff)
euro_price = np.exp(-r * T) * np.mean(euro_payoff)

cui_price_rounded = round(cui_price, 3)
euro_price_rounded = round(euro_price, 3)

print(f"European Up-and-In Call (Barrier $95) Price: ${cui_price_rounded:.3f}")
print(f"Standard European Call (Strike $95) Price: ${euro_price_rounded:.3f}")
print(f"Difference: ${euro_price_rounded - cui_price_rounded:.3f}")

European Up-and-In Call (Barrier $95) Price: $0.030
Standard European Call (Strike $95) Price: $0.030
Difference: $0.000

The up-and-in call and the standard European call both price at $0.030, meaning every path that ends in-the-money must have crossed the $95 barrier at some point. Under the Heston model’s continuous dynamics, reaching $S_T > 95$ guarantees $\max_{t\le T}S_t \ge 95$, so the barrier never excludes any payoff. Consequently, the two contracts share identical payoffs and identical prices. Even with strong negative correlation ($\rho = -0.70$) driving volatility higher when returns fall, there are no price “gaps” that evade the barrier, so the deep out‐of‐the‐money call remains cheap. This exact equivalence confirms that, for continuous processes, an up-and-in call with barrier equal to strike is mathematically the same as a plain European call.

Question 15

To solve this question we used the following formulas: Payoff: $\text{Payoff} = \ \text{Payoff} = \begin{cases} \max(K - S_T, 0) & \text{if } \min_{t \leq T} S_t \leq H \\ 0 & \text{otherwise} \end{cases}$
(Hull 530), and pricing as discounted expected payoff with barrier condition.

import numpy as np

# Parameters from Question 8 (Merton model with λ=0.75)
S0 = 80.0
K = 65.0
H = 65.0  # Barrier
r = 0.055
T = 0.25
sigma = 0.35
lam = 0.75
mu_j = -0.5
delta_j = 0.22
n_simulations = 500000
n_steps = 100
dt = T / n_steps

# Jump compensator
kappa = np.exp(mu_j + 0.5 * delta_j**2) - 1
drift_adj = r - 0.5 * sigma**2 - lam * kappa

np.random.seed(42)


lnS = np.log(S0) * np.ones(n_simulations)
min_price = S0 * np.ones(n_simulations)  # Track minimum price per path

# Simulate paths
for t in range(n_steps):
    # Diffusion component
    Z_diff = np.random.normal(0, 1, n_simulations)
    diffusion = drift_adj * dt + sigma * np.sqrt(dt) * Z_diff

    # Jump component
    n_jumps = np.random.poisson(lam * dt, n_simulations)
    jump_sizes = np.zeros(n_simulations)
    for i in range(n_simulations):
        if n_jumps[i] > 0:
            jumps = np.random.normal(mu_j, delta_j, n_jumps[i])
            jump_sizes[i] = np.sum(jumps)

    # Update log-price and minimum price
    lnS += diffusion + jump_sizes
    S_t = np.exp(lnS)
    min_price = np.minimum(min_price, S_t)

# Terminal prices
S_T = np.exp(lnS)

# Payoffs
barrier_breached = min_price <= H
pdi_payoff = np.where(barrier_breached, np.maximum(K - S_T, 0), 0)
euro_put_payoff = np.maximum(K - S_T, 0)

# Discounted prices
pdi_price = np.exp(-r * T) * np.mean(pdi_payoff)
euro_put_price = np.exp(-r * T) * np.mean(euro_put_payoff)

print(f"European Down-and-In Put (Barrier $65) Price: ${pdi_price:.3f}")
print(f"Standard European Put (Strike $65) Price: ${euro_put_price:.3f}")
print(f"Difference: ${euro_put_price - pdi_price:.3f}")

European Down-and-In Put (Barrier $65) Price: $2.779
Standard European Put (Strike $65) Price: $2.779
Difference: $0.000

The down-and-in put option with barrier equal to the strike price is equivalent to a standard European put option. This is because the barrier condition (stock price falling to or below the barrier) is automatically satisfied whenever the option finishes in the money (i.e., the terminal stock price is below the strike). Moreover, if the option finishes out of the money, the barrier condition being met does not result in a payoff. Therefore, the two options have identical payoffs and thus the same price. The Merton model with jumps does not alter this equivalence because the barrier is monitored continuously and the terminal condition is included. The result of $2.779 for both confirms this equivalence.

Conclusion

Monte Carlo implementations of the Heston and Merton models produce option prices and sensitivities that behave as theory predicts: stronger negative stock‐volatility correlation under Heston slightly raises calls and puts and tightens skew, while jump intensity in Merton meaningfully raises OTM call values and puts’ curvature. Put–call parity holds within simulation tolerances in all cases. Early exercise value for an American Heston call on a non‐dividend stock is negligible, matching its European price, and barrier options with barriers at the strike reduce to their vanilla counterparts under continuous monitoring, regardless of model features. Overall, Heston captures smooth volatility clustering and leverage effects, whereas Merton’s jumps better reflect fat tails and crash risk guiding model choice based on the need for skew versus tail accuracy.

References

Andersen, Leif B. G., and Brotherton-Ratcliffe, Mark L. B. "Efficient Simulation of the Heston Stochastic Volatility Model." Journal of Computational Finance, 2007, pp. 1–38.

Glasserman, Paul. Monte Carlo Methods in Financial Engineering. Springer, 2004.

Heston, Steven L. "A Closed-Form Solution for Options with Stochastic Volatility." Review of Financial Studies, vol. 6, no. 2, 1993, pp. 327–343, https://doi.org/10.1093/rfs/6.2.327.

Hull, John C. Options, Futures, and Other Derivatives. 10th ed., Pearson, 2017.

Kahl, Clemens, and Jäckel, Peter. "Not-So-Complex Logarithms in the Heston Model." Wilmott Magazine, Mar. 2006, pp. 110–119.

Longstaff, Francis A., and Schwartz, Eduardo S. "Valuing American Options by Simulation." Review of Financial Studies, vol. 14, no. 1, 2001, pp. 113–147, https://doi.org/10.1093/rfs/14.1.113.

Merton, Robert C. "Option Pricing When Underlying Stock Returns Are Discontinuous." Journal of Financial Economics, vol. 3, no. 1–2, 1976, pp. 125–144, https://doi.org/10.1016/0304-405X(76)90022-2.