Volatility Modeling: GARCH, Realized Volatility, and Implied Vol Surface

1import numpy as np
2import pandas as pd
3from scipy.optimize import minimize
4from numba import jit
5
6class GARCH11:
7    """
8    GARCH(1,1) model: σ²ₜ = ω + α·r²ₜ₋₁ + β·σ²ₜ₋₁
9    
10    where:
11        ω: long-run variance
12        α: reaction to new information
13        β: persistence of variance
14        α + β < 1 for stationarity
15    """
16    
17    def __init__(self):
18        self.omega = None
19        self.alpha = None
20        self.beta = None
21        self.fitted = False
22        
23    def fit(self, returns: np.ndarray, method='MLE') -> dict:
24        """
25        Fit GARCH(1,1) parameters using Maximum Likelihood Estimation.
26        
27        Args:
28            returns: Array of log returns
29            method: 'MLE' or 'QML' (Quasi-Maximum Likelihood)
30            
31        Returns:
32            Dictionary with parameters and fit statistics
33        """
34        # Initial parameter guess
35        # ω, α, β
36        x0 = np.array([
37            returns.var() * 0.1,  # ω
38            0.1,                   # α
39            0.8                    # β
40        ])
41        
42        # Bounds: all positive, α + β < 1
43        bounds = (
44            (1e-6, returns.var()),  # ω
45            (1e-6, 0.3),            # α
46            (1e-6, 0.99)            # β
47        )
48        
49        # Optimize
50        result = minimize(
51            lambda params: -self._log_likelihood(returns, params),
52            x0,
53            method='L-BFGS-B',
54            bounds=bounds,
55            options={'maxiter': 1000}
56        )
57        
58        if not result.success:
59            raise RuntimeError(f"GARCH optimization failed: {result.message}")
60        
61        self.omega, self.alpha, self.beta = result.x
62        self.fitted = True
63        
64        # Calculate fit statistics
65        log_likelihood = -result.fun
66        n = len(returns)
67        k = 3  # number of parameters
68        
69        aic = 2 * k - 2 * log_likelihood
70        bic = k * np.log(n) - 2 * log_likelihood
71        
72        return {
73            'omega': self.omega,
74            'alpha': self.alpha,
75            'beta': self.beta,
76            'persistence': self.alpha + self.beta,
77            'log_likelihood': log_likelihood,
78            'aic': aic,
79            'bic': bic,
80            'converged': result.success
81        }
82    
83    @staticmethod
84    @jit(nopython=True)
85    def _compute_variance_path(returns, omega, alpha, beta):
86        """Fast variance path computation using Numba."""
87        n = len(returns)
88        variance = np.zeros(n)
89        
90        # Initialize with sample variance
91        variance[0] = returns[:10].var()
92        
93        for t in range(1, n):
94            variance[t] = omega + alpha * returns[t-1]**2 + beta * variance[t-1]
95        
96        return variance
97    
98    def _log_likelihood(self, returns, params):
99        """Compute log-likelihood for parameter estimation."""
100        omega, alpha, beta = params
101        
102        # Constraint: α + β < 1
103        if alpha + beta >= 1.0:
104            return -1e10
105        
106        variance = self._compute_variance_path(returns, omega, alpha, beta)
107        
108        # Avoid log(0)
109        variance = np.maximum(variance, 1e-8)
110        
111        # Log-likelihood (assuming normal distribution)
112        ll = -0.5 * (np.log(2 * np.pi) + np.log(variance) + returns**2 / variance)
113        
114        return np.sum(ll)
115    
116    def forecast(self, returns: np.ndarray, horizon: int = 1) -> np.ndarray:
117        """
118        Forecast variance for multiple horizons.
119        
120        Args:
121            returns: Historical returns
122            horizon: Number of steps ahead to forecast
123            
124        Returns:
125            Array of variance forecasts
126        """
127        if not self.fitted:
128            raise RuntimeError("Model not fitted. Call fit() first.")
129        
130        # Current variance
131        variance = self._compute_variance_path(
132            returns, self.omega, self.alpha, self.beta
133        )
134        current_var = variance[-1]
135        
136        # Long-run variance
137        long_run_var = self.omega / (1 - self.alpha - self.beta)
138        
139        # Forecast recursively
140        forecasts = np.zeros(horizon)
141        
142        for h in range(horizon):
143            if h == 0:
144                # One-step ahead
145                forecasts[h] = (self.omega + 
146                               self.alpha * returns[-1]**2 + 
147                               self.beta * current_var)
148            else:
149                # Multi-step: mean-reverting to long-run variance
150                weight = (self.alpha + self.beta) ** h
151                forecasts[h] = (weight * current_var + 
152                               (1 - weight) * long_run_var)
153        
154        return forecasts
155    
156    def volatility_forecast(self, returns: np.ndarray, 
157                           horizon: int = 1,
158                           annualize: bool = True,
159                           periods_per_year: int = 252) -> np.ndarray:
160        """
161        Forecast volatility (std dev) instead of variance.
162        """
163        var_forecast = self.forecast(returns, horizon)
164        vol_forecast = np.sqrt(var_forecast)
165        
166        if annualize:
167            vol_forecast *= np.sqrt(periods_per_year)
168        
169        return vol_forecast
170
171# Example usage
172if __name__ == "__main__":
173    # Load returns
174    prices = pd.read_csv('spy_prices.csv')['close']
175    returns = np.log(prices / prices.shift(1)).dropna().values
176    
177    # Fit GARCH
178    garch = GARCH11()
179    params = garch.fit(returns)
180    
181    print(f"GARCH(1,1) Parameters:")
182    print(f"ω = {params['omega']:.6f}")
183    print(f"α = {params['alpha']:.6f}")
184    print(f"β = {params['beta']:.6f}")
185    print(f"Persistence (α+β) = {params['persistence']:.6f}")
186    print(f"Log-likelihood = {params['log_likelihood']:.2f}")
187    
188    # Forecast 10 days ahead
189    vol_forecast = garch.volatility_forecast(returns, horizon=10)
190    print(f"\nVolatility Forecast (annualized):")
191    for h, vol in enumerate(vol_forecast, 1):
192        print(f"  {h} day: {vol*100:.2f}%")
193

