<jupyter_start><jupyter_text>There are now 20 targets! There are 10 different types of targets constructed and 2 versions of each constructed with a 20 day window or a 60 day window You are only scored on 'target' which currently corresponds to 'target_nomi_20' but could change in the future The other  auxillary targets can be very useful for training good models. Indeed you may find that a model trained on some of the auxillary targets generalize out-of-sample to 'target' better than a model trained on 'target'<jupyter_code>target = "target"
targets = [c for c in df if c.startswith("target")]
len(targets)
# The targets have a wide range of correlations with each other from ~0.3 to ~0.9
# This should allow the construction of many diverse models which ensemble nicely
df[targets].corr()<jupyter_output><empty_output><jupyter_text>1050 features, no feature groups!<jupyter_code># The features are no longer explicit feature groups, but you can see some pretty obvious patterns if you
# visualize the feature correlation matrix. Feel free to construct your own groupings!
plt.figure(figsize = (8,8))
plt.imshow(df[df.era=='0001'][features].corr())
# The models should be scored based on the rank-correlation (spearman) with the target
def numerai_score(y_true, y_pred):
    rank_pred = y_pred.groupby(eras).apply(lambda x: x.rank(pct=True, method="first"))
    return np.corrcoef(y_true, rank_pred)[0,1]

# It can also be convenient while working to evaluate based on the regular (pearson) correlation
def correlation_score(y_true, y_pred):
    return np.corrcoef(y_true, y_pred)[0,1]
# There are 574 eras numbered from 1 to 574
eras.describe()
# The earlier eras are smaller, but generally each era is 4000-5000 rows
df.groupby(eras).size().plot()
# The target is discrete and takes on 5 different values with about 5% in 0 & 1, 20% in 0.75 & 0.25, and 50% in 0.50
df.groupby(target).size()/len(df)
# Some new targets have a different distribution and the bins are less rigid in the percent of samples contained
print(df.groupby('target_william_20').size()/len(df))
print()
print(df.groupby('target_arthur_20').size()/len(df))<jupyter_output>target_william_20
0.000000    0.029735
0.166667    0.059873
0.333333    0.211150
0.500000    0.391798
0.666667    0.203550
0.833333    0.065247
1.000000    0.038511
dtype: float64

