PCA_September_2025.py

from systems.accounts.account_forecast import pandl_for_instrument_forecast

import datetime
from copy import copy

import matplotlib
import scipy

matplotlib.rcParams.update({"font.size": 16})


matplotlib.use("TkAgg")

from typing import Tuple, List




from syscore.interactive.display import set_pd_print_options
set_pd_print_options()

import pandas as pd
from sklearn.decomposition import PCA as sklearnPCA
import statsmodels.formula.api as sm


from systems.provided.basic.system import basic_db_futures_system
from syscore.interactive.progress_bar import progressBar
from sysquant.fitting_dates import generate_fitting_dates, fitDates

## REPLACE THIS WITH YOUR OWN DATA IF YOU WISH
system = basic_db_futures_system()

MASTER_INSTRUMENT_LIST = ['AEX', 'ALUMINIUM', 'AUD', 'BBCOMM', 'BITCOIN', 'BOBL', 'BONO', 'BRE', 'BRENT-LAST', 'BTP', 'BTP3', 'BUND', 'BUXL', 'CAC', 'CAD', 'CHEESE', 'CHF', 'CNH', 'COPPER', 'CORN', 'CRUDE_W_mini', 'DAX', 'DJSTX-SMALL', 'DOW',  'ETHEREUM', 'EU-AUTO', 'EU-BASIC', 'EU-DIV30', 'EU-HEALTH', 'EU-INSURE', 'EU-OIL', 'EU-TECH', 'EU-TRAVEL', 'EU-UTILS', 'EUR', 'EURO600', 'EUROSTX', 'FEEDCOW', 'FTSECHINAA', 'FTSECHINAH', 'FTSETAIWAN', 'GAS-LAST', 'GASOILINE', 'GAS_US_mini', 'GBP', 'GBPEUR', 'GOLD_micro', 'HEATOIL', 'INR', 'IRON', 'JGB', 'JPY', 'KOSDAQ', 'KOSPI', 'KR10', 'KR3', 'LEANHOG', 'LIVECOW', 'MSCIASIA', 'MSCISING', 'MUMMY', 'MXP', 'NASDAQ_micro',  'NIKKEI', 'NIKKEI400', 'NOK', 'NZD', 'OAT', 'OATIES', 'PALLAD', 'PLAT', 'R1000', 'REDWHEAT', 'RICE', 'RUR', 'RUSSELL', 'SEK', 'SGD', 'SHATZ', 'SILVER', 'SMI', 'SOYBEAN', 'SOYMEAL', 'SOYOIL', 'SP400', 'SP500_micro', 'TOPIX', 'US10', 'US10U', 'US2', 'US20', 'US3', 'US30', 'US5',  'V2X', 'VIX', 'WHEAT', 'YENEUR']

prices = pd.DataFrame(dict([(instrument, system.data.daily_prices(instrument)) for instrument in MASTER_INSTRUMENT_LIST]))
returns = pd.DataFrame(dict([(instrument, system.rawdata.get_daily_percentage_returns(instrument)) for instrument in MASTER_INSTRUMENT_LIST]))
returns = returns.cumsum().ffill()
returns = returns.diff()

## setup
interval_frequency="1M"
rollyears = 12 ## 'years' is actually multiples of interval frequency
MAX_COMPONENTS = 1


def generate_all_regr_coeff():
    fit_dates = generate_fitting_dates(data=returns, date_method='rolling', rollyears=rollyears, interval_frequency=interval_frequency)
    p = progressBar(len(fit_dates))
    std_returns = returns.rolling(30).std().ffill()

    results = []
    for fit_date in fit_dates:
        res_this_period = calculate_regression_alpha_and_betas_and_pca_in_window(
            fit_date=fit_date,
            returns=returns,
            std_returns=std_returns
        )
        results.append(res_this_period)
        p.iterate()

    return results