1class EGARCH:
2    """
3    EGARCH(1,1): log(σ²ₜ) = ω + β·log(σ²ₜ₋₁) + α·|zₜ₋₁| + γ·zₜ₋₁
4    
5    where zₜ = rₜ/σₜ (standardized return)
6    γ < 0 captures leverage effect
7    """
8    
9    def __init__(self):
10        self.omega = None
11        self.alpha = None
12        self.beta = None
13        self.gamma = None
14        self.fitted = False
15    
16    def fit(self, returns: np.ndarray) -> dict:
17        x0 = np.array([
18            -0.1,  # ω
19            0.1,   # α
20            0.95,  # β
21            -0.05  # γ (negative for leverage effect)
22        ])
23        
24        bounds = (
25            (-10, 0),      # ω
26            (0, 1),        # α
27            (0.8, 0.999),  # β
28            (-1, 0.1)      # γ
29        )
30        
31        result = minimize(
32            lambda params: -self._log_likelihood(returns, params),
33            x0,
34            method='L-BFGS-B',
35            bounds=bounds
36        )
37        
38        self.omega, self.alpha, self.beta, self.gamma = result.x
39        self.fitted = True
40        
41        return {
42            'omega': self.omega,
43            'alpha': self.alpha,
44            'beta': self.beta,
45            'gamma': self.gamma,
46            'log_likelihood': -result.fun
47        }
48    
49    @staticmethod
50    @jit(nopython=True)
51    def _compute_log_variance_path(returns, omega, alpha, beta, gamma):
52        n = len(returns)
53        log_var = np.zeros(n)
54        
55        # Initialize
56        log_var[0] = np.log(returns[:10].var())
57        
58        for t in range(1, n):
59            var_prev = np.exp(log_var[t-1])
60            z = returns[t-1] / np.sqrt(var_prev) if var_prev > 0 else 0
61            
62            log_var[t] = (omega + 
63                         beta * log_var[t-1] + 
64                         alpha * abs(z) + 
65                         gamma * z)
66        
67        return np.exp(log_var)
68    
69    def _log_likelihood(self, returns, params):
70        omega, alpha, beta, gamma = params
71        
72        variance = self._compute_log_variance_path(
73            returns, omega, alpha, beta, gamma
74        )
75        variance = np.maximum(variance, 1e-8)
76        
77        ll = -0.5 * (np.log(2 * np.pi) + np.log(variance) + 
78                     returns**2 / variance)
79        
80        return np.sum(ll)
81

