From 8db76568e4101b0fe7426cd95e6279244edf84af Mon Sep 17 00:00:00 2001 From: pepicello Date: Sun, 9 Aug 2020 17:34:50 +0100 Subject: [PATCH] ENH: add expected max drawdown function (#90) --- empyrical/__init__.py | 1 + empyrical/stats.py | 115 +++++++++++++++++++++++++++++++++- empyrical/tests/test_stats.py | 13 ++++ 3 files changed, 128 insertions(+), 1 deletion(-) diff --git a/empyrical/__init__.py b/empyrical/__init__.py index cf0168b..7268d9a 100644 --- a/empyrical/__init__.py +++ b/empyrical/__init__.py @@ -43,6 +43,7 @@ downside_risk, excess_sharpe, max_drawdown, + expected_max_drawdown, omega_ratio, roll_alpha, roll_alpha_aligned, diff --git a/empyrical/stats.py b/empyrical/stats.py index 45e7d69..d3a5c52 100644 --- a/empyrical/stats.py +++ b/empyrical/stats.py @@ -19,7 +19,7 @@ import pandas as pd import numpy as np from math import pow -from scipy import stats, optimize +from scipy import stats, optimize, interpolate from six import iteritems from sys import float_info @@ -405,6 +405,119 @@ def max_drawdown(returns, out=None): roll_max_drawdown = _create_unary_vectorized_roll_function(max_drawdown) +def expected_max_drawdown(mu, sigma, t, gbm=False): + """ + Determines the expected maximum drawdown of a brownian motion, + given drift and diffusion + + If a geometric Brownian motion with stochastic differential equation + dS(t) = Mu0 * S(t) * dt + Sigma0 * S(t) * dW(t) , + it converts to the form here by Ito's Lemma with X(t) = log(S(t)) such that + Mu = Mu0 - 0.5 * Sigma0^2 + Sigma = Sigma0 . + + Parameters + ---------- + mu : float + The drift term of a Brownian motion with drift. + sigma : float + The diffusion term of a Brownian motion with drift. + t : float + A time period of interest + gbp : bool + If true, compute for geometric brownian motion + + Returns + ------- + expected_max_drawdown : float + + Note + ----- + See http://www.cs.rpi.edu/~magdon/ps/journal/drawdown_journal.pdf + for more details. + """ + + if gbm: + new_mu = mu - 0.5 * sigma**2 + else: + new_mu = mu + + def emdd_qp(x): + """ Q function based on lookup table, for positive drifts """ + A = [0.0005, 0.001, 0.0015, 0.002, 0.0025, 0.005, + 0.0075, 0.01, 0.0125, 0.015, 0.0175, 0.02, + 0.0225, 0.025, 0.0275, 0.03, 0.0325, 0.035, + 0.0375, 0.04, 0.0425, 0.045, 0.0475, 0.05, + 0.055, 0.06, 0.065, 0.07, 0.075, 0.08, + 0.085, 0.09, 0.095, 0.1, 0.15, 0.2, + 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, + 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, + 4.0, 4.5, 5.0, 10.0, 15.0, 20.0, + 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, + 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, + 400.0, 450.0, 500.0, 1000.0, 1500.0, 2000.0, + 2500.0, 3000.0, 3500.0, 4000.0, 4500.0, 5000.0] + B = [0.01969, 0.027694, 0.033789, 0.038896, 0.043372, 0.060721, + 0.073808, 0.084693, 0.094171, 0.102651, 0.110375, 0.117503, + 0.124142, 0.130374, 0.136259, 0.141842, 0.147162, 0.152249, + 0.157127, 0.161817, 0.166337, 0.170702, 0.174924, 0.179015, + 0.186842, 0.194248, 0.201287, 0.207999, 0.214421, 0.220581, + 0.226505, 0.232212, 0.237722, 0.24305, 0.288719, 0.325071, + 0.355581, 0.382016, 0.405415, 0.426452, 0.445588, 0.463159, + 0.588642, 0.668992, 0.72854, 0.775976, 0.815456, 0.849298, + 0.878933, 0.905305, 0.92907, 1.088998, 1.184918, 1.253794, + 1.307607, 1.351794, 1.389289, 1.42186, 1.450654, 1.476457, + 1.647113, 1.747485, 1.818873, 1.874323, 1.919671, 1.958037, + 1.991288, 2.02063, 2.046885, 2.219765, 2.320983, 2.392826, + 2.448562, 2.494109, 2.532622, 2.565985, 2.595416, 2.621743] + if x > 5000: + return 0.25 * np.log(x) + 0.49088 + elif x < 0.0005: + return 0.5 * np.sqrt(np.pi * x) + else: + return interpolate.interp1d(A, B)([x])[0] + + def emdd_qn(x): + """ Q function based on lookup table, for negative drifts """ + A = [0.0005, 0.001, 0.0015, 0.002, 0.0025, 0.005, + 0.0075, 0.01, 0.0125, 0.015, 0.0175, 0.02, + 0.0225, 0.025, 0.0275, 0.03, 0.0325, 0.035, + 0.0375, 0.04, 0.0425, 0.045, 0.0475, 0.05, + 0.055, 0.06, 0.065, 0.07, 0.075, 0.08, + 0.085, 0.09, 0.095, 0.1, 0.15, 0.2, + 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, + 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, + 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, + 3.75, 4.0, 4.25, 4.5, 4.75, 5.0] + B = [0.019965, 0.028394, 0.034874, 0.040369, 0.045256, 0.064633, + 0.079746, 0.092708, 0.104259, 0.114814, 0.124608, 0.133772, + 0.142429, 0.150739, 0.158565, 0.166229, 0.173756, 0.180793, + 0.187739, 0.194489, 0.201094, 0.207572, 0.213877, 0.220056, + 0.231797, 0.243374, 0.254585, 0.265472, 0.27607, 0.286406, + 0.296507, 0.306393, 0.316066, 0.325586, 0.413136, 0.491599, + 0.564333, 0.633007, 0.698849, 0.762455, 0.824484, 0.884593, + 1.17202, 1.44552, 1.70936, 1.97074, 2.22742, 2.48396, + 2.73676, 2.99094, 3.24354, 3.49252, 3.74294, 3.99519, + 4.24274, 4.49238, 4.73859, 4.99043, 5.24083, 5.49882] + if x > 5: + return x + 0.5 + elif x < 0.0005: + return 0.5 * np.sqrt(np.pi * x) + else: + return interpolate.interp1d(A, B)([x])[0] + + if ((not np.isfinite(new_mu)) | (not np.isfinite(sigma)) | (sigma <= 0)): + return np.nan + else: + alpha = new_mu / (2 * sigma**2) + if new_mu == 0: + return - np.sqrt(np.pi / 2) * sigma * np.sqrt(t) + elif new_mu > 0: + return - emdd_qp(alpha * new_mu * t) / alpha + else: + return emdd_qn(alpha * new_mu * t) / alpha + + def annual_return(returns, period=DAILY, annualization=None): """ Determines the mean annual growth rate of returns. This is equivilent diff --git a/empyrical/tests/test_stats.py b/empyrical/tests/test_stats.py index 7dd5e47..8b42213 100644 --- a/empyrical/tests/test_stats.py +++ b/empyrical/tests/test_stats.py @@ -282,6 +282,19 @@ def test_max_drawdown_translation(self, returns, constant): assert max_dd <= raised_dd assert depressed_dd <= max_dd + @parameterized.expand([ + (0, 1, 1, False, -1.2533141373155001), + (1, 1, 1, False, -0.926318), + (-1, 1, 1, False, -1.769186), + (0, 1, 1, True, -1.477444), + ]) + def test_expected_max_drawdown(self, mu, sigma, t, gbm, expected): + assert_almost_equal( + self.empyrical.expected_max_drawdown(mu, sigma, t, gbm), + expected, + DECIMAL_PLACES, + ) + @parameterized.expand([ (mixed_returns, empyrical.DAILY, 1.9135925373194231), (weekly_returns, empyrical.WEEKLY, 0.24690830513998208),