def calculate_regression_alpha_and_betas_and_pca_in_window(returns: pd.DataFrame,std_returns: pd.DataFrame, fit_date: fitDates):
    norm_returns =calculate_norm_returns_in_sub_period(returns=returns, std_returns=std_returns, fit_date=fit_date, fitting=True)

    if len(norm_returns.columns)<2:
        return dict(fit_date=fit_date,regr_coef ={}, residuals_following_period={})

    instruments = norm_returns.columns
    sklearn_pca = sklearnPCA(n_components=max(1, MAX_COMPONENTS))
    try:
        sklearn_pca.fit_transform(norm_returns)
    except:
        print("SK learn failed with %d rows of data and %d columns use longer time period to fit or fewer components" % (len(norm_returns), len(norm_returns.columns)))
        return dict(fit_date=fit_date, regr_coef={}, residuals_following_period={})

    comp_count = len(sklearn_pca.components_)
    comp_count = min([MAX_COMPONENTS, comp_count])

    weights_by_component = [pd.DataFrame(sklearn_pca.components_[component,:], index=instruments) for
                            component in range(comp_count)]

    pca_returns = get_all_pca_returns(norm_returns=norm_returns, weights_by_component=weights_by_component)

    regr_coef = dict([(instrument, get_regression_coeff_for_instrument(instrument, pca_returns=pca_returns, norm_returns=norm_returns))
                 for instrument in instruments])

    norm_returns_following_period = calculate_norm_returns_in_sub_period(
        returns=returns, std_returns=std_returns, fit_date=fit_date,
        fitting=False
    )

    pca_returns_following_period = get_all_pca_returns(norm_returns=norm_returns_following_period,
                                                       weights_by_component=weights_by_component)

    residuals_following_period_as_dict =dict([
        (instrument_code, get_residuals(
        instrument_code=instrument_code,
        pca_returns_following_period=pca_returns_following_period,
        norm_returns_following_period=norm_returns_following_period,
        regr_coef=regr_coef
    )) for instrument_code in instruments])

    residuals_following_period_as_dict = dict([code, r] for code,r in residuals_following_period_as_dict.items() if len(r)>0)
    residuals_following_period = pd.DataFrame(residuals_following_period_as_dict)

    return dict(fit_date=fit_date,regr_coef =regr_coef, residuals_following_period=residuals_following_period)

def get_all_pca_returns( norm_returns: pd.DataFrame, weights_by_component: List[pd.DataFrame]):
    pca_returns = [get_pca_return_for_component(
        weight_for_component, norm_returns
    ) for weight_for_component in weights_by_component]
    pca_returns = pd.DataFrame(pca_returns).transpose()
    pca_names = ['P%d' % i for i in range(len(pca_returns.columns))]
    pca_returns.columns = pca_names

    return pca_returns

def get_pca_return_for_component(weight_for_component: pd.DataFrame, norm_returns: pd.DataFrame):
    weight_for_component_scaled = pd.concat([weight_for_component.transpose()]*len(norm_returns), axis=0)
    weight_for_component_scaled.index = norm_returns.index
    p_returns = weight_for_component_scaled *norm_returns

    return p_returns.sum(axis=1)

def get_regression_coeff_for_instrument(instrument_code, pca_returns: pd.DataFrame, norm_returns: pd.DataFrame):
    returns_for_instrument = norm_returns[instrument_code]
    both = pd.concat([returns_for_instrument, pca_returns], axis=1)
    pca_names = list(pca_returns.columns)
    both.columns = ['Y']+pca_names
    formula = "Y ~ "+"+ ".join(pca_names)

    result = sm.ols(formula=formula, data=both).fit()

    return result.params.to_dict()

