sharperatioandrsquared.py

import datetime
import random
import pandas as pd
import numpy as np

from syscore.interactive.progress_bar import progressBar
from systems.accounts.account_forecast import pandl_for_instrument_forecast
import statsmodels.formula.api as sm
import matplotlib

matplotlib.use("TkAgg")
matplotlib.rcParams.update({'font.size': 22})

#### FROM FORECAST ACCURACY TO SR

# generate a price series / series of returns at horizon


def avg_over_many_runs(forecast_horizon: int, noise_coef: float, length_days: int,
                       p: progressBar,
                       monte: int=100):
    many_runs=[]
    for __ in range(monte):
        p.iterate()
        many_runs.append(one_instant_of_regr_results_and_sr(
        forecast_horizon=forecast_horizon,
        noise_coef=noise_coef,
        length_days=length_days
    ))

    pd_all = pd.DataFrame(many_runs)

    return pd_all.mean()

def one_instant_of_regr_results_and_sr(forecast_horizon: int, noise_coef: float, length_days: int):
    price = generate_random_price_series(length_days)
    risk_adj_return = future_risk_adjusted_returns(price, forecast_horizon)
    forecast = calculate_forecast(risk_adj_return, noise_coef)

    try:
        r_squared, tstat = regression_results(risk_adj_return, forecast)
        SR = sharpe_ratio(price, forecast)
    except:
        return dict(SR=np.nan, r_squared=np.nan, tstat=np.nan)

    return dict(SR=SR, r_squared=r_squared, tstat=tstat)

def regression_results(risk_adj_change: pd.Series, forecast: pd.Series):
    both = pd.concat([risk_adj_change, forecast], axis=1)
    both.columns = ['price', 'forecast']
    both = both.dropna()
    result = sm.ols(formula="price ~ forecast", data=both).fit()

    return result.rsquared, result.tvalues.forecast

def sharpe_ratio( price: pd.Series, forecast: pd.Series):
    abs_median_forecast = forecast.abs().median()
    pandl = pandl_for_instrument_forecast(price=price, forecast=forecast, target_abs_forecast=abs_median_forecast)

    return pandl.sharpe()

def calculate_forecast( risk_adj_change: pd.Series, noise_coeff: float=10):
    noise = [random.gauss(0, noise_coeff) for __ in range(len(risk_adj_change))]

    return risk_adj_change+noise

def future_risk_adjusted_returns( price: pd.Series, forecast_horizon: int):
    stdev_estimate = price.diff().rolling(30, min_periods=3).std()
    future_price = price.shift(-forecast_horizon)
    px_change = future_price - price
    risk_adj_change = px_change / stdev_estimate

    return risk_adj_change

def generate_random_price_series(length_days: int):
    daily_returns = pd.Series(generate_random_returns(length_days), index=arbitrary_timeindex(length_days))
    price = daily_returns.cumsum()

    return price

def generate_random_returns(length_days: int):
    ## mean / stdev don't matter
    return [random.gauss(0,1) for __ in range(length_days)]

def arbitrary_timeindex(Nperiods, index_start=datetime.date(2000, 1, 1)):
    """
    For nice plotting, convert a list of prices or returns into an arbitrary pandas time series
    """

    ans = pd.bdate_range(start=index_start, periods=Nperiods)

    return ans


