ENH: add expected max drawdown function (#90)

Author: pepicello

empyrical/__init__.py

@@ -43,6 +43,7 @@
     downside_risk,
     excess_sharpe,
     max_drawdown,
+    expected_max_drawdown,
     omega_ratio,
     roll_alpha,
     roll_alpha_aligned,

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,118 @@ 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):
+        """ Return value for Q functions 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):
+        """ Return value for Q functions 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

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),