def get_residuals(instrument_code, pca_returns_following_period: pd.DataFrame,
                  norm_returns_following_period: pd.DataFrame, regr_coef: dict):

    try:
        norm_returns_this_instrument = norm_returns_following_period[instrument_code]
    except KeyError:
        ## gone from data
        return pd.DataFrame()

    regr_coef_this_instrument = copy(regr_coef[instrument_code])
    regr_coef_intercept = regr_coef_this_instrument.pop("Intercept")
    beta_times_returns = [pca_returns_following_period[pca_name]*regr_coef_this_instrument[pca_name]
                          for pca_name in regr_coef_this_instrument.keys()]
    beta_times_returns_as_df = pd.DataFrame(beta_times_returns).transpose()
    explained_return_from_beta = beta_times_returns_as_df.sum(axis=1)
    explained_return = regr_coef_intercept + explained_return_from_beta

    both = pd.concat([norm_returns_this_instrument, - explained_return], axis=1)
    residual = both.sum(axis=1)

    return residual

def calculate_norm_returns_in_sub_period(returns: pd.DataFrame, std_returns: pd.DataFrame,
                                         fit_date: fitDates, fitting: bool = False):
    if fitting:
        returns_subperiod = returns[fit_date.fit_start:fit_date.fit_end]
        std_returns_subperiod = std_returns[fit_date.fit_start:fit_date.fit_end]
    else:
        returns_subperiod = returns[fit_date.period_start:fit_date.period_end]
        std_returns_subperiod = std_returns[fit_date.period_start:fit_date.period_end]

    norm_returns = returns_subperiod/std_returns_subperiod
    norm_returns = norm_returns.dropna(axis=1)

    return norm_returns



### alphas
def compare_alpha_vs_following_period_return( all_regr_coeff):
    std_returns = returns.rolling(30).std().ffill()
    comp_all = []
    for regr_coeff_et_al_in_period in all_regr_coeff:
        fit_date = regr_coeff_et_al_in_period['fit_date']
        regr_coef_this_period = regr_coeff_et_al_in_period['regr_coef']
        if len(regr_coef_this_period)==0:
            continue
        norm_returns_following_period = calculate_norm_returns_in_sub_period(
            returns=returns, std_returns=std_returns, fit_date=fit_date,
            fitting=False
        )
        alphas = dict([(instrument_code, regr_coef_this_period_and_instrument['Intercept']) for instrument_code, regr_coef_this_period_and_instrument in regr_coef_this_period.items()])
        avg_norm_returns = norm_returns_following_period.mean().to_dict()
        both = list(set(list(alphas.keys())).intersection(list(avg_norm_returns.keys())))
        compare = pd.DataFrame([(alphas.get(instrument_code), avg_norm_returns.get(instrument_code)) for instrument_code in both])
        compare.columns = ['alpha', 'ex_post_return']
        comp_all.append(compare)

    comp_all_as_single_df = pd.concat(comp_all, axis=0)

    return comp_all_as_single_df

all_regr_coeff = generate_all_regr_coeff()

comp_all_as_single_df = compare_alpha_vs_following_period_return(all_regr_coeff)

#comp_all_as_single_df.plot.scatter(x='alpha', y='ex_post_return')
matplotlib.pyplot.show(block=True)

formula = "ex_post_return ~ alpha"
result = sm.ols(formula=formula, data=comp_all_as_single_df).fit()
print(result.summary())

scipy.stats.ttest_ind(comp_all_as_single_df[comp_all_as_single_df.alpha>=0].ex_post_return,
                      comp_all_as_single_df[comp_all_as_single_df.alpha<0].ex_post_return)

print(comp_all_as_single_df[comp_all_as_single_df.alpha>=0].ex_post_return.mean())
print(comp_all_as_single_df[comp_all_as_single_df.alpha<0].ex_post_return.mean())


def get_alphas_as_df( all_regr_coeff):
    all_alphas = []
    for regr_coeff_et_al_in_period in all_regr_coeff:
        fit_date = regr_coeff_et_al_in_period['fit_date']
        regr_coef_this_period = regr_coeff_et_al_in_period['regr_coef']
        if len(regr_coef_this_period)==0:
            continue
        alphas = dict([(instrument_code, regr_coef_this_period_and_instrument['Intercept']) for instrument_code, regr_coef_this_period_and_instrument in regr_coef_this_period.items()])
        alphas_as_df_this_period = pd.DataFrame(alphas, index=[fit_date.period_start])
        all_alphas.append(alphas_as_df_this_period)

    all_alphas = pd.concat(all_alphas)

    return all_alphas

