volatility_functions.py

# -*- coding: utf-8 -*-
"""Volatility Functions.ipynb

Automatically generated by Colab.

Original file is located at
    https://colab.research.google.com/drive/1V-ta8SQ3N-cp9q5GkkK3fDlLBJ2pE30v
"""

# Garch (1,1) Implementation

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from arch import arch_model
import yfinance as yf

# Download Apple stock price data from 2022 to 2023
apple_data = yf.download('AAPL', start='2022-01-01', end='2023-12-31')
prices = apple_data['Close']

# Calculate log returns and scale them
returns = np.log(prices / prices.shift(1)).dropna() * 100  # Scaling to avoid optimizer warning

def fit_garch_model(returns, p=1, q=1):
    """Fits a GARCH(p,q) model to the return series and returns the fitted model."""
    model = arch_model(returns, vol='Garch', p=p, q=q, dist='normal')
    fitted_model = model.fit(disp='off')
    return fitted_model

# Fit GARCH(1,1) model
garch_model = fit_garch_model(returns)

# Forecast volatility at next time step
forecasts = garch_model.forecast(horizon=1)
volatility_t = np.sqrt(forecasts.variance.iloc[-1, 0]) * np.sqrt(360)  # Apply annualization

print(f"Annualized Estimated Volatility at time t: {volatility_t:.4f}")

# Plot conditional volatility with annualization
plt.figure(figsize=(10, 5))
plt.plot(garch_model.conditional_volatility * np.sqrt(360), label='Annualized Conditional Volatility')
plt.xlabel('Time')
plt.ylabel('Volatility')
plt.title('Estimated Volatility using GARCH(1,1) for AAPL (2022-2023)')
plt.legend()
plt.show()



# Historical Volatility

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import yfinance as yf

# Download Apple stock price data from 2022 to 2023
apple_data = yf.download('AAPL', start='2022-01-01', end='2023-12-31')
prices = apple_data['Close']

# Calculate log returns
returns = np.log(prices / prices.shift(1)).dropna() * 100  # Scaling to avoid optimizer warning

# Function to calculate historical volatility
def historical_volatility(returns, window):
    """Computes historical volatility using a given window size."""
    rolling_std = returns.rolling(window=window).std()
    annualized_vol = rolling_std * np.sqrt(360)
    return annualized_vol  # No need to drop NaN values explicitly

# Compute historical volatility for different window sizes
window_sizes = [10, 30, 360, 21, 252]
hv_dict = {f"HV_{window}": historical_volatility(returns, window).values.flatten() for window in window_sizes}
hv_df = pd.DataFrame(hv_dict, index=returns.index)

# Plot historical volatility
plt.figure(figsize=(10, 5))
for window in window_sizes:
    plt.plot(hv_df.index, hv_df[f"HV_{window}"], label=f'Historical Vol {window} days')
plt.xlabel('Time')
plt.ylabel('Volatility')
plt.title('Historical Volatility Estimates')
plt.legend()
plt.show()

# American Option Binomial Asset Pricing Model

import numpy as np

def binomial_american_option(S, K, T, r, sigma, N, option_type="call"):
    dt = T / N  # Time step
    u = np.exp(sigma * np.sqrt(dt))  # Up factor
    d = 1 / u  # Down factor
    p = (np.exp(r * dt) - d) / (u - d)  # Risk-neutral probability

    # Price tree
    price_tree = np.zeros((N+1, N+1))
    for j in range(N+1):
        for i in range(j+1):
            price_tree[i, j] = S * (u ** (j - i)) * (d ** i)

    # Option value tree
    option_tree = np.zeros((N+1, N+1))

    # Compute terminal values
    if option_type == "call":
        option_tree[:, N] = np.maximum(price_tree[:, N] - K, 0)
    elif option_type == "put":
        option_tree[:, N] = np.maximum(K - price_tree[:, N], 0)

    # Backward induction
    for j in range(N-1, -1, -1):
        for i in range(j+1):
            expected_value = np.exp(-r * dt) * (p * option_tree[i, j+1] + (1 - p) * option_tree[i+1, j+1])

            if option_type == "call":
                option_tree[i, j] = max(price_tree[i, j] - K, expected_value)
            elif option_type == "put":
                option_tree[i, j] = max(K - price_tree[i, j], expected_value)

    return option_tree[0, 0]

# Example Usage
S = 100    # Initial stock price
K = 100    # Strike price
T = 1      # Time to maturity in years
r = 0.05   # Risk-free rate
sigma = 0.2 # Volatility
N = 100    # Number of steps

american_call_price = binomial_american_option(S, K, T, r, sigma, N, "call")
american_put_price = binomial_american_option(S, K, T, r, sigma, N, "put")

print(f"American Call Option Price: {american_call_price:.4f}")
print(f"American Put Option Price: {american_put_price:.4f}")

# American Options Pricing Using Least Squares Monte Carlo (Longstaff and Schwarz 2001)

import numpy as np
import scipy.stats as sp

# Simulating Geometric Brownian Motion paths
def generate_asset_paths(S0, r, sigma, T, M, I):
    dt = T / M
    S = np.zeros((M + 1, I))
    S[0] = S0
    for t in range(1, M + 1):
        z = np.random.standard_normal(I)
        S[t] = S[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * z)
    return S