1class RealizedVolatility:
2    """
3    Various realized volatility estimators.
4    """
5    
6    @staticmethod
7    def parkinson(high: np.ndarray, low: np.ndarray, 
8                  window: int = 20) -> np.ndarray:
9        """
10        Parkinson high-low volatility estimator.
11        More efficient than close-to-close (uses intraday range).
12        
13        σ² = (1/(4·ln(2))) · E[(ln(H/L))²]
14        """
15        hl_ratio = np.log(high / low)
16        factor = 1.0 / (4 * np.log(2))
17        
18        # Rolling variance
19        hl_var = pd.Series(hl_ratio).rolling(window).var().values
20        
21        return np.sqrt(factor * hl_var)
22    
23    @staticmethod
24    def garman_klass(open_price: np.ndarray, 
25                     high: np.ndarray,
26                     low: np.ndarray,
27                     close: np.ndarray,
28                     window: int = 20) -> np.ndarray:
29        """
30        Garman-Klass OHLC volatility estimator.
31        
32        More efficient than Parkinson, uses all OHLC data.
33        """
34        log_hl = np.log(high / low)
35        log_co = np.log(close / open_price)
36        
37        gk = 0.5 * log_hl**2 - (2*np.log(2) - 1) * log_co**2
38        
39        return pd.Series(gk).rolling(window).mean().pipe(np.sqrt).values
40    
41    @staticmethod
42    def rogers_satchell(open_price: np.ndarray,
43                       high: np.ndarray, 
44                       low: np.ndarray,
45                       close: np.ndarray,
46                       window: int = 20) -> np.ndarray:
47        """
48        Rogers-Satchell volatility estimator.
49        Handles drift (trending markets) better than Garman-Klass.
50        """
51        log_ho = np.log(high / open_price)
52        log_lo = np.log(low / open_price)
53        log_hc = np.log(high / close)
54        log_lc = np.log(low / close)
55        
56        rs = log_ho * log_hc + log_lo * log_lc
57        
58        return pd.Series(rs).rolling(window).mean().pipe(np.sqrt).values
59    
60    @staticmethod
61    def yang_zhang(open_price: np.ndarray,
62                  high: np.ndarray,
63                  low: np.ndarray, 
64                  close: np.ndarray,
65                  window: int = 20) -> np.ndarray:
66        """
67        Yang-Zhang volatility estimator.
68        Combines overnight, opening, and intraday volatility.
69        Most accurate for markets with gaps.
70        """
71        # Overnight volatility
72        prev_close = np.roll(close, 1)
73        log_oc = np.log(open_price[1:] / prev_close[1:])
74        overnight_var = np.var(log_oc)
75        
76        # Opening volatility  
77        log_co = np.log(close / open_price)
78        open_var = np.var(log_co)
79        
80        # Rogers-Satchell (intraday)
81        rs = RealizedVolatility.rogers_satchell(
82            open_price, high, low, close, window
83        )**2
84        
85        # Combine (simplified Yang-Zhang)
86        k = 0.34 / (1 + (window + 1) / (window - 1))
87        
88        yz_var = overnight_var + k * open_var + (1 - k) * rs
89        
90        return np.sqrt(yz_var)
91
92# Example: Compare estimators
93ohlc = pd.read_csv('spy_ohlc.csv')
94
95rv_parkinson = RealizedVolatility.parkinson(
96    ohlc['high'].values, ohlc['low'].values
97)
98
99rv_gk = RealizedVolatility.garman_klass(
100    ohlc['open'].values,
101    ohlc['high'].values, 
102    ohlc['low'].values,
103    ohlc['close'].values
104)
105
106rv_yz = RealizedVolatility.yang_zhang(
107    ohlc['open'].values,
108    ohlc['high'].values,
109    ohlc['low'].values, 
110    ohlc['close'].values
111)
112
113print(f"Parkinson Vol: {rv_parkinson[-1]*100:.2f}%")
114print(f"Garman-Klass Vol: {rv_gk[-1]*100:.2f}%")
115print(f"Yang-Zhang Vol: {rv_yz[-1]*100:.2f}%")
116

