Model Gold Volatility Surfaces with GS-Quant in 20 Minutes

# Personal note: Learned this after hitting rate limits without proper auth
# Install gs-quant
!pip install gs-quant==3.1.2

from gs_quant.session import GsSession, Environment
from gs_quant.markets import PricingContext
from gs_quant.instrument import FXOption
from gs_quant.risk import MarketDataShockBasedScenario, MarketDataPattern
import pandas as pd
import numpy as np
from datetime import datetime, timedelta

# Authenticate (get credentials from developer.gs.com)
client_id = 'YOUR_CLIENT_ID'  # Replace with your ID
client_secret = 'YOUR_SECRET'  # Replace with your secret

# Watch out: Use Environment.PROD for real data, BETA for testing
GsSession.use(Environment.PROD, client_id=client_id, client_secret=client_secret)
print(f"✓ Connected to GS Markets API at {datetime.now().strftime('%H:%M:%S')}")

# Personal note: Gold options trade in $10 strike increments on COMEX
# Current gold spot price (example: $1,950/oz)
gold_spot = 1950.0

# Define strike range: ±20% around spot in $50 increments
strike_range = np.arange(
    gold_spot * 0.80,  # 80% of spot = $1,560
    gold_spot * 1.20,  # 120% of spot = $2,340
    50  # $50 strike increments
)

# Define tenors: 1W, 1M, 2M, 3M, 6M, 1Y
tenor_days = [7, 30, 60, 90, 180, 365]
expiration_dates = [
    (datetime.now() + timedelta(days=d)).strftime('%Y-%m-%d')
    for d in tenor_days
]

print(f"✓ Strike grid: {len(strike_range)} strikes from ${strike_range[0]:.0f} to ${strike_range[-1]:.0f}")
print(f"✓ Tenor grid: {len(tenor_days)} expirations up to {tenor_days[-1]} days")

# Watch out: GS-Quant uses different ticker formats for different exchanges
gold_ticker = 'GC'  # COMEX gold futures symbol

# Personal note: This took me 6 hours to figure out - GS uses CrossAsset class for commodities
from gs_quant.markets.securities import SecurityMaster, AssetIdentifier
from gs_quant.markets import HistoricalPricingContext
from gs_quant.data import Dataset

# Build the volatility surface dataframe
vol_surface_data = []

with PricingContext(pricing_date=datetime.now().date()):
    for expiry_date in expiration_dates:
        for strike in strike_range:
            try:
                # Create option object (American-style for gold)
                option = FXOption(
                    buy_sell='Buy',
                    option_type='Call',
                    pair='XAUUSD',  # Gold vs USD
                    strike_price=strike,
                    expiration_date=expiry_date,
                    notional_amount=100  # 100 oz
                )
                
                # Fetch implied vol (this is the magic line)
                implied_vol = option.calc_implied_volatility()
                
                vol_surface_data.append({
                    'Strike': strike,
                    'Expiry': expiry_date,
                    'ImpliedVol': implied_vol,
                    'Moneyness': strike / gold_spot,
                    'DaysToExpiry': (datetime.strptime(expiry_date, '%Y-%m-%d') - datetime.now()).days
                })
                
            except Exception as e:
                print(f"✗ Failed for strike ${strike}, expiry {expiry_date}: {e}")
                continue

# Convert to DataFrame
vol_df = pd.DataFrame(vol_surface_data)
print(f"✓ Fetched {len(vol_df)} implied vol quotes in {datetime.now().strftime('%H:%M:%S')}")
print(f"\nSample data:\n{vol_df.head(3)}")

from scipy.interpolate import griddata
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Create meshgrid for smooth surface
strike_grid = np.linspace(strike_range.min(), strike_range.max(), 50)
days_grid = np.linspace(min(tenor_days), max(tenor_days), 50)
strike_mesh, days_mesh = np.meshgrid(strike_grid, days_grid)

# Interpolate volatility surface using cubic method
vol_mesh = griddata(
    points=(vol_df['Strike'], vol_df['DaysToExpiry']),
    values=vol_df['ImpliedVol'],
    xi=(strike_mesh, days_mesh),
    method='cubic'  # Smooth surface, handles volatility smile
)

# Watch out: Check for NaN values at surface edges
valid_points = ~np.isnan(vol_mesh)
print(f"✓ Surface coverage: {valid_points.sum() / vol_mesh.size * 100:.1f}% valid points")

# Calculate surface statistics
print(f"\nVolatility Surface Stats:")
print(f"  Min IV: {np.nanmin(vol_mesh):.2%}")
print(f"  Max IV: {np.nanmax(vol_mesh):.2%}")
print(f"  ATM Vol (30d): {vol_df[vol_df['DaysToExpiry'] == 30].iloc[len(strike_range)//2]['ImpliedVol']:.2%}")

# Create 3D surface plot
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

# Plot surface
surf = ax.plot_surface(
    strike_mesh, days_mesh, vol_mesh * 100,  # Convert to percentage
    cmap='viridis',
    alpha=0.8,
    edgecolor='none'
)

# Overlay actual data points
ax.scatter(
    vol_df['Strike'], 
    vol_df['DaysToExpiry'], 
    vol_df['ImpliedVol'] * 100,
    c='red', 
    s=20, 
    label='Market quotes'
)

ax.set_xlabel('Strike Price ($)')
ax.set_ylabel('Days to Expiry')
ax.set_zlabel('Implied Volatility (%)')
ax.set_title(f'Gold Implied Volatility Surface\n{datetime.now().strftime("%Y-%m-%d %H:%M")} | Spot: ${gold_spot:.2f}')
ax.legend()

plt.colorbar(surf, shrink=0.5, aspect=5)
plt.savefig('gold_vol_surface.png', dpi=300, bbox_inches='tight')
print("✓ Saved visualization to gold_vol_surface.png")

# Export to CSV for production systems
export_df = vol_df.copy()
export_df['ImpliedVol_Pct'] = export_df['ImpliedVol'] * 100
export_df.to_csv('gold_ivol_surface.csv', index=False)
print(f"✓ Exported {len(export_df)} data points to gold_ivol_surface.csv")

# Export surface grid for interpolation
surface_export = pd.DataFrame({
    'Strike': strike_mesh.flatten(),
    'DaysToExpiry': days_mesh.flatten(),
    'ImpliedVol': vol_mesh.flatten()
}).dropna()
surface_export.to_csv('gold_ivol_surface_grid.csv', index=False)
print(f"✓ Exported {len(surface_export)} surface points to gold_ivol_surface_grid.csv")