sim.py

# simple trading strategy simulator

import pandas as pd
from pandas.tools.plotting import autocorrelation_plot
from pandas.tools.plotting import scatter_matrix

import numpy as np
from scipy import stats

import sklearn
from sklearn import preprocessing as pp

import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import interactive
interactive(True)

import sys
import time

import logging as log
log.basicConfig(level=log.DEBUG)

import glob
import os.path
import pickle

import logging as log
log.basicConfig(level=log.DEBUG)

import random
import pdb

pd.set_option('display.width',500)

# define constant friction function
DefaultBPS = 10
def FrictionInBps(U, cfg, kvargs):
    """ default FrictionInBps function just returns default,
        but the interface receives all strategy info after 
        strategy is run, so one can create more realistic 
        impact models """
    return DefaultBPS

""" default simulator cfg dictionary.
    default keys/values:
    FrictionInBps - function that takes same args as strategy. 
                    by default, returns DefaultBps.
    InitBal - in $s
    Reinvest - should we reinvest our winnings or constantly assume we have InitBal?
    Verbose 
"""
DEF_SIM_CFG= { 'FrictionInBps': FrictionInBps,
               'Verbose'      : True,
               'InitBal'      : 1e7,
               'Reinvest'     : True  }

# columns in prepped univ
SIM_COLS = ["Sym","Product","Instrument",
            "Multiplier","Expiry","Strike",
            "Open","High","Low","Close","Volume"]

SIM_COLS_OUT = ["Prev_Weight", "Weight", "Prev_Qty", "Qty",
                "Trade_Qty",  "Trade_Fric", "PNL", "NET_PNL"]

SIM_COL_BALS =[ "NAV","Friction","PNL","NET_PNL", "Longs","Shorts",
                "Long_Dlrs","Short_Dlrs","Num_Trades","Turnover","NET_Return"]

  
def squarem( df, sym='Sym', min_pct=.9 ) :
# sim_squarem solves the common problem in which you have a large table of
#  data grouped by symbols, some of which have missing data.  You want to
#  'square' the data such that any symbol which is missing 'too much' data
#  is expunged and the remaining data is filled appropriately, leaving you
#  with a dataset which has the same # of observations for each symbol.
#
    bysyms = df.groupby(sym).size()
    idx = df.index.unique()
    onumsyms = len(bysyms)

    minlen = int(round(len(idx) * .9 ))
    keep = bysyms[bysyms > minlen]
    u = df[ df[sym].isin(keep.index) ]
    numsyms = len(keep)
    log.info('Got rid of %d/%d symbols',(numsyms-onumsyms),onumsyms)

    u.replace(0,np.nan,inplace=True)
    u.replace([np.inf, -np.inf], np.nan,inplace=True)

    u.sort_index(inplace=True)

    uidx = u.index.unique()
    # groupby and reindex magic
    z = u.groupby(sym).apply(
        lambda x:  x.reindex(uidx).ffill()).reset_index(0,drop=True)

    # badz = z[z.isnull().any(axis=1)]
    # if len(badz.index) > 0 :
    #     badtimes = badz.index.unique().values
    #     z.drop( badtimes, inplace=True )
    #     for dt in badtimes:
    #         log.info('removed %s for NaNs',pd.to_datetime(str(dt)).strftime(
    #             '%Y-%m-%d'))
    return z