target_arthur_20
0.000000    0.028859
0.166667    0.073267
0.333333    0.258668
0.500000    0.384931
0.666667    0.173287
0.833333    0.052486
1.000000    0.028367
dtype: float64<jupyter_text>Some of the features are very correlatedSome even perfectly correlated within a particular era<jupyter_code>feature_corrs = df[df.era=='0001'][features].corr()
feature_corrs.stack().head()
tdf = feature_corrs.stack()
tdf = tdf[tdf.index.get_level_values(0) < tdf.index.get_level_values(1)]
tdf.sort_values()<jupyter_output><empty_output><jupyter_text>The correlation can change over timeYou can see this by comparing feature correlations on the first half and second half on the training set<jupyter_code>df1 = df[eras<=eras.median()]
df2 = df[eras>eras.median()]
corr1 = df1[df1.era=='0001'][features].corr().unstack()
corr1 = corr1[corr1.index.get_level_values(0) < corr1.index.get_level_values(1)]
corr2 = df2[df2.era=='0574'][features].corr().unstack()
corr2 = corr2[corr2.index.get_level_values(0) < corr2.index.get_level_values(1)]
tdf = pd.DataFrame({
    "corr1": corr1,
    "corr2": corr2,
})
tdf["corr_diff"] = tdf.corr2 - tdf.corr1
tdf.sort_values(by="corr_diff")<jupyter_output><empty_output><jupyter_text>Some features are predictive on their own<jupyter_code>feature_scores = {feature: score for feature, score in zip(features, np.corrcoef(df[df.era=='0001'][[target]+features].T)[1:,0])}
pd.Series(feature_scores).sort_values()
# Single features do not work consistently though
by_era_correlation = pd.Series({
    era: np.corrcoef(tdf[target], tdf["feature_grandmotherly_circumnavigable_homonymity"])[0,1]
    for era, tdf in df.groupby(eras)
})
by_era_correlation.plot()
# With a rolling 10 era average you can see some trends
by_era_correlation.rolling(10).mean().plot()<jupyter_output><empty_output><jupyter_text>The feature exposures of different targets can be very different, even orthogonal<jupyter_code>feature_scores_nomi20 = {feature: score for feature, score in zip(features, np.corrcoef(df[df.era=='0001'][['target_nomi_20']+features].T)[1:,0])}
feature_scores_nomi60 = {feature: score for feature, score in zip(features, np.corrcoef(df[df.era=='0001'][['target_nomi_60']+features].fillna(0.5).T)[1:,0])}
plt.scatter(feature_scores_nomi20.values(), feature_scores_nomi60.values())
print(np.corrcoef(list(feature_scores_nomi20.values()), list(feature_scores_nomi60.values()))[0,1])
feature_scores_jerome = {feature: score for feature, score in zip(features, np.corrcoef(df[df.era=='0001'][['target_jerome_20']+features].T)[1:,0])}
feature_scores_janet = {feature: score for feature, score in zip(features, np.corrcoef(df[df.era=='0001'][['target_janet_20']+features].T)[1:,0])}
plt.scatter(feature_scores_jerome.values(), feature_scores_janet.values())
print(np.corrcoef(list(feature_scores_jerome.values()), list(feature_scores_janet.values()))[0,1])<jupyter_output>0.05603328047427541<jupyter_text>Gotcha: MSE looks worse than correlation out of sampleModels will generally be overconfident, so even if they are good at ranking rows, the Mean-Squared-Error of the residuals could be larger than even the Mean-Squared-Error of the target (r-squared<0)<jupyter_code>df1 = df[eras<=eras.median()]
df2 = df[eras>eras.median()]
linear1 = linear_model.LinearRegression()
linear1.fit(df1[features], df1[target])
linear2 = linear_model.LinearRegression()
linear2.fit(df2[features], df2[target])
# Note in particular that the R-squared of (train_on_1, eval_on_2) is slightly negative!
r2 = [
    [
        model.score(dfX[features], dfX[target])
        for dfX in [df1, df2]
    ]
    for model in [linear1, linear2]
]
pd.DataFrame(r2, columns=["eval_on_1", "eval_on_2"], index=["train_on_1", "train_on_2"])
# Note in particular that the correlation of (train_on_1, eval_on_2) is quite decent
corrs = [
    [
        numerai_score(dfX[target], pd.Series(model.predict(dfX[features]), index=dfX.index))
        for dfX in [df1, df2]
    ]
    for model in [linear1, linear2]
]
pd.DataFrame(corrs, columns=["eval_on_1", "eval_on_2"], index=["train_on_1", "train_on_2"])
# This can be be run with LightGBM as well
lgb1 = lightgbm.LGBMRegressor()
lgb1.fit(df1[features], df1[target])
lgb2 = lightgbm.LGBMRegressor()
lgb2.fit(df2[features], df2[target])
r2 = [
    [
        model.score(dfX[features], dfX[target])
        for dfX in [df1, df2]
    ]
    for model in [lgb1, lgb2]
]
pd.DataFrame(r2, columns=["eval_on_1", "eval_on_2"], index=["train_on_1", "train_on_2"])
corrs = [
    [
        numerai_score(dfX[target], pd.Series(model.predict(dfX[features]), index=dfX.index))
        for dfX in [df1, df2]
    ]
    for model in [lgb1, lgb2]
]
pd.DataFrame(corrs, columns=["eval_on_1", "eval_on_2"], index=["train_on_1", "train_on_2"])<jupyter_output><empty_output><jupyter_text>Gotcha: Eras are now overlapping!Eras are weekly (5 days), but the target is four weeks out (20 days). You can subsample the eras (every 4th era) to prevent models from overfitting due to overlapping eras. You can then ensemble the models trained on each of the 4 different ways to subsample.<jupyter_code># train model on all eras in df1
lgb_all = lightgbm.LGBMRegressor()
lgb_all.fit(df1[features], df1[target])
lgb_all.feature_name_
# train models on subsamples eras in df1
lgb1 = lightgbm.LGBMRegressor()
lgb1.fit(df1[eras.isin(np.arange(1, 304, 4))][features], df1[eras.isin(np.arange(1, 304, 4))][target])

