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