# Computing the cumulative distribution of max drawdown of random walk with drift

I am interested in the distribution of the maximum drawdown of a random walk: Let $X_0 = 0, X_{i+1} = X_i + Y_{i+1}$ where $Y_i \sim \mathcal{N}(\mu,1)$. The maximum drawdown after $n$ periods is $\max_{0 \le i \le j \le n} (X_i - X_j)$. A paper by Magdon-Ismail et. al. gives the distribution for maximum drawdown of a Brownian motion with drift. The expression involves an infinite sum which includes some terms defined only implicitly. I am having problems writing an implementation which converges. Is anyone aware of an alternative expression of the CDF or a reference implementation in code?

• How accurate do you need it? Could you just simulate the walk and avoid fully-functional solutions?
– kyle
Sep 20, 2010 at 23:20
• good point. I don't need atomic-physics level accuracy. in fact, 3 sigfigs is probably just fine.... Sep 21, 2010 at 0:11
• That would require about a million simulated random walks...
– whuber
Oct 27, 2010 at 22:40

This is an alternating sum. Each successive pair almost cancels; such pair-sums eventually decrease monotonically.

One approach, then, is to compute the sum in pairs where $n$ = {1,2}, {3,4}, {5,6}, etc. (Doing so eliminates a lot of floating point error, too.) Some more tricks can help:

(1) To solve $\tan(t) = t / \alpha$ for a positive constant $\alpha$, a good starting value to search--and an excellent approximation for the $n^\text{th}$ largest root-- is $t = (n + 1/2)\pi - \frac{\alpha}{(n + 1/2)\pi}$. I suspect Newton-Raphson should work really well.

(2) After a small number of initial terms, the sums of pairs start decreasing in size very, very consistently. The logarithms of the absolute values of exponentially-spaced pairs quickly decrease almost linearly. This means you can interpolate among a very small number of calculated pair-sums to estimate all the pair-sums you did not compute. For example, by computing the values for only pairs (2,3), (4,5), (8,9), (16,17), ..., (16384, 16385) and constructing the interpolating polynomial for these (thought of as the values of a function at 1, 2, ..., 14) and using the arguments $h = \mu = \sigma = 1$, I was able to achieve six-figure precision for the worst-case errors. (Even nicer, the errors oscillate in sign, suggesting the precision in the summed interpolated values might be quite a bit better than six figures.) You could probably estimate the limiting sum to good precision by extrapolating linearly off the end of these values (which translates to a power law) and integrating the extrapolating function out to infinity. To complete this example calculation you also need the first term. That gives six-figure precision by means of only 29 computed terms in the summation.

(3) Note that the function really depends on $h/\sigma$ and $\mu/\sigma$, not on all three of these variables independently. The dependence on $T$ is weak (as it should be); you might be content to fix its value throughout all your calculations.

(4) On top of all this, consider using some series-acceleration methods, like Aitken's method. A good accounting of this appears in Numerical Recipes.

(5) You can estimate the tail of the sum with an integral. Upon writing $\theta_n = (n + 1/2)\pi - 1/t_n$, the equation $\tan(\theta_n) = \theta_n / \alpha$ (with $\alpha = \mu h / \sigma^2$) can be solved for $t_n$, which is small, and then for $\theta_n$ by substituting back. Expanding the tangent in a Taylor series in $t_n$ gives the approximate solution

$$\theta_n = z - \frac{\alpha}{ z} - \frac{\alpha^2 - \alpha^3/3 }{z^3} + O\left((\frac{\alpha}{n})^5\right)$$

where $z = (n + 1/2)\pi$.

Provided $n$ is sufficiently large, the exponential factors of the form $1 - \exp(\frac{-\sigma^2 \theta_n^2 T}{2 h^2}) \exp(\frac{-\mu^2 T}{2 \sigma^2})$ become extremely close to 1 so you can neglect them. Typically these terms can be neglected even for small $n$ because $\theta_n^2$ is $\Theta\left(n^2\right)$, making the first exponential go to zero extremely quickly. (This happens once $n$ substantially exceeds $\alpha / T^{1/2}$. Do your calculations for large $T$ if you can!)

Using this expression for $\theta_n$ to sum the terms for $n$ and $n+1$ lets us approximate them (once all the smoke clears) as

$$\frac{2}{\pi n^2}-\frac{4}{\pi n^3}+\frac{13 \pi ^2+6 (4-3 \alpha ) \alpha }{2 \pi ^3 n^4}+O\left(\frac{1}{n^5}\right) \text{.}$$

Replacing the sum starting at $n = 2N$ by an integral over $N$ starting at $N - 1/4$ approximates the tail. (The integral has to be multiplied by a common factor of $\exp(-\alpha)$.) The error in the integral is $O(1/n^4)$. Thus, to achieve three significant figures you will typically need to compute around eight or so of the terms in the sum and then add this tail approximation.

• this is really great, and should go a long way towards the CDF. Above and Beyond badge material. Oct 29, 2010 at 2:21

You might start by looking at the drawdown distribution functions in fBasics. So you could easily simulate the brownian motion with drift and apply these functions as a start.

• +1 That's a pretty direct answer, considering these functions implement the formulas in the paper!
– whuber
Oct 28, 2010 at 15:43
• It seems like this package computes the expected max drawdown based on the paper, but does not compute the CDF. The paper gives 'shortcut' results, IIRC, to compute that expectation. Oct 28, 2010 at 22:39
• @shabbychef Sorry, I missed that nicety. I see how obtaining the entire CDF can be more useful than just knowing the expectation. (Financial risk is much more than just expected losses...) But now I feel a little bit better about the work I did to approximate the CDF!
– whuber
Oct 29, 2010 at 1:17