2015-10-09 3 views
10

Я новичок в python, поэтому надеюсь, что мои два вопроса ясны и полны. Я разместил фактический код и набор тестовых данных в формате csv ниже.Создайте несколько столбцов в Pandas Dataframe из одной функции.

Я смог построить следующий код (в основном при помощи вкладчиков StackOverflow) для расчета подразумеваемой волатильности опционного контракта с использованием метода Ньютона-Рафсона. Процесс вычисляет Вега при определении подразумеваемой волатильности. Хотя я могу создать новый столбец DataFrame для Implied Volatility с использованием метода применения Pandas DataFrame, я не могу создать второй столбец для Vega. Есть ли способ создать два отдельных столбца DataFrame, когда функция возвращает IV & Vega вместе?

Я пробовал:

  • return iv, vega из функции
  • df[['myIV', 'Vega']] = df.apply(newtonRap, axis=1)
  • Got ValueError: Shape of passed values is (56, 2), indices imply (56, 13)

Также пробовал:

  • return iv, vega из функции
  • df['myIV'], df['Vega'] = df.apply(newtonRap, axis=1)
  • Got ValueError: Shape of passed values is (56, 2), indices imply (56, 13)

Кроме того, процесс вычисления медленно. Я импортировал numba и реализовал декоратор @jit (nogil = True), но я вижу только улучшение производительности на 25%. Набор тестовых данных - это тест производительности, который содержит почти 900 000 записей. Время работы составляет 2 часа и 9 минут без numba или с numba, но witout nogil = True. Время выполнения при использовании numba и @jit (nogil = True) составляет 1 час 32 минуты. Могу ли я сделать лучше?

from datetime import datetime 
from math import sqrt, pi, log, exp, isnan 
from scipy.stats import norm 
from numba import jit 


# dff = Daily Fed Funds (Posted rate is usually one day behind) 
dff = pd.read_csv('https://research.stlouisfed.org/fred2/data/DFF.csv', parse_dates=[0], index_col='DATE') 
rf = float('%.4f' % (dff['VALUE'][-1:][0]/100)) 
# rf = .0015      # Get Fed Funds Rate https://research.stlouisfed.org/fred2/data/DFF.csv 
tradingMinutesDay = 450    # 7.5 hours per day * 60 minutes per hour 
tradingMinutesAnnum = 113400  # trading minutes per day * 252 trading days per year 
cal = USFederalHolidayCalendar() # Load US Federal holiday calendar 


@jit(nogil=True)        # nogil=True arg improves performance by 25% 
def newtonRap(row): 
    """Estimate Implied Volatility (IV) using Newton-Raphson method 

    :param row (dataframe): Options contract params for function 
     TimeStamp (datetime): Close date 
     Expiry (datetime): Option contract expiration date 
     Strike (float): Option strike 
     OptType (object): 'C' for call; 'P' for put 
     RootPrice (float): Underlying close price 
     Bid (float): Option contact closing bid 
     Ask (float): Option contact closing ask 

    :return: 
     float: Estimated implied volatility 
    """ 
    if row['Bid'] == 0.0 or row['Ask'] == 0.0 or row['RootPrice'] == 0.0 or row['Strike'] == 0.0 or \ 
     row['TimeStamp'] == row['Expiry']: 
     iv, vega = 0.0, 0.0   # Set iv and vega to zero if option contract is invalid or expired 
    else: 
     # dte (Days to expiration) uses pandas bdate_range method to determine the number of business days to expiration 
     # minus USFederalHolidays minus constant of 1 for the TimeStamp date 
     dte = float(len(pd.bdate_range(row['TimeStamp'], row['Expiry'])) - 
        len(cal.holidays(row['TimeStamp'], row['Expiry']).to_pydatetime()) - 1) 
     mark = (row['Bid'] + row['Ask'])/2 
     cp = 1 if row['OptType'] == 'C' else -1 
     S = row['RootPrice'] 
     K = row['Strike'] 
     # T = the number of trading minutes to expiration divided by the number of trading minutes in year 
     T = (dte * tradingMinutesDay)/tradingMinutesAnnum 
     # TODO get dividend value 
     d = 0.00 
     iv = sqrt(2 * pi/T) * mark/S  # Closed form estimate of IV Brenner and Subrahmanyam (1988) 
     vega = 0.0 
     for i in range(1, 100): 
      d1 = (log(S/K) + T * (rf - d + iv ** 2/2))/(iv * sqrt(T)) 
      d2 = d1 - iv * sqrt(T) 
      vega = S * norm.pdf(d1) * sqrt(T) 
      model = cp * S * norm.cdf(cp * d1) - cp * K * exp(-rf * T) * norm.cdf(cp * d2) 
      iv -= (model - mark)/vega 
      if abs(model - mark) < 1.0e-9: 
       break 
     if isnan(iv) or isnan(vega): 
      iv, vega = 0.0, 0.0 
    # TODO Return vega with iv if add'l pandas column possible 
    # return iv, vega 
    return iv 