# Least Squares Monte Carlo for American Option Pricing
def least_squares_mc(S0, K, r, sigma, T, M, I, option_type="put"):
    dt = T / M
    discount = np.exp(-r * dt)
    S = generate_asset_paths(S0, r, sigma, T, M, I)

    if option_type == "put":
        payoff = np.maximum(K - S, 0)
    else:
        payoff = np.maximum(S - K, 0)

    V = np.copy(payoff)

    for t in range(M - 1, 0, -1):
        itm = np.where(payoff[t] > 0)[0]
        X = S[t, itm]
        Y = V[t + 1, itm] * discount

        if len(X) > 0:
            A = np.vstack([np.ones_like(X), X, X ** 2]).T
            beta = np.linalg.lstsq(A, Y, rcond=None)[0]
            continuation_value = A @ beta
            exercise = payoff[t, itm] > continuation_value
            V[t, itm] = np.where(exercise, payoff[t, itm], V[t + 1, itm] * discount)

    option_price = np.mean(V[1] * discount)
    return option_price

# Example Usage
S0 = 100    # Initial stock price
K = 100     # Strike price
r = 0.05    # Risk-free rate
sigma = 0.2 # Volatility
T = 1       # Time to maturity (1 year)
M = 50      # Number of time steps
I = 10000   # Number of paths


price = least_squares_mc(S0, K, r, sigma, T, M, I, option_type="call")
print(f"American Call Option Price: {price:.4f}")
price = least_squares_mc(S0, K, r, sigma, T, M, I, option_type="put")
print(f"American Put Option Price: {price:.4f}")

# NOT NECESSARILY CORRECT

import numpy as np
import scipy.stats as sp
from numpy.polynomial.laguerre import lagval

# Binomial Tree Method for American Options
def binomial_american_option(S, K, T, r, sigma, N, option_type="call"):
    dt = T / N  # Time step
    u = np.exp(sigma * np.sqrt(dt))  # Up factor
    d = 1 / u  # Down factor
    p = (np.exp(r * dt) - d) / (u - d)  # Risk-neutral probability
    discount = np.exp(-r * dt)

    # Initialize price tree
    price_tree = np.zeros((N+1, N+1))
    for j in range(N+1):
        for i in range(j+1):
            price_tree[i, j] = S * (u**(j - i)) * (d**i)

    # Option value tree
    option_tree = np.zeros_like(price_tree)

    # Compute terminal values
    if option_type == "call":
        option_tree[:, -1] = np.maximum(price_tree[:, -1] - K, 0)
    else:
        option_tree[:, -1] = np.maximum(K - price_tree[:, -1], 0)

    # Backward induction
    for j in range(N-1, -1, -1):
        for i in range(j+1):
            continuation = discount * (p * option_tree[i, j+1] + (1 - p) * option_tree[i+1, j+1])
            if option_type == "call":
                option_tree[i, j] = max(price_tree[i, j] - K, continuation)
            else:
                option_tree[i, j] = max(K - price_tree[i, j], continuation)

    return option_tree[0, 0]

# Least Squares Monte Carlo for American Options
def generate_asset_paths(S0, r, sigma, T, M, I):
    dt = T / M
    S = np.zeros((M + 1, I))
    S[0] = S0
    for t in range(1, M + 1):
        z = np.random.standard_normal(I)
        S[t] = S[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * z)
    return S

def least_squares_mc(S0, K, r, sigma, T, M, I, option_type="put"):
    dt = T / M
    discount = np.exp(-r * dt)
    S = generate_asset_paths(S0, r, sigma, T, M, I)

    payoff = np.maximum(K - S, 0) if option_type == "put" else np.maximum(S - K, 0)
    V = np.copy(payoff)

    for t in range(M - 1, 0, -1):
        itm = np.where(payoff[t] > 0)[0]
        X = S[t, itm] / S0  # Normalize for stability
        Y = V[t + 1, itm] * np.exp(-r * dt * (M - t))  # Proper discounting

        if len(X) > 0:
            A = np.vstack([np.ones_like(X), X, 0.5 * (X**2 - 2*X + 1)]).T  # Improved basis
            beta = np.linalg.lstsq(A, Y, rcond=None)[0]
            continuation_value = A @ beta
            exercise = payoff[t, itm] > continuation_value
            V[t, itm] = np.where(exercise, payoff[t, itm], V[t + 1, itm] * discount)

    return np.mean(V[1] * discount)

# Parameters
S0 = 100    # Initial stock price
K = 100     # Strike price
r = 0.05    # Risk-free rate
sigma = 0.2 # Volatility
T = 1       # Time to maturity (1 year)
M = 100     # Increased time steps for LSM
I = 100000  # Increased number of paths for LSM
N = 100     # Steps in binomial tree

# Compute American Option Prices
binomial_call = binomial_american_option(S0, K, T, r, sigma, N, "call")
binomial_put = binomial_american_option(S0, K, T, r, sigma, N, "put")
lsm_call = least_squares_mc(S0, K, r, sigma, T, M, I, "call")
lsm_put = least_squares_mc(S0, K, r, sigma, T, M, I, "put")

# Output results
print(f"Binomial American Call Option Price: {binomial_call:.4f}")
print(f"Binomial American Put Option Price: {binomial_put:.4f}")
print(f"LSM American Call Option Price: {lsm_call:.4f}")
print(f"LSM American Put Option Price: {lsm_put:.4f}")