from scipy import *
from scipy.integrate import quad

#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)

#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)

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fmin

#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)

#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

#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

#calibration of parameters
kappa, theta, sigma, rho, v0 = calibrate(init_val, market_datas)

#
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)

#plot result
plt.plot(K, market_prices, 'g*',K, heston_prices, 'b')
plt.xlabel('Strike (K)')
plt.ylabel('Price')
plt.show()

#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

print call_price(kappa, theta, sigma, rho, v0, r, T, s0, 2205)

print call_price(kappa, theta, sigma, rho, v0, r, T, s0, 1.0)

print call_price(kappa, theta, sigma, rho, v0, r, T, s0, 1.5)   

print put_price(kappa, theta, sigma, rho, v0, r, T, s0, 1.0)

print put_price(kappa, theta, sigma, rho, v0, r, T, s0, 2205)

from statsmodels.nonparametric.kernel_regression import KernelReg
import numpy as np
import matplotlib.pyplot as plt

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

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)

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)

import pandas as pd
df = pd.read_csv("my_csv4.csv", parse_dates=['DateType'])
feed = PandasBarFeed(df, "SPX", barfeed.Frequency.DAY)

bars = feed.getNextBars()

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)

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, '+')

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

init_val = [1.1, 0.1, 0.4, -0.0, 0.1]

#calibration of parameters
kappa, theta, sigma, rho, v0 = calibrate(init_val, items)

#
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)

#plot result
plt.plot(K, market_prices, 'g*',K, heston_prices, 'b')
plt.xlabel('Strike (K)')
plt.ylabel('Price')
plt.show()