luphord / longstaff_schwartz Goto Github PK

from longstaff_schwartz.binomial import create_binomial_model
from longstaff_schwartz.stochastic_process import GeometricBrownianMotion
from longstaff_schwartz.algorithm import longstaff_schwartz
from numpy.polynomial import Polynomial
import numpy as np


def call_payoff(spot):
    return np.maximum(spot - strike, 0.0)

# Discount factor function
def constant_rate_df(t_from, t_to):
    return np.exp(-riskless_rate * (t_to - t_from))

# Approximation of continuation value
def fit_quadratic(x, y):
    return np.polynomial.Polynomial.fit(x, y, 2, rcond=None)

# Selection of paths to consider for exercise
# (and continuation value approxmation)
def itm(payoff, spot):
    return payoff > 0

# Model parameters
timegrid = np.linspace(0, 5, 100)  # timegrid for simulation
riskless_rate = 0.01  # riskless rate
sigma = 0.15  # annual volatility of underlying
simnum = 100  # number of simulated paths

# Simulate the underlying
gbm = GeometricBrownianMotion(mu=riskless_rate, sigma=sigma)
rnd = np.random.RandomState(1234)
underlying = gbm.simulate(timegrid, simnum, rnd)  # x.shape == (t.size, n)

# Payoff (exercise) function
strike = 0.95

def call_payoff_func(strike):
    return lambda S: np.maximum(S - strike, 0)

def american_call_exercise_barrier(mdl, strike):
    exercises = []
    payoff = call_payoff_func(strike)
    for cnt, s, ex, opt in mdl.evaluate_american_exercisable_iter(payoff):
        ex_idx = (ex >= cnt) & (ex > 0)
        ex_spots = s[ex_idx]
        exercises.append(ex_spots.max() if ex_idx.any() else np.nan)
    exercises.reverse()
    return np.array(exercises)

def american_call_exercise_barrier_fitted(mdl, strike, degree=4):
    barrier = american_call_exercise_barrier(mdl, strike)
    t = np.linspace(0, mdl.T, mdl.n)
    ex_exists = ~np.isnan(barrier)
    return Polynomial.fit(t[ex_exists], barrier[ex_exists], degree)

npv_american = longstaff_schwartz(underlying, timegrid, constant_rate_df,
                                  fit_quadratic, call_payoff, itm)
npv_european = constant_rate_df(timegrid[0], timegrid[-1]) * call_payoff(underlying[-1]).mean()

# Exercise Boundary
mdl = create_binomial_model(sigma=sigma, r=riskless_rate, S0=strike, T=timegrid[-1], n=len(timegrid))
exercise_barrier = american_call_exercise_barrier_fitted(mdl, strike, 3)