if __name__ == "__main__": 
    # test function from baseline data 
    get_csv = True 

    if get_csv: 
     csvHeaderList = ['TimeStamp', 'OpraSymbol', 'RootSymbol', 'Expiry', 'Strike', 'OptType', 'RootPrice', 'Last', 
         'Bid', 'Ask', 'Volume', 'OpenInt', 'IV'] 
     fileName = 'C:/tmp/test-20150930-56records.csv' 
     df = pd.read_csv(fileName, parse_dates=[0, 3], names=csvHeaderList) 
    else: 
     pass 

    start = datetime.now() 
    # TODO Create add'l pandas dataframe column, if possible, for vega 
    # df[['myIV', 'Vega']] = df.apply(newtonRap, axis=1) 
    # df['myIV'], df['Vega'] = df.apply(newtonRap, axis=1) 
    df['myIV'] = df.apply(newtonRap, axis=1) 
    end = datetime.now() 
    print end - start 

Данные испытаний: C /tmp/test-20150930-56records.csv

2015-09-30 16: 00: 00, AAPL151016C00109000, AAPL, 2015-10-16 16:00: 00,109, C, 109,95,3,46,3,6,3,7,1565,1290,0.3497 2015-09-30 16: 00: 00, AAPL151016P00109000, AAPL, 2015-10-16 16: 00: 00,109, P, 109,95,2,4, 2.34.2.42,3790,3087,0.3146 2015-09-30 16: 00: 00, AAPL151016C00110000, AAPL, 2015-10-16 16: 00: 00,110, C, 109,95,3,2,86,3,10217,28850, 0,3288 2015-09-30 16: 00: 00, AAPL151016P00110000, AAPL, 2015-10-16 16: 00: 00,110, P, 109,95,2,81,2,74,2,8,12113,44427,0.3029 2015-09-30 16 : 00: 00, AAPL151016C00111000, AAPL, 2015-10-16 16 : 00: 00,111, C, 109,95,2,35,2,44,2,45,6674,2318,0,3187 2015-09-30 16: 00: 00, AAPL151016P00111000, AAPL, 2015-10-16 16: 00: 00,111, P, 109,95 , 3,2,3,1,3,25,2031,3773,0,2926 2015-09-30 16: 00: 00, AAPL151120C00110000, AAPL, 2015-11-20 16: 00: 00,110, C, 109,95,5,9,5,7,5,95,5330 , 17112,0.3635 2015-09-30 16: 00: 00, AAPL151120P00110000, AAPL, 2015-11-20 16: 00: 00,110, P, 109,95,6,15,6,1,6,3,3724,15704,0,3842

ответ

13

Если я правильно понимаю, что вы должны делать, это вернуть Серию из вашей функции.Что-то вроде:

return pandas.Series({"IV": iv, "Vega": vega}) 

Если вы хотите поместить результат в новые столбцы одного и того же входного DataFrame, то просто сделать:

df[["IV", "Vega"]] = df.apply(newtonRap, axis=1) 
+0

Спасибо, BrenBarn. Это сработало отлично. – vlmercado

2

Что касается производительности с Numba обеспокоен, Numba Безразлично» t знать что-либо о pandas dataframes и не может скомпилировать операции с ними до быстрого машинного кода. Лучше всего определить, какая часть вашего метода медленная (например, line_profiler), а затем разгрузить эту часть другому методу, который вы создаете входы, используя атрибуты .values столбцов данных, которые дают вам доступ к базовому numpy массив. В противном случае numba просто будет работать в основном в «объектном режиме» (см. numba glossary) и не будет резко улучшать производительность.

+0

Хорошая точка. Я исследую line_profiler. Благодарю. – vlmercado

1

Трюк для векторизации кода - не думать в терминах строк, а вместо этого думать в терминах столбцов ,

Я почти иметь эту работу (я постараюсь закончить его позже), но вы хотите сделать что-то вдоль линий этого:

