期权波动率交易实战：隐含波动率曲面与Delta中性策略

import numpy as np
from scipy.optimize import brentq
from scipy.stats import norm

def implied_volatility(option_price, S, K, T, r, q=0, option_type='call'):
    """
    通过Newton-Raphson方法计算隐含波动率
    """
    def bs_price(sigma):
        d1 = (np.log(S/K) + (r - q + 0.5*sigma**2)*T) / (sigma*np.sqrt(T))
        d2 = d1 - sigma*np.sqrt(T)
        if option_type == 'call':
            return S*np.exp(-q*T)*norm.cdf(d1) - K*np.exp(-r*T)*norm.cdf(d2)
        else:
            return K*np.exp(-r*T)*norm.cdf(-d2) - S*np.exp(-q*T)*norm.cdf(-d1)
    
    # 使用Brent方法求解
    try:
        iv = brentq(lambda x: bs_price(x) - option_price, 0.01, 5.0)
        return iv
    except:
        return np.nan

# 示例：计算茅台期权的隐含波动率
option_market_price = 15.5  # 期权市场价格
underlying_price = 1800  # 标的股价
strike = 1850  # 行权价
time_to_expiry = 30/365  # 30天到期
risk_free_rate = 0.025  # 无风险利率2.5%

iv = implied_volatility(option_market_price, underlying_price, strike, 
                        time_to_expiry, risk_free_rate, option_type='call')
print(f"隐含波动率: {iv:.2%}")

import plotly.graph_objects as go
import pandas as pd

def plot_volatility_surface(options_data):
    """
    绘制波动率曲面
    options_data: DataFrame with columns [strike, expiry, implied_vol]
    """
    # 构建网格
    K_range = np.linspace(options_data['strike'].min(), 
                          options_data['strike'].max(), 50)
    T_range = np.linspace(options_data['expiry'].min(), 
                          options_data['expiry'].max(), 50)
    K_grid, T_grid = np.meshgrid(K_range, T_range)
    
    # 插值得到曲面
    from scipy.interpolate import griddata
    points = options_data[['strike', 'expiry']].values
    values = options_data['implied_vol'].values
    IV_grid = griddata(points, values, (K_grid, T_grid), method='cubic')
    
    # 绘制3D曲面
    fig = go.Figure(data=[go.Surface(z=IV_grid, x=K_grid, y=T_grid)])
    fig.update_layout(
        title='期权隐含波动率曲面',
        scene=dict(
            xaxis_title='行权价 (Strike)',
            yaxis_title='到期时间 (Time to Expiry)',
            zaxis_title='隐含波动率 (IV)'
        )
    )
    fig.show()
    
    return fig

# 示例：绘制上证50ETF期权波动率曲面
options_chain = pd.read_csv('sse50_etf_options.csv')
fig = plot_volatility_surface(options_chain)
fig.write_html('vol_surface.html')

def smile_arbitrage(options_chain, S, T):
    """
    波动率微笑套利策略
    """
    # 计算每个行权价的理论IV（基于BSM模型）
    options_chain['theoretical_iv'] = options_chain.apply(
        lambda row: theoretical_smile_iv(S, row['strike'], T), axis=1
    )
    
    # 找出IV偏差最大的期权
    options_chain['iv_mispricing'] = options_chain['implied_vol'] - options_chain['theoretical_iv']
    
    # 交易信号
    long_iv = options_chain[options_chain['iv_mispricing'] < -0.05]  # IV低估
    short_iv = options_chain[options_chain['iv_mispricing'] > 0.05]  # IV高估
    
    return long_iv, short_iv

# 示例：找出沪深300期权波动率微笑套利机会
underlying_price = 4100
days_to_expiry = 30
options = load_options_chain('hs300_etf', days_to_expiry)
underpriced, overpriced = smile_arbitrage(options, underlying_price, days_to_expiry/365)

print("低估期权（买入）：")
print(underpriced[['strike', 'implied_vol', 'iv_mispricing']])
print("\n高估期权（卖出）：")
print(overpriced[['strike', 'implied_vol', 'iv_mispricing']])

def calendar_spread_arbitrage(options_chain):
    """
    跨期波动率套利：卖出高IV期权，买入低IV期权
    """
    # 按到期日分组计算平均IV
    iv_term_structure = options_chain.groupby('expiry')['implied_vol'].mean()
    
    # 找出期限结构异常点
    iv_diff = iv_term_structure.diff().abs()
    anomaly_expiries = iv_diff[iv_diff > 0.1]  # IV跳变超过10%
    
    trades = []
    for exp in anomaly_expiries.index:
        near_term = iv_term_structure[iv_term_structure.index < exp].iloc[-1]
        far_term = iv_term_structure[iv_term_structure.index > exp].iloc[0]
        
        if near_term < far_term:
            # 做多近期IV，做空远期IV
            trades.append({
                'buy_expiry': near_term.name,
                'sell_expiry': far_term.name,
                'spread': far_term - near_term
            })
    
    return trades