"""
all_horizons = [5,10,21, 43, 62, 125,250] ## 1 week, 2 weeks, 1 month, 2m, 3m, 6m, 1 year
all_noises = [1, 3, 5, 7, 10,20,30, 40,60,70, 80, 100, 150, 200]

length_days=2500 ## 10 years

monte=500
all_results = {}

p = progressBar(monte*len(all_horizons)*len(all_noises))
for forecast_horizon in all_horizons:
    all_results[forecast_horizon] = {}
    for noise_coef in all_noises:
        run_results = avg_over_many_runs(forecast_horizon=forecast_horizon,
                                         noise_coef=noise_coef,
                                         length_days=length_days,
                                         monte=monte,
                                         p=p)

        all_results[forecast_horizon][noise_coef] = run_results


def extract_data_to_plot(forecast_horizon: int, all_results: dict):
    all_this_item=[]
    for noise_coef in all_noises:
        SR = getattr(all_results[forecast_horizon][noise_coef], 'SR')
        r2 = getattr(all_results[forecast_horizon][noise_coef], 'r_squared')
        all_this_item.append(dict(SR=SR, r2=r2))

    return pd.DataFrame(all_this_item)


forecast_horizon=250
extract_data_to_plot(forecast_horizon, all_results).plot( 'r2','SR', logx=True)
matplotlib.pyplot.title("%d days" % forecast_horizon)

extract_data_to_plot(5, all_results).plot('SR', 'r2', logy=True)

import pickle
with open("/home/rob/temp/chapter16.pck", 'wb') as f:
    pickle.dump(file=f, obj=all_results)

"""

from systems.provided.basic.system import basic_db_futures_system
from systems.provided.rules.carry import carry
from systems.provided.rules.ewmac import ewmac
from systems.trading_rules import TradingRule, create_variations


ewmacrule = TradingRule(ewmac, ['rawdata.get_daily_prices', 'rawdata.daily_returns_volatility'],  dict(Lfast=4,Lslow=16))
carryrule = TradingRule(carry, ['rawdata.raw_carry'])
mom_rules = create_variations(ewmacrule, [dict(Lfast=2, Lslow=8), dict(Lfast=4,Lslow=16),
                              dict(Lfast=8, Lslow=32), dict(Lfast=16, Lslow=64),
                              dict(Lfast=32, Lslow=128),
                              dict(Lfast=64, Lslow=256)], key_argname="Lfast", nameformat="Mom_%s_%s")
mom_rules['carry'] = carryrule

system=basic_db_futures_system(trading_rules=mom_rules)

from syscore.pandas.strategy_functions import turnover
from private.projects.futures_book.generate_instrument_list import MASTER_INSTRUMENT_LIST
from systems.accounts.account_forecast import pandl_for_instrument_forecast


def get_results_for_instrument(instrument_code: str, rule_name:str, system):
    forecast = system.rules.get_raw_forecast(instrument_code, rule_name)
    forecast = 10*forecast/forecast.abs().median()
    price = system.rawdata.get_daily_prices(instrument_code)

    annual_turnover = turnover(forecast, 10)
    forecast_horizon = int((365/annual_turnover)/2.0)

    risk_adj_return = future_risk_adjusted_returns(price, forecast_horizon)

    r_squared, tstat = regression_results(risk_adj_return, forecast)
    SR = sharpe_ratio(price, forecast)

    return dict(horizon=forecast_horizon, r2=r_squared, SR=SR)

all_results = []
p=progressBar(len(MASTER_INSTRUMENT_LIST)*len(system.rules.trading_rules()))
for instrument in MASTER_INSTRUMENT_LIST:
    for rule in list(system.rules.trading_rules().keys()):
        p.iterate()
        all_results.append(get_results_for_instrument(instrument, rule, system))

all_results = pd.DataFrame(all_results)
all_results = all_results[all_results.SR>0]

#all_horizons = [5,10,21, 43, 62, 125,250] ## 1 week, 2 weeks, 1 month, 2m, 3m, 6m, 1 year
hedges = [1,8,15,30,50,100,175,300]

hindex=7 ## can't be zero
hmin = hedges[hindex-1]
hmax = hedges[hindex]-1 ##ensure no overlap
subset_of_results = all_results[all_results['horizon'].between(hmin,hmax, inclusive='both')]
subset_of_results.plot.scatter( 'r2','SR', logx=True)
matplotlib.pyplot.title("Horizons between %d and %d days" % (hmin, hmax))
matplotlib.pyplot.show(block=True)