from datetime import datetime 
from math import sqrt, pi, log, exp, isnan 
from numpy import inf, nan 
from scipy.stats import norm 
import pandas as pd 
from pandas import Timestamp 
from pandas.tseries.holiday import USFederalHolidayCalendar 

# Initial parameters 
rf = .0015       # Get Fed Funds Rate https://research.stlouisfed.org/fred2/data/DFF.csv 
tradingMinutesDay = 450    # 7.5 hours per day * 60 minutes per hour 
tradingMinutesAnnum = 113400  # trading minutes per day * 252 trading days per year 
cal = USFederalHolidayCalendar() # Load US Federal holiday calendar 
two_pi = 2 * pi      # 2 * Pi (to reduce computations) 
threshold = 1.0e-9     # convergence threshold. 

# Create sample data: 
col_order = ['TimeStamp', 'OpraSymbol', 'RootSymbol', 'Expiry', 'Strike', 'OptType', 'RootPrice', 'Last', 'Bid', 'Ask', 'Volume', 'OpenInt', 'IV'] 
df = pd.DataFrame({'Ask': {0: 3.7000000000000002, 1: 2.4199999999999999, 2: 3.0, 3: 2.7999999999999998, 4: 2.4500000000000002, 5: 3.25, 6: 5.9500000000000002, 7: 6.2999999999999998}, 
        'Bid': {0: 3.6000000000000001, 1: 2.3399999999999999, 2: 2.8599999999999999, 3: 2.7400000000000002, 4: 2.4399999999999999, 5: 3.1000000000000001, 6: 5.7000000000000002, 7: 6.0999999999999996}, 
        'Expiry': {0: Timestamp('2015-10-16 16:00:00'), 1: Timestamp('2015-10-16 16:00:00'), 2: Timestamp('2015-10-16 16:00:00'), 3: Timestamp('2015-10-16 16:00:00'), 4: Timestamp('2015-10-16 16:00:00'), 5: Timestamp('2015-10-16 16:00:00'), 6: Timestamp('2015-11-20 16:00:00'), 7: Timestamp('2015-11-20 16:00:00')}, 
        'IV': {0: 0.3497, 1: 0.3146, 2: 0.3288, 3: 0.3029, 4: 0.3187, 5: 0.2926, 6: 0.3635, 7: 0.3842}, 
        'Last': {0: 3.46, 1: 2.34, 2: 3.0, 3: 2.81, 4: 2.35, 5: 3.20, 6: 5.90, 7: 6.15}, 
        'OpenInt': {0: 1290.0, 1: 3087.0, 2: 28850.0, 3: 44427.0, 4: 2318.0, 5: 3773.0, 6: 17112.0, 7: 15704.0}, 
        'OpraSymbol': {0: 'AAPL151016C00109000', 1: 'AAPL151016P00109000', 2: 'AAPL151016C00110000', 3: 'AAPL151016P00110000', 4: 'AAPL151016C00111000', 5: 'AAPL151016P00111000', 6: 'AAPL151120C00110000', 7: 'AAPL151120P00110000'}, 
        'OptType': {0: 'C', 1: 'P', 2: 'C', 3: 'P', 4: 'C', 5: 'P', 6: 'C', 7: 'P'}, 
        'RootPrice': {0: 109.95, 1: 109.95, 2: 109.95, 3: 109.95, 4: 109.95, 5: 109.95, 6: 109.95, 7: 109.95}, 
        'RootSymbol': {0: 'AAPL', 1: 'AAPL', 2: 'AAPL', 3: 'AAPL', 4: 'AAPL', 5: 'AAPL', 6: 'AAPL', 7: 'AAPL'}, 
        'Strike': {0: 109.0, 1: 109.0, 2: 110.0, 3: 110.0, 4: 111.0, 5: 111.0, 6: 110.0, 7: 110.0}, 
        'TimeStamp': {0: Timestamp('2015-09-30 16:00:00'), 1: Timestamp('2015-09-30 16:00:00'), 2: Timestamp('2015-09-30 16:00:00'), 3: Timestamp('2015-09-30 16:00:00'), 4: Timestamp('2015-09-30 16:00:00'), 5: Timestamp('2015-09-30 16:00:00'), 6: Timestamp('2015-09-30 16:00:00'), 7: Timestamp('2015-09-30 16:00:00')}, 
        'Volume': {0: 1565.0, 1: 3790.0, 2: 10217.0, 3: 12113.0, 4: 6674.0, 5: 2031.0, 6: 5330.0, 7: 3724.0}}) 
