CoCalc Public FilesHeston_Model_Test.ipynbOpen with one click!
Authors: Erik Brandberg, Paul Cao
Views : 40

Heston Model Calibration

In [ ]:
from scipy import * from scipy.integrate import quad
In [ ]:
#public def call_price(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K): p1 = __p1(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K) p2 = __p2(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K) return (s0 * p1 - K * exp(-r * T) * p2) #calculated using put-call parity def put_price(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K): return K*exp(-r*T) - s0 + call_price(kappa, theta, sigma, rho, v0, r, T, s0, K)
In [ ]:
#private def __p(kappa, theta, sigma, rho, v0 ,r ,T ,s0 , K, status): integrand = lambda phi: (exp(-1j * phi * log(K)) * __f(phi, kappa, theta, sigma, rho, v0, r, T, s0, status) / (1j * phi)).real return (0.5 + (1 / pi) * quad(integrand, 0, 100)[0]) def __p1(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K): return __p(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K, 1) def __p2(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K): return __p(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K, 2) def __f(phi, kappa, theta, sigma, rho, v0, r, T, s0, status): if status == 1: u = 0.5 b = kappa - rho * sigma else: u = -0.5 b = kappa a = kappa * theta x = log(s0) d = sqrt((rho * sigma * phi * 1j - b)**2 - sigma**2 * (2 * u * phi * 1j - phi**2)) g = (b - rho * sigma * phi * 1j + d) / (b - rho * sigma * phi * 1j - d) C = r * phi * 1j * T + (a / sigma**2)*((b - rho * sigma * phi * 1j + d) * T - 2 * log((1 - g * exp(d * T))/(1 - g))) D = (b - rho * sigma * phi * 1j + d) / sigma**2 * ((1 - exp(d * T)) / (1 - g * exp(d * T))) return exp(C + D * v0 + 1j * phi * x)
In [ ]:
import numpy as np import matplotlib.pyplot as plt from scipy.optimize import fmin
In [ ]:
#sample market data def sample_data(): x = [x.split() for x in open('marketdata.txt')] header = x[0] market_datas = [] for market_data in x[1:]: market_datas += [map(lambda z:float(z), market_data)] return (header, market_datas)
In [ ]:
#parameter calibration(kappa, theta, sigma, rho, v0) def calibrate(init_val, market_datas): def error(x, market_datas): kappa, theta, sigma, rho, v0 = x print kappa, theta, sigma, rho, v0 result = 0.0 for market_data in market_datas: s0, k, market_price, r, T = market_data #print s0, k, market_price, r, T heston_price = put_price(kappa, theta, sigma, rho, v0, r, T, s0, k) result += (heston_price - market_price)**2 return result opt = fmin(error, init_val, args = (market_datas,), maxiter = 20) return opt
In [ ]:
#load market data header, market_datas = sample_data() #Initialize kappa, theta, sigma, rho, v0 init_val = [1.1, 0.1, 0.4, -0.0, 0.1] print market_datas
In [ ]:
#calibration of parameters kappa, theta, sigma, rho, v0 = calibrate(init_val, market_datas)
In [ ]:
# market_prices = np.array([]) heston_prices = np.array([]) K = np.array([]) for market_data in market_datas: s0, k, market_price, r, T = market_data heston_prices = np.append(heston_prices, call_price(kappa, theta, sigma, rho, v0, r, T, s0, k)) market_prices = np.append(market_prices, market_price) K = np.append(K,k)
In [ ]:
#plot result plt.plot(K, market_prices, 'g*',K, heston_prices, 'b') plt.xlabel('Strike (K)') plt.ylabel('Price') plt.show()

Test Put Pricing

In [ ]:
#maturity T = 0.0833 #risk free rate r = 0.05 #long term volatility(equiribrium level) theta = 0.1 #Mean reversion speed of volatility kappa = 1.1 #sigma(volatility of Volatility) sigma = 0.4 #rho rho = -0.6 #Initial stock price s0 = 2205 #Initial volatility v0 = 0.12
In [ ]:
print call_price(kappa, theta, sigma, rho, v0, r, T, s0, 2205)
In [ ]:
print call_price(kappa, theta, sigma, rho, v0, r, T, s0, 1.0)
In [ ]:
print call_price(kappa, theta, sigma, rho, v0, r, T, s0, 1.5)
In [21]:
print put_price(kappa, theta, sigma, rho, v0, r, T, s0, 1.0)
0.162755390807
In [ ]:
print put_price(kappa, theta, sigma, rho, v0, r, T, s0, 2205)
In [ ]:
from statsmodels.nonparametric.kernel_regression import KernelReg import numpy as np import matplotlib.pyplot as plt