def backtest_straddle_strategy(underlying, options_chain, entry_date, expiry_date):
    """
    回测跨式组合策略
    """
    # 1. 选择ATM期权（行权价最接近当前股价）
    S = underlying.loc[entry_date, 'close']
    atm_strike = options_chain.loc[entry_date]['strike'].iloc[
        (options_chain.loc[entry_date]['strike'] - S).abs().argsort()[0]
    ]
    
    # 2. 计算初始成本
    call_price = options_chain.loc[entry_date][
        (options_chain['strike'] == atm_strike) & 
        (options_chain['type'] == 'call')
    ]['price'].values[0]
    
    put_price = options_chain.loc[entry_date][
        (options_chain['strike'] == atm_strike) & 
        (options_chain['type'] == 'put')
    ]['price'].values[0]
    
    initial_cost = call_price + put_price
    
    # 3. 持有至到期 or 动态对冲
    results = []
    for date in pd.date_range(entry_date, expiry_date, freq='D'):
        if date in underlying.index:
            S_t = underlying.loc[date, 'close']
            
            # 计算当前组合价值
            call_value = black_scholes_call(S_t, atm_strike, 
                                           (expiry_date - date).days/365, 
                                           risk_free_rate, implied_vol)
            put_value = black_scholes_put(S_t, atm_strike, 
                                          (expiry_date - date).days/365, 
                                          risk_free_rate, implied_vol)
            
            portfolio_value = call_value + put_value
            pnl = portfolio_value - initial_cost
            
            results.append({
                'date': date,
                'underlying_price': S_t,
                'portfolio_value': portfolio_value,
                'pnl': pnl
            })
    
    return pd.DataFrame(results)

# 示例：回测沪深300ETF跨式组合
underlying_data = pd.read_csv('hs300_etf_daily.csv', index_col='date', parse_dates=True)
options_data = pd.read_csv('hs300_options_chain.csv', index_col='date', parse_dates=True)

results = backtest_straddle_strategy(
    underlying_data, options_data, 
    entry_date='2026-01-15',
    expiry_date='2026-02-15'
)

print(f"策略总收益: {results['pnl'].iloc[-1]:.2f}")
print(f"最大回撤: {results['pnl'].min():.2f}")

def delta_hedging_simulation(S0, K, T, r, sigma, option_type='call', 
                            hedge_freq='daily', transaction_cost=0.001):
    """
    Delta动态对冲模拟
    """
    dt = {'daily': 1/252, 'weekly': 5/252, 'monthly': 22/252}[hedge_freq]
    n_steps = int(T / dt)
    
    # 模拟标的资产路径（几何布朗运动）
    np.random.seed(42)
    S_path = [S0]
    for t in range(1, n_steps + 1):
        dS = S_path[-1] * (r*dt + sigma*np.sqrt(dt)*np.random.randn())
        S_path.append(S_path[-1] + dS)
    
    # 动态对冲
    portfolio_value = []
    hedge_error = []
    transaction_costs = []
    
    delta_prev = 0
    for t, S in enumerate(S_path):
        # 计算当前期权Delta
        d1 = (np.log(S/K) + (r + 0.5*sigma**2)*(T - t*dt)) / (sigma*np.sqrt(T - t*dt))
        delta = norm.cdf(d1) if option_type == 'call' else norm.cdf(d1) - 1
        
        # 调整标的持仓（买入/卖出股票）
        shares_to_trade = delta - delta_prev
        cost = abs(shares_to_trade * S * transaction_cost)
        transaction_costs.append(cost)
        
        # 计算组合价值
        option_value = black_scholes_option(S, K, T - t*dt, r, sigma, option_type)
        portfolio = option_value - delta * S  # 空头对冲组合
        portfolio_value.append(portfolio)
        
        # 对冲误差
        if t > 0:
            hedge_error.append(abs(portfolio - portfolio_value[0]))
        
        delta_prev = delta
    
    return {
        'final_pnl': portfolio_value[-1] - portfolio_value[0],
        'mean_hedge_error': np.mean(hedge_error),
        'total_transaction_cost': sum(transaction_costs),
        'sharpe_ratio': np.mean(portfolio_value) / np.std(portfolio_value) * np.sqrt(252)
    }

# 示例：比较不同对冲频率的效果
results = {}
for freq in ['daily', 'weekly', 'monthly']:
    stats = delta_hedging_simulation(
        S0=100, K=100, T=1, r=0.025, sigma=0.25,
        option_type='call', hedge_freq=freq, transaction_cost=0.001
    )
    results[freq] = stats

print("Delta对冲效果对比：")
print(pd.DataFrame(results).T)

