optimsation_dec_2023.py

"""
This piece of code gets fake data with the appropriate characteristics from my trading backtesting system:

import datetime
from systems.provided.basic.system import basic_db_futures_system
system = basic_db_futures_system()
sp500 = system.rawdata.get_daily_percentage_returns('SP500')
us10 = system.rawdata.get_daily_percentage_returns('US10')
us10 = us10.reindex(sp500.index)

import pandas as pd
both = pd.concat([sp500, us10], axis=1)
both.columns = ['equity', 'bonds']
both = both.dropna()
us10 = both.bonds
equity = both.equity

adj = 0.05 / (us10.std()*16)
us10 = adj*us10
adj=  (0.035/256) - us10.mean()
us10 = us10+adj

adj = 0.17 / (sp500.std()*16)
sp500 = adj * sp500
adj=  (0.0575/256) - sp500.mean()
sp500 = sp500+adj


both = pd.concat([sp500, us10], axis=1)
both.columns = ['equity', 'bonds']

both.to_csv('/home/rob/temp/optimisation_data.csv')


from syscore.interactive.display import set_pd_print_options
set_pd_print_options()

"""
import pickle

import matplotlib

matplotlib.use("TkAgg")
import pandas as pd

def pd_readcsv(
    filename: str,
    date_index_name: str = "index",
    date_format: str = "%Y-%m-%d",
) -> pd.DataFrame:

    df = pd.read_csv(filename)

    ## Add time index as index
    df = add_datetime_index(
        df=df, date_index_name=date_index_name, date_format=date_format
    )

    return df

def add_datetime_index(
    df: pd.DataFrame,
    date_index_name: str,
    date_format: str = "%Y-%m-%d",

) -> pd.DataFrame:
    date_index = df[date_index_name]
    date_index = date_index.astype(str)

    df.index = pd.to_datetime(date_index, format=date_format).values
    del df[date_index_name]
    df.index.name = None

    return df

data = pd_readcsv('/home/rob/temp/optimisation_data.csv')
import random

risk_free =.025
def single_bootstrap(data: pd.DataFrame) -> pd.DataFrame:
    len_data= len(data)
    list_of_rows = [data.iloc[random.randint(0,len_data-1)].values for __ in range(len_data)]
    list_of_rows = pd.DataFrame(list_of_rows)
    list_of_rows.index = data.index
    list_of_rows.columns = data.columns

    return list_of_rows

#single_bootstrap(data).cumsum().plot()

allocation = dict(equity=.6, bonds=.4)

def apply_portfolio_weights(data: pd.DataFrame, allocation: dict) -> pd.Series:
    weighted_returns = [data[instrument]*allocation[instrument] for instrument in allocation.keys()]
    weighted_returns = pd.concat(weighted_returns, axis=1)
    return weighted_returns.sum(axis=1)

def apply_leverage_and_portfolio_weights(data: pd.DataFrame, leverage: float, allocation: dict) -> pd.Series:
    portfolio_returns = apply_portfolio_weights(data, allocation=allocation)
    leveraged_returns = portfolio_returns * leverage
    annual_borrowing_cost = risk_free*(leverage-1)
    daily_borrowing_cost = annual_borrowing_cost / 256
    return leveraged_returns - daily_borrowing_cost

def single_bootstrap_applying_leverage_and_portfolio_weights(data: pd.DataFrame, leverage: float, allocation: dict):
    returns_for_bootstrap = single_bootstrap(data)

    return apply_leverage_and_portfolio_weights(data=returns_for_bootstrap, leverage=leverage, allocation=allocation)


def multiple_geometric_mean_bootstrap(data: pd.DataFrame, leverage: float, allocation: dict, number_of_runs: int = 50) -> pd.Series:

    values = [annual_geometric_return_given_account_curve(single_bootstrap_applying_leverage_and_portfolio_weights(
        data=data, leverage=leverage, allocation=allocation
    ))
    for __ in range(number_of_runs)]

    return pd.Series(values)

def annual_geometric_return_given_account_curve(account_curve: pd.Series) -> float:
    terminal_value = terminal_value_given_account_curve(account_curve)

    return annual_geometric_return_given_terminal_value(terminal_value, len(account_curve))