1#include <cmath>
2#include <algorithm>
3
4class BlackScholes {
5public:
6    // Standard normal CDF
7    static double norm_cdf(double x) {
8        return 0.5 * std::erfc(-x / std::sqrt(2.0));
9    }
10    
11    // Standard normal PDF
12    static double norm_pdf(double x) {
13        return std::exp(-0.5 * x * x) / std::sqrt(2.0 * M_PI);
14    }
15    
16    // Black-Scholes option price
17    static double call_price(double S, double K, double r, 
18                            double sigma, double T) {
19        if (T <= 0) return std::max(S - K, 0.0);
20        
21        double d1 = (std::log(S/K) + (r + 0.5*sigma*sigma)*T) 
22                    / (sigma * std::sqrt(T));
23        double d2 = d1 - sigma * std::sqrt(T);
24        
25        return S * norm_cdf(d1) - K * std::exp(-r*T) * norm_cdf(d2);
26    }
27    
28    static double put_price(double S, double K, double r,
29                           double sigma, double T) {
30        if (T <= 0) return std::max(K - S, 0.0);
31        
32        double d1 = (std::log(S/K) + (r + 0.5*sigma*sigma)*T)
33                    / (sigma * std::sqrt(T));
34        double d2 = d1 - sigma * std::sqrt(T);
35        
36        return K * std::exp(-r*T) * norm_cdf(-d2) - S * norm_cdf(-d1);
37    }
38    
39    // Vega (derivative w.r.t. volatility)
40    static double vega(double S, double K, double r,
41                      double sigma, double T) {
42        if (T <= 0) return 0.0;
43        
44        double d1 = (std::log(S/K) + (r + 0.5*sigma*sigma)*T)
45                    / (sigma * std::sqrt(T));
46        
47        return S * norm_pdf(d1) * std::sqrt(T);
48    }
49    
50    // Implied volatility using Newton-Raphson
51    static double implied_volatility(double price, double S, double K,
52                                    double r, double T, bool is_call,
53                                    double tol = 1e-6, int max_iter = 100) {
54        // Initial guess using Brenner-Subrahmanyam approximation
55        double moneyness = std::log(S / K);
56        double sigma = std::sqrt(2.0 * M_PI / T) * price / S;
57        
58        // Bound sigma
59        sigma = std::max(0.001, std::min(sigma, 5.0));
60        
61        for (int i = 0; i < max_iter; ++i) {
62            double theoretical_price = is_call 
63                ? call_price(S, K, r, sigma, T)
64                : put_price(S, K, r, sigma, T);
65            
66            double diff = theoretical_price - price;
67            
68            if (std::abs(diff) < tol) {
69                return sigma;
70            }
71            
72            double vega_val = vega(S, K, r, sigma, T);
73            
74            if (vega_val < 1e-10) {
75                break;  // Vega too small, can't converge
76            }
77            
78            // Newton-Raphson update
79            sigma = sigma - diff / vega_val;
80            
81            // Keep sigma in reasonable bounds
82            sigma = std::max(0.001, std::min(sigma, 5.0));
83        }
84        
85        return sigma;  // Return best estimate even if not converged
86    }
87};
88
89// Example usage
90int main() {
91    double S = 100.0;   // Spot price
92    double K = 105.0;   // Strike
93    double r = 0.05;    // Risk-free rate
94    double T = 0.25;    // Time to maturity (3 months)
95    
96    // Market price
97    double market_price = 2.50;
98    
99    // Calculate implied vol
100    double iv = BlackScholes::implied_volatility(
101        market_price, S, K, r, T, true
102    );
103    
104    std::cout << "Implied Volatility: " << iv * 100 << "%\n";
105    
106    return 0;
107}
108