lgb2 = lightgbm.LGBMRegressor()
lgb2.fit(df1[eras.isin(np.arange(2, 304, 4))][features], df1[eras.isin(np.arange(2, 304, 4))][target])

lgb3 = lightgbm.LGBMRegressor()
lgb3.fit(df1[eras.isin(np.arange(3, 304, 4))][features], df1[eras.isin(np.arange(3, 304, 4))][target])

lgb4 = lightgbm.LGBMRegressor()
lgb4.fit(df1[eras.isin(np.arange(4, 304, 4))][features], df1[eras.isin(np.arange(4, 304, 4))][target])
# calculate predictions for each model
preds = [
        pd.Series(model.predict(df2[features]), index=df2.index)
    for model in [lgb_all, lgb1, lgb2, lgb3, lgb4]
]
# numerai score of model trained on all eras of df1
numerai_score(df2[target], preds[0])
# numerai score of ensemble of models trained on all subsampled eras of df1
numerai_score(df2[target], pd.concat(preds, axis=1)[[1,2,3,4]].mean(axis=1))<jupyter_output><empty_output><jupyter_text>Gotcha: eras are homogenous, but different from each other  Random cross-validation will look much better than cross-validating by eraEven for a simple linear model, taking a random shuffle reports a correlation of ~5.1%, but a time series split reports a lower score of ~2.6%  Gotcha: Eras are now overlapping!  Even era-wise cross-validation will now have leakage!Eras are weekly (5 days), but the target is four weeks out (20 days). You can either use blocks of eras and purge the overlapping sections or subsample the eras (every 4th era) to make them not overlapping and prevent leakage<jupyter_code># Subsample data to every 4th era to prevent overlapping eras
df = df[eras.isin(np.arange(1, 575, 4))]
df1 = df[eras<=eras.median()]
df2 = df[eras>eras.median()]
eras = eras[eras.isin(np.arange(1, 575, 4))]
# Because the TimeSeriesSplit class in sklearn does not use groups and won't respect era boundries, we implement
# a version that will

from sklearn.model_selection._split import _BaseKFold, indexable, _num_samples