def terminal_value_given_account_curve(account_curve: pd.Series) -> float:
    one_plus = account_curve+1
    cumulative_returns = one_plus.cumprod()
    print(cumulative_returns.iloc[-1])
    return cumulative_returns.iloc[-1]

import math

def annual_geometric_return_given_terminal_value(terminal_value:float, number_of_points: int) -> float:
    return (1+daily_geometric_return_given_terminal_value(terminal_value, number_of_points))**256-1

import numpy as np
def daily_geometric_return_given_terminal_value(terminal_value:float, number_of_points: int) -> float:
    try:
        return math.exp((1/number_of_points)*math.log(terminal_value))-1
    except:
        return np.nan


x1=multiple_geometric_mean_bootstrap(data, 1.0, allocation, 1000)
x2=multiple_geometric_mean_bootstrap(data, 2.0, allocation, 1000)


leverage_options = [1, 1.25, 1.5, 1.75, 2.0, 2.173, 2.5, 3.5, 4.346, 5.0]
all_results = {}
for l in leverage_options:
    all_results[l] = multiple_geometric_mean_bootstrap(data, l, allocation, 250)

perc90 = [x.quantile(.9) for x in all_results.values()]
perc75 = [x.quantile(.75) for x in all_results.values()]
medians = [x.median() for x in all_results.values()]
perc25 = [x.quantile(.25) for x in all_results.values()]
perc10 = [x.quantile(.1) for x in all_results.values()]

all = pd.DataFrame({90: perc90, 75: perc75, 50: medians, 25: perc25, 10: perc10}, index =leverage_options)
all.plot()


from scipy.optimize import minimize
import numpy as np

def optimise_with_corr_and_std(mean_list, stdev_list, corrmatrix):
    sigma = sigma_from_corr_and_std(stdev_list, corrmatrix)

    weights = optimise_with_sigma(sigma, mean_list)
    return weights


def sigma_from_corr_and_std(stdev_list, corrmatrix):
    stdev = np.array(stdev_list, ndmin=2).transpose()
    sigma = stdev * corrmatrix * stdev

    return sigma

def optimise_with_sigma(sigma, mean_list):
    mus = np.array(mean_list, ndmin=2).transpose()
    number_assets = sigma.shape[1]
    start_weights = [1.0 / number_assets] * number_assets
    # Constraints - positive weights, adding to 1.0
    bounds = [(0.0, 1.0)] * number_assets
    cdict = [{'type': 'eq', 'fun': addem}]
    ans = minimize(
        neg_sharpe_ratio,
        start_weights, (sigma, mus),
        method='SLSQP',
        bounds=bounds,
        constraints=cdict,
        tol=0.00001)
    weights = ans['x']
    return weights

def neg_sharpe_ratio(weights, sigma, mus):
    # Returns minus the portfolio sharpe ratio (as we're minimising)
    weights = np.matrix(weights)
    estreturn = (weights * mus)[0, 0]
    std_dev = (variance(weights, sigma)**.5)

    return -(estreturn - risk_free) / std_dev

def portfolio_stdev(weights, sigma):
    std_dev = (variance(weights, sigma)**.5)

    return std_dev


def variance(weights, sigma):
    # returns the variance (NOT standard deviation) given weights and sigma
    return (weights * sigma * weights.transpose())[0, 0]


def addem(weights):
    # Used for constraints
    gap = 1.0 - sum(weights)
    return gap




optimise_with_corr_and_std([.035,.0575], [0.05,0.17], np.array([[1,0.0],[0.0,1]]))

optimise_with_corr_and_std([.035,.0575], [0.05,0.17], np.array([[1,0.65],[0.65,1]]))

def risk_weights_from_cas_weights(cash_weights: np.array, stdev: np.array) -> np.array:
    approx_weights = cash_weights * stdev
    sum_approx = np.sum(approx_weights)

    return approx_weights / sum_approx

risk_weights_from_cash_weights(np.array([0.77857023, 0.22142977]), np.array([0.05,0.17]))

#### BOOTSTRAPPING

