import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm, weibull_min, lognorm, genextreme
from scipy.interpolate import CubicSpline, interp1d, UnivariateSpline
from scipy.optimize import minimize, minimize_scalar, root_scalar, brentq
from scipy.special import gamma, gammainc, gammaincc
from scipy.ndimage import gaussian_filter1d
import warnings
from dataclasses import dataclass
from typing import Dict, List, Tuple, Optional, Union
import os
from datetime import datetime, timedelta
warnings.filterwarnings('ignore')
# Set matplotlib style to match paper
plt.rcParams.update({
'font.family': ['DejaVu Sans', 'SimHei'],
'font.size': 10,
'axes.labelsize': 10,
'axes.titlesize': 11,
'legend.fontsize': 9,
'xtick.labelsize': 9,
'ytick.labelsize': 9,
'figure.dpi': 100,
'savefig.dpi': 300,
'axes.linewidth': 0.8,
'grid.alpha': 0.3,
'grid.linewidth': 0.5,
'axes.grid': True,
'figure.figsize': (12, 8)
})
@dataclass
class OptionData:
"""Option data structure matching paper requirements"""
strikes: np.ndarray
call_prices: np.ndarray
put_prices: np.ndarray
implied_vols: np.ndarray
call_deltas: np.ndarray
vegas: np.ndarray
F: float # Futures/forward price
r: float # Risk-free rate
T: float # Time to expiry
sigma_atm: float = 0.25
# Additional fields for paper replication
all_strikes: np.ndarray = None
all_call_prices: np.ndarray = None
all_put_prices: np.ndarray = None
in_sample_mask: np.ndarray = None
class BlackScholes:
"""Black-Scholes pricing engine"""
@staticmethod
def call_price(S, K, r, T, sigma):
"""Black (1976) futures option pricing"""
if T <= 0 or sigma <= 0:
return max(S - K * np.exp(-r * T), 0)
d1 = (np.log(S/K) + 0.5 * sigma**2 * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
price = np.exp(-r * T) * (S * norm.cdf(d1) - K * norm.cdf(d2))
return max(price, 1e-6)
@staticmethod
def put_price(S, K, r, T, sigma):
"""Put option pricing"""
if T <= 0 or sigma <= 0:
return max(K * np.exp(-r * T) - S, 0)
d1 = (np.log(S/K) + 0.5 * sigma**2 * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
price = np.exp(-r * T) * (K * norm.cdf(-d2) - S * norm.cdf(-d1))
return max(price, 1e-6)
@staticmethod
def implied_vol(price, S, K, r, T, option_type='call'):
"""Calculate implied volatility"""
if price <= 0 or T <= 0:
return 0.25
def objective(sigma):
if option_type == 'call':
model_price = BlackScholes.call_price(S, K, r, T, sigma)
else:
model_price = BlackScholes.put_price(S, K, r, T, sigma)
return abs(model_price - price)
try:
result = minimize_scalar(objective, bounds=(0.01, 5.0), method='bounded')
return max(min(result.x, 5.0), 0.01)
except:
return 0.25
@staticmethod
def delta(S, K, r, T, sigma, option_type='call'):
"""Calculate option delta"""
if T <= 0:
return 1.0 if (option_type == 'call' and S > K) else 0.0
d1 = (np.log(S/K) + 0.5 * sigma**2 * T) / (sigma * np.sqrt(T))
if option_type == 'call':
return np.exp(-r * T) * norm.cdf(d1)
else:
return np.exp(-r * T) * (norm.cdf(d1) - 1)
@staticmethod
def vega(S, K, r, T, sigma):
"""Calculate option vega"""
if T <= 0:
return 0.0
d1 = (np.log(S/K) + 0.5 * sigma**2 * T) / (sigma * np.sqrt(T))
return np.exp(-r * T) * S * norm.pdf(d1) * np.sqrt(T)
def load_50etf_data(filename='50ETF期权2025.csv'):
"""Load and process 50ETF option data"""
print("\n" + "="*80)
print("Loading 50ETF Option Data")
print("="*80)
# Try different encodings
for encoding in ['utf-8', 'gbk', 'gb18030', 'utf-8-sig']:
try:
df = pd.read_csv(filename, encoding=encoding)
print(f"Successfully loaded data with encoding: {encoding}")
break
except:
continue
print(f"Data shape: {df.shape}")
print(f"Columns: {list(df.columns)[:10]}...") # Show first 10 columns
# Process numeric columns
numeric_cols = ['行权价', '收盘价', '结算价', 'Delta', 'Gamma', 'Vega', 'Theta']
for col in numeric_cols:
if col in df.columns:
df[col] = pd.to_numeric(df[col].astype(str).str.replace(',', ''), errors='coerce')
# Identify option type from name
df['option_type'] = df['期权简称'].apply(lambda x: 'call' if '购' in str(x) else 'put')
# Filter valid data
df = df[(df['行权价'] > 0) & (df['收盘价'] > 0.0001) & df['option_type'].notna()]
# Calculate underlying price from put-call parity
F = estimate_underlying_price(df)
print(f"Estimated underlying price: {F:.4f}")
print(f"Strike range: {df['行权价'].min():.4f} - {df['行权价'].max():.4f}")
print(f"Valid options: {len(df)}")
return df, F
def estimate_underlying_price(df):
"""Estimate underlying price using put-call parity"""
r = 0.025
T = 5/365 # Approximate 5 days to expiry
estimates = []
strikes = df['行权价'].unique()
for K in strikes[:20]: # Use first 20 strikes
calls = df[(df['行权价'] == K) & (df['option_type'] == 'call')]
puts = df[(df['行权价'] == K) & (df['option_type'] == 'put')]
if len(calls) > 0 and len(puts) > 0:
C = calls['收盘价'].mean()
P = puts['收盘价'].mean()
if C > 0.001 and P > 0.001:
F_est = (C - P) * np.exp(r * T) + K
if 4.0 < F_est < 10.0: # Reasonable range
estimates.append(F_est)
if estimates:
return np.median(estimates)
else:
# Fallback: use ATM strikes
return df['行权价'].median()
def create_option_data(df, F):
"""Create OptionData object from dataframe"""
print("\nCreating OptionData object...")
bs = BlackScholes()
r = 0.025
T = 5/365 # 5 days to expiry
# Get unique strikes and sort
strikes = sorted(df['行权价'].unique())
# Process each strike
processed_data = []
for K in strikes:
call_data = df[(df['行权价'] == K) & (df['option_type'] == 'call')]
put_data = df[(df['行权价'] == K) & (df['option_type'] == 'put')]
# Get prices
if len(call_data) > 0:
call_price = call_data['收盘价'].mean()
else:
if len(put_data) > 0:
put_price = put_data['收盘价'].mean()
call_price = put_price + F - K * np.exp(-r * T)
else:
continue
if len(put_data) > 0:
put_price = put_data['收盘价'].mean()
else:
put_price = call_price - F + K * np.exp(-r * T)
# Ensure positive prices
call_price = max(call_price, 1e-4)
put_price = max(put_price, 1e-4)
# Calculate implied vol and Greeks
iv = bs.implied_vol(call_price, F, K, r, T, 'call')
delta = bs.delta(F, K, r, T, iv, 'call')
vega = bs.vega(F, K, r, T, iv)
processed_data.append({
'strike': K,
'call_price': call_price,
'put_price': put_price,
'iv': iv,
'delta': delta,
'vega': max(vega, 1e-6)
})
# Convert to arrays
processed_df = pd.DataFrame(processed_data)
# Select in-sample options (relaxed delta range for 50ETF)
in_sample = (processed_df['delta'] >= 0.10) & (processed_df['delta'] <= 0.90)
if in_sample.sum() < 5:
# If too few options, use all with reasonable delta
in_sample = (processed_df['delta'] >= 0.01) & (processed_df['delta'] <= 0.99)
sample_data = processed_df[in_sample]
# Find ATM volatility
atm_idx = np.argmin(np.abs(processed_df['strike'] - F))
sigma_atm = processed_df.iloc[atm_idx]['iv']
option_data = OptionData(
strikes=sample_data['strike'].values,
call_prices=sample_data['call_price'].values,
put_prices=sample_data['put_price'].values,
implied_vols=sample_data['iv'].values,
call_deltas=sample_data['delta'].values,
vegas=sample_data['vega'].values,
F=F,
r=r,
T=T,
sigma_atm=sigma_atm,
all_strikes=processed_df['strike'].values,
all_call_prices=processed_df['call_price'].values,
all_put_prices=processed_df['put_price'].values,
in_sample_mask=in_sample.values
)
print(f"In-sample options: {len(option_data.strikes)}")
print(f"Strike range: [{option_data.strikes.min():.3f}, {option_data.strikes.max():.3f}]")
print(f"Delta range: [{option_data.call_deltas.min():.3f}, {option_data.call_deltas.max():.3f}]")
print(f"ATM volatility: {sigma_atm:.3f}")
return option_data
def recover_interior_rnd(option_data, n_points=1000):
"""Recover interior RND using Breeden-Litzenberger method"""
print("\nRecovering Interior RND...")
bs = BlackScholes()
# Step 1: Create spline in (delta, IV) space
# Sort by delta
sort_idx = np.argsort(option_data.call_deltas)
sorted_deltas = option_data.call_deltas[sort_idx]
sorted_ivs = option_data.implied_vols[sort_idx]
sorted_vegas = option_data.vegas[sort_idx]
# Weight by vega
weights = sorted_vegas / np.max(sorted_vegas)
weights = np.clip(weights, 0.1, 1.0)
# Fit cubic spline
try:
spline = UnivariateSpline(sorted_deltas, sorted_ivs, w=weights, s=0.01, k=3)
except:
spline = interp1d(sorted_deltas, sorted_ivs, kind='linear', fill_value='extrapolate')
# Step 2: Generate dense strike grid
K_min = option_data.strikes.min() * 0.95
K_max = option_data.strikes.max() * 1.05
dense_strikes = np.linspace(K_min, K_max, n_points)
# Step 3: Calculate call prices on dense grid
dense_calls = []
for K in dense_strikes:
# Convert strike to delta
sigma_atm = option_data.sigma_atm
d1 = (np.log(option_data.F/K) + 0.5*sigma_atm**2*option_data.T) / (sigma_atm*np.sqrt(option_data.T))
point_delta = np.exp(-option_data.r * option_data.T) * norm.cdf(d1)
# Get IV from spline
if sorted_deltas.min() <= point_delta <= sorted_deltas.max():
if callable(spline):
iv = float(spline(point_delta))
else:
iv = float(spline(point_delta))
else:
# Extrapolate carefully
if point_delta < sorted_deltas.min():
iv = sorted_ivs[0]
else:
iv = sorted_ivs[-1]
iv = np.clip(iv, 0.01, 5.0)
# Calculate call price
call_price = bs.call_price(option_data.F, K, option_data.r, option_data.T, iv)
dense_calls.append(call_price)
dense_calls = np.array(dense_calls)
# Step 4: Calculate RND using numerical derivatives
rnd = np.zeros_like(dense_strikes)
rncdf = np.zeros_like(dense_strikes)
# Use finite differences for second derivative
h = dense_strikes[1] - dense_strikes[0]
for i in range(1, len(dense_strikes)-1):
# Second derivative for RND
rnd[i] = np.exp(option_data.r * option_data.T) * \
(dense_calls[i+1] - 2*dense_calls[i] + dense_calls[i-1]) / (h**2)
# First derivative for CDF
rncdf[i] = 1.0 + np.exp(option_data.r * option_data.T) * \
(dense_calls[i+1] - dense_calls[i-1]) / (2*h)
# Handle boundaries
rnd[0] = rnd[1]
rnd[-1] = rnd[-2]
rncdf[0] = 0.0
rncdf[-1] = 1.0
# Clean up
rnd = np.maximum(rnd, 0)
rncdf = np.clip(rncdf, 0, 1)
# Normalize
integral = np.trapz(rnd, dense_strikes)
if integral > 0:
rnd = rnd / integral
print(f"Interior RND points: {len(dense_strikes)}")
print(f"RND integral: {np.trapz(rnd, dense_strikes):.4f}")
return dense_strikes, dense_calls, rncdf, rnd
class TailHAP:
"""TailHAP method with Weibull tails"""
def attach_tails(self, strikes, calls, rncdf, rnd, option_data):
"""Attach Weibull tails using HAP constraints"""
print("\nApplying TailHAP method...")
# Find attachment points
valid_idx = np.where((rnd > 1e-8) & (rncdf > 0.01) & (rncdf < 0.99))[0]
if len(valid_idx) < 4:
print("Warning: Few valid points for tail attachment")
return strikes, rnd, rncdf
left_idx = valid_idx[1] # Skip first point
right_idx = valid_idx[-2] # Skip last point
X_left = strikes[left_idx]
X_right = strikes[right_idx]
rnd_left = rnd[left_idx]
rnd_right = rnd[right_idx]
cdf_left = rncdf[left_idx]
cdf_right = rncdf[right_idx]
print(f"Left attachment: X={X_left:.3f}, RND={rnd_left:.6f}, CDF={cdf_left:.3f}")
print(f"Right attachment: X={X_right:.3f}, RND={rnd_right:.6f}, CDF={cdf_right:.3f}")
# Extended strike range
extended_strikes = np.linspace(0.5*option_data.F, 1.5*option_data.F, 2000)
extended_rnd = np.zeros_like(extended_strikes)
# Simple Weibull tails
# Left tail
for i, x in enumerate(extended_strikes):
if x < X_left:
# Exponential decay
if X_left > 0 and x > 0:
decay = 2.0 / option_data.F
extended_rnd[i] = rnd_left * np.exp(-decay * (X_left - x))
elif x > X_right:
# Right tail exponential decay
decay = 1.5 / option_data.F
extended_rnd[i] = rnd_right * np.exp(-decay * (x - X_right))
else:
# Interior: interpolate
idx = np.searchsorted(strikes, x)
if 0 < idx < len(strikes):
w = (x - strikes[idx-1]) / (strikes[idx] - strikes[idx-1])
extended_rnd[i] = rnd[idx-1] * (1-w) + rnd[idx] * w
# Normalize
integral = np.trapz(extended_rnd, extended_strikes)
if integral > 0:
extended_rnd = extended_rnd / integral
# Recalculate CDF
extended_cdf = np.zeros_like(extended_strikes)
for i in range(1, len(extended_strikes)):
extended_cdf[i] = extended_cdf[i-1] + \
0.5 * (extended_rnd[i] + extended_rnd[i-1]) * \
(extended_strikes[i] - extended_strikes[i-1])
print(f"TailHAP complete. RND integral: {integral:.4f}")
return extended_strikes, extended_rnd, extended_cdf
class BPMethod:
"""BP method with lognormal tails"""
def attach_tails(self, strikes, calls, rncdf, rnd, option_data):
"""Attach horizontal volatility extension (lognormal tails)"""
print("\nApplying BP method...")
# Use first and last implied volatilities
sigma_left = option_data.implied_vols[0]
sigma_right = option_data.implied_vols[-1]
# Extended strike range
extended_strikes = np.linspace(0.5*option_data.F, 1.5*option_data.F, 2000)
extended_rnd = np.zeros_like(extended_strikes)
# Find attachment points
valid_idx = np.where((rnd > 1e-8))[0]
if len(valid_idx) > 0:
left_idx = valid_idx[0]
right_idx = valid_idx[-1]
X_left = strikes[left_idx]
X_right = strikes[right_idx]
rnd_left = rnd[left_idx]
rnd_right = rnd[right_idx]
else:
X_left = strikes[0]
X_right = strikes[-1]
rnd_left = rnd_right = 1e-6
print(f"Left: X={X_left:.3f}, sigma={sigma_left:.3f}")
print(f"Right: X={X_right:.3f}, sigma={sigma_right:.3f}")
for i, x in enumerate(extended_strikes):
if x < X_left:
# Left lognormal tail
if x > 0:
mu = np.log(option_data.F) - 0.5 * sigma_left**2 * option_data.T
extended_rnd[i] = (1/(x * sigma_left * np.sqrt(2*np.pi*option_data.T))) * \
np.exp(-0.5*((np.log(x) - mu)/(sigma_left*np.sqrt(option_data.T)))**2)
extended_rnd[i] *= np.exp(-2*(X_left - x)/option_data.F) # Decay factor
elif x > X_right:
# Right lognormal tail
mu = np.log(option_data.F) - 0.5 * sigma_right**2 * option_data.T
extended_rnd[i] = (1/(x * sigma_right * np.sqrt(2*np.pi*option_data.T))) * \
np.exp(-0.5*((np.log(x) - mu)/(sigma_right*np.sqrt(option_data.T)))**2)
extended_rnd[i] *= np.exp(-2*(x - X_right)/option_data.F) # Decay factor
else:
# Interior
idx = np.searchsorted(strikes, x)
if 0 < idx < len(strikes):
w = (x - strikes[idx-1]) / (strikes[idx] - strikes[idx-1])
extended_rnd[i] = rnd[idx-1] * (1-w) + rnd[idx] * w
# Normalize
integral = np.trapz(extended_rnd, extended_strikes)
if integral > 0:
extended_rnd = extended_rnd / integral
print(f"BP complete. RND integral: {integral:.4f}")
return extended_strikes, extended_rnd, extended_rnd
class FBMethod:
"""Figlewski-Birru method with GEV tails"""
def attach_tails(self, strikes, calls, rncdf, rnd, option_data):
"""Attach GEV tails"""
print("\nApplying FB method...")
# Extended strike range
extended_strikes = np.linspace(0.5*option_data.F, 1.5*option_data.F, 2000)
extended_rnd = np.zeros_like(extended_strikes)
# Find percentile points
valid_idx = np.where((rnd > 1e-8) & (rncdf > 0.02) & (rncdf < 0.98))[0]
if len(valid_idx) > 4:
# Find points close to 5% and 95% percentiles
idx_5 = valid_idx[np.argmin(np.abs(rncdf[valid_idx] - 0.05))]
idx_95 = valid_idx[np.argmin(np.abs(rncdf[valid_idx] - 0.95))]
X_5 = strikes[idx_5]
X_95 = strikes[idx_95]
rnd_5 = rnd[idx_5]
rnd_95 = rnd[idx_95]
else:
X_5 = strikes[0]
X_95 = strikes[-1]
rnd_5 = rnd[0]
rnd_95 = rnd[-1]
print(f"Attachment points: X_5%={X_5:.3f}, X_95%={X_95:.3f}")
# Simplified GEV-like tails
for i, x in enumerate(extended_strikes):
if x < X_5:
# Left tail with power law decay
if x > 0:
xi = 0.2 # Shape parameter
scale = X_5 / 5
extended_rnd[i] = rnd_5 * (x / X_5)**(-1-xi) * np.exp(-(X_5-x)/scale)
elif x > X_95:
# Right tail
xi = 0.1
scale = X_95 / 3
extended_rnd[i] = rnd_95 * (X_95 / x)**(1+xi) * np.exp(-(x-X_95)/scale)
else:
# Interior
idx = np.searchsorted(strikes, x)
if 0 < idx < len(strikes):
w = (x - strikes[idx-1]) / (strikes[idx] - strikes[idx-1])
extended_rnd[i] = rnd[idx-1] * (1-w) + rnd[idx] * w
# Normalize
integral = np.trapz(extended_rnd, extended_strikes)
if integral > 0:
extended_rnd = extended_rnd / integral
print(f"FB complete. RND integral: {integral:.4f}")
return extended_strikes, extended_rnd, extended_rnd
def create_figure1_tailhap(option_data, strikes, rnd):
"""Create Figure 1 - TailHAP visualization"""
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))
# Top panel: RND and option prices
ax1_twin = ax1.twinx()
# Plot RND
ax1.plot(strikes, rnd, 'r-', linewidth=2, label='Risk Neutral Density')
# Mark attachment points (simplified)
valid_idx = np.where(rnd > 1e-8)[0]
if len(valid_idx) > 4:
left_attach = strikes[valid_idx[2]]
right_attach = strikes[valid_idx[-3]]
ax1.axvline(left_attach, color='g', linestyle='--', alpha=0.7, label='Left Attachment')
ax1.axvline(right_attach, color='g', linestyle='--', alpha=0.7, label='Right Attachment')
# Plot option prices
ax1_twin.scatter(option_data.strikes, option_data.call_prices,
marker='+', s=100, c='b', label='In-sample Options')
# Out-of-sample options
if hasattr(option_data, 'all_strikes'):
out_mask = ~option_data.in_sample_mask
if out_mask.any():
ax1_twin.scatter(option_data.all_strikes[out_mask],
option_data.all_call_prices[out_mask],
marker='o', s=30, c='k', alpha=0.5, label='Out-of-sample Options')
ax1.set_xlabel('Strike Price')
ax1.set_ylabel('Density', color='r')
ax1_twin.set_ylabel('Option Price', color='b')
integral = np.trapz(rnd, strikes)
mean = np.trapz(strikes * rnd, strikes) if integral > 0 else np.nan
ax1.set_title(f'TailHAP Risk Neutral Density\nIntegral={integral:.3f}, Mean={mean:.3f}')
ax1.legend(loc='upper left')
ax1_twin.legend(loc='upper right')
ax1.grid(True, alpha=0.3)
# Bottom panel: IV smile
vegas = option_data.vegas
sizes = 50 * vegas / np.max(vegas) + 10
scatter = ax2.scatter(option_data.call_deltas, option_data.implied_vols,
s=sizes, c=vegas, cmap='viridis', alpha=0.7,
edgecolors='black', linewidth=0.5)
# Fit spline
try:
sort_idx = np.argsort(option_data.call_deltas)
spline = UnivariateSpline(option_data.call_deltas[sort_idx],
option_data.implied_vols[sort_idx],
w=vegas[sort_idx]/np.max(vegas), s=0.01, k=3)
delta_grid = np.linspace(option_data.call_deltas.min(),
option_data.call_deltas.max(), 100)
vol_fit = spline(delta_grid)
ax2.plot(delta_grid, vol_fit, 'r-', linewidth=2, label='Cubic Spline Fit')
except:
pass
ax2.set_xlabel('Call Delta')
ax2.set_ylabel('Implied Volatility')
ax2.set_title('Spline Fit in (Call Delta, Implied Volatility) Space')
ax2.legend()
ax2.grid(True, alpha=0.3)
plt.colorbar(scatter, ax=ax2, label='Vega')
plt.tight_layout()
return fig
def create_figure2_bp(option_data, strikes, rnd):
"""Create Figure 2 - BP visualization"""
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))
# Top panel: RND
ax1_twin = ax1.twinx()
ax1.plot(strikes, rnd, 'r-', linewidth=2, label='Risk Neutral Density')
ax1_twin.scatter(option_data.strikes, option_data.call_prices,
marker='+', s=100, c='b', label='Options Used')
ax1.set_xlabel('Strike Price')
ax1.set_ylabel('Density', color='r')
ax1_twin.set_ylabel('Option Price', color='b')
integral = np.trapz(rnd, strikes)
mean = np.trapz(strikes * rnd, strikes) if integral > 0 else np.nan
ax1.set_title(f'BP Risk Neutral Density\nIntegral={integral:.3f}, Mean={mean:.3f}')
ax1.legend(loc='upper left')
ax1_twin.legend(loc='upper right')
ax1.grid(True, alpha=0.3)
# Bottom panel: IV in delta space
ax2.scatter(option_data.call_deltas, option_data.implied_vols,
s=50, c='blue', alpha=0.7, edgecolors='black', linewidth=0.5)
# Show horizontal extension
if len(option_data.implied_vols) > 0:
ax2.axhline(option_data.implied_vols[0], color='r', linestyle='--',
alpha=0.5, label='Left Extension')
ax2.axhline(option_data.implied_vols[-1], color='r', linestyle='--',
alpha=0.5, label='Right Extension')
ax2.set_xlabel('Call Delta')
ax2.set_ylabel('Implied Volatility')
ax2.set_title('Implied Volatility with Horizontal Extensions')
ax2.legend()
ax2.grid(True, alpha=0.3)
plt.tight_layout()
return fig
def create_figure3_fb(option_data, strikes, rnd):
"""Create Figure 3 - FB visualization"""
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))
# Top panel: RND
ax1_twin = ax1.twinx()
ax1.plot(strikes, rnd, 'r-', linewidth=2, label='Risk Neutral Density')
ax1_twin.scatter(option_data.strikes, option_data.call_prices,
marker='+', s=100, c='b', label='Options Used')
ax1.set_xlabel('Strike Price')
ax1.set_ylabel('Density', color='r')
ax1_twin.set_ylabel('Option Price', color='b')
integral = np.trapz(rnd, strikes)
mean = np.trapz(strikes * rnd, strikes) if integral > 0 else np.nan
ax1.set_title(f'FB Risk Neutral Density\nIntegral={integral:.3f}, Mean={mean:.3f}')
ax1.legend(loc='upper left')
ax1_twin.legend(loc='upper right')
ax1.grid(True, alpha=0.3)
# Bottom panel: IV smile in strike space
ax2.scatter(option_data.strikes, option_data.implied_vols,
s=50, c='blue', alpha=0.7, edgecolors='black', linewidth=0.5)
# Fourth order spline fit
try:
from scipy.interpolate import splrep, splev
tck = splrep(option_data.strikes, option_data.implied_vols, k=min(4, len(option_data.strikes)-1))
strike_grid = np.linspace(option_data.strikes.min(), option_data.strikes.max(), 100)
vol_fit = splev(strike_grid, tck)
ax2.plot(strike_grid, vol_fit, 'r-', linewidth=2, label='Fourth Order Spline')
except:
pass
ax2.set_xlabel('Strike Price')
ax2.set_ylabel('Implied Volatility')
ax2.set_title('Spline Fit in (Strike Price, Implied Volatility) Space')
ax2.legend()
ax2.grid(True, alpha=0.3)
plt.tight_layout()
return fig
def create_comparison_plot(results):
"""Create method comparison plot"""
fig, axes = plt.subplots(2, 2, figsize=(15, 12))
axes = axes.ravel()
colors = {'TailHAP': 'red', 'BP': 'blue', 'FB': 'green'}
for method, data in results.items():
color = colors.get(method, 'black')
strikes, rnd = data['strikes'], data['rnd']
# Plot on all subplots
for ax in axes:
ax.plot(strikes, rnd, color=color, linewidth=2, alpha=0.7, label=method)
# Customize each subplot
axes[0].set_title('RND Comparison - Full View')
axes[1].set_title('RND Comparison - Zoomed Center')
axes[2].set_title('RND Comparison - Log Scale')
axes[2].set_yscale('log')
axes[3].set_title('RND Comparison - Left Tail')
# Adjust x-limits for specific views
if results:
sample_strikes = list(results.values())[0]['strikes']
F = sample_strikes.mean()
axes[1].set_xlim(0.9*F, 1.1*F)
axes[3].set_xlim(sample_strikes.min(), F)
for ax in axes:
ax.set_xlabel('Strike Price')
ax.set_ylabel('Density')
ax.legend()
ax.grid(True, alpha=0.3)
plt.suptitle('Risk Neutral Density Methods Comparison', fontsize=16)
plt.tight_layout()
return fig
def calculate_statistics(results, option_data):
"""Calculate performance statistics for Table 5"""
stats = []
for method, data in results.items():
strikes, rnd = data['strikes'], data['rnd']
# Basic statistics
integral = np.trapz(rnd, strikes)
mean = np.trapz(strikes * rnd, strikes) if integral > 0 else np.nan
# Pricing errors (simplified)
futures_diff = mean - option_data.F if not np.isnan(mean) else np.nan
stats.append({
'Method': method,
'Integral': integral,
'Mean': mean,
'Futures_Price': option_data.F,
'Futures_Diff': futures_diff,
'Futures_Diff_Pct': 100*futures_diff/option_data.F if not np.isnan(futures_diff) else np.nan
})
return pd.DataFrame(stats)
def main():
"""Main execution function"""
print("\n" + "="*80)
print("RISK NEUTRAL DENSITY RECOVERY ANALYSIS")
print("Paper: Principled pasting - Bollinger, Melick & Thomas (2023)")
print("="*80)
# Load data
df, F = load_50etf_data('50ETF期权2025.csv')
# Create option data object
option_data = create_option_data(df, F)
# Recover interior RND
interior_strikes, interior_calls, interior_cdf, interior_rnd = recover_interior_rnd(option_data)
# Apply different tail methods
results = {}
# TailHAP method
tailhap = TailHAP()
th_strikes, th_rnd, th_cdf = tailhap.attach_tails(
interior_strikes, interior_calls, interior_cdf, interior_rnd, option_data)
results['TailHAP'] = {'strikes': th_strikes, 'rnd': th_rnd, 'cdf': th_cdf}
# BP method
bp = BPMethod()
bp_strikes, bp_rnd, bp_cdf = bp.attach_tails(
interior_strikes, interior_calls, interior_cdf, interior_rnd, option_data)
results['BP'] = {'strikes': bp_strikes, 'rnd': bp_rnd, 'cdf': bp_cdf}
# FB method
fb = FBMethod()
fb_strikes, fb_rnd, fb_cdf = fb.attach_tails(
interior_strikes, interior_calls, interior_cdf, interior_rnd, option_data)
results['FB'] = {'strikes': fb_strikes, 'rnd': fb_rnd, 'cdf': fb_cdf}
# Create figures
print("\n" + "="*80)
print("GENERATING FIGURES")
print("="*80)
# Figure 1: TailHAP
fig1 = create_figure1_tailhap(option_data, th_strikes, th_rnd)
plt.savefig('Figure1_TailHAP.png', dpi=300, bbox_inches='tight')
print("Figure 1 saved: Figure1_TailHAP.png")
# Figure 2: BP
fig2 = create_figure2_bp(option_data, bp_strikes, bp_rnd)
plt.savefig('Figure2_BP.png', dpi=300, bbox_inches='tight')
print("Figure 2 saved: Figure2_BP.png")
# Figure 3: FB
fig3 = create_figure3_fb(option_data, fb_strikes, fb_rnd)
plt.savefig('Figure3_FB.png', dpi=300, bbox_inches='tight')
print("Figure 3 saved: Figure3_FB.png")
# Comparison figure
fig_comp = create_comparison_plot(results)
plt.savefig('Methods_Comparison.png', dpi=300, bbox_inches='tight')
print("Comparison saved: Methods_Comparison.png")
# Generate statistics table
print("\n" + "="*80)
print("PERFORMANCE STATISTICS")
print("="*80)
stats_df = calculate_statistics(results, option_data)
print("\n", stats_df.round(4))
stats_df.to_csv('RND_Statistics.csv', index=False)
print("\nStatistics saved: RND_Statistics.csv")
print("\n" + "="*80)
print("ANALYSIS COMPLETE")
print("="*80)
# Show all figures
plt.show()
return results, option_data
if __name__ == "__main__":
results, option_data = main()请分析上述代码是否复现附件英文文献，并详细总结TailHAP方法、BP方法、FB方法模型的构建
最新发布
函数方法与jQuery在JavaScript中的应用
本文探讨了函数方法和jQuery在JavaScript编程中的高效应用，重点介绍了如何使用$.each方法处理数组，并通过条件判断返回true或false。讨论了如何在循环内部避免使用break和continue，提供了替代的实现策略。
扫一扫
30
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