class TimeSeriesSplitGroups(_BaseKFold):
    def __init__(self, n_splits=5):
        super().__init__(n_splits, shuffle=False, random_state=None)

    def split(self, X, y=None, groups=None):
        X, y, groups = indexable(X, y, groups)
        n_samples = _num_samples(X)
        n_splits = self.n_splits
        n_folds = n_splits + 1
        group_list = np.unique(groups)
        n_groups = len(group_list)
        if n_folds > n_groups:
            raise ValueError(
                ("Cannot have number of folds ={0} greater"
                 " than the number of samples: {1}.").format(n_folds,
                                                             n_groups))
        indices = np.arange(n_samples)
        test_size = (n_groups // n_folds)
        test_starts = range(test_size + n_groups % n_folds,
                            n_groups, test_size)
        test_starts = list(test_starts)[::-1]
        for test_start in test_starts:
            
            yield (indices[groups.isin(group_list[:test_start])],
                   indices[groups.isin(group_list[test_start:test_start + test_size])])

crossvalidators = [
    model_selection.KFold(5),
    model_selection.KFold(5, shuffle=True),
    model_selection.GroupKFold(5),
    TimeSeriesSplitGroups(5)
]
def correlation_score(y_true, y_pred):
    return np.corrcoef(y_true, y_pred)[0,1]
for cv in crossvalidators:
    print(cv)
    print(np.mean(
            model_selection.cross_val_score(
            linear_model.LinearRegression(),
            df[features],
            df[target],
            cv=cv,
            n_jobs=1,
            groups=eras,
            scoring=metrics.make_scorer(correlation_score, greater_is_better=True)
        )))
    print()<jupyter_output>KFold(n_splits=5, random_state=None, shuffle=False)
0.04081357263615691

KFold(n_splits=5, random_state=None, shuffle=True)
0.05205511367759082

GroupKFold(n_splits=5)
0.04341132164952086

TimeSeriesSplitGroups(n_splits=5)
0.02598605396562758<jupyter_text>Gotcha:  {0, 1} are noticeably different from {0.25, 0.75}This makes training a classifier one-versus-rest behave counterintuitively.Specifically, a multinomial classifier seem to learn how to pick out extreme targets, and their predictions are the most correlated<jupyter_code># Train a standard logistic regression as a classifier
logistic = linear_model.LogisticRegression(multi_class='multinomial', solver='saga')
logistic.fit(df1[features], (df1[target]*4).astype(int))
logistic.score(df1[features], (df1[target]*4).astype(int))
# The first and last class are highly correlated
corrs=np.corrcoef(logistic.predict_proba(df2[features]).T)
plt.imshow(corrs, vmin=-1, vmax=1, cmap="RdYlGn")
corrs
# Out-of-sample correlation is 3.5%
preds = pd.Series(logistic.predict_proba(df2[features]).dot(logistic.classes_), index=df2.index)
numerai_score(df2[target], preds)
# A standard linear model has a slightly lower correlation out-of-sample
linear = linear_model.LinearRegression()
linear.fit(df1[features], df1[target])
linear.score(df2[features], df2[target])
preds = pd.Series(linear.predict(df2[features]), index=df2.index)
numerai_score(df2[target], preds)<jupyter_output><empty_output><jupyter_text>Eras can be more or less applicable to other erasYou can test this be splitting the eras into blocks of 15, training on each block, and evaluating on each other block.<jupyter_code># convert the subsampled era labels to continuous integers
eras_sub = (eras-1)//4
eras15 = (eras_sub // 15) * 15
eras15.value_counts()
results15 = []
for train_era, tdf in df.groupby(eras15):
    print(train_era)
    model = linear_model.LinearRegression()
    model.fit(tdf[features], tdf[target])
    for test_era, tdf in df.groupby(eras15):
        results15.append([
            train_era,
            test_era,
            correlation_score(tdf[target], model.predict(tdf[features]))
        ])
results_df = pd.DataFrame(
    results15,
    columns=["train_era", "test_era", "score"]
).pivot(index="train_era", columns="test_era", values="score")
results_df
# Each row here is the training block of eras, each column is a testing block of eras.
# Note that there is a period in the early groups that does not seem to be relevant to other eras, and the
# overall performance seems to decrease a bit over time.
plt.figure(figsize=(15,15))
plt.imshow(results_df, vmin=-0.04, vmax=0.04, cmap="RdYlGn")<jupyter_output><empty_output><jupyter_text>Here is an advanced paper that talks about generalization.Eras can be thought about in the same way that "distributions" or "environments" are talked about herehttps://arxiv.org/pdf/1907.02893.pdf  Gotcha: Since the signal-to-noise ratio is so low, models can take many more iterations than expected, and have scarily high in-sample performance<jupyter_code>def our_score(preds, dtrain):
    return "score", -np.corrcoef(preds, dtrain.get_label())[0,1], False

df1 = df[eras<=eras.median()]
df2 = df[eras>eras.median()]
dtrain = lightgbm.Dataset(df1[features], df1[target])
dtest = lightgbm.Dataset(df2[features], df2[target])
dall = lightgbm.Dataset(df[features], df[target])
param = {
    'max_depth':3,
    'eta':0.1,
    'objective':'regression',
    'metric':'mse',
    'nthread': -1,
}
evals_result = {}
bst = lightgbm.train(
    params=param,
    train_set=dtrain,
    feval=our_score,
    num_boost_round=500,
    valid_sets=[dtrain, dtest],
    valid_names=['train', 'test'],
    evals_result=evals_result,
    verbose_eval=20,
)

(-pd.DataFrame({k: v['score'] for k,v in evals_result.items()})).plot(ylim=[0,0.25])<jupyter_output>[LightGBM] [Warning] Auto-choosing col-wise multi-threading, the overhead of testing was 0.395800 seconds.
You can set `force_col_wise=true` to remove the overhead.
[LightGBM] [Info] Total Bins 5250
[LightGBM] [Info] Number of data points in the train set: 305247, number of used features: 1050
[LightGBM] [Info] Start training from score 0.499998
[20]	train's l2: 0.0498386	train's score: -0.0696572	test's l2: 0.0499211	test's score: -0.042272
[40]	train's l2: 0.0497466	train's score: -0.0845498	test's l2: 0.0499021	test's score: -0.0449487
[60]	train's l2: 0.0496736	train's score: -0.0950992	test's l2: 0.049894	test's score: -0.0462429
[80]	train's l2: 0.0496094	train's score: -0.103748	test's l2: 0.0498883	test's score: -0.0474028
[100]	train's l2: 0.0495503	train's score: -0.111415	test's l2: 0.049883	test's score: -0.0485958
[120]	train's l2: 0.0494963	train's score: -0.118318	test's l2: 0.0498832	test's score: -0.0487773
[140]	train's l2: 0.0494427	train's score: -0.124419	test's l2[...]<jupyter_text>The results are sensitive to the choice of parameters, which should be picked through cross-validation<jupyter_code>df1 = df[eras<=eras.median()]
df2 = df[eras>eras.median()]
models = [
    linear_model.LinearRegression(),
] + [
    linear_model.ElasticNet(alpha=alpha)
    for alpha in [0.01, 0.005, 0.002, 0.001, 0.0005, 0.0002, 0.0001, 0.00005, 0.00002, 0.00001]
] + [
    lightgbm.LGBMRegressor(n_jobs=-1),
    lightgbm.LGBMRegressor(n_jobs=-1, learning_rate=0.01, n_estimators=1000),
    lightgbm.LGBMRegressor(n_jobs=-1, colsample_bytree=0.1, learning_rate=0.01, n_estimators=1000),
    lightgbm.LGBMRegressor(n_jobs=-1, colsample_bytree=0.1, learning_rate=0.01, n_estimators=1000, max_depth=5),
    lightgbm.LGBMRegressor(n_jobs=-1, colsample_bytree=0.1, learning_rate=0.001, n_estimators=10000, max_depth=5),
]
for model in models:
    print(" -- ", model)   
    model.fit(df1[features], df1[target])
    outsample = numerai_score(df2[target], pd.Series(model.predict(df2[features]), index=df2.index))
    insample = numerai_score(df1[target], pd.Series(model.predict(df1[features]), index=df1.index))
    print(
        f"outsample: {outsample}, insample: {insample}"
    )
    print()<jupyter_output>--  LinearRegression()
outsample: 0.03449532642406725, insample: 0.08811577962731876

 --  ElasticNet(alpha=0.01)
outsample: -0.0005817296572705758, insample: 0.0011047675627246087

 --  ElasticNet(alpha=0.005)
outsample: -0.0005817296572705758, insample: 0.0011047675627246087

 --  ElasticNet(alpha=0.002)
outsample: 0.028356554996391864, insample: 0.036860395887039446

 --  ElasticNet(alpha=0.001)
outsample: 0.03991189570044001, insample: 0.05365003421490132

 --  ElasticNet(alpha=0.0005)
outsample: 0.04275518297522411, insample: 0.06223910243312713

 --  ElasticNet(alpha=0.0002)
outsample: 0.04408886141392253, insample: 0.07101617648855875

 --  ElasticNet(alpha=0.0001)
outsample: 0.042497050882822206, insample: 0.07683662591705001

 --  ElasticNet(alpha=5e-05)
outsample: 0.04031370910547076, insample: 0.0826422949678376

 --  ElasticNet(alpha=2e-05)<jupyter_text>Gotcha: Models with large exposures to individual features tend to perform poorly or inconsistently out of sample<jupyter_code>import numpy as np
import scipy
# Train a standard xgboost on half the train eras
lgb = lightgbm.LGBMRegressor(n_estimators=1000, max_depth=5, learning_rate=0.01, n_jobs=-1)
lgb.fit(df1[features], df1[target])
lgb_preds = lgb.predict(df2[features])<jupyter_output><jupyter_text>Our predictions have correlation  0.20 for some single features!Sure hope those features continue to act as they have in the past!<jupyter_code>corr_list = []
for feature in features:
    corr_list.append(np.corrcoef(df2[feature], lgb_preds)[0,1])
corr_series = pd.Series(corr_list, index=features)
corr_series.describe()
from sklearn.preprocessing import MinMaxScaler
import scipy

def neutralize(df,
               columns,
               neutralizers=None,
               proportion=1.0,
               normalize=True,
               era_col="era"):
    if neutralizers is None:
        neutralizers = []
    unique_eras = df[era_col].unique()
    computed = []
    for u in unique_eras:
        print(u, end="\r")
        df_era = df[df[era_col] == u]
        scores = df_era[columns].values
        if normalize:
            scores2 = []
            for x in scores.T:
                x = (scipy.stats.rankdata(x, method='ordinal') - .5) / len(x)
                x = scipy.stats.norm.ppf(x)
                scores2.append(x)
            scores = np.array(scores2).T
        exposures = df_era[neutralizers].values

        scores -= proportion * exposures.dot(
            np.linalg.pinv(exposures.astype(np.float32)).dot(scores.astype(np.float32)))

        scores /= scores.std(ddof=0)

        computed.append(scores)

    return pd.DataFrame(np.concatenate(computed),
                        columns=columns,
                        index=df.index)
df2["preds"] = lgb_preds
df2["preds_neutralized"] = neutralize(df2, 
                                      columns=["preds"], 
                                      neutralizers=features, 
                                      proportion=0.5, # neutralize by 50% within each era
                                      normalize=True,
                                      era_col='era') 

scaler = MinMaxScaler()
df2["preds_neutralized"] = scaler.fit_transform(df2[["preds_neutralized"]]) # transform back to 0-1<jupyter_output><ipython-input-69-4aa87b479082>:1: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  df2["preds"] = lgb_preds<jupyter_text>Now our single feature exposures are much smaller<jupyter_code>corr_list2 = []
for feature in features:
    corr_list2.append(np.corrcoef(df2[feature], df2["preds_neutralized"])[0,1])
corr_series2 = pd.Series(corr_list2, index=features)
corr_series2.describe()<jupyter_output><empty_output><jupyter_text>Our overall score goes down, but the scores are more consistent than before. This leads to a higher sharpe<jupyter_code>unbalanced_scores_per_era = df2.groupby("era").apply(lambda d: np.corrcoef(d["preds"], d[target])[0,1])
balanced_scores_per_era = df2.groupby("era").apply(lambda d: np.corrcoef(d["preds_neutralized"], d[target])[0,1])

print(f"score for high feature exposure: {unbalanced_scores_per_era.mean()}")
print(f"score for balanced feature expo: {balanced_scores_per_era.mean()}")

print(f"std for high feature exposure: {unbalanced_scores_per_era.std(ddof=0)}")
print(f"std for balanced feature expo: {balanced_scores_per_era.std(ddof=0)}")

print(f"sharpe for high feature exposure: {unbalanced_scores_per_era.mean()/unbalanced_scores_per_era.std(ddof=0)}")
print(f"sharpe for balanced feature expo: {balanced_scores_per_era.mean()/balanced_scores_per_era.std(ddof=0)}")
balanced_scores_per_era.describe()
unbalanced_scores_per_era.describe()<jupyter_output><empty_output><jupyter_text>We can also try to be more clever about what features we neutralize The "riskiest" features are the features the predictions are most exposed to that also have the highest volatility<jupyter_code># compute feature correlations with target on the first half data
all_feature_corrs = df1.groupby('era').apply(lambda d: d[features].corrwith(d[target]))
# compute the volatility of the feature correlations
feature_corr_volatility = all_feature_corrs.std()
# calculate the feature exposures of the predictions
feature_exposure_list = []
for feature in features:
    feature_exposure_list.append(np.corrcoef(df2[feature], lgb_preds)[0,1])
feature_exposure_list = pd.Series(feature_exposure_list, index=features)
# get list of 100 riskiest features
riskiest_features = (feature_exposure_list.abs()*feature_corr_volatility).sort_values()[-100:].index.tolist()
# 50% neutralize to the riskiest features
df2["preds_neutralized_riskiest_100"] = neutralize(df2, 
                                                  columns=["preds"], 
                                                  neutralizers=riskiest_features, 
                                                  proportion=0.5, # 50% neutralize only those 100 features
                                                  normalize=True, 
                                                  era_col='era')

scaler = MinMaxScaler()
df2["preds_neutralized_riskiest_100"] = scaler.fit_transform(df2[["preds_neutralized_riskiest_100"]]) # transform back to 0-1<jupyter_output>0573<jupyter_text>The score goes down less and the sharpe goes up more than 50% neutralizing all the features<jupyter_code>unbalanced_scores_per_era = df2.groupby("era").apply(lambda d: np.corrcoef(d["preds"], d[target])[0,1])
balanced_scores_per_era = df2.groupby("era").apply(lambda d: np.corrcoef(d["preds_neutralized_riskiest_100"], d[target])[0,1])

print(f"score for high feature exposure: {unbalanced_scores_per_era.mean()}")
print(f"score for balanced feature expo: {balanced_scores_per_era.mean()}")

print(f"std for high feature exposure: {unbalanced_scores_per_era.std(ddof=0)}")
print(f"std for balanced feature expo: {balanced_scores_per_era.std(ddof=0)}")

print(f"sharpe for high feature exposure: {unbalanced_scores_per_era.mean()/unbalanced_scores_per_era.std(ddof=0)}")
print(f"sharpe for balanced feature expo: {balanced_scores_per_era.mean()/balanced_scores_per_era.std(ddof=0)}")<jupyter_output>score for high feature exposure: 0.05233661570835508
score for balanced feature expo: 0.04784620200240574
std for high feature exposure: 0.022859528419203096
std for balanced feature expo: 0.01865889557683214
sharpe for high feature exposure: 2.2894879871795526
sharpe for balanced feature expo: 2.5642569146382965<jupyter_text>Interestingly, the maximal feature exposure magnitude is actually lower (~0.28 vs ~0.31) than when 50% neutralizing all features<jupyter_code>corr_list3 = []
for feature in features:
    corr_list3.append(np.corrcoef(df2[feature], df2["preds_neutralized_riskiest_100"])[0,1])
corr_series3 = pd.Series(corr_list3, index=features)
corr_series3.describe()
corr_series3.describe()<jupyter_output><empty_output><jupyter_text>Training on the alternative targets can result in models with different patterns of feature exposures This is a great way to go after MMC and very useful for ensembling<jupyter_code>lgb = lightgbm.LGBMRegressor(n_jobs=-1, learning_rate=0.01, n_estimators=1000, max_depth=5)
lgb.fit(df1[features], df1['target'])
lgb_jerome = lightgbm.LGBMRegressor(n_jobs=-1, learning_rate=0.01, n_estimators=1000, max_depth=5)
lgb_jerome.fit(df1[features], df1['target_jerome_20'])
# The feature importances for each model are about 63% correlated
plt.figure(figsize = (16,6))
plt.subplot(1,2,1)
plt.plot(lgb.feature_importances_, alpha=.5)
plt.plot(lgb_jerome.feature_importances_, alpha=.5)
plt.subplot(1,2,2)
plt.scatter(lgb.feature_importances_, lgb_jerome.feature_importances_)

print(np.corrcoef([lgb.feature_importances_, lgb_jerome.feature_importances_])[0,1])
lgb_preds = lgb.predict(df2[features])
lgb_jerome_preds = lgb_jerome.predict(df2[features])
# The predictions are about 75% correlated
np.corrcoef([lgb_preds, lgb_jerome_preds])[0,1]
df2["preds"] = lgb_preds
df2["preds_jerome"] = lgb_jerome_preds
# 50% neutralize to riskest features for each model as above
for preds in ["preds", "preds_jerome"]:
    feature_exposure_list = []
    for feature in features:
        feature_exposure_list.append(np.corrcoef(df2[feature], df2[preds])[0,1])
    feature_exposure_list = pd.Series(feature_exposure_list, index=features)

    riskiest_features = (feature_exposure_list.abs()*feature_corr_volatility).sort_values()[-100:].index.tolist()
    
    df2[f"{preds}_neutralized_riskiest_100"] = neutralize(df2, 
                                                  columns=[f"{preds}"], 
                                                  neutralizers=riskiest_features, 
                                                  proportion=0.5, # 50% neutralize only those 100 features
                                                  normalize=True, 
                                                  era_col='era')
# create ensemble
df2["preds_ens"] = (df2["preds"] + df2["preds_jerome"])/2
scores = df2.groupby("era").apply(lambda d: np.corrcoef(d["preds"], d[target])[0,1])
scores_jerome = df2.groupby("era").apply(lambda d: np.corrcoef(d["preds_jerome"], d[target])[0,1])

scores_derisked = df2.groupby("era").apply(lambda d: np.corrcoef(d["preds_neutralized_riskiest_100"], d[target])[0,1])
scores_jerome_derisked = df2.groupby("era").apply(lambda d: np.corrcoef(d["preds_jerome_neutralized_riskiest_100"], d[target])[0,1])

scores_ens = df2.groupby("era").apply(lambda d: np.corrcoef(d["preds_ens"], d[target])[0,1])

print(f"score for target: {scores.mean()}")
print(f"score for target_jerome: {scores_jerome.mean()}")
print()
print(f"sharpe for target: {scores.mean()/scores.std(ddof=0)}")
print(f"sharpe for target_jerome: {scores_jerome.mean()/scores_jerome.std(ddof=0)}")
print()
print(f"score for target after neutralization: {scores_derisked.mean()}")
print(f"score for target_jerome after neutralization: {scores_jerome_derisked.mean()}")
print()
print(f"sharpe for target after neutralization: {scores_derisked.mean()/scores_derisked.std(ddof=0)}")
print(f"sharpe for target_jerome after neutralization: {scores_jerome_derisked.mean()/scores_jerome_derisked.std(ddof=0)}")
print()
print(f"score for target + target_jerome ensemble: {scores_ens.mean()}")
print(f"sharpe for target + target_jerome ensemble: {scores_ens.mean()/scores_ens.std(ddof=0)}")<jupyter_output>score for target: 0.05233661570835508
score for target_jerome: 0.05293083310575564

sharpe for target: 2.2894879871795526
sharpe for target_jerome: 1.9501765858045648

score for target after neutralization: 0.04784620200240574
score for target_jerome after neutralization: 0.050838690600204337

sharpe for target after neutralization: 2.5642569146382965
sharpe for target_jerome after neutralization: 2.2734387549255732

score for target + target_jerome ensemble: 0.05627261521740282
sharpe for target + target_jerome ensemble: 2.1496760967783635<jupyter_text>We see here that the model trained on `target_jerome` actually predicts `target` slightly better than the model trained on `target` (albeit with a lower Sharpe) both before and after neutralization of the riskiest features<jupyter_code># calculate correlation of neutralized predictions
df2[[f"preds_neutralized_riskiest_100", f"preds_jerome_neutralized_riskiest_100"]].corr()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1197984449561184, 'fraction_numerical': 0.06657323055360898, 'mean_word_length': 4.686527514231499, 'pattern_counts': {'":': 2, '<': 89, '<?xml version=': 0, '>': 86, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '262383', 'license': 'MIT', 'n_tokens_mistral': 11709, 'n_tokens_neox': 10314, 'n_words': 2393, 'path': 'docs/analysis_and_tips.ipynb', 'repo_name': 'goldinger/numerai-tournament', 'size': '29967'}