1import numpy as np
2from scipy.optimize import minimize
3
4class SVIVolatilitySurface:
5    """
6    SVI (Stochastic Volatility Inspired) parameterization.
7    
8    Total variance: w(k) = a + b(ρ(k-m) + sqrt((k-m)² + σ²))
9    
10    where:
11        k = log(K/F) (log-moneyness)
12        a = ATM total variance
13        b = slope of vol smile
14        ρ = skew (-1 to 1)
15        m = horizontal shift
16        σ = smoothness
17    """
18    
19    def __init__(self):
20        self.params = {}
21        
22    def fit(self, log_moneyness: np.ndarray,
23            total_variance: np.ndarray,
24            T: float) -> dict:
25        """
26        Fit SVI parameters to observed implied variances.
27        
28        Args:
29            log_moneyness: log(K/F) for each strike
30            total_variance: σ²·T for each strike
31            T: Time to maturity
32            
33        Returns:
34            Fitted parameters
35        """
36        # Initial guess
37        atm_var = total_variance[len(total_variance)//2]
38        
39        x0 = np.array([
40            atm_var,      # a
41            0.1,          # b
42            -0.3,         # ρ
43            0.0,          # m
44            0.1           # σ
45        ])
46        
47        # Constraints for arbitrage-free surface
48        def constraints(params):
49            a, b, rho, m, sigma = params
50            
51            # No calendar spread arbitrage: a + b*σ*sqrt(1-ρ²) >= 0
52            c1 = a + b * sigma * np.sqrt(1 - rho**2)
53            
54            # No butterfly arbitrage (simplified)
55            c2 = b * (1 + abs(rho))
56            
57            return [c1, c2]
58        
59        cons = {'type': 'ineq', 'fun': constraints}
60        
61        bounds = (
62            (0, None),      # a >= 0
63            (0, 1),         # b >= 0
64            (-0.999, 0.999),  # |ρ| < 1
65            (-1, 1),        # m
66            (0.01, 1)       # σ > 0
67        )
68        
69        result = minimize(
70            lambda p: self._objective(p, log_moneyness, total_variance),
71            x0,
72            method='SLSQP',
73            bounds=bounds,
74            constraints=cons
75        )
76        
77        self.params[T] = {
78            'a': result.x[0],
79            'b': result.x[1],
80            'rho': result.x[2],
81            'm': result.x[3],
82            'sigma': result.x[4],
83            'rmse': np.sqrt(result.fun / len(log_moneyness))
84        }
85        
86        return self.params[T]
87    
88    def _objective(self, params, k, w_market):
89        """Sum of squared errors."""
90        w_model = self._total_variance(k, params)
91        return np.sum((w_model - w_market)**2)
92    
93    @staticmethod
94    def _total_variance(k, params):
95        """SVI total variance formula."""
96        a, b, rho, m, sigma = params
97        
98        sqrt_term = np.sqrt((k - m)**2 + sigma**2)
99        
100        return a + b * (rho * (k - m) + sqrt_term)
101    
102    def interpolate(self, log_moneyness: float, T: float) -> float:
103        """Get implied variance for given strike and maturity."""
104        if T not in self.params:
105            raise ValueError(f"No parameters fitted for T={T}")
106        
107        params = self.params[T]
108        w = self._total_variance(
109            log_moneyness,
110            [params['a'], params['b'], params['rho'], 
111             params['m'], params['sigma']]
112        )
113        
114        return w
115    
116    def implied_vol(self, log_moneyness: float, T: float) -> float:
117        """Get implied volatility (not variance)."""
118        w = self.interpolate(log_moneyness, T)
119        return np.sqrt(w / T)
120
121# Example usage
122if __name__ == "__main__":
123    # Market data
124    strikes = np.array([90, 95, 100, 105, 110])
125    spot = 100.0
126    T = 0.25
127    
128    # Observed implied vols
129    ivs = np.array([0.25, 0.22, 0.20, 0.21, 0.23])
130    
131    # Convert to log-moneyness and total variance
132    forward = spot  # Simplified, assume r=q=0
133    log_moneyness = np.log(strikes / forward)
134    total_variance = ivs**2 * T
135    
136    # Fit SVI
137    svi = SVIVolatilitySurface()
138    params = svi.fit(log_moneyness, total_variance, T)
139    
140    print(f"SVI Parameters:")
141    print(f"  a = {params['a']:.6f}")
142    print(f"  b = {params['b']:.6f}")
143    print(f"  ρ = {params['rho']:.6f}")
144    print(f"  m = {params['m']:.6f}")
145    print(f"  σ = {params['sigma']:.6f}")
146    print(f"  RMSE = {params['rmse']:.6f}")
147    
148    # Interpolate at new strike
149    new_strike = 97.5
150    new_k = np.log(new_strike / forward)
151    interp_vol = svi.implied_vol(new_k, T)
152    print(f"\nInterpolated IV at K={new_strike}: {interp_vol*100:.2f}%")
153

1Metric                   GARCH(1,1)   Historical   Implied Vol
2──────────────────────────────────────────────────────────────
31-day forecast RMSE      0.031        0.045        0.028
45-day forecast RMSE      0.042        0.058        0.039
520-day forecast RMSE     0.055        0.071        0.052
6
7Correlation with realized: 0.68      0.52         0.71
8

1Estimator           Efficiency   Best For
2────────────────────────────────────────────
3Close-to-close      1.00×        Baseline
4Parkinson (HL)      4.90×        Clean markets
5Garman-Klass (OHLC) 7.40×        Low noise
6Rogers-Satchell     6.20×        Trending
7Yang-Zhang          8.30×        Gaps/jumps
8