1
2# This library is free software; you can redistribute it and/or
3# modify it under the terms of the GNU Library General Public
4# License as published by the Free Software Foundation; either
5# version 2 of the License, or (at your option) any later version.
6#
7# This library is distributed in the hope that it will be useful,
8# but WITHOUT ANY WARRANTY; without even the implied warranty of
9# MERCHANTABILITY or FITNESS FOR A PARTICULAR Description. See the
10# GNU Library General Public License for more details.
11#
12# You should have received a copy of the GNU Library General
13# Public License along with this library; if not, write to the
14# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
15# MA 02111-1307 USA
16
17
18################################################################################
19# FUNCTION:                 DESCRIPTION:
20#  maxddStats                Expectation of Drawdowns for BM with drift
21#  .maxddStats               Utility function called by "maxddStats"
22# FUNCTION:                 DISTRIBUTION AND RANDOM VARIATES:
23#  dmaxdd                    Density function of mean Max-Drawdowns
24#  pmaxdd                    Probability function of mean Max-Drawdowns
25#  rmaxdd                    Random Variates of mean Max-Drawdowns
26################################################################################
27
28
29maxddStats <-
30function(mean = 0, sd = 1, horizon = 1000)
31{
32    # A function implemented by Diethelm Wuertz
33
34    # Description:
35    #   Calculates Expectation Value E[D] of maximum Drawdowns of
36    #   Brownian Motion for a given drift "mu", variance "sigma",
37    #   and runtime "horizon" of the Brownian Motion process.
38
39    # Arguments:
40    #   mu - Drift of Brownian Motion, a numeric vector.
41    #   sigma - Standard Deviation of Brownian Motion, a numeric value.
42    #   horizon - Runtime of the process, the time horizon of the
43    #     investor, a numeric value.
44
45    # Details:
46    #   Interpolates the data given in the table of Appendix B
47    #   in Magdon-Ismail et al. [2003], "On the Maximum Drawdown
48    #   of Brownian Motion". The interpolation is done with R's
49    #   function spline().
50
51    # Notes.
52    #   This computes for a vector of horizon values.
53
54    # FUNCTION:
55
56    # Statistics:
57    ans = NULL
58    for (i in 1:length(horizon)) {
59        ans = c(ans, .maxddStats(mu = mean, sigma = sd, horizon = horizon[i]))
60    }
61
62    # Return Value:
63    ans
64}
65
66
67# ------------------------------------------------------------------------------
68
69
70.maxddStats <-
71function(mu = 0, sigma = 1, horizon = 1000)
72{
73    # A function implemented by Diethelm Wuertz
74
75    # Description:
76    #   Utility function called by "maxddStats"
77
78    # Arguments:
79    #   see function "maxddStats
80
81    # FUNCTION:
82
83    # Internal Function - POSITIVE CASE: mu > 0
84    # Left Table from Appendix B:
85    QP = function(x) {
86        gamma = sqrt(pi/8)
87        vqn = t(matrix(c(
88            0.0005,0.019690, 0.0010,0.027694, 0.0015,0.033789, 0.0020,0.038896,
89            0.0025,0.043372, 0.0050,0.060721, 0.0075,0.073808, 0.0100,0.084693,
90            0.0125,0.094171, 0.0150,0.102651, 0.0175,0.110375, 0.0200,0.117503,
91            0.0225,0.124142, 0.0250,0.130374, 0.0275,0.136259, 0.0300,0.141842,
92            0.0325,0.147162, 0.0350,0.152249, 0.0375,0.157127, 0.0400,0.161817,
93            0.0425,0.166337, 0.0450,0.170702, 0.0500,0.179015, 0.0600,0.194248,
94            0.0700,0.207999, 0.0800,0.220581, 0.0900,0.232212, 0.1000,0.243050,
95            0.2000,0.325071, 0.3000,0.382016, 0.4000,0.426452, 0.5000,0.463159,
96            1.5000,0.668992, 2.5000,0.775976, 3.5000,0.849298, 4.5000,0.905305,
97            10.000,1.088998, 20.000,1.253794, 30.000,1.351794, 40.000,1.421860,
98            50.000,1.476457, 150.00,1.747485, 250.00,1.874323, 350.00,1.958037,
99            450.00,2.020630, 1000.0,2.219765, 2000.0,2.392826, 3000.0,2.494109,
100            4000.0,2.565985, 5000.0,2.621743), 2))
101        # Interpolation:
102        if (x < 0.0005) { y = gamma*sqrt(2*x) }
103        if (x >= 0.0005 & x <= 5000) {
104            y = spline(log(vqn[,1]), vqn[,2], n = 1, xmin = log(x),
105                xmax = log(x))$y
106        }
107        if (x > 5000) { y = 0.25*log(x) + 0.49088}
108        # Return Value:
109        y
110    }
111
112    # Internal Function - NEGATIVE CASE: mu < 0
113    # Right Table from Appendix B:
114    QN = function(x) {
115        gamma = sqrt(pi/8)
116        vqn = t(matrix(c(
117            0.0005,0.019965, 0.0010,0.028394, 0.0015,0.034874, 0.0020,0.040369,
118            0.0025,0.045256, 0.0050,0.064633, 0.0075,0.079746, 0.0100,0.092708,
119            0.0125,0.104259, 0.0150,0.114814, 0.0175,0.124608, 0.0200,0.133772,
120            0.0225,0.142429, 0.0250,0.150739, 0.0275,0.158565, 0.0300,0.166229,
121            0.0325,0.173756, 0.0350,0.180793, 0.0375,0.187739, 0.0400,0.194489,
122            0.0425,0.201094, 0.0450,0.207572, 0.0475,0.213877, 0.0500,0.220056,
123            0.0550,0.231797, 0.0600,0.243374, 0.0650,0.254585, 0.0700,0.265472,
124            0.0750,0.276070, 0.0800,0.286406, 0.0850,0.296507, 0.0900,0.306393,
125            0.0950,0.316066, 0.1000,0.325586, 0.1500,0.413136, 0.2000,0.491599,
126            0.2500,0.564333, 0.3000,0.633007, 0.3500,0.698849, 0.4000,0.762455,
127            0.5000,0.884593, 1.0000,1.445520, 1.5000,1.970740, 2.0000,2.483960,
128            2.5000,2.990940, 3.0000,3.492520, 3.5000,3.995190, 4.0000,4.492380,
129            4.5000,4.990430, 5.0000,5.498820), 2))
130        # Interpolation:
131        if (x < 0.0005) { y = gamma*sqrt(2*x) }
132        if (x >= 0.0005 & x <= 5) {
133            y = spline(vqn[,1], vqn[,2], n = 1, xmin = x, xmax = x)$y }
134        if (x > 5) { y = x + 1/2 }
135        # Return Value:
136        y
137    }
138
139    # Result:
140    ED = NULL
141    for (i in 1:length(mu)) {
142        if (mu[i] == 0) {
143            gamma = sqrt(pi/8)
144            ED[i] = 2 * gamma * sigma * sqrt(horizon)
145        } else {
146            x = mu[i]^2 * horizon / ( 2 * sigma^2 )
147            if (mu[i] > 0) { ED[i] = ( 2* sigma^2 / mu[i] ) * QP(x) }
148            if (mu[i] < 0) { ED[i] = - ( 2* sigma^2 / mu[i] ) * QN(x) }
149        }
150    }
151
152    # Return Value:
153    ED
154}
155
156
157################################################################################
158
159
160dmaxdd <-
161function(x, sd = 1, horizon = 100, N = 1000)
162{
163    # Description:
164    #   Calculates for a trendless Brownian process (mean=0) and
165    #   standard deviation "sd" the density from the probability
166    #   that the maximum drawdown "D" is larger or equal to "h"
167    #   in the interval [0,T], where "T" denotes the time horizon
168    #   of the investor.
169
170    # Arguments:
171    #   x - Vector of Drawdowns
172    #   sd - Standard Deviation
173    #   horizon - Time horizon of the investor
174    #   N - number of summands
175
176    # Notes:
177    #   The drift dependent case is not yet implemented!
178
179    # FUNCTION:
180
181    # Use "h", "sigma" to be conform with Magdon-Ismael et al. [2003]
182    h = x
183    sigma = sd
184
185    # Settings:
186    n = 1:N
187    pn = 2 * sin((n-0.5)*pi) * sigma^2 * (n-0.5) * pi * horizon
188    en = sigma^2 * (n-0.5)^2 * pi^2 * horizon / 2
189
190    # Loop over all Drawdowns:
191    result = rep(NA, times = length(h))
192    for (i in 1:length(h)) {
193        if (h[i] == 0) {
194            result[i] = 0
195        } else {
196            g = pn *  exp(-en/(h[i]^2)) / h[i]^3
197            result[i] = sum(g)
198        }
199    }
200
201    # Return Value:
202    result
203}
204
205
206# ------------------------------------------------------------------------------
207
208
209pmaxdd <-
210function(q, sd = 1, horizon = 100, N = 1000)
211{
212    # Description:
213    #   Calculates for a trendless Brownian process (mean=0)
214    #   with standard deviation "sd" the probability that the
215    #   maximum drawdown "D" is larger or equal to "h" in the
216    #   interval [0,T], where "T" denotes the time horizon of
217    #   the investor.
218
219    # Arguments:
220    #   q - Vector of Drawdowns
221    #   sd - Standard Deviation
222    #   horizon - Time horizon of the investor
223    #   N - number of summands
224
225    # Value:
226    #   GD(h) - eqn(14) In Magdon-Ismail et al.
227
228    # Notes:
229    #   The drift dependent case is not yet implemented!
230
231    # FUNCTION:
232
233    # Use "h", "sigma" to be conform with Magdon-Ismael et al. [2003]
234    h = q
235    sigma = sd
236
237    # Settings:
238    n = 1:N
239    pn = 2 * ( sin((n-0.5)*pi) / ((n-0.5)*pi) )
240    en = sigma^2 * (n-0.5)^2 * pi^2 * horizon / 2
241
242    # Loop over all Drawdowns:
243    result = rep(NA, times = length(h))
244    for (i in 1:length(h)) {
245        g = pn * ( 1 - exp(-en/(h[i]^2)) )
246        result[i] = sum(g)
247    }
248
249    # Return Value:
250    result
251}
252
253
254# ------------------------------------------------------------------------------
255
256
257rmaxdd <-
258function(n, mean = 0, sd = 1, horizon = 100)
259{
260    # Description:
261    #   Generates for a Brownian process with mean "mean" and
262    #   standard deviation "sd" random variates of maximum
263    #   Drawdowns.
264
265    # Arguments:
266    #   n - Number of Drawdown rvs
267    #   mean - Drift
268    #   sd - Standard Deviation
269    #   horizon - Time horizon of the investor
270
271    # FUNCTION:
272
273    # Simulation of "n" max Drawdowns "h":
274    result = NULL
275    for (i in 1:n) {
276        D = cumsum(rnorm(horizon, mean = mean, sd = sd))
277        result[i] = max(cummax(D)-D)
278    }
279
280    # Return Value:
281    result
282}
283
284
285################################################################################
286
287