def multiple_sharpe_ratio_bootstrap_allocation_to_equities(data: pd.DataFrame, leverage: float, allocation_to_equities: float, number_of_runs: int = 50) -> pd.Series:
    allocation={}
    allocation['equity'] = allocation_to_equities
    allocation['bonds'] = 1-allocation_to_equities

    values = [annual_sharpe_ratio_given_account_curve(single_bootstrap_applying_leverage_and_portfolio_weights(
        data=data, leverage=leverage, allocation=allocation
    ))
    for __ in range(number_of_runs)]

    return pd.Series(values)

def annual_sharpe_ratio_given_account_curve(account_curve: pd.Series)-> float:
    sr= 16*account_curve.mean()/account_curve.std()
    print(sr)
    return sr

allocation_options = [0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, .5, .6, .7, .8,.9, 1.0]
all_results = {}
for a in allocation_options:
    all_results[a] = multiple_sharpe_ratio_bootstrap_allocation_to_equities(data, leverage=1.0, allocation_to_equities=a, number_of_runs=250)

import pickle
with open('/home/rob/temp/all_results.ck', 'wb') as f:
    pickle.dump(all_results, f)

perc90 = [x.quantile(.9) for x in all_results.values()]
perc75 = [x.quantile(.75) for x in all_results.values()]
medians = [x.median() for x in all_results.values()]
perc25 = [x.quantile(.25) for x in all_results.values()]
perc10 = [x.quantile(.1) for x in all_results.values()]

all = pd.DataFrame({90: perc90, 75: perc75, 50: medians, 25: perc25, 10: perc10}, index =allocation_options)
all.plot()


def multiple_geometric_mean_bootstrap_allocation_to_equities(data: pd.DataFrame, leverage: float, allocation_to_equities: float, number_of_runs: int = 50) -> pd.Series:
    allocation={}
    allocation['equity'] = allocation_to_equities
    allocation['bonds'] = 1-allocation_to_equities

    values = [annual_geometric_return_given_account_curve(single_bootstrap_applying_leverage_and_portfolio_weights(
        data=data, leverage=leverage, allocation=allocation
    ))
    for __ in range(number_of_runs)]

    return pd.Series(values)


leverage_options = np.arange(start=1.0, stop=10.0, step=0.25)
allocation_options = [0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, .5, .6, .7, .8,.9, 1.0]

joint_results = {}
for a in allocation_options:
    joint_results[a] = {}
    for l in leverage_options:
        joint_results[a][l]= multiple_geometric_mean_bootstrap_allocation_to_equities(data, leverage=l, allocation_to_equities=a, number_of_runs=250)


with open('/home/rob/temp/join_results.ck', 'wb') as f:
    pickle.dump(joint_results, f)


leverage_options = np.arange(start=1.0, stop=10.0, step=0.25)
allocation_options = [0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, .5, .6, .7, .8,.9, 1.0]


with open('/home/rob/temp/join_results.ck', 'rb') as f:
    joint_results = pickle.load(f)


medians = [[x.median() for x in results_for_allocation.values()] for results_for_allocation in joint_results.values()]
df_medians = pd.DataFrame(medians)
df_medians.index = allocation_options
df_medians.columns = leverage_options

import matplotlib.pyplot as plt
def heatmap(df: pd.DataFrame):
    plt.imshow(df, cmap="RdYlBu")
    plt.colorbar()
    plt.xticks(range(len(df.columns)), df.columns)
    plt.yticks(range(len(df.index)), df.index)
    plt.show()

heatmap(df_medians)
df_medians[df_medians<0.03] = np.nan
heatmap(df_medians)


perc_25 = [[x.quantile(.25) for x in results_for_allocation.values()] for results_for_allocation in joint_results.values()]
perc_25 = pd.DataFrame(perc_25)
perc_25.index = allocation_options
perc_25.columns = leverage_options

perc_25[perc_25<0.01] = np.nan


def max_df_point(df):
    return df.stack().index[np.argmax(df.values)]

max_df_point(df_medians)
(0.3, 4.5)
max_df_point(perc_25)
(0.2, 3.25)