feat(Metrics): added probabilistic sharpe ratio (#82)

* feat(Metrics): added probabilistic sharpe ratio

* Apply suggestions from code review

run_pipeline.py

@@ -82,9 +82,16 @@ def pipeline(project_name:str, wandb, sweep:bool, model_config:dict, training_co
     level1_columns = results[[column for column in results.columns if 'Ensemble' not in column]]
     ensemble_columns = results[[column for column in results.columns if 'Ensemble' in column]]
 
-    print("Mean no of samples: ", results.loc['no_of_samples'].mean())
+    print("\n--------\n")
+    print("Benchmark buy-and-hold sharpe: ", round(results.loc['benchmark_sharpe'].mean(), 3))
+
+    print("Level-1: Number of samples evaluated: ", level1_columns.loc['no_of_samples'].sum())
     print("Mean Sharpe ratio for Level-1 models: ", round(level1_columns.loc['sharpe'].mean(), 3))
+    print("Mean Probabilistic Sharpe ratio for Level-1 models: ", round(level1_columns.loc['prob_sharpe'].mean(), 3))
+
+    print("Level-2 (Ensemble): Number of samples evaluated: ", ensemble_columns.loc['no_of_samples'].sum())
     print("Mean Sharpe ratio for Level-2 (Ensemble) models: ", round(ensemble_columns.loc['sharpe'].mean(), 3))
+    print("Mean Probabilistic Sharpe ratio for Level-2 (Ensemble) models: ", round(ensemble_columns.loc['prob_sharpe'].mean(), 3))
 
     if sweep:
         if wandb.run is not None:

utils/evaluate.py

@@ -1,6 +1,7 @@
 from typing import Literal
 from sklearn.metrics import mean_absolute_error, accuracy_score, r2_score, f1_score, precision_score, recall_score
-from quantstats.stats import expected_return, sharpe, skew, sortino
+from quantstats.stats import skew, sortino
+from utils.metrics import probabilistic_sharpe_ratio, sharpe_ratio, average_holding_period
 from utils.helpers import get_first_valid_return_index
 import pandas as pd
 import numpy as np
@@ -54,44 +55,48 @@ def evaluate_predictions(
     df = __preprocess(target_returns, y_pred, y_true, method, no_of_classes)
     
     scorecard = pd.Series()
-    if method == 'regression':
+    # we probably will not need regression models at all
-        scorecard.loc['RSQ'] = r2_score(df.target_returns, df.y_pred)
+    # if method == 'regression':
-        scorecard.loc['MAE'] = mean_absolute_error(df.target_returns, df.y_pred)
+    #     scorecard.loc['RSQ'] = r2_score(df.target_returns, df.y_pred)
-    elif method == 'classification':
+    #     scorecard.loc['MAE'] = mean_absolute_error(df.target_returns, df.y_pred)
-        scorecard.loc['RSQ'] = 0.
+    # elif method == 'classification':
-        scorecard.loc['MAE Matrix'] = 0.
+    #     scorecard.loc['RSQ'] = 0.
-    sign_true = df.sign_true.astype(int)
+    #     scorecard.loc['MAE Matrix'] = 0.
-    sign_pred = df.sign_pred.astype(int)
+    
     
     def count_non_zero(series: pd.Series) -> int:
         return len(series[series != 0])
-    scorecard.loc['no_of_samples'] = count_non_zero(y_pred[evaluate_from:])
+    no_of_samples = count_non_zero(df.y_pred)
-    scorecard.loc['sharpe'] = sharpe(df.result)
+    scorecard.loc['no_of_samples'] = no_of_samples
+    sharpe = sharpe_ratio(df.result)
+    scorecard.loc['sharpe'] = sharpe
+    benchmark_sharpe = sharpe_ratio(df.target_returns)
+    scorecard.loc['benchmark_sharpe'] = benchmark_sharpe
+    scorecard.loc['prob_sharpe'] = probabilistic_sharpe_ratio(sharpe, benchmark_sharpe, no_of_samples)
     scorecard.loc['sortino'] = sortino(df.result)
     scorecard.loc['skew'] = skew(df.result)
 
     labels = [1, -1] if no_of_classes == 'two' else [1, -1, 0]
     avg_type = 'weighted' if no_of_classes == 'two' else 'macro'
 
-    scorecard.loc['accuracy'] = accuracy_score(sign_true, sign_pred) * 100
+    scorecard.loc['accuracy'] = accuracy_score(df.sign_true, df.sign_pred) * 100
-    scorecard.loc['recall'] = recall_score(sign_true, sign_pred, labels = labels, average=avg_type)
+    scorecard.loc['recall'] = recall_score(df.sign_true, df.sign_pred, labels = labels, average=avg_type)
-    scorecard.loc['precision'] = precision_score(sign_true, sign_pred, labels = labels, average=avg_type)
+    scorecard.loc['precision'] = precision_score(df.sign_true, df.sign_pred, labels = labels, average=avg_type)
-    scorecard.loc['f1_score'] = f1_score(sign_true, sign_pred, labels = labels, average=avg_type)
+    scorecard.loc['f1_score'] = f1_score(df.sign_true, df.sign_pred, labels = labels, average=avg_type)
     scorecard.loc['edge'] = df.result.mean()
     scorecard.loc['noise'] = df.y_pred.diff().abs().mean()
     scorecard.loc['edge_to_noise'] = scorecard.loc['edge'] / scorecard.loc['noise']
     
-    for index, row in sign_true.value_counts().iteritems():
+    for index, row in df.sign_true.value_counts().iteritems():
-        scorecard.loc['sign_true_ratio_' + str(index)] = row / len(sign_true)
+        scorecard.loc['sign_true_ratio_' + str(index)] = row / len(df.sign_true)
-    
-    for index, row in sign_pred.value_counts().iteritems():
-        scorecard.loc['sign_pred_ratio_' + str(index)] = row / len(sign_pred)
     
+    for index, row in df.sign_pred.value_counts().iteritems():
+        scorecard.loc['sign_pred_ratio_' + str(index)] = row / len(df.sign_pred)
 
-    if method == 'regression':
+    # if method == 'regression':
-        scorecard.loc['edge_to_mae'] = scorecard.loc['edge'] / scorecard.loc['MAE']
+    #     scorecard.loc['edge_to_mae'] = scorecard.loc['edge'] / scorecard.loc['MAE']
-    elif method == 'classification':
+    # elif method == 'classification':
-        scorecard.loc['edge_to_mae'] = 0.
+    #     scorecard.loc['edge_to_mae'] = 0.
 
     scorecard = scorecard.round(3)
     print("Model name: ", model_name)

utils/metrics.py

@@ -0,0 +1,336 @@
+import warnings
+import pandas as pd
+import scipy.stats as ss
+import numpy as np
+
+def timing_of_flattening_and_flips(target_positions: pd.Series) -> pd.DatetimeIndex:
+    """
+    Advances in Financial Machine Learning, Snippet 14.1, page 197
+    Derives the timestamps of flattening or flipping trades from a pandas series
+    of target positions. Can be used for position changes analysis, such as
+    frequency and balance of position changes.
+    Flattenings - times when open position is bing closed (final target position is 0).
+    Flips - times when positive position is reversed to negative and vice versa.
+    :param target_positions: (pd.Series) Target position series with timestamps as indices
+    :return: (pd.DatetimeIndex) Timestamps of trades flattening, flipping and last bet
+    """
+
+    empty_positions = target_positions[(target_positions == 0)].index  # Empty positions index
+    previous_positions = target_positions.shift(1)  # Timestamps pointing at previous positions
+
+    # Index of positions where previous one wasn't empty
+    previous_positions = previous_positions[(previous_positions != 0)].index
+
+    # FLATTENING - if previous position was open, but current is empty
+    flattening = empty_positions.intersection(previous_positions)
+
+    # Multiplies current position with value of next one
+    multiplied_posions = target_positions.iloc[1:] * target_positions.iloc[:-1].values
+
+    # FLIPS - if current position has another direction compared to the next
+    flips = multiplied_posions[(multiplied_posions < 0)].index
+    flips_and_flattenings = flattening.union(flips).sort_values()
+    if target_positions.index[-1] not in flips_and_flattenings:  # Appending with last bet
+        flips_and_flattenings = flips_and_flattenings.append(target_positions.index[-1:])
+
+    return flips_and_flattenings
+
+
+def average_holding_period(target_positions: pd.Series) -> float:
+    """
+    Advances in Financial Machine Learning, Snippet 14.2, page 197
+    Estimates the average holding period (in days) of a strategy, given a pandas series
+    of target positions using average entry time pairing algorithm.
+    Idea of an algorithm:
+    * entry_time = (previous_time * weight_of_previous_position + time_since_beginning_of_trade * increase_in_position )
+      / weight_of_current_position
+    * holding_period ['holding_time' = time a position was held, 'weight' = weight of position closed]
+    * res = weighted average time a trade was held
+    :param target_positions: (pd.Series) Target position series with timestamps as indices
+    :return: (float) Estimated average holding period, NaN if zero or unpredicted
+    """
+
+    holding_period = pd.DataFrame(columns=['holding_time', 'weight'])
+    entry_time = 0
+    position_difference = target_positions.diff()
+
+    # Time elapsed from the starting time for each position
+    time_difference = (target_positions.index - target_positions.index[0]) / np.timedelta64(1, 'D')
+    for i in range(1, target_positions.size):
+
+        # Increased or unchanged position
+        if float(position_difference.iloc[i] * target_positions.iloc[i - 1]) >= 0:
+            if float(target_positions.iloc[i]) != 0:  # And not an empty position
+                entry_time = (entry_time * target_positions.iloc[i - 1] +
+                              time_difference[i] * position_difference.iloc[i]) / target_positions.iloc[i]
+
+        # Decreased
+        if float(position_difference.iloc[i] * target_positions.iloc[i - 1]) < 0:
+            hold_time = time_difference[i] - entry_time
+
+            # Flip of a position
+            if float(target_positions.iloc[i] * target_positions.iloc[i - 1]) < 0:
+                weight = abs(target_positions.iloc[i - 1])
+                holding_period.loc[target_positions.index[i], ['holding_time', 'weight']] = (hold_time, weight)
+                entry_time = time_difference[i]  # Reset entry time
+
+            # Only a part of position is closed
+            else:
+                weight = abs(position_difference.iloc[i])
+                holding_period.loc[target_positions.index[i], ['holding_time', 'weight']] = (hold_time, weight)
+
+    if float(holding_period['weight'].sum()) > 0:  # If there were closed trades at all
+        avg_holding_period = float((holding_period['holding_time'] * \
+                                    holding_period['weight']).sum() / holding_period['weight'].sum())
+    else:
+        avg_holding_period = float('nan')
+
+    return avg_holding_period
+
+
+def bets_concentration(returns: pd.Series) -> float:
+    """
+    Advances in Financial Machine Learning, Snippet 14.3, page 201
+    Derives the concentration of returns from given pd.Series of returns.
+    Algorithm is based on Herfindahl-Hirschman Index where return weights
+    are taken as an input.
+    :param returns: (pd.Series) Returns from bets
+    :return: (float) Concentration of returns (nan if less than 3 returns)
+    """
+
+    if returns.size <= 2:
+        return float('nan')  # If less than 3 bets
+    weights = returns / returns.sum()  # Weights of each bet
+    hhi = (weights ** 2).sum()  # Herfindahl-Hirschman Index for weights
+    hhi = float((hhi - returns.size ** (-1)) / (1 - returns.size ** (-1)))
+
+    return hhi
+
+
+def all_bets_concentration(returns: pd.Series, frequency: str = 'M') -> tuple:
+    """
+    Advances in Financial Machine Learning, Snippet 14.3, page 201
+    Given a pd.Series of returns, derives concentration of positive returns, negative returns
+    and concentration of bets grouped by time intervals (daily, monthly etc.).
+    If after time grouping less than 3 observations, returns nan.
+    Properties or results:
+    * low positive_concentration ⇒ no right fat-tail of returns (desirable)
+    * low negative_concentration ⇒ no left fat-tail of returns (desirable)
+    * low time_concentration ⇒ bets are not concentrated in time, or are evenly concentrated (desirable)
+    * positive_concentration == 0 ⇔ returns are uniform
+    * positive_concentration == 1 ⇔ only one non-zero return exists
+    :param returns: (pd.Series) Returns from bets
+    :param frequency: (str) Desired time grouping frequency from pd.Grouper
+    :return: (tuple of floats) Concentration of positive, negative and time grouped concentrations
+    """
+
+    # Concentration of positive returns per bet
+    positive_concentration = bets_concentration(returns[returns >= 0])
+
+    # Concentration of negative returns per bet
+    negative_concentration = bets_concentration(returns[returns < 0])
+
+    # Concentration of bets/time period (month by default)
+    time_concentration = bets_concentration(returns.groupby(pd.Grouper(freq=frequency)).count())
+
+    return (positive_concentration, negative_concentration, time_concentration)
+
+
+def drawdown_and_time_under_water(returns: pd.Series, dollars: bool = False) -> tuple:
+    """
+    Advances in Financial Machine Learning, Snippet 14.4, page 201
+    Calculates drawdowns and time under water for pd.Series of either relative price of a
+    portfolio or dollar price of a portfolio.
+    Intuitively, a drawdown is the maximum loss suffered by an investment between two consecutive high-watermarks.
+    The time under water is the time elapsed between an high watermark and the moment the PnL (profit and loss)
+    exceeds the previous maximum PnL. We also append the Time under water series with period from the last
+    high-watermark to the last return observed.
+    Return details:
+    * Drawdown series index is the time of a high watermark and the value of a
+      drawdown after it.
+    * Time under water index is the time of a high watermark and how much time
+      passed till the next high watermark in years. Also includes time between
+      the last high watermark and last observation in returns as the last element.
+    :param returns: (pd.Series) Returns from bets
+    :param dollars: (bool) Flag if given dollar performance and not returns.
+                    If dollars, then drawdowns are in dollars, else as a %.
+    :return: (tuple of pd.Series) Series of drawdowns and time under water
+    """
+
+    frame = returns.to_frame('pnl')
+    frame['hwm'] = returns.expanding().max()  # Adding high watermarks as column
+
+    # Grouped as min returns by high watermarks
+    high_watermarks = frame.groupby('hwm').min().reset_index()
+    high_watermarks.columns = ['hwm', 'min']
+
+    # Time high watermark occurred
+    high_watermarks.index = frame['hwm'].drop_duplicates(keep='first').index
+
+    # Picking ones that had a drawdown after high watermark
+    high_watermarks = high_watermarks[high_watermarks['hwm'] > high_watermarks['min']]
+    if dollars:
+        drawdown = high_watermarks['hwm'] - high_watermarks['min']
+    else:
+        drawdown = 1 - high_watermarks['min'] / high_watermarks['hwm']
+
+    time_under_water = ((high_watermarks.index[1:] - high_watermarks.index[:-1]) / np.timedelta64(1, 'Y')).values
+
+    # Adding also period from last High watermark to last return observed.
+    time_under_water = np.append(time_under_water,
+                                 (returns.index[-1] - high_watermarks.index[-1]) / np.timedelta64(1, 'Y'))
+
+    time_under_water = pd.Series(time_under_water, index=high_watermarks.index)
+
+    return drawdown, time_under_water
+
+
+def sharpe_ratio(returns: pd.Series, entries_per_year: int = 252, risk_free_rate: float = 0) -> float:
+    """
+    Calculates annualized Sharpe ratio for pd.Series of normal or log returns.
+    Risk_free_rate should be given for the same period the returns are given.
+    For example, if the input returns are observed in 3 months, the risk-free
+    rate given should be the 3-month risk-free rate.
+    :param returns: (pd.Series) Returns - normal or log
+    :param entries_per_year: (int) Times returns are recorded per year (252 by default)
+    :param risk_free_rate: (float) Risk-free rate (0 by default)
+    :return: (float) Annualized Sharpe ratio
+    """
+
+    sharpe_r = (returns.mean() - risk_free_rate) / returns.std() * (entries_per_year) ** (1 / 2)
+
+    return sharpe_r
+
+
+def information_ratio(returns: pd.Series, benchmark: float = 0, entries_per_year: int = 252) -> float:
+    """
+    Calculates annualized information ratio for pd.Series of normal or log returns.
+    Benchmark should be provided as a return for the same time period as that between
+    input returns. For example, for the daily observations it should be the
+    benchmark of daily returns.
+    It is the annualized ratio between the average excess return and the tracking error.
+    The excess return is measured as the portfolio’s return in excess of the benchmark’s
+    return. The tracking error is estimated as the standard deviation of the excess returns.
+    :param returns: (pd.Series) Returns - normal or log
+    :param benchmark: (float) Benchmark for performance comparison (0 by default)
+    :param entries_per_year: (int) Times returns are recorded per year (252 by default)
+    :return: (float) Annualized information ratio
+    """
+
+    excess_returns = returns - benchmark
+    information_r = sharpe_ratio(excess_returns, entries_per_year)
+
+    return information_r
+
+
+def probabilistic_sharpe_ratio(observed_sr: float, benchmark_sr: float, number_of_returns: int,
+                               skewness_of_returns: float = 0, kurtosis_of_returns: float = 3) -> float:
+    """
+    Calculates the probabilistic Sharpe ratio (PSR) that provides an adjusted estimate of SR,
+    by removing the inflationary effect caused by short series with skewed and/or
+    fat-tailed returns.
+    Given a user-defined benchmark Sharpe ratio and an observed Sharpe ratio,
+    PSR estimates the probability that SR ̂is greater than a hypothetical SR.
+    - It should exceed 0.95, for the standard significance level of 5%.
+    - It can be computed on absolute or relative returns.
+    :param observed_sr: (float) Sharpe ratio that is observed
+    :param benchmark_sr: (float) Sharpe ratio to which observed_SR is tested against
+    :param number_of_returns: (int) Times returns are recorded for observed_SR
+    :param skewness_of_returns: (float) Skewness of returns (0 by default)
+    :param kurtosis_of_returns: (float) Kurtosis of returns (3 by default)
+    :return: (float) Probabilistic Sharpe ratio
+    """
+
+    test_value = ((observed_sr - benchmark_sr) * np.sqrt(number_of_returns - 1)) / \
+                  ((1 - skewness_of_returns * observed_sr + (kurtosis_of_returns - 1) / \
+                    4 * observed_sr ** 2)**(1 / 2))
+
+    if np.isnan(test_value):
+        warnings.warn('Test value is nan. Please check the input values.', UserWarning)
+        return test_value
+
+    if isinstance(test_value, complex):
+        warnings.warn('Output is a complex number. You may want to check the input skewness (too high), '
+                      'kurtosis (too low), or observed_sr values.', UserWarning)
+
+    if np.isinf(test_value):
+        warnings.warn('Test value is infinite. You may want to check the input skewness, '
+                      'kurtosis, or observed_sr values.', UserWarning)
+
+    probab_sr = ss.norm.cdf(test_value)
+
+    return probab_sr
+
+
+def deflated_sharpe_ratio(observed_sr: float, sr_estimates: list, number_of_returns: int,
+                          skewness_of_returns: float = 0, kurtosis_of_returns: float = 3,
+                          estimates_param: bool = False, benchmark_out: bool = False) -> float:
+    """
+    Calculates the deflated Sharpe ratio (DSR) - a PSR where the rejection threshold is
+    adjusted to reflect the multiplicity of trials. DSR is estimated as PSR[SR∗], where
+    the benchmark Sharpe ratio, SR∗, is no longer user-defined, but calculated from
+    SR estimate trails.
+    DSR corrects SR for inflationary effects caused by non-Normal returns, track record
+    length, and multiple testing/selection bias.
+    - It should exceed 0.95, for the standard significance level of 5%.
+    - It can be computed on absolute or relative returns.
+    Function allows the calculated SR benchmark output and usage of only
+    standard deviation and number of SR trails instead of full list of trails.
+    :param observed_sr: (float) Sharpe ratio that is being tested
+    :param sr_estimates: (list) Sharpe ratios estimates trials list or
+        properties list: [Standard deviation of estimates, Number of estimates]
+        if estimates_param flag is set to True.
+    :param  number_of_returns: (int) Times returns are recorded for observed_SR
+    :param skewness_of_returns: (float) Skewness of returns (0 by default)
+    :param kurtosis_of_returns: (float) Kurtosis of returns (3 by default)
+    :param estimates_param: (bool) Flag to use properties of estimates instead of full list
+    :param benchmark_out: (bool) Flag to output the calculated benchmark instead of DSR
+    :return: (float) Deflated Sharpe ratio or Benchmark SR (if benchmark_out)
+    """
+
+    # Calculating benchmark_SR from the parameters of estimates
+    if estimates_param:
+        benchmark_sr = sr_estimates[0] * \
+                       ((1 - np.euler_gamma) * ss.norm.ppf(1 - 1 / sr_estimates[1]) +
+                        np.euler_gamma * ss.norm.ppf(1 - 1 / sr_estimates[1] * np.e ** (-1)))
+
+    # Calculating benchmark_SR from a list of estimates
+    else:
+        benchmark_sr = np.array(sr_estimates).std() * \
+                       ((1 - np.euler_gamma) * ss.norm.ppf(1 - 1 / len(sr_estimates)) +
+                        np.euler_gamma * ss.norm.ppf(1 - 1 / len(sr_estimates) * np.e ** (-1)))
+
+    deflated_sr = probabilistic_sharpe_ratio(observed_sr, benchmark_sr, number_of_returns,
+                                             skewness_of_returns, kurtosis_of_returns)
+
+    if benchmark_out:
+        return benchmark_sr
+
+    return deflated_sr
+
+
+def minimum_track_record_length(observed_sr: float, benchmark_sr: float,
+                                skewness_of_returns: float = 0,
+                                kurtosis_of_returns: float = 3,
+                                alpha: float = 0.05) -> float:
+    """
+    Calculates the minimum track record length (MinTRL) - "How long should a track
+    record be in order to have statistical confidence that its Sharpe ratio is above
+    a given threshold?”
+    If a track record is shorter than MinTRL, we do not  have  enough  confidence
+    that  the  observed Sharpe ratio ̂is above the designated Sharpe ratio threshold.
+    MinTRLis expressed in terms of number of observations, not annual or calendar terms.
+    :param observed_sr: (float) Sharpe ratio that is being tested
+    :param benchmark_sr: (float) Sharpe ratio to which observed_SR is tested against
+    :param  number_of_returns: (int) Times returns are recorded for observed_SR
+    :param skewness_of_returns: (float) Skewness of returns (0 by default)
+    :param kurtosis_of_returns: (float) Kurtosis of returns (3 by default)
+    :param alpha: (float) Desired significance level (0.05 by default)
+    :return: (float) Minimum number of track records
+    """
+
+    track_rec_length = 1 + (1 - skewness_of_returns * observed_sr +
+                            (kurtosis_of_returns - 1) / 4 * observed_sr ** 2) * \
+                       (ss.norm.ppf(1 - alpha) / (observed_sr - benchmark_sr)) ** (2)
+
+    return track_rec_length