def prep_univ( dateTime, symbol, 
               open, high, low, close, volume,
               product, instrument='STK', multiplier=1.0,expiry=None,
               strike=None,adv_days=20,sd_days=20, open2close_returns=True,
               scaleAndCenter=False,  **more_cols) :
    # constructs universe appropriate for use with simulator; any additional columns
    #  passed-in via ellipsis will be added to table as named
    #
    
    U = pd.DataFrame({'Sym': symbol, 
                      'Product' : product, 'Instrument':instrument,
                      'Multiplier': 1.0, 'Expiry': None, 'Strike':None,
                      'Open':open,'High':high, 'Low':low, 'Close':close,
                      'Volume':volume }, index=dateTime )
    U = U[ SIM_COLS ]
    if len(more_cols) > 0:
        U = pd.concat( [U, pd.DataFrame(more_cols)], axis=1 )

    U.reset_index( inplace=True)
    U.sort_values(['Sym','Date'],inplace=True)
    U.Date = pd.to_datetime(U.Date)
    U.set_index('Date',inplace=True)

    if scaleAndCenter :
        log.debug('prep_univ: scaling & centering')
        raw_scaled = U.groupby('Sym').transform(
            lambda x : (x - x.mean())/x.std())
        U = pd.concat([ u.Sym, raw_scaled], axis=1)
    
    # calculate adv, returns, fwd_returns & change in volume
    U['ADV'] = U.groupby('Sym')['Volume'].apply(
        pd.rolling_mean, adv_days, 1).shift()
    U['DeltaV'] = U.groupby('Sym')['Volume'].transform(
        lambda x : np.log(x / x.shift()) ) 
    U['Return'] = U.groupby('Sym')['Close'].transform(
        lambda x : np.log(x / x.shift()) ) 
    U['Fwd_Close'] = U.groupby('Sym')['Close'].shift(-1)
    U['Fwd_Return'] = U.groupby('Sym')['Close'].transform(
        lambda x : np.log(x / x.shift()).shift(-1) ) # fwd.returns
    U['SD'] = U.groupby('Sym')['Return'].apply(
        pd.rolling_std, sd_days, 1).shift()
    
    if open2close_returns:
        U['Fwd_Open'] = U.groupby('Sym')['Open'].shift(-1)
        U['Fwd_COReturn'] = np.divide(np.add( U.Fwd_Open, -U.Close ),U.Close)

    U.ffill(inplace=True)
    U.sort_index(inplace=True)
        
    return U


# simple, default strategy: equal weight universe on daily basis
def eq_wt( U, cfg, kvargs ) :
    #pdb.set_trace()
    U.Weight = 1/float(len(U.index))
    return U

# given today's Universe U and Yesterday's Y, set U's
#  Prev_Weight and Prev_Qty to Y's Weight & Qty
#  TODO: clean-up
def _getprevs( U, Y ) :
    #  TODO: surely there's a cleaner way to do this...
    wts = Y.reset_index()[['Sym','Weight']]
    wts.columns = ['Sym','Prev_Weight']
    pwts = U[['Sym']].merge( wts, on = 'Sym' )['Prev_Weight']
    U.Prev_Weight=pwts.values
    qts = Y.reset_index()[['Sym','Qty']]
    qts.columns = ['Sym','Prev_Qty']
    pqts = U[['Sym']].merge( qts, on = 'Sym' )['Prev_Qty']
    U.Prev_Qty=pqts.values

# functor to run strategy each day and update tbls ...
# TODO: clean-up
def __sim ( U, FUN, cfg, B, kvargs) :
    # run sim to set weights
    U = FUN( U, cfg, kvargs)

    # set prev values for weight & qty...
    Y = kvargs.pop('_Y', None)
    if Y is not None and not np.all(Y.index==U.index):
        _getprevs(U,Y)
        loop = 1 + int(kvargs.pop('_L'))
    else:
        loop = 0
    kvargs['_L'] = loop
    kvargs['_Y'] = U

    bb = B.iloc[loop]
    #  fill-out trade details 
    NAV = bb.NAV
    tospend = NAV/U.Weight
    U.Qty = np.round((NAV*U.Weight) / (U.Multiplier*U.Close))
    U.Trade_Qty  = U.Qty - U.Prev_Qty
    fbps = 1e-4 * cfg['FrictionInBps'](U,cfg,kvargs)
    U.Trade_Fric = U.Trade_Qty * U.Close * U.Multiplier * fbps 
    U.PNL = (U.Fwd_Close - U.Close) * U.Qty * U.Multiplier
    U.NET_PNL = U.PNL - U.Trade_Fric

    # today's balances are based on yesterday's posns...
    longs = U[U.Qty > 0]
    shorts = U[U.Qty < 0]
    trades = U[U.Trade_Qty != 0]
    
    bb.Friction = U.Trade_Fric.sum()
    bb.PNL = U.PNL.sum()
    bb.NET_PNL = U.NET_PNL.sum()
    bb.Longs = len(longs.index)
    bb.Shorts = len(shorts.index)
    bb.Long_Dlrs = (longs.Close * longs.Multiplier * longs.Qty).sum()
    bb.Short_Dlrs = (shorts.Close * shorts.Multiplier * shorts.Qty).sum()
    bb.Num_Trades = len(trades.index)
    bb.Turnover = (trades.Close * trades.Multiplier
                   * trades.Trade_Qty.abs()).sum()/NAV
    
    if loop > 0 :
        yb = B.iloc[loop-1]
        ynav = yb.NAV
        tnav = ynav + yb.NET_PNL
        bb.NAV = tnav
        bb.NET_Return = (tnav-ynav)/ynav
        
    B.iloc[loop] = bb

    # pdb.set_trace()
    return U