df = df[col_order] 

# Vectorize columns 
df['mark'] = (df.Bid + df.Ask)/2 
df['cp'] = df.OptType.map({'C': 1, 'P': -1}) 
df['Log_S_K'] = (df.RootPrice/df.Strike).apply(log) 
df['divs'] = 0 # TODO: Get dividend value. 
df['vega'] = 0. 
df['converged'] = False 

# Vectorized datetime calculations 
date_pairs = set(zip(df.TimeStamp, df.Expiry)) 
total_days = {(t1, t2): len(pd.bdate_range(t1, t2)) 
         for t1, t2 in date_pairs} 
hols = {(t1, t2): len(cal.holidays(t1, t2).to_pydatetime()) 
        for t1, t2 in date_pairs} 
del date_pairs 

df['total_days'] = [total_days.get((t1, t2)) 
        for t1, t2 in zip(df.TimeStamp, df.Expiry)] 
df['hols'] = [hols.get((t1, t2)) 
       for t1, t2 in zip(df.TimeStamp, df.Expiry)] 
df['days_to_exp'] = df.total_days - df.hols - 1 
df.loc[df.days_to_exp < 0, 'days_to_exp'] = 0 # Min zero. 
df.drop(['total_days', 'hols'], axis='columns', inplace=True) 
df['years_to_expiry'] = (df.days_to_exp * tradingMinutesDay/tradingMinutesAnnum) 

# Initial implied vol 'guess' 
df['implied_vol'] = (two_pi/df.years_to_expiry) ** 0.5 * df.mark/df.RootPrice 

for i in xrange(100): # range(100) in Python 3.x 
    # Create mask of options where the vol has not converged. 
    mask = [not c for c in df.converged.values] 
    if df.converged.all(): 
     break 

    # Aliases. 
    data = df.loc[mask, :] 
    cp = data.cp 
    mark = data.mark 
    S = data.RootPrice 
    K = data.Strike 
    d = data.divs 
    T = data.years_to_expiry 
    log_S_K = data.Log_S_K 
    iv = data.implied_vol 

    # Calcs. 
    d1 = (log_S_K + T * (rf - d + .5 * iv ** 2))/(iv * T ** 0.5) 
    d2 = d1 - iv * T ** 0.5 
    df.loc[mask, 'vega'] = vega = S * d1.apply(norm.pdf) * T ** 0.5 
    model = cp * (S * (cp * d1).apply(norm.cdf) 
        - K * (-rf * T).apply(exp) * (cp * d2).apply(norm.cdf)) 
    iv_delta = (model - mark)/vega 
    df.loc[mask, 'implied_vol'] = iv - iv_delta 

    # Clean-up and check for convergence. 
    df.loc[df.implied_vol < 0, 'implied_vol'] = 0 
    idx = model[(model - mark).abs() < threshold].index 
    df.ix[idx, 'converged'] = True 
    df.loc[:, 'implied_vol'].fillna(0, inplace=True) 
    df.loc[:, 'implied_vol'].replace([inf, -inf], nan, inplace=True) 
    df.loc[:, 'vega'].fillna(0, inplace=True) 
    df.loc[:, 'vega'].replace([inf, -inf], nan, inplace=True) 
+0

Александр, я вижу, где ты собираешься с этим. Я с нетерпением жду конечного результата. В строке 38 вашего кода он читает «df ['implied_vol'] = df.IV # Первоначальное« угадание ». IV в наборе данных, который я предоставил, является фактическим IV, основанным на количестве календарных дней до истечения срока действия и 365 дней за год. Фактическое первоначальное предположение IV в моем коде: «iv = sqrt (2 * pi/T) * mark/S # Закрытая форма оценки IV Brenner и Subrahmanyam (1988)». – vlmercado

+0

S - цена акций, которая в основном равна знаку, поэтому отметка/S составляет около 1,0, правильно? – Alexander

+0

Alexander, S - базовая цена акций. Знак - это середина между Bid и Ask of Option Contract. Используя AMZN в качестве примера, AMZN закрылся на отметке 539.80 в пятницу, 10 октября. Однако, ноябрь 20 540 Call закрылся с предложением 32.30. Задать вопрос на 32.60. Знак, следовательно, составляет 32.45. – vlmercado

Смежные вопросы