def gamma_scalping_strategy(underlying, options_chain, entry_date, holding_days=22):
    """
    Gamma Scalping策略
    """
    # 1. 选择高Gamma期权（短期ATM）
    S = underlying.loc[entry_date, 'close']
    atm_option = options_chain.loc[entry_date][
        (options_chain['strike'] - S).abs().argsort()[0]
    ]
    
    # 2. 计算初始Greek
    initial_delta = atm_option['delta']
    initial_gamma = atm_option['gamma']
    
    # 3. 动态再平衡
    portfolio = []
    shares_held = 0
    
    for i, date in enumerate(pd.date_range(entry_date, periods=holding_days, freq='D')):
        if date not in underlying.index:
            continue
            
        S_t = underlying.loc[date, 'close']
        
        # 计算当前Delta
        delta_t = calculate_delta(S_t, atm_option['strike'], 
                                  (atm_option['expiry'] - date).days/365, 
                                  risk_free_rate, atm_option['implied_vol'])
        
        # 再平衡：买入/卖出标的使Delta中性
        target_shares = -delta_t * 100  # 假设1手期权对应100股
        shares_to_trade = target_shares - shares_held
        
        if abs(shares_to_trade) > 0:
            # 记录交易
            portfolio.append({
                'date': date,
                'action': 'buy' if shares_to_trade > 0 else 'sell',
                'shares': abs(shares_to_trade),
                'price': S_t,
                'cost': shares_to_trade * S_t
            })
            shares_held = target_shares
        
        # 计算当日组合价值
        option_value = calculate_option_value(S_t, atm_option)
        portfolio_value = option_value + shares_held * S_t
        
        portfolio.append({
            'date': date,
            'option_value': option_value,
            'shares_held': shares_held,
            'portfolio_value': portfolio_value
        })
    
    return pd.DataFrame(portfolio)

# 示例：在波动率飙升期间运行Gamma Scalping
high_vol_period_start = '2026-02-01'  # 假设这段时间波动率飙升
results = gamma_scalping_strategy(
    underlying_data, options_data, 
    entry_date=high_vol_period_start, 
    holding_days=22
)

print("Gamma Scalping策略表现：")
print(results[['date', 'option_value', 'portfolio_value']].tail(10))

def vega_risk_analysis(options_portfolio, vol_scenarios):
    """
    Vega风险分析：模拟不同波动率情景下的组合价值
    """
    results = []
    
    for vol_change in vol_scenarios:
        portfolio_value = 0
        for option in options_portfolio:
            # 计算Vega暴露
            vega = option['vega']
            price_change = vega * vol_change
            new_option_value = option['current_price'] + price_change
            portfolio_value += new_option_value * option['position']
        
        results.append({
            'vol_change': vol_change,
            'portfolio_value': portfolio_value,
            'pnl': portfolio_value - initial_portfolio_value
        })
    
    return pd.DataFrame(results)

# 示例：分析波动率±20%情景下的组合风险
portfolio = [
    {'type': 'call', 'strike': 100, 'vega': 0.35, 'current_price': 8.5, 'position': 10},
    {'type': 'put', 'strike': 95, 'vega': 0.28, 'current_price': 5.2, 'position': -10}  # 做空
]

vol_scenarios = np.arange(-0.5, 0.51, 0.05)  # IV变化±50%
risk_report = vega_risk_analysis(portfolio, vol_scenarios)

print("Vega风险分析：")
print(risk_report[['vol_change', 'pnl']].to_string(index=False))

def theta_decay_curve(S, K, T_range, r, sigma, option_type='call'):
    """
    绘制Theta衰减曲线
    """
    theta_values = []
    option_values = []
    
    for T in T_range:
        # 计算期权价格
        option_price = black_scholes_option(S, K, T, r, sigma, option_type)
        
        # 计算Theta（数值微分）
        dT = 1/365  # 1天
        option_price_next = black_scholes_option(S, K, T - dT, r, sigma, option_type)
        theta = (option_price_next - option_price) / dT
        
        theta_values.append(theta)
        option_values.append(option_price)
    
    # 绘图
    fig, ax1 = plt.subplots(figsize=(10, 6))
    
    ax1.plot(T_range, option_values, 'b-', label='Option Price')
    ax1.set_xlabel('Time to Expiry (Years)')
    ax1.set_ylabel('Option Price', color='b')
    ax1.tick_params(axis='y', labelcolor='b')
    
    ax2 = ax1.twinx()
    ax2.plot(T_range, theta_values, 'r--', label='Theta')
    ax2.set_ylabel('Theta (Decay per Day)', color='r')
    ax2.tick_params(axis='y', labelcolor='r')
    
    plt.title('Option Theta Decay Curve')
    fig.tight_layout()
    plt.savefig('theta_decay.png', dpi=300, bbox_inches='tight')
    
    return fig

# 示例：绘制ATM期权的Theta衰减曲线
T_range = np.linspace(1, 0, 252)  # 从1年到到期
fig = theta_decay_curve(S=100, K=100, T_range=T_range, r=0.025, sigma=0.25)