PAM-Driven Volatility and Option Pricing¶
This notebook explores a toy model connecting the parabolic Anderson model / stochastic heat equation to volatility modelling in quantitative finance.
We begin with the Black--Scholes model and Monte Carlo pricing of European call options. We then replace constant volatility with time-dependent volatility paths. Finally, we generate a toy PAM/SHE field and use one spatial slice of that field to drive the volatility of a stock price process.
The main question is:
How does intermittent volatility change endpoint distributions and option prices?
This notebook is intended as a mathematical and computational experiment; AI is used throughout this document.
Imports¶
import math
import random
import matplotlib.pyplot as plt
2. Black--Scholes terminal stock simulation¶
In the Black--Scholes model,
$$ dS_t = rS_t\,dt + \sigma S_t\,dB_t. $$
The terminal stock price is
$$ S_T = S_0 \exp\left((r-\tfrac12\sigma^2)T + \sigma\sqrt{T}Z\right), $$
where $$Z \sim N(0,1)$$.
def simulate_stock_price(S0, T, r, sigma):
"""
Simulate one terminal stock price under the Black-Scholes model.
"""
Z = random.gauss(0, 1)
ST = S0 * math.exp(
(r - 0.5 * sigma ** 2) * T
+ sigma * math.sqrt(T) * Z
)
return ST
3. European call payoff and Monte Carlo pricing¶
A European call option with strike (K) has payoff
$$ (S_T-K)^+ = \max(S_T-K,0). $$
We estimate the option price by simulating many terminal stock prices, averaging the payoffs, and discounting back to time zero:
$$ C_0 \approx e^{-rT}\frac{1}{N}\sum_{j=1}^N (S_T^{(j)}-K)^+. $$
def monte_carlo_call_price(S0, K, T, r, sigma, n_paths):
"""
Estimate the price of a European call option by Monte Carlo simulation.
"""
payoffs = []
for _ in range(n_paths):
ST = simulate_stock_price(S0, T, r, sigma)
payoff = max(ST - K, 0.0)
payoffs.append(payoff)
average_payoff = sum(payoffs) / n_paths
price = math.exp(-r * T) * average_payoff
return price
The Monte Carlo price should be close to the exact Black-Scholes price. Since the method uses random simulation, the estimate will fluctuate slightly from run to run. Increasing n_paths reduces the Monte Carlo error.
4. Exact Black-Scholes formula¶
For a European call option, the exact Black-Scholes price is
$$ C_0 = S_0 \Phi(d_1) - K e^{-rT}\Phi(d_2), $$
where
$$ d_1 = \frac{\log(S_0/K) + (r+\frac12\sigma^2)T}{\sigma\sqrt{T}}, $$
and
$$ d_2 = d_1 - \sigma\sqrt{T}. $$
Here $\Phi$ denotes the standard normal cumulative distribution function.
def normal_cdf(x):
"""
Standard normal cumulative distribution function.
"""
return 0.5 * (1.0 + math.erf(x / math.sqrt(2.0)))
def black_scholes_call_price(S0, K, T, r, sigma):
"""
Compute the exact Black-Scholes price of a European call option.
"""
if T <= 0:
return max(S0 - K, 0.0)
d1 = (
math.log(S0 / K)
+ (r + 0.5 * sigma ** 2) * T
) / (sigma * math.sqrt(T))
d2 = d1 - sigma * math.sqrt(T)
price = (
S0 * normal_cdf(d1)
- K * math.exp(-r * T) * normal_cdf(d2)
)
return price
5. Simulating full stock paths¶
So far, we have simulated only the terminal stock price (S_T). To study path-dependent behaviour, such as barrier crossings or time-dependent volatility, we need to simulate the full stock path.
We divide the time interval ([0,T]) into n_steps small time intervals of size
$$ \Delta t = \frac{T}{n_{\text{steps}}}. $$
The Black-Scholes update over one small step is
$$ S_{t+\Delta t} = S_t \exp\left((r-\frac12\sigma^2)\Delta t + \sigma\sqrt{\Delta t}Z \right), $$
where $Z \sim N(0,1)$.
def simulate_stock_path(S0, T, r, sigma, n_steps):
"""
Simulate a full stock price path under the Black-Scholes model.
Returns a list:
[S_0, S_dt, S_2dt, ..., S_T]
"""
dt = T / n_steps
path = [S0]
S = S0
for _ in range(n_steps):
Z = random.gauss(0, 1)
S = S * math.exp(
(r - 0.5 * sigma ** 2) * dt
+ sigma * math.sqrt(dt) * Z
)
path.append(S)
return path
This plot shows one possible stock price path under the Black-Scholes model. Each run produces a different path because the model is driven by random Gaussian increments.
6. Time-dependent volatility paths¶
The Black-Scholes model assumes constant volatility. To move toward PAM-driven volatility, we first allow the variance to change over time.
We write
$$ v_t = \sigma_t^2. $$
Over a small time interval, the stock price update becomes
$$ S_{t+\Delta t} = S_t \exp\left((r-\frac12 v_t)\Delta t+\sqrt{v_t\Delta t}Z\right), $$
where $Z \sim N(0,1)$.
In the code below, volatility_path is actually a path of variance values:
$$ v_0, v_1, \dots, v_{n-1}. $$
def simulate_stock_path_with_volatility_path(S0, T, r, volatility_path):¶
"""
Simulate a stock price path using a time-dependent variance path.
volatility_path contains variance values v_m, not volatility values sigma_m.
"""
n_steps = len(volatility_path)
dt = T / n_steps
stock_path = [S0]
S = S0
for v in volatility_path:
Z = random.gauss(0, 1)
S = S * math.exp(
(r - 0.5 * v) * dt
+ math.sqrt(v * dt) * Z
)
stock_path.append(S)
return stock_path 6. Time-dependent volatility paths¶
The Black-Scholes model assumes constant volatility. To move toward PAM-driven volatility, we first allow the variance to change over time.
Instead of using one constant volatility parameter $\sigma$, we use a variance path
$$ v_0, v_1, \dots, v_{n-1}. $$
Over one small time step, the stock price update becomes
$$ S_{t+\Delta t} = S_t \exp\left((r-\frac12 v_t)\Delta t + \sqrt{v_t \Delta t}Z \right), $$
where $Z \sim N(0,1)$.
This allows us to study models where volatility is calm during some periods and high during others.
def simulate_stock_path_with_volatility_path(S0, T, r, volatility_path):
"""
Simulate a stock price path using a time-dependent variance path.
volatility_path contains variance values v_m, not volatility values sigma_m.
"""
n_steps = len(volatility_path)
dt = T / n_steps
stock_path = [S0]
S = S0
for v in volatility_path:
Z = random.gauss(0, 1)
S = S * math.exp(
(r - 0.5 * v) * dt
+ math.sqrt(v * dt) * Z
)
stock_path.append(S)
return stock_path
random.seed(42)
S0 = 100
T = 1
r = 0.05
n_steps = 100
constant_variance_path = [0.2 ** 2] * n_steps
stock_path = simulate_stock_path_with_volatility_path(
S0=S0,
T=T,
r=r,
volatility_path=constant_variance_path
)
times = [i * (T / n_steps) for i in range(n_steps + 1)]
plt.figure(figsize=(8, 5))
plt.plot(times, stock_path)
plt.xlabel("Time")
plt.ylabel("Stock price")
plt.title("Stock path with constant variance path")
plt.show()
When the variance path is constant, this reproduces the usual Black-Scholes path simulation. The advantage of this formulation is that we can later replace the constant variance path with a volatility-spike path or a PAM-driven variance path.
7. Volatility spikes¶
Before introducing PAM-driven volatility, we first create a simple hand-made volatility spike model.
The idea is that markets may have calm periods and crisis periods. During a crisis window, the variance becomes much larger.
We start with a baseline variance
$$ v_{\text{base}} = 0.2^2, $$
and then manually insert periods of higher variance.
n_steps = 100
base_variance = 0.2 ** 2
high_variance = 0.6 ** 2
spike_variance_path = [base_variance] * n_steps
# Insert a volatility spike between time steps 40 and 60
for i in range(40, 60):
spike_variance_path[i] = high_variance
variance_times = [i * (T / n_steps) for i in range(n_steps)]
plt.figure(figsize=(8, 5))
plt.plot(variance_times, spike_variance_path)
plt.xlabel("Time")
plt.ylabel("Variance")
plt.title("Hand-made volatility spike path")
plt.show()
random.seed(42)
stock_path_spike = simulate_stock_path_with_volatility_path(
S0=100,
T=1,
r=0.05,
volatility_path=spike_variance_path
)
times = [i * (T / n_steps) for i in range(n_steps + 1)]
plt.figure(figsize=(8, 5))
plt.plot(times, stock_path_spike)
plt.xlabel("Time")
plt.ylabel("Stock price")
plt.title("Stock path with volatility spike")
plt.show()
The volatility spike increases the randomness of the stock path during the crisis window. It does not force the stock upward or downward. Instead, it makes large movements in either direction more likely.
8. Endpoint distributions¶
We now compare the terminal stock price distribution under two models:
- constant volatility,
- hand-made spike volatility.
For each model, we simulate many stock paths and record only the final value (S_T). This lets us see how volatility spikes change the endpoint distribution.
def simulate_endpoint_distribution_with_volatility_path(S0, T, r, volatility_path, n_paths):
"""
Simulate many stock paths and return their terminal values.
"""
endpoints = []
for _ in range(n_paths):
stock_path = simulate_stock_path_with_volatility_path(
S0=S0,
T=T,
r=r,
volatility_path=volatility_path
)
ST = stock_path[-1]
endpoints.append(ST)
return endpoints
random.seed(42)
constant_variance_path = [base_variance] * n_steps
constant_endpoints = simulate_endpoint_distribution_with_volatility_path(
S0=100,
T=1,
r=0.05,
volatility_path=constant_variance_path,
n_paths=5000
)
spike_endpoints = simulate_endpoint_distribution_with_volatility_path(
S0=100,
T=1,
r=0.05,
volatility_path=spike_variance_path,
n_paths=5000
)
plt.figure(figsize=(8, 5))
plt.hist(constant_endpoints, bins=50, alpha=0.5, label="Constant volatility")
plt.hist(spike_endpoints, bins=50, alpha=0.5, label="Spike volatility")
plt.xlabel("Terminal stock price $S_T$")
plt.ylabel("Frequency")
plt.title("Endpoint distributions: constant vs spike volatility")
plt.legend()
plt.show()
The spike-volatility model produces a wider endpoint distribution. Some paths finish much lower, while others finish much higher. This reflects the fact that volatility spikes increase randomness rather than simply pushing the stock in one direction.
9. Option price comparison¶
We now compare European call option prices under constant volatility and spike volatility.
For a call option with strike (K), the payoff is
$$ (S_T-K)^+ = \max(S_T-K,0). $$
Using the simulated endpoint distributions, we approximate the option price by
$$ C_0 \approx e^{-rT}\frac{1}{N}\sum_{j=1}^N (S_T^{(j)}-K)^+. $$
def call_price_from_endpoints(endpoints, K, r, T):
"""
Estimate a European call price from simulated terminal stock prices.
"""
payoffs = []
for ST in endpoints:
payoff = max(ST - K, 0.0)
payoffs.append(payoff)
average_payoff = sum(payoffs) / len(payoffs)
price = math.exp(-r * T) * average_payoff
return price
strikes = [80, 100, 120, 150, 200]
print(f"{'K':>10} {'Constant Vol Price':>20} {'Spike Vol Price':>20} {'Difference':>15}")
print("-" * 75)
for K in strikes:
constant_price = call_price_from_endpoints(
endpoints=constant_endpoints,
K=K,
r=0.05,
T=1
)
spike_price = call_price_from_endpoints(
endpoints=spike_endpoints,
K=K,
r=0.05,
T=1
)
difference = spike_price - constant_price
print(f"{K:>10} {constant_price:>20.4f} {spike_price:>20.4f} {difference:>15.4f}")
K Constant Vol Price Spike Vol Price Difference
---------------------------------------------------------------------------
80 24.7031 25.7745 1.0714
100 10.6071 14.2162 3.6091
120 3.2847 7.1233 3.8387
150 0.3685 2.2930 1.9245
200 0.0035 0.2405 0.2370
The spike volatility model tends to increase call option prices, especially for higher strikes. This is because volatility spikes make the endpoint distribution wider and create more extreme right-tail outcomes. Since call option payoffs are capped below at zero but uncapped above, rare large upward moves can significantly increase the average payoff.
Cell In[16], line 1 The spike volatility model tends to increase call option prices, especially for higher strikes. This is because volatility spikes make the endpoint distribution wider and create more extreme right-tail outcomes. Since call option payoffs are capped below at zero but uncapped above, rare large upward moves can significantly increase the average payoff. ^ SyntaxError: invalid syntax
10. Toy PAM/SHE field¶
We now introduce a toy discrete version of the parabolic Anderson model / stochastic heat equation.
The continuum equation is formally
$$ \partial_t u = \frac12 \Delta u + \beta u \xi, $$
where (\xi) is space-time white noise.
The heat equation term smooths the field, while the multiplicative noise term can create random peaks. This random peak formation is related to the idea of intermittency.
We use a simple explicit finite-difference scheme. This is a toy numerical model rather than a fully rigorous SPDE solver.
10. Toy PAM/SHE field¶
We now introduce a toy discrete version of the parabolic Anderson model / stochastic heat equation.
The continuum equation is formally
$$ \partial_t u = \frac12 \Delta u + \beta u \xi, $$
where (\xi) is space-time white noise.
The heat equation term smooths the field, while the multiplicative noise term can create random peaks. This random peak formation is related to the idea of intermittency.
We use a simple explicit finite-difference scheme. This is a toy numerical model rather than a fully rigorous SPDE solver.
def simulate_toy_pam(T, n_time, n_space, beta):
"""
Simulate a toy discrete PAM/SHE field.
pam_path[m][i] is the value u_m(i),
where m is the time index and i is the space index.
"""
dt = T / n_time
dx = 1 / n_space
alpha = 0.5 * dt / (dx ** 2)
# Initial condition: flat positive field
u = [1.0] * n_space
# Store the whole evolution
pam_path = [u.copy()]
for _ in range(n_time):
new_u = u.copy()
for i in range(1, n_space - 1):
laplacian = u[i + 1] - 2 * u[i] + u[i - 1]
noise = random.gauss(0, 1)
new_u[i] = (
u[i]
+ alpha * laplacian
+ beta * u[i] * math.sqrt(dt / dx) * noise
)
# Toy positivity fix
if new_u[i] < 0:
new_u[i] = 0.0
# Simple boundary condition
new_u[0] = new_u[1]
new_u[-1] = new_u[-2]
u = new_u
pam_path.append(u.copy())
return pam_path
random.seed(42)
T = 1
n_time = 10000
n_space = 20
beta = 0.05
pam_path = simulate_toy_pam(
T=T,
n_time=n_time,
n_space=n_space,
beta=beta
)
len(pam_path)
10001
The object pam_path stores the full simulated field. The value pam_path[m][i] represents the field value at time index (m) and space index (i).
11. Extracting a PAM time series¶
The toy PAM simulation produces a space-time field (u(t,x)). To use this field as a volatility input, we extract the value of the field at one fixed spatial point (x_0) through time.
This gives a time series
$$ u(t,x_0), $$
which we will later transform into a variance path for the stock price model.
space_index = n_space // 2
pam_time_series = []
for m in range(len(pam_path)):
value = pam_path[m][space_index]
pam_time_series.append(value)
len(pam_time_series)
10001
pam_times = [m * (T / n_time) for m in range(n_time + 1)]
plt.figure(figsize=(8, 5))
plt.plot(pam_times, pam_time_series)
plt.xlabel("Time")
plt.ylabel("PAM value at fixed space point")
plt.title("Toy PAM time series at one spatial point")
plt.show()
The plot shows the evolution of the PAM/SHE field at one fixed spatial point. In this mild parameter regime, the field fluctuates around (1). These fluctuations will be transformed into a time-dependent variance path in the next section.
12. PAM-driven volatility¶
We now use the extracted PAM time series to define a time-dependent variance path.
A mild choice would be
$$ v_t = v_{\text{base}}u(t,x_0), $$
where
$$ v_{\text{base}} = 0.2^2. $$
However, this often produces only small fluctuations. To create a more visibly intermittent volatility path, we use an amplified transformation:
$$ v_t = v_{\text{base}}\exp\left(\gamma(u(t,x_0)-1)\right). $$
The parameter (\gamma) controls how strongly PAM fluctuations are amplified into volatility fluctuations.
base_variance = 0.2 ** 2
gamma = 20
pam_volatility_path = []
for value in pam_time_series[:-1]:
variance = base_variance * math.exp(gamma * (value - 1.0))
pam_volatility_path.append(variance)
len(pam_volatility_path)
10000
variance_times = [m * (T / n_time) for m in range(n_time)]
plt.figure(figsize=(8, 5))
plt.plot(variance_times, pam_volatility_path)
plt.xlabel("Time")
plt.ylabel("Variance")
plt.title("Amplified PAM-driven variance path")
plt.show()
The amplified PAM-driven variance path transforms small fluctuations in the PAM field into larger volatility fluctuations. This gives a toy mechanism for intermittent volatility: calm periods are interrupted by periods of higher variance.
13. Stock paths with PAM-driven volatility¶
We now plug the PAM-driven variance path into the stock price simulation.
The stock price model is
$$ dS_t = rS_t\,dt + \sqrt{v_t}S_t\,dB_t, $$
where the variance path (v_t) is generated from the toy PAM/SHE field.
This creates a stock model where volatility is random and time-dependent.
random.seed(123)
stock_path_pam = simulate_stock_path_with_volatility_path(
S0=100,
T=T,
r=0.05,
volatility_path=pam_volatility_path
)
stock_times = [m * (T / n_time) for m in range(n_time + 1)]
plt.figure(figsize=(8, 5))
plt.plot(stock_times, stock_path_pam)
plt.xlabel("Time")
plt.ylabel("Stock price")
plt.title("Stock path with PAM-driven volatility")
plt.show()
This plot shows one possible stock path driven by the PAM-generated variance path. The PAM field does not directly push the stock up or down. Instead, it changes the local volatility, making the stock path more or less random at different times.
14. Endpoint distributions: constant vs PAM-driven volatility¶
We now compare the terminal stock price distribution under two models:
- constant volatility,
- PAM-driven volatility.
For each model, we simulate many stock paths and record the final value (S_T). This lets us see whether PAM-driven volatility changes the shape of the endpoint distribution.
random.seed(42)
constant_variance_path = [0.2 ** 2] * n_time
constant_endpoints = simulate_endpoint_distribution_with_volatility_path(
S0=100,
T=T,
r=0.05,
volatility_path=constant_variance_path,
n_paths=500
)
random.seed(42)
pam_endpoints = simulate_endpoint_distribution_with_volatility_path(
S0=100,
T=T,
r=0.05,
volatility_path=pam_volatility_path,
n_paths=500
)
plt.figure(figsize=(8, 5))
plt.hist(constant_endpoints, bins=40, alpha=0.5, label="Constant volatility")
plt.hist(pam_endpoints, bins=40, alpha=0.5, label="PAM-driven volatility")
plt.xlabel("Terminal stock price $S_T$")
plt.ylabel("Frequency")
plt.title("Endpoint distributions: constant vs PAM-driven volatility")
plt.legend()
plt.show()
The PAM-driven volatility model changes the shape of the endpoint distribution. In particular, intermittent volatility can create more extreme endpoints: some paths finish lower, while others finish much higher. This affects option prices because call options are sensitive to the right tail of the terminal stock price distribution.
15. Option price comparison: constant vs PAM-driven volatility¶
We now compare European call option prices under constant volatility and PAM-driven volatility.
For each strike (K), we estimate
$$ C_0 \approx e^{-rT}\frac{1}{N}\sum_{j=1}^N (S_T^{(j)}-K)^+. $$
This lets us see how the PAM-driven endpoint distribution changes option prices across different strikes.
strikes = [80, 100, 120, 150, 200]
print(f"{'K':>10} {'Constant Vol Price':>20} {'PAM Vol Price':>20} {'Difference':>15}")
print("-" * 75)
for K in strikes:
constant_price = call_price_from_endpoints(
endpoints=constant_endpoints,
K=K,
r=0.05,
T=1
)
pam_price = call_price_from_endpoints(
endpoints=pam_endpoints,
K=K,
r=0.05,
T=1
)
difference = pam_price - constant_price
print(f"{K:>10} {constant_price:>20.4f} {pam_price:>20.4f} {difference:>15.4f}")
K Constant Vol Price PAM Vol Price Difference
---------------------------------------------------------------------------
80 23.3911 24.2122 0.8211
100 9.6102 11.5104 1.9002
120 2.8200 4.4477 1.6278
150 0.2596 0.7991 0.5395
200 0.0000 0.0400 0.0400
print(f"{'K':>10} {'P_const(ST>K)':>20} {'P_PAM(ST>K)':>20}")
print("-" * 55)
for K in strikes:
p_constant = sum(ST > K for ST in constant_endpoints) / len(constant_endpoints)
p_pam = sum(ST > K for ST in pam_endpoints) / len(pam_endpoints)
print(f"{K:>10} {p_constant:>20.4f} {p_pam:>20.4f}")
The PAM-driven volatility model changes option prices by changing the shape of the terminal stock price distribution. A useful decomposition is
$$ \mathbb E[(S_T-K)^+] = \mathbb P(S_T>K) \mathbb E[S_T-K \mid S_T>K]. $$
Thus, a call price depends both on how often paths finish in the money and on how large the payoff is conditional on finishing in the money. PAM-driven volatility can reduce the number of ordinary in-the-money paths while increasing the size of rare extreme payoffs.
16. Interpretation and next steps¶
This notebook studied a toy connection between PAM/SHE-type intermittency and volatility modelling in option pricing.
We began with the Black-Scholes model, where volatility is constant:
$$ dS_t = rS_t\,dt + \sigma S_t\,dB_t. $$
We then replaced constant volatility with a time-dependent variance path:
$$ dS_t = rS_t\,dt + \sqrt{v_t}S_t\,dB_t. $$
First, we considered hand-made volatility spikes. These showed that volatility spikes widen the terminal stock price distribution and can increase call option prices, especially for higher strikes.
We then generated a toy PAM/SHE field
$$ \partial_t u = \frac12 \Delta u + \beta u \xi, $$
and extracted a time series (u(t,x_0)) at one fixed spatial point. This was transformed into a variance path using
$$ v_t = v_{\text{base}}\exp\left(\gamma(u(t,x_0)-1)\right). $$
This produced a toy PAM-driven volatility model.
The main qualitative conclusion is that PAM-driven volatility changes the shape of the endpoint distribution. It can create a more intermittent distribution: more extreme low outcomes and more extreme high outcomes.
For call options, the key identity is
$$ \mathbb E[(S_T-K)^+] = \mathbb P(S_T>K) \mathbb E[S_T-K \mid S_T>K]. $$
This shows that option prices depend on two effects:
- how often paths finish in the money,
- how large the payoff is conditional on finishing in the money.
In the simulations, PAM-driven volatility could reduce the number of ordinary in-the-money paths while increasing the size of rare extreme payoffs. This gives a possible mechanism by which intermittent volatility affects option prices differently across strikes.
Next steps¶
Several improvements would make the project more realistic:
- Use a more stable and principled numerical scheme for the PAM/SHE.
- Compare different values of $\beta$ and $\gamma$.
- Study barrier options, where volatility spikes may reduce prices by knocking paths out.
- Compare the simulated endpoint distributions with real stock return data.
- Calibrate the baseline volatility using historical data.
- Compare constant-volatility, spike-volatility, and PAM-driven-volatility option prices across many strikes.
The current notebook is therefore a proof of concept: it shows how stochastic PDE ideas can be used to generate intermittent volatility and study its effect on option prices.