all_alphas = get_alphas_as_df(all_regr_coeff)
all_sr = {}
for instrument_code in system.get_instrument_list():
    p = progressBar(len(system.get_instrument_list()))
    try:
        pandl = pandl_for_instrument_forecast(
            forecast=all_alphas[instrument_code],
            price = prices[instrument_code]
        )
        p.iterate()
    except:
        continue

    sr = pandl.sharpe()
    all_sr[instrument_code]= sr

print("Average SR across instruments trading alpha")
print(pd.Series(all_sr).median())


def get_residuals_as_single_df( all_regr_coeff):
    all_residuals = []
    for regr_coeff_et_al_in_period in all_regr_coeff:
        residuals_this_period = regr_coeff_et_al_in_period['residuals_following_period']
        if len(residuals_this_period)==0:
            continue
        all_residuals.append(residuals_this_period)

    all_residuals_as_df = pd.concat(all_residuals, axis=0)
    all_residuals_as_df = all_residuals_as_df.resample("1B").last()

    return all_residuals_as_df

all_residuals_as_df = get_residuals_as_single_df(all_regr_coeff)


CUM_COUNT=22
print("Calculating rolling residuals...")
roll_residual = all_residuals_as_df.rolling(CUM_COUNT).sum()

all_sr = {}
p = progressBar(len(all_residuals_as_df.columns))
for instrument_code in all_residuals_as_df.columns:
    pandl = pandl_for_instrument_forecast(
        forecast=-roll_residual[instrument_code],
        price = prices[instrument_code]
    )
    sr = pandl.sharpe()
    all_sr[instrument_code]= sr
    p.iterate()

print("Median SR for residuals")
print(pd.Series(all_sr).mean())


##### Use alphas in trading model


from systems.provided.basic.system import *

class AlphaRawData(RawData):
    def get_alpha(self, code):
        alpha = all_alphas[code]
        daily_prices = self.get_daily_prices(code)

        return alpha.reindex(daily_prices.index)

def basic_futures_system(
    trading_rules
):
    config = Config()

    rules = Rules(trading_rules)
    data = csvFuturesSimData()
    system = System(
        [
            Account(),
            Portfolios(),
            PositionSizing(),
            AlphaRawData(),
            ForecastCombine(),
            ForecastScaleCap(),
            rules,
        ],
        data,
        config,
    )
    system.get_instrument_list(force_to_passed_list=MASTER_INSTRUMENT_LIST)
    system.config.notional_trading_capital = 500000000
    system.config.use_instrument_div_mult_estimates = True
    system.config.use_instrument_weight_estimates = False
    system.config.use_forecast_weight_estimates = False
    system.config.use_forecast_scale_estimates = True
    system.config.use_forecast_div_mult_estimates = False

    return system

def pass_thru_rawdata(item_from_raw_data):
    return item_from_raw_data

from systems.trading_rules import TradingRule


rule = TradingRule(pass_thru_rawdata, ['rawdata.get_alpha'])

system = basic_futures_system(dict(alpha=rule))
ac = system.accounts.portfolio()

from systems.provided.rules.ewmac import ewmac
rule2 = TradingRule(ewmac, ['rawdata.get_daily_prices', 'rawdata.daily_returns_volatility'], dict(Lfast=16, Lslow=64))

system2 = basic_futures_system(dict(mom=rule2))
ac2 = system2.accounts.portfolio()

both = pd.concat([ac, ac2], axis=1)
both.columns = ['alpha', 'ewmac']

system3 = basic_futures_system(dict(mom=rule2, alpha=rule))
system3.config.forecast_weights = dict(mom=.9, alpha=.1)
ac3 = system3.accounts.portfolio()

triple_thread = pd.concat([ac,ac2, ac3], axis=1)
triple_thread.columns = ['alpha', 'ewmac', 'both']