def sim( univ, sim_FUN=eq_wt, cfg=DEF_SIM_CFG.copy(), kvargs={} ) :
    """ simulator: runs simulation and returns a table of activity and balances.
        args:
          univ - historical data that's been produced by prep_univ
          sim_FUN - strategy function.  by default, equal weights univ.
          cfg - cfg info.  by default
          kvargs - strat-specific extra data in a dict

    """
#
    t0 = time.time()
    all_times = univ.index.unique().values

    # prepare writable/output side of universe
    W = pd.DataFrame( columns=SIM_COLS_OUT, index = univ.index).fillna(0.0)
    U = pd.concat( [univ, W], axis=1 )
    
    # create balances table: one per day
    B = pd.DataFrame( columns = SIM_COL_BALS, index = all_times ).fillna(0.0)
    B.NAV = cfg['InitBal']

    # 'daily' loop
    Z = U.groupby(U.index).apply( __sim, FUN=sim_FUN,
                                  cfg=cfg, B=B, kvargs=kvargs )
    
    log.info('ran over %d days and %d rows in %d secs', len(all_times),
             len(U.index),time.time()-t0)
        
    # summarize results a bit more...?
    #ts=xts(B$Net.Return,order.by=B$DateTime)
    
    # return universe and balances
    #list(U=U,B=B, ts=ts)
    return Z, B

def sharpe(Returns) :
    return np.sqrt(252) * np.mean(Returns)/np.std(Returns)

def random_strat( U, cfg, kvargs ) :
    # random portfolio strategy: picks 'num_names' randomly
    nnames = kvargs.get('num_names',10)
    names = random.sample(U.Sym, nnames )
    U.Weight = np.where( U.Sym.isin( names ), 1/float(nnames), 0 )
                 
    return U

def best_strat( U, cfg, kvargs ) :
    # portfolio strategy: picks 'num_names' based on trailing return
    nnames = kvargs.get('num_names',10)
    #pdb.set_trace()
    best = U.sort_values('Return',ascending=False,
                         na_position='last')['Sym'].head(10).values
    U.Weight = np.where( U.Sym.isin( best ), 1/float(nnames), 0 )                 
    return U

def worst_strat( U, cfg, kvargs ) :
    # portfolio strategy: picks 'num_names' based on trailing return
    nnames = kvargs.get('num_names',10)
    #pdb.set_trace()
    worst = U.sort_values('Return',ascending=True,
                          na_position='last')['Sym'].head(10).values
    U.Weight = np.where( U.Sym.isin( worst ), 1/float(nnames), 0 )                 
    return U


def rtest(U,FUN=random_strat, runs=10):
# run given strat repeatedly, plotting NAVs and Returning them
#   nb: this only makes sense if the strategy is random...
    # run random_strat 'runs' times and plot NAVs
    N = None
    for i in range(runs) :
        _,b = sim( U, sim_FUN=FUN )
        n = pd.DataFrame(b.NAV)
        N = n if N is None else pd.concat([N,n],axis=1)
    N.plot(legend=False)
    return N

def sim_test():
# dev driver
    f = 'U.pkl'
    P = pickle.load(open(f))
    log.info('loaded <%s>',f)
    P.describe()
    U = P[P.index >= '2005-01-01']
    U.describe()
    import sim
    _,B = sim.sim(U)
    #plot NAV
    B.NAV.plot(title='Equal Weight Everyone')
    return B