Pandas Bar, Bars & Feed

In [ ]:
from pyalgotrade.bar import Bar from pyalgotrade import barfeed class PandasBar(Bar): def __init__(self, df, ix, bars): self.df = df self.ix = ix self.bars = bars def getRow(self): rowkey = self.df.index[self.ix] return self.df.ix[rowkey] def getDateTime(self): """Returns the :class:`datetime.datetime`.""" return self.bars.date def getOpen(self, adjusted=False): """Returns the opening price.""" ret = ( float(self.getRow()["Ask"]) + float(self.getRow()["Bid"]) ) / 2 if ret > 0: return ret else: return 0.01 #fix for worthless options def getHigh(self, adjusted=False): """Returns the highest price.""" return float(self.getRow()["Ask"]) def getLow(self, adjusted=False): """Returns the lowest price.""" return float(self.getRow()["Bid"]) def getClose(self, adjusted=False): """Returns the closing price.""" return self.getOpen() def getVolume(self): """Returns the volume.""" return float(self.getRow()["Volume"]) def getOpenInterest(self): """Returns the open interest.""" if "OpenInterest" in self.getRow(): return float(self.getRow()["OpenInterest"]) else: return float(0) def getAdjClose(self): """Returns the adjusted closing price.""" return self.getOpen() def getFrequency(self): """The bar's period.""" return barfeed.Frequency.DAY def getPrice(self): """Returns the closing or adjusted closing price.""" return self.getOpen() def setUseAdjustedValue(self, useAdjusted): pass def getUseAdjValue(self): return False
In [ ]:
class PandasBars(object): """A group of :class:`Bar` objects. :param barDict: A map of instrument to :class:`Bar` objects. :type barDict: map. .. note:: All bars must have the same datetime. """ def __init__(self, df, date, arr_ix=[]): self.date = date self.df = df self.arr_ix = arr_ix #key_val_arr = map(lambda ix: (df.ix[df.index[ix]]["Symbol"], PandasBar(df, ix, self)), self.arr_ix) #self.bar_dict = dict(key_val_arr) def set_bars_dict(self, bar_dict): self.bar_dict = bar_dict def keys(self): return self.bar_dict.keys() def __getitem__(self, instrument): """Returns the :class:`pyalgotrade.bar.Bar` for the given instrument. If the instrument is not found an exception is raised.""" return self.bar_dict[instrument] def __contains__(self, instrument): """Returns True if a :class:`pyalgotrade.bar.Bar` for the given instrument is available.""" return instrument in self.bar_dict def items(self): return self.bar_dict.items() def keys(self): return self.bar_dict.keys() def getInstruments(self): """Returns the instrument symbols.""" return self.bar_dict.keys() def getDateTime(self): """Returns the :class:`datetime.datetime` for this set of bars.""" return self.date def getBar(self, instrument): """Returns the :class:`pyalgotrade.bar.Bar` for the given instrument or None if the instrument is not found.""" return self.bar_dict.get(instrument, None)
In [ ]:
from more_itertools import peekable from pyalgotrade.bar import Bar from pyalgotrade import barfeed class PandasBarFeed(barfeed.BaseBarFeed): def __init__(self, dataframe, instrument, frequency): super(PandasBarFeed, self).__init__(frequency) self.registerInstrument(instrument) self.__df = dataframe self.__instrument = instrument self.__next = 0 self.started = False groups = self.__df.groupby(["DateType","Symbol"]).groups timestamp_dict = {} bars_dict = {} def loadBars(timestamp, symbol, index): if timestamp not in timestamp_dict: timestamp_dict[timestamp] = {} bars_dict[timestamp] = PandasBars(self.__df, timestamp) timestamp_dict[timestamp][symbol] = PandasBar(self.__df, index, bars_dict[timestamp]) map(lambda ((timestamp, symbol),index):loadBars(timestamp, symbol, index), groups.iteritems()) map(lambda (timestamp, bars):bars.set_bars_dict(timestamp_dict[timestamp]), bars_dict.iteritems()) dates = sorted(timestamp_dict.keys()) self.bars = map(lambda d:bars_dict[d], dates) self.bar_iter = peekable(self.bars) def reset(self): super(PandasBarFeed, self).reset() self.bar_iter = peekable(self.bars) self.started = False def peekDateTime(self): return self.getCurrentDateTime() def getCurrentDateTime(self): try: return self.bar_iter.peek().date except StopIteration: return None def barsHaveAdjClose(self): return True def getNextBars(self): try: return self.bar_iter.next() except StopIteration: return None def start(self): #super(PandasBarFeed, self).start() self.__started = True def stop(self): pass def join(self): pass def eof(self): try: if self.bar_iter.peek() and not self.bar_iter.peek() == self.bars[-1]: return False else: return True except StopIteration: return True def getCurrentBars(self): try: return self.bar_iter.peek() except StopIteration: return None def getNextValuesAndUpdateDS(self): dateTime, values = self.getNextValues() if dateTime is not None: for key, value in values.items(): pass return (dateTime, values)

Load a Single Volatility Simile

In [ ]:
import pandas as pd df = pd.read_csv("my_csv4.csv", parse_dates=['DateType']) feed = PandasBarFeed(df, "SPX", barfeed.Frequency.DAY) bars = feed.getNextBars()
In [ ]:
import datetime def get_chain(bars): def parse_option_symbol(symbol): try: opra_symbol = symbol[:-15] opra_expiry = datetime.datetime.strptime(symbol[-15:-9], '%y%m%d').date() opra_cp = symbol[-9] opra_price = int(symbol[-8:]) * .001 return {"Symbol":opra_symbol, "Expiration":opra_expiry, "Type":opra_cp, "Price":opra_price, "OptionSymbol":symbol} except Exception as e: print e chain = [] for symbol in bars.keys(): if symbol == "SPX" or symbol == "SPXW": continue chain.append(parse_option_symbol(symbol)) return chain def get_closet_expiration(n_days, chain, current_datetime): expiration_datetime = current_datetime + datetime.timedelta(days=n_days) option_chain = filter(lambda option:option['Expiration'] > expiration_datetime, chain) expiration_dates = map(lambda option:option['Expiration'], option_chain) return min(expiration_dates)
In [ ]:
import math from vollib.black_scholes.implied_volatility import implied_volatility chain = get_chain(bars) expiration_date = get_closet_expiration(45, chain, bars.getDateTime().date()) t = float(((expiration_date - bars.getDateTime().date()).total_seconds()) / (60 * 60 * 24 * 365)) def calc_iv(filtered_option_chain, bars, t): chains = [] for option in filtered_option_chain: price = float(bars[option['OptionSymbol']].getClose()) cp = option['Type'].lower() K = float(option['Price']) S = float((bars['SPX']).getClose()) r = 0.028 iv = implied_volatility(price, S, K, t, r, str(cp)) option["IV"] = iv chains.append(option) return chains spx_price = (bars['SPX']).getClose() option_chain = filter(lambda option:option['Expiration'] == expiration_date and option['Type'] == 'P' and math.fabs(option['Price']-spx_price) < 100, chain) #filter by 45-Day Put Simile #option_chain = calc_iv(option_chain, bars, t) put_iv = map(lambda x:(x['Price'], float(bars[x['OptionSymbol']].getClose())), option_chain) put_iv = sorted(put_iv, key=lambda x:x[0]) strikes = map(lambda x:x[0], put_iv) option_prices = map(lambda x:x[1], put_iv) plt.plot(strikes, option_prices, '+')

Load into Heston Model's Sample Data Format

In [ ]:
s0 = spx_price T = float(45.0/365.0) r = 0.028 i = 0 items = [] while i < len(strikes): K = strikes[i] option_price = option_prices[i] current_item = [s0, K, option_price, r, T] items.append(current_item) i = i + 1 print items
In [ ]:
init_val = [1.1, 0.1, 0.4, -0.0, 0.1] #calibration of parameters kappa, theta, sigma, rho, v0 = calibrate(init_val, items)
In [65]:
# market_prices = np.array([]) heston_prices = np.array([]) K = np.array([]) for market_data in items: s0, k, market_price, r, T = market_data heston_prices = np.append(heston_prices, put_price(kappa, theta, sigma, rho, v0, r, T, s0, k)) market_prices = np.append(market_prices, market_price) K = np.append(K,k)
In [66]:
#plot result plt.plot(K, market_prices, 'g*',K, heston_prices, 'b') plt.xlabel('Strike (K)') plt.ylabel('Price') plt.show()

FairPricingTool for IB. Given any option series; Drop any bids/asks where the sizes are significantly smaller; caliberate the option chain using Heston model; Compute the fair price of the option series.

FairPricingTool, BflyEntry, BflyLiquidate

In [ ]: