Compilation: prefer glib functions over goffice equivalents
[gnumeric.git] / plugins / fn-derivatives / options.c
blob299776d597a7c37459ab01cef14e7250e21e262b
1 /* vim: set sw=8: -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
2 /*
3 * Options pricing
5 * Authors:
6 * Elliot Lee <sopwith@redhat.com> All initial work.
7 * Morten Welinder <terra@gnome.org> Port to new plugin framework.
8 * Cleanup.
9 * Hal Ashburner <hal_ashburner@yahoo.co.uk>
10 * Black Scholes Code re-structure, optional asset leakage paramaters,
11 * American approximations, alternative models to Black-Scholes
12 * and All exotic Options Functions.
14 * This program is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU General Public License as published by
16 * the Free Software Foundation; either version 2 of the License, or
17 * (at your option) any later version.
19 * This program is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
24 * You should have received a copy of the GNU General Public License
25 * along with this program; if not, see <https://www.gnu.org/licenses/>.
27 #include <gnumeric-config.h>
28 #include <gnumeric.h>
29 #include <gnm-i18n.h>
30 #include <goffice/goffice.h>
31 #include <gnm-plugin.h>
33 #include <func.h>
34 #include <mathfunc.h>
35 #include <sf-dpq.h>
36 #include <sf-gamma.h>
37 #include <value.h>
39 #include <math.h>
40 #include <string.h>
42 GNM_PLUGIN_MODULE_HEADER;
44 /* Some common decriptors */
45 #define DEF_ARG_CALL_PUT_FLAG { GNM_FUNC_HELP_ARG, F_("call_put_flag:'c' for a call and 'p' for a put") }
46 #define DEF_ARG_SPOT { GNM_FUNC_HELP_ARG, F_("spot:spot price") }
47 #define DEF_ARG_STRIKE { GNM_FUNC_HELP_ARG, F_("strike:strike price") }
48 #define DEF_ARG_TIME_MATURITY_Y { GNM_FUNC_HELP_ARG, F_("time:time to maturity in years") }
49 #define DEF_ARG_TIME_MATURITY_D { GNM_FUNC_HELP_ARG, F_("time:time to maturity in days") }
50 #define DEF_ARG_TIME_DIVIDEND { GNM_FUNC_HELP_ARG, F_("time_payout:time to dividend payout") }
51 #define DEF_ARG_TIME_EXPIRATION { GNM_FUNC_HELP_ARG, F_("time_exp:time to expiration") }
52 #define DEF_ARG_RATE_RISKFREE { GNM_FUNC_HELP_ARG, F_("rate:risk-free interest rate to the exercise date in percent") }
53 #define DEF_ARG_RATE_ANNUALIZED { GNM_FUNC_HELP_ARG, F_("rate:annualized interest rate") }
54 #define DEF_ARG_RATE_RISKFREE_ANN { GNM_FUNC_HELP_ARG, F_("rate:annualized risk-free interest rate") }
55 #define DEF_ARG_VOLATILITY { GNM_FUNC_HELP_ARG, F_("volatility:annualized volatility of the asset in percent for the period through to the exercise date") }
56 #define DEF_ARG_VOLATILITY_SHORT { GNM_FUNC_HELP_ARG, F_("volatility:annualized volatility of the asset") }
57 #define DEF_ARG_AMOUNT { GNM_FUNC_HELP_ARG, F_("d:amount of the dividend to be paid expressed in currency") }
58 #define DEF_ARG_CC_OPT { GNM_FUNC_HELP_ARG, F_("cost_of_carry:net cost of holding the underlying asset (for common stocks, the risk free rate less the dividend yield), defaults to 0") }
59 #define DEF_ARG_CC { GNM_FUNC_HELP_ARG, F_("cost_of_carry:net cost of holding the underlying asset") }
61 #define DEF_NOTE_UNITS { GNM_FUNC_HELP_NOTE, F_("The returned value will be expressed in the same units as @{strike} and @{spot}.")}
64 typedef enum {
65 OS_Call,
66 OS_Put,
67 OS_Error
68 } OptionSide;
70 typedef enum{
71 OT_Euro,
72 OT_Amer,
73 OT_Error
74 } OptionType;
76 static gnm_float opt_baw_call (gnm_float s, gnm_float x, gnm_float t, gnm_float r, gnm_float b, gnm_float v);
77 static gnm_float opt_baw_put (gnm_float s, gnm_float x, gnm_float t, gnm_float r, gnm_float b, gnm_float v);
78 static gnm_float NRA_c (gnm_float x, gnm_float t, gnm_float r, gnm_float b, gnm_float v);
79 static gnm_float NRA_p (gnm_float x, gnm_float t, gnm_float r, gnm_float b, gnm_float v);
80 static gnm_float opt_bjer_stens1_c (gnm_float s, gnm_float x, gnm_float t, gnm_float r, gnm_float b, gnm_float v);
81 /* static gnm_float opt_bjer_stens1_p (gnm_float s, gnm_float x, gnm_float t, gnm_float r, gnm_float b, gnm_float v); */
82 static gnm_float phi (gnm_float s, gnm_float t, gnm_float gamma, gnm_float H, gnm_float I, gnm_float r, gnm_float b, gnm_float v);
83 static gnm_float CriticalValueOptionsOnOptions (OptionSide side, gnm_float x1, gnm_float x2, gnm_float t,
84 gnm_float r, gnm_float b, gnm_float v);
85 static gnm_float opt_crit_val_chooser (gnm_float s,gnm_float xc,gnm_float xp,gnm_float t,
86 gnm_float tc, gnm_float tp, gnm_float r, gnm_float b, gnm_float v);
89 static OptionSide
90 option_side (char const *s)
92 if (s[0] == 'p' || s[0] == 'P')
93 return OS_Put;
94 else if (s[0] == 'c' || s[0] == 'C')
95 return OS_Call;
96 else
97 return OS_Error;
100 static OptionType
101 option_type (char const *s)
103 if (s[0] == 'a' || s[0] == 'A')
104 return OT_Amer;
105 else if (s[0] == 'e' || s[0] == 'E')
106 return OT_Euro;
107 else
108 return OT_Error;
111 /* The normal distribution function */
112 static gnm_float
113 ncdf (gnm_float x)
115 return pnorm (x, 0.0, 1.0, TRUE, FALSE);
118 static gnm_float
119 npdf (gnm_float x)
121 return dnorm (x, 0.0, 1.0, FALSE);
124 static int
125 Sgn (gnm_float a)
127 if ( a >0)
128 return 1;
129 else if (a < 0)
130 return -1;
131 else
132 return 0;
135 /* The cumulative bivariate normal distribution function */
136 static gnm_float
137 cum_biv_norm_dist1 (gnm_float a, gnm_float b, gnm_float rho)
139 gnm_float rho1, rho2, delta;
140 gnm_float a1, b1, sum = 0.0;
141 int i, j;
143 static const gnm_float x[] = {0.24840615, 0.39233107, 0.21141819, 0.03324666, 0.00082485334};
144 static const gnm_float y[] = {0.10024215, 0.48281397, 1.0609498, 1.7797294, 2.6697604};
145 a1 = a / gnm_sqrt (2.0 * (1 - (rho * rho)));
146 b1 = b / gnm_sqrt (2.0 * (1 - (rho * rho)));
148 if (a <= 0 && b <= 0 && rho <= 0) {
149 for (i = 0; i != 5; ++i) {
150 for (j = 0; j != 5; ++j) {
151 sum = sum + x[i] * x[j] * gnm_exp (a1 * (2.0 * y[i] - a1) + b1 * (2.0 *
152 y[j] - b1) + 2 * rho * (y[i] - a1) * (y[j] - b1));
155 return gnm_sqrt (1.0 - (rho * rho)) / M_PIgnum * sum;
156 } else if (a <= 0 && b >= 0 && rho >= 0)
157 return ncdf (a) - cum_biv_norm_dist1 (a,-b,-rho);
158 else if (a >= 0 && b <= 0 && rho >= 0)
159 return ncdf (b) - cum_biv_norm_dist1 (-a,b,-rho);
160 else if (a >= 0 && b >= 0 && rho <= 0)
161 return ncdf (a) + ncdf (b) - 1.0 + cum_biv_norm_dist1 (-a,-b,rho);
162 else if ((a * b * rho) > 0) {
163 rho1 = (rho * a - b) * Sgn (a) / gnm_sqrt ((a * a) - 2 * rho * a
164 * b + (b * b));
165 rho2 = (rho * b - a) * Sgn (b) / gnm_sqrt ((a * a) - 2 * rho * a
166 * b + (b * b));
167 delta = (1.0 - Sgn (a) * Sgn (b)) / 4.0;
168 return (cum_biv_norm_dist1 (a,0.0,rho1) +
169 cum_biv_norm_dist1 (b,0.0,rho2) -
170 delta);
172 return gnm_nan;
176 static GnmValue *
177 cum_biv_norm_dist(GnmFuncEvalInfo *ei, GnmValue const * const *argv)
179 gnm_float a = value_get_as_float (argv[0]);
180 gnm_float b = value_get_as_float (argv[1]);
181 gnm_float rho = value_get_as_float (argv[2]);
182 gnm_float result = cum_biv_norm_dist1 (a,b,rho);
184 if (gnm_isnan (result))
185 return value_new_error_NUM (ei->pos);
186 else
187 return value_new_float (result);
190 static GnmFuncHelp const help_cum_biv_norm_dist[] = {
191 { GNM_FUNC_HELP_NAME, F_("CUM_BIV_NORM_DIST:cumulative bivariate normal distribution")},
192 { GNM_FUNC_HELP_ARG, F_("a:limit for first random variable")},
193 { GNM_FUNC_HELP_ARG, F_("b:limit for second random variable")},
194 { GNM_FUNC_HELP_ARG, F_("rho:correlation of the two random variables")},
195 { GNM_FUNC_HELP_DESCRIPTION, F_("CUM_BIV_NORM_DIST calculates the probability that two standard "
196 "normal distributed random variables with correlation @{rho} are "
197 "respectively each less than @{a} and @{b}.")},
198 { GNM_FUNC_HELP_EXAMPLES, "=CUM_BIV_NORM_DIST(0,0,0.5)" },
199 { GNM_FUNC_HELP_END}
204 /* the generalized Black and Scholes formula*/
205 static gnm_float
206 opt_bs1 (OptionSide side,
207 gnm_float s, gnm_float x, gnm_float t, gnm_float r, gnm_float v,
208 gnm_float b)
210 gnm_float d1 = (gnm_log (s / x) + (b + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
211 gnm_float d2 = d1 - v * gnm_sqrt (t);
213 switch (side) {
214 case OS_Call:
215 return (s * gnm_exp ((b - r) * t) * ncdf (d1) -
216 x * gnm_exp (-r * t) * ncdf (d2));
217 case OS_Put:
218 return (x * gnm_exp (-r * t) * ncdf (-d2) -
219 s * gnm_exp ((b - r) * t) * ncdf (-d1));
220 default:
221 return gnm_nan;
226 static GnmValue *
227 opt_bs (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
229 OptionSide call_put = option_side (value_peek_string (argv[0]));
230 gnm_float s = value_get_as_float (argv[1]);
231 gnm_float x = value_get_as_float (argv[2]);
232 gnm_float t = value_get_as_float (argv[3]);
233 gnm_float r = value_get_as_float (argv[4]);
234 gnm_float v = value_get_as_float (argv[5]);
235 gnm_float b = argv[6] ? value_get_as_float (argv[6]) : 0;
236 gnm_float gfresult = opt_bs1 (call_put, s, x, t, r, v, b);
238 if (gnm_isnan (gfresult))
239 return value_new_error_NUM (ei->pos);
240 return value_new_float (gfresult);
243 static GnmFuncHelp const help_opt_bs[] = {
244 { GNM_FUNC_HELP_NAME, F_("OPT_BS:price of a European option")},
245 DEF_ARG_CALL_PUT_FLAG,
246 DEF_ARG_SPOT,
247 DEF_ARG_STRIKE,
248 DEF_ARG_TIME_MATURITY_Y,
249 DEF_ARG_RATE_RISKFREE,
250 DEF_ARG_VOLATILITY,
251 DEF_ARG_CC_OPT,
252 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_BS uses the Black-Scholes model to calculate "
253 "the price of a European option struck at @{strike} "
254 "on an asset with spot price @{spot}.")},
255 DEF_NOTE_UNITS,
256 { GNM_FUNC_HELP_SEEALSO, "OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_VEGA,OPT_BS_GAMMA"},
257 { GNM_FUNC_HELP_END}
260 /* Delta for the generalized Black and Scholes formula */
261 static gnm_float
262 opt_bs_delta1 (OptionSide side,
263 gnm_float s, gnm_float x, gnm_float t, gnm_float r, gnm_float v, gnm_float b)
265 gnm_float d1 =
266 (gnm_log (s / x) + (b + (v * v) / 2.0) * t) /
267 (v * gnm_sqrt (t));
269 switch (side) {
270 case OS_Call:
271 return gnm_exp ((b - r) * t) * ncdf (d1);
273 case OS_Put:
274 return gnm_exp ((b - r) * t) * (ncdf (d1) - 1.0);
276 default:
277 return gnm_nan;
282 static GnmValue *
283 opt_bs_delta (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
285 OptionSide call_put = option_side (value_peek_string (argv[0]));
286 gnm_float s = value_get_as_float (argv[1]);
287 gnm_float x = value_get_as_float (argv[2]);
288 gnm_float t = value_get_as_float (argv[3]);
289 gnm_float r = value_get_as_float (argv[4]);
290 gnm_float v = value_get_as_float (argv[5]);
291 gnm_float b = argv[6] ? value_get_as_float (argv[6]) : 0.0;
292 gnm_float gfresult = opt_bs_delta1 (call_put, s, x, t, r, v, b);
294 if (gnm_isnan (gfresult))
295 return value_new_error_NUM (ei->pos);
297 return value_new_float (gfresult);
300 static GnmFuncHelp const help_opt_bs_delta[] = {
301 { GNM_FUNC_HELP_NAME, F_("OPT_BS_DELTA:delta of a European option")},
302 DEF_ARG_CALL_PUT_FLAG,
303 DEF_ARG_SPOT,
304 DEF_ARG_STRIKE,
305 DEF_ARG_TIME_MATURITY_Y,
306 DEF_ARG_RATE_RISKFREE,
307 DEF_ARG_VOLATILITY,
308 DEF_ARG_CC_OPT,
309 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_BS_DELTA uses the Black-Scholes model to calculate "
310 "the 'delta' of a European option struck at @{strike} "
311 "on an asset with spot price @{spot}.")},
312 DEF_NOTE_UNITS,
313 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_VEGA,OPT_BS_GAMMA"},
314 { GNM_FUNC_HELP_END}
318 /* Gamma for the generalized Black and Scholes formula */
319 static gnm_float
320 opt_bs_gamma1 (gnm_float s,gnm_float x,gnm_float t,gnm_float r,gnm_float v,gnm_float b)
322 gnm_float d1;
324 d1 = (gnm_log (s / x) + (b + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
325 return gnm_exp ((b - r) * t) * npdf (d1) / (s * v * gnm_sqrt (t));
329 static GnmValue *
330 opt_bs_gamma (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
332 gnm_float s = value_get_as_float (argv[0]);
333 gnm_float x = value_get_as_float (argv[1]);
334 gnm_float t = value_get_as_float (argv[2]);
335 gnm_float r = value_get_as_float (argv[3]);
336 gnm_float v = value_get_as_float (argv[4]);
337 gnm_float b = argv[5] ? value_get_as_float (argv[5]) : 0.0;
338 gnm_float gfresult = opt_bs_gamma1 (s,x,t,r,v,b);
339 return value_new_float (gfresult);
342 static GnmFuncHelp const help_opt_bs_gamma[] = {
343 { GNM_FUNC_HELP_NAME, F_("OPT_BS_GAMMA:gamma of a European option")},
344 DEF_ARG_SPOT,
345 DEF_ARG_STRIKE,
346 DEF_ARG_TIME_MATURITY_Y,
347 DEF_ARG_RATE_RISKFREE,
348 DEF_ARG_VOLATILITY,
349 DEF_ARG_CC_OPT,
350 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_BS_GAMMA uses the Black-Scholes model to calculate "
351 "the 'gamma' of a European option struck at @{strike} "
352 "on an asset with spot price @{spot}. The gamma of an "
353 "option is the second derivative of its price "
354 "with respect to the price of the underlying asset.")},
355 { GNM_FUNC_HELP_NOTE, F_("Gamma is expressed as the rate of change "
356 "of delta per unit change in @{spot}.")},
357 { GNM_FUNC_HELP_NOTE, F_("Gamma is the same for calls and puts.")},
358 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_VEGA"},
359 { GNM_FUNC_HELP_END}
362 /* theta for the generalized Black and Scholes formula */
363 static gnm_float
364 opt_bs_theta1 (OptionSide side,
365 gnm_float s,gnm_float x,gnm_float t,gnm_float r,gnm_float v,gnm_float b)
367 gnm_float d1 = (gnm_log (s / x) + (b + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
368 gnm_float d2 = d1 - v * gnm_sqrt (t);
370 switch (side) {
371 case OS_Call:
372 return -s * gnm_exp ((b - r) * t) * npdf (d1) * v / (2.0 * gnm_sqrt (t)) -
373 (b - r) * s * gnm_exp ((b - r) * t) * ncdf (d1) - r * x * gnm_exp (-r * t) * ncdf (d2);
374 case OS_Put:
375 return -s * gnm_exp ((b - r) * t) * npdf (d1) * v / (2.0 * gnm_sqrt (t)) +
376 (b - r) * s * gnm_exp ((b - r) * t) * ncdf (-d1) + r * x * gnm_exp (-r * t) * ncdf (-d2);
377 default:
378 return gnm_nan;
382 static GnmValue *
383 opt_bs_theta (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
385 OptionSide call_put = option_side (value_peek_string (argv[0]));
386 gnm_float s = value_get_as_float (argv[1]);
387 gnm_float x = value_get_as_float (argv[2]);
388 gnm_float t = value_get_as_float (argv[3]);
389 gnm_float r = value_get_as_float (argv[4]);
390 gnm_float v = value_get_as_float (argv[5]);
391 gnm_float b = argv[6] ? value_get_as_float (argv[6]) : 0.0;
392 gnm_float gfresult = opt_bs_theta1 (call_put, s, x, t, r, v, b);
394 if (gnm_isnan (gfresult))
395 return value_new_error_NUM (ei->pos);
396 return value_new_float (gfresult);
399 static GnmFuncHelp const help_opt_bs_theta[] = {
400 { GNM_FUNC_HELP_NAME, F_("OPT_BS_THETA:theta of a European option")},
401 DEF_ARG_CALL_PUT_FLAG,
402 DEF_ARG_SPOT,
403 DEF_ARG_STRIKE,
404 DEF_ARG_TIME_MATURITY_Y,
405 DEF_ARG_RATE_RISKFREE,
406 DEF_ARG_VOLATILITY,
407 DEF_ARG_CC_OPT,
408 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_BS_THETA uses the Black-Scholes model to calculate "
409 "the 'theta' of a European option struck at @{strike} "
410 "on an asset with spot price @{spot}. The theta of an "
411 "option is the rate of change of its price with "
412 "respect to time to expiry.")},
413 { GNM_FUNC_HELP_NOTE, F_("Theta is expressed as the negative of the rate of change "
414 "of the option value, per 365.25 days.")},
415 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_VEGA,OPT_BS_GAMMA"},
416 { GNM_FUNC_HELP_END}
420 /* Vega for the generalized Black and Scholes formula */
421 static gnm_float
422 opt_bs_vega1 (gnm_float s, gnm_float x, gnm_float t,
423 gnm_float r, gnm_float v, gnm_float b)
425 gnm_float d1 = (gnm_log (s / x) + (b + (v * v) / 2.0) * t) /
426 (v * gnm_sqrt (t));
427 return s * gnm_exp ((b - r) * t) * npdf (d1) * gnm_sqrt (t);
430 static GnmValue *
431 opt_bs_vega (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
433 gnm_float s = value_get_as_float (argv[0]);
434 gnm_float x = value_get_as_float (argv[1]);
435 gnm_float t = value_get_as_float (argv[2]);
436 gnm_float r = value_get_as_float (argv[3]);
437 gnm_float v = value_get_as_float (argv[4]);
438 gnm_float b = argv[5] ? value_get_as_float (argv[5]) : 0.0;
440 return value_new_float (opt_bs_vega1 (s, x, t, r, v, b));
443 static GnmFuncHelp const help_opt_bs_vega[] = {
444 { GNM_FUNC_HELP_NAME, F_("OPT_BS_VEGA:vega of a European option")},
445 DEF_ARG_SPOT,
446 DEF_ARG_STRIKE,
447 DEF_ARG_TIME_MATURITY_Y,
448 DEF_ARG_RATE_RISKFREE,
449 DEF_ARG_VOLATILITY,
450 DEF_ARG_CC_OPT,
451 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_BS_VEGA uses the Black-Scholes model to calculate "
452 "the 'vega' of a European option struck at @{strike} "
453 "on an asset with spot price @{spot}. The vega of an "
454 "option is the rate of change of its price with "
455 "respect to volatility.")},
456 { GNM_FUNC_HELP_NOTE, F_("Vega is the same for calls and puts.")},
457 /* xgettext:no-c-format */
458 { GNM_FUNC_HELP_NOTE, F_("Vega is expressed as the rate of change "
459 "of option value, per 100% volatility.")},
460 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
461 { GNM_FUNC_HELP_END}
465 /* Rho for the generalized Black and Scholes formula */
466 static gnm_float
467 opt_bs_rho1 (OptionSide side, gnm_float s, gnm_float x,
468 gnm_float t, gnm_float r, gnm_float v, gnm_float b)
470 gnm_float d1 = (gnm_log (s / x) + (b + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
471 gnm_float d2 = d1 - v * gnm_sqrt (t);
472 switch (side) {
473 case OS_Call:
474 if (b != 0)
475 return t * x * gnm_exp (-r * t) * ncdf (d2);
476 else
477 return -t * opt_bs1 (side, s, x, t, r, v, b);
479 case OS_Put:
480 if (b != 0)
481 return -t * x * gnm_exp (-r * t) * ncdf (-d2);
482 else
483 return -t * opt_bs1 (side, s, x, t, r, v, b);
485 default:
486 return gnm_nan;
491 static GnmValue *
492 opt_bs_rho (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
494 OptionSide call_put = option_side (value_peek_string (argv[0]));
495 gnm_float s = value_get_as_float (argv[1]);
496 gnm_float x = value_get_as_float (argv[2]);
497 gnm_float t = value_get_as_float (argv[3]);
498 gnm_float r = value_get_as_float (argv[4]);
499 gnm_float v = value_get_as_float (argv[5]);
500 gnm_float b = argv[6] ? value_get_as_float (argv[6]) : 0.0;
501 gnm_float gfresult = opt_bs_rho1 (call_put, s, x, t, r, v, b);
503 if (gnm_isnan (gfresult))
504 return value_new_error_NUM (ei->pos);
505 return value_new_float (gfresult);
508 static GnmFuncHelp const help_opt_bs_rho[] = {
509 { GNM_FUNC_HELP_NAME, F_("OPT_BS_RHO:rho of a European option")},
510 DEF_ARG_CALL_PUT_FLAG,
511 DEF_ARG_SPOT,
512 DEF_ARG_STRIKE,
513 DEF_ARG_TIME_MATURITY_Y,
514 DEF_ARG_RATE_RISKFREE,
515 DEF_ARG_VOLATILITY,
516 DEF_ARG_CC_OPT,
517 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_BS_RHO uses the Black-Scholes model to calculate "
518 "the 'rho' of a European option struck at @{strike} "
519 "on an asset with spot price @{spot}. The rho of an "
520 "option is the rate of change of its price with "
521 "respect to the risk free interest rate.")},
522 /* xgettext:no-c-format */
523 { GNM_FUNC_HELP_NOTE, F_("Rho is expressed as the rate of change "
524 "of the option value, per 100% change in @{rate}.")},
525 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_THETA,OPT_BS_VEGA,OPT_BS_GAMMA"},
526 { GNM_FUNC_HELP_END}
529 /* Carry for the generalized Black and Scholes formula */
530 static gnm_float
531 opt_bs_carrycost1 (OptionSide side, gnm_float s, gnm_float x, gnm_float t, gnm_float r, gnm_float v, gnm_float b)
533 gnm_float d1 = (gnm_log (s / x) + (b + (v * v) / 2.0) * t) /
534 (v * gnm_sqrt (t));
536 switch (side) {
537 case OS_Call:
538 return t * s * gnm_exp ((b - r) * t) * ncdf (d1);
539 case OS_Put:
540 return -t * s * gnm_exp ((b - r) * t) * ncdf (-d1);
541 default:
542 return gnm_nan; /*should never get to here*/
546 static GnmValue *
547 opt_bs_carrycost (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
549 OptionSide call_put = option_side (value_peek_string (argv[0]));
550 gnm_float s = value_get_as_float (argv[1]);
551 gnm_float x = value_get_as_float (argv[2]);
552 gnm_float t = value_get_as_float (argv[3]);
553 gnm_float r = value_get_as_float (argv[4]);
554 gnm_float v = value_get_as_float (argv[5]);
555 gnm_float b = argv[6] ? value_get_as_float (argv[6]) : 0.0;
556 gnm_float gfresult = opt_bs_carrycost1 (call_put, s, x, t, r, v, b);
558 if (gnm_isnan (gfresult))
559 return value_new_error_NUM (ei->pos);
560 return value_new_float (gfresult);
564 static GnmFuncHelp const help_opt_bs_carrycost[] = {
565 { GNM_FUNC_HELP_NAME, F_("OPT_BS_CARRYCOST:elasticity of a European option")},
566 DEF_ARG_CALL_PUT_FLAG,
567 DEF_ARG_SPOT,
568 DEF_ARG_STRIKE,
569 DEF_ARG_TIME_MATURITY_Y,
570 DEF_ARG_RATE_RISKFREE,
571 DEF_ARG_VOLATILITY,
572 DEF_ARG_CC_OPT,
573 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_BS_CARRYCOST uses the Black-Scholes model to calculate "
574 "the 'elasticity' of a European option struck at @{strike} "
575 "on an asset with spot price @{spot}. The elasticity of an option "
576 "is the rate of change of its price "
577 "with respect to its @{cost_of_carry}.")},
578 /* xgettext:no-c-format */
579 { GNM_FUNC_HELP_NOTE, F_("Elasticity is expressed as the rate of change "
580 "of the option value, per 100% volatility.")},
581 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
582 { GNM_FUNC_HELP_END}
586 /* Currency Options - Garman and Kohlhagen */
587 static gnm_float
588 opt_garman_kohlhagen1 (OptionSide side,
589 gnm_float s, gnm_float x, gnm_float t,
590 gnm_float r, gnm_float rf, gnm_float v)
592 gnm_float d1 = (gnm_log (s / x) + (r - rf + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
593 gnm_float d2 = d1 - v * gnm_sqrt (t);
594 switch (side) {
595 case OS_Call:
596 return s * gnm_exp (-rf * t) * ncdf (d1) - x * gnm_exp (-r * t) * ncdf (d2);
597 case OS_Put:
598 return x * gnm_exp (-r * t) * ncdf (-d2) - s * gnm_exp (-rf * t) * ncdf (-d1);
599 default:
600 return gnm_nan; /*should never get to here*/
604 static GnmValue *
605 opt_garman_kohlhagen (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
607 OptionSide call_put = option_side (value_peek_string (argv[0]));
608 gnm_float s = value_get_as_float (argv[1]);
609 gnm_float x = value_get_as_float (argv[2]);
610 gnm_float t = value_get_as_float (argv[3]);
611 gnm_float r = value_get_as_float (argv[4]);
612 gnm_float rf = value_get_as_float (argv[5]);
613 gnm_float v = value_get_as_float (argv[6]);
614 gnm_float gfresult = opt_garman_kohlhagen1 (call_put, s, x, t, r, rf, v);
616 if (gnm_isnan (gfresult))
617 return value_new_error_NUM (ei->pos);
618 else
619 return value_new_float (gfresult);
622 static GnmFuncHelp const help_opt_garman_kohlhagen[] = {
623 { GNM_FUNC_HELP_NAME, F_("OPT_GARMAN_KOHLHAGEN:theoretical price of a European currency option")},
624 DEF_ARG_CALL_PUT_FLAG,
625 DEF_ARG_SPOT,
626 DEF_ARG_STRIKE,
627 { GNM_FUNC_HELP_ARG, F_("time:number of days to exercise")},
628 { GNM_FUNC_HELP_ARG, F_("domestic_rate:domestic risk-free interest rate to the exercise date in percent")},
629 { GNM_FUNC_HELP_ARG, F_("foreign_rate:foreign risk-free interest rate to the exercise date in percent")},
630 DEF_ARG_VOLATILITY,
631 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_GARMAN_KOHLHAGEN values the theoretical price of a European "
632 "currency option struck at @{strike} on an asset with spot price @{spot}.")},
633 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
634 { GNM_FUNC_HELP_END}
638 /* French (1984) adjusted Black and scholes model for trading day volatility */
639 static gnm_float
640 opt_french1 (OptionSide side, gnm_float s, gnm_float x, gnm_float tradingt, gnm_float calendart,
641 gnm_float r, gnm_float v, gnm_float b)
643 gnm_float d1 = (gnm_log (s / x) + b * calendart + ((v * v) / 2.0) * tradingt) / (v * gnm_sqrt (tradingt));
644 gnm_float d2 = d1 - v * gnm_sqrt (tradingt);
646 switch (side) {
647 case OS_Call:
648 return s * gnm_exp ((b - r) * calendart) * ncdf (d1) - x * gnm_exp (-r * calendart) * ncdf (d2);
649 case OS_Put:
650 return x * gnm_exp (-r * calendart) * ncdf (-d2) - s * gnm_exp ((b - r) * calendart) * ncdf (-d1);
651 default:
652 return gnm_nan;
657 static GnmValue *
658 opt_french (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
660 OptionSide call_put = option_side (value_peek_string (argv[0]));
661 gnm_float s = value_get_as_float (argv[1]);
662 gnm_float x = value_get_as_float (argv[2]);
663 gnm_float t = value_get_as_float (argv[3]);
664 gnm_float t1 = value_get_as_float (argv[4]);
665 gnm_float r = value_get_as_float (argv[5]);
666 gnm_float v = value_get_as_float (argv[6]);
667 gnm_float b = value_get_as_float (argv[7]);
668 gnm_float gfresult = opt_french1 (call_put, s, x, t, t1, r, v, b);
670 if (gnm_isnan (gfresult))
671 return value_new_error_NUM (ei->pos);
672 else
673 return value_new_float (gfresult);
676 static GnmFuncHelp const help_opt_french[] = {
677 { GNM_FUNC_HELP_NAME, F_("OPT_FRENCH:theoretical price of a European option adjusted for trading day volatility")},
678 DEF_ARG_CALL_PUT_FLAG,
679 DEF_ARG_SPOT,
680 DEF_ARG_STRIKE,
681 { GNM_FUNC_HELP_ARG, F_("time:ratio of the number of calendar days to exercise and the number of calendar days in the year")},
682 { GNM_FUNC_HELP_ARG, F_("ttime:ratio of the number of trading days to exercise and the number of trading days in the year")},
683 DEF_ARG_RATE_RISKFREE,
684 DEF_ARG_VOLATILITY,
685 DEF_ARG_CC_OPT,
686 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_FRENCH values the theoretical price of a "
687 "European option adjusted for trading day volatility, struck at "
688 "@{strike} on an asset with spot price @{spot}.")},
689 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
690 { GNM_FUNC_HELP_END}
693 /* Merton jump diffusion model*/
694 static gnm_float
695 opt_jump_diff1 (OptionSide side, gnm_float s, gnm_float x, gnm_float t, gnm_float r, gnm_float v,
696 gnm_float lambda, gnm_float gamma)
698 gnm_float delta, sum;
699 gnm_float Z, vi;
700 int i;
702 delta = gnm_sqrt (gamma * (v * v) / lambda);
703 Z = gnm_sqrt ((v * v) - lambda * (delta * delta));
704 sum = 0.0;
705 for (i = 0; i != 11; ++i) {
706 vi = gnm_sqrt ((Z * Z) + (delta * delta) * (i / t));
707 sum = sum + gnm_exp (-lambda * t) * gnm_pow (lambda * t, i) / gnm_fact(i) *
708 opt_bs1 (side, s, x, t, r, vi, r);
710 return sum;
713 static GnmValue *
714 opt_jump_diff (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
716 OptionSide call_put = option_side (value_peek_string (argv[0]));
717 gnm_float s = value_get_as_float (argv[1]);
718 gnm_float x = value_get_as_float (argv[2]);
719 gnm_float t = value_get_as_float (argv[3]);
720 gnm_float r = value_get_as_float (argv[4]);
721 gnm_float v = value_get_as_float (argv[5]);
722 gnm_float lambda = value_get_as_float (argv[6]);
723 gnm_float gamma = value_get_as_float (argv[7]);
724 gnm_float gfresult =
725 opt_jump_diff1 (call_put, s, x, t, r, v, lambda, gamma);
726 return value_new_float (gfresult);
729 static GnmFuncHelp const help_opt_jump_diff[] = {
730 { GNM_FUNC_HELP_NAME, F_("OPT_JUMP_DIFF:theoretical price of an option according to the Jump Diffusion process")},
731 DEF_ARG_CALL_PUT_FLAG,
732 DEF_ARG_SPOT,
733 DEF_ARG_STRIKE,
734 DEF_ARG_TIME_MATURITY_Y,
735 { GNM_FUNC_HELP_ARG, F_("rate:the annualized rate of interest")},
736 DEF_ARG_VOLATILITY,
737 { GNM_FUNC_HELP_ARG, F_("lambda:expected number of 'jumps' per year")},
738 { GNM_FUNC_HELP_ARG, F_("gamma:proportion of volatility explained by the 'jumps'")},
739 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_JUMP_DIFF models the theoretical price of an option according "
740 "to the Jump Diffusion process (Merton).")},
741 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
742 { GNM_FUNC_HELP_END}
747 /* Miltersen schwartz (1997) commodity option model */
748 static gnm_float
749 opt_miltersen_schwartz1 (OptionSide side, gnm_float p_t, gnm_float f_t, gnm_float x, gnm_float t1,
750 gnm_float t2, gnm_float v_s, gnm_float v_e, gnm_float v_f, gnm_float rho_se,
751 gnm_float rho_sf, gnm_float rho_ef, gnm_float kappa_e, gnm_float kappa_f)
753 gnm_float vz, vxz;
754 gnm_float d1, d2;
756 vz = (v_s * v_s) * t1 + 2.0 * v_s * (v_f * rho_sf * 1.0/ kappa_f * (t1 - 1.0/ kappa_f * gnm_exp (-kappa_f * t2) * (gnm_exp (kappa_f * t1) - 1.0))
757 - v_e * rho_se * 1.0/ kappa_e * (t1 - 1.0/ kappa_e * gnm_exp (-kappa_e * t2) * (gnm_exp (kappa_e * t1) - 1.0)))
758 + (v_e * v_e) * 1.0/ (kappa_e * kappa_e) * (t1 + 1.0/ (2.0 * kappa_e) * gnm_exp (-2 * kappa_e * t2) * (gnm_exp (2.0 * kappa_e * t1) - 1.0)
759 - 2.0 * 1.0/ kappa_e * gnm_exp (-kappa_e * t2) * (gnm_exp (kappa_e * t1) - 1.0))
760 + (v_f * v_f) * 1.0/ (kappa_f * kappa_f) * (t1 + 1.0/ (2.0 * kappa_f) * gnm_exp (-2.0 * kappa_f * t2) * (gnm_exp (2.0 * kappa_f * t1) - 1.0)
761 - 2.0 * 1.0/ kappa_f * gnm_exp (-kappa_f * t2) * (gnm_exp (kappa_f * t1) - 1.0))
762 - 2.0 * v_e * v_f * rho_ef * 1.0/ kappa_e * 1.0/ kappa_f * (t1 - 1.0/ kappa_e * gnm_exp (-kappa_e * t2) * (gnm_exp (kappa_e * t1) - 1.0)
763 - 1.0/ kappa_f * gnm_exp (-kappa_f * t2) * (gnm_exp (kappa_f * t1) - 1.0)
764 + 1.0/ (kappa_e + kappa_f) * gnm_exp (-(kappa_e + kappa_f) * t2) * (gnm_exp ((kappa_e + kappa_f) * t1) - 1.0));
766 vxz = v_f * 1.0/ kappa_f * (v_s * rho_sf * (t1 - 1.0/ kappa_f * (1.0 - gnm_exp (-kappa_f * t1)))
767 + v_f * 1.0/ kappa_f * (t1 - 1.0/ kappa_f * gnm_exp (-kappa_f * t2) * (gnm_exp (kappa_f * t1) - 1.0) - 1.0/ kappa_f * (1 - gnm_exp (-kappa_f * t1))
768 + 1.0/ (2.0 * kappa_f) * gnm_exp (-kappa_f * t2) * (gnm_exp (kappa_f * t1) - gnm_exp (-kappa_f * t1)))
769 - v_e * rho_ef * 1.0/ kappa_e * (t1 - 1.0/ kappa_e * gnm_exp (-kappa_e * t2) * (gnm_exp (kappa_e * t1) - 1.0) - 1.0/ kappa_f * (1.0 - gnm_exp (-kappa_f * t1))
770 + 1.0/ (kappa_e + kappa_f) * gnm_exp (-kappa_e * t2) * (gnm_exp (kappa_e * t1) - gnm_exp (-kappa_f * t1))));
772 vz = gnm_sqrt (vz);
774 d1 = (gnm_log (f_t / x) - vxz + (vz * vz) / 2.0) / vz;
775 d2 = (gnm_log (f_t / x) - vxz - (vz * vz) / 2.0) / vz;
777 switch (side) {
778 case OS_Call:
779 return p_t * (f_t * gnm_exp (-vxz) * ncdf (d1) - x * ncdf (d2));
780 case OS_Put:
781 return p_t * (x * ncdf (-d2) - f_t * gnm_exp (-vxz) * ncdf (-d1));
782 default:
783 return gnm_nan;
787 static GnmValue *
788 opt_miltersen_schwartz (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
790 OptionSide call_put = option_side (value_peek_string (argv[0]));
791 gnm_float p_t = value_get_as_float (argv[1]);
792 gnm_float f_t = value_get_as_float (argv[2]);
793 gnm_float x = value_get_as_float (argv[3]);
794 gnm_float t1 = value_get_as_float (argv[4]);
795 gnm_float t2 = value_get_as_float (argv[5]);
796 gnm_float v_s = value_get_as_float (argv[6]);
797 gnm_float v_e = value_get_as_float (argv[7]);
798 gnm_float v_f = value_get_as_float (argv[8]);
799 gnm_float rho_se = value_get_as_float (argv[9]);
800 gnm_float rho_sf = value_get_as_float (argv[10]);
801 gnm_float rho_ef = value_get_as_float (argv[11]);
802 gnm_float kappa_e = value_get_as_float (argv[12]);
803 gnm_float kappa_f = value_get_as_float (argv[13]);
805 gnm_float gfresult =
806 opt_miltersen_schwartz1 (call_put, p_t, f_t, x, t1, t2,
807 v_s, v_e, v_f,
808 rho_se, rho_sf, rho_ef, kappa_e, kappa_f);
810 if (gnm_isnan (gfresult))
811 return value_new_error_NUM (ei->pos);
812 return value_new_float (gfresult);
816 static GnmFuncHelp const help_opt_miltersen_schwartz[] = {
817 { GNM_FUNC_HELP_NAME, F_("OPT_MILTERSEN_SCHWARTZ:theoretical price of options on commodities futures according to Miltersen & Schwartz")},
818 DEF_ARG_CALL_PUT_FLAG,
819 { GNM_FUNC_HELP_ARG, F_("p_t:zero coupon bond with expiry at option maturity")},
820 { GNM_FUNC_HELP_ARG, F_("f_t:futures price")},
821 DEF_ARG_STRIKE,
822 { GNM_FUNC_HELP_ARG, F_("t1:time to maturity of the option")},
823 { GNM_FUNC_HELP_ARG, F_("t2:time to maturity of the underlying commodity futures contract")},
824 { GNM_FUNC_HELP_ARG, F_("v_s:volatility of the spot commodity price")},
825 { GNM_FUNC_HELP_ARG, F_("v_e:volatility of the future convenience yield")},
826 { GNM_FUNC_HELP_ARG, F_("v_f:volatility of the forward rate of interest")},
827 { GNM_FUNC_HELP_ARG, F_("rho_se:correlation between the spot commodity price and the convenience yield")},
828 { GNM_FUNC_HELP_ARG, F_("rho_sf:correlation between the spot commodity price and the forward interest rate")},
829 { GNM_FUNC_HELP_ARG, F_("rho_ef:correlation between the forward interest rate and the convenience yield")},
830 { GNM_FUNC_HELP_ARG, F_("kappa_e:speed of mean reversion of the convenience yield")},
831 { GNM_FUNC_HELP_ARG, F_("kappa_f:speed of mean reversion of the forward interest rate")},
832 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
833 { GNM_FUNC_HELP_END}
836 /* American options */
839 /* American Calls on stocks with known dividends, Roll-Geske-Whaley */
840 static gnm_float opt_rgw1 (gnm_float s, gnm_float x, gnm_float t1, gnm_float t2, gnm_float r, gnm_float d, gnm_float v)
841 /*t1 time to dividend payout
842 t2 time to option expiration */
844 gnm_float sx, i;
845 gnm_float a1, a2, b1, b2;
846 gnm_float HighS, LowS, epsilon;
847 gnm_float ci, infinity;
848 gnm_float gfresult;
850 if (!(s > 0))
851 return gnm_nan;
853 infinity = 100000000;
854 epsilon = 0.00001;
855 sx = s - d * gnm_exp (-r * t1);
856 if (d <= (x * (1.0 - gnm_exp (-r * (t2 - t1)))))
857 /* Not optimal to exercise */
858 return opt_bs1 (OS_Call, sx, x, t2, r, v,0.0);
860 ci = opt_bs1 (OS_Call, s, x, t2 - t1, r, v,0.0);
861 HighS = s;
862 while ((ci - HighS - d + x) > 0.0 && HighS < infinity) {
864 HighS *= 2.0;
865 ci = opt_bs1 (OS_Call, HighS, x, t2 - t1, r, v,0.0);
867 if (HighS > infinity)
868 return opt_bs1 (OS_Call, sx, x, t2, r, v,0.0);
870 LowS = 0.0;
871 i = HighS * 0.5;
872 ci = opt_bs1 (OS_Call, i, x, t2 - t1, r, v, 0.0);
874 /* search algorithm to find the critical stock price i */
875 while (gnm_abs (ci - i - d + x) > epsilon && HighS - LowS > epsilon) {
876 if ((ci - i - d + x) < 0)
877 HighS = i;
878 else
879 LowS = i;
880 i = (HighS + LowS) / 2.0;
881 ci = opt_bs1 (OS_Call, i, x, (t2 - t1), r, v, 0.0);
884 a1 = (gnm_log (sx / x) + (r + (v * v) / 2.0) * t2) / (v * gnm_sqrt (t2));
885 a2 = a1 - v * gnm_sqrt (t2);
886 b1 = (gnm_log (sx / i) + (r + (v * v) / 2.0) * t1) / (v * gnm_sqrt (t1));
887 b2 = b1 - v * gnm_sqrt (t1);
889 gfresult = sx * ncdf (b1) + sx * cum_biv_norm_dist1 (a1, -b1, -gnm_sqrt (t1 / t2))
890 - x * gnm_exp (-r * t2) * cum_biv_norm_dist1 (a2, -b2, -gnm_sqrt (t1 / t2)) - (x - d)
891 * gnm_exp (-r * t1) * ncdf (b2);
892 return gfresult;
896 static GnmValue *
897 opt_rgw(GnmFuncEvalInfo *ei, GnmValue const * const *argv)
899 gnm_float s = value_get_as_float (argv[0]);
900 gnm_float x = value_get_as_float (argv[1]);
901 gnm_float t1 = value_get_as_float (argv[2]);
902 gnm_float t2 = value_get_as_float (argv[3]);
903 gnm_float r = value_get_as_float (argv[4]);
904 gnm_float d = value_get_as_float (argv[5]);
905 gnm_float v = value_get_as_float (argv[6]);
906 gnm_float gfresult = 0.0;
908 gfresult = opt_rgw1 (s, x, t1, t2, r, d, v);
910 return value_new_float (gfresult);
913 static GnmFuncHelp const help_opt_rgw[] = {
914 { GNM_FUNC_HELP_NAME, F_("OPT_RGW:theoretical price of an American option according to the Roll-Geske-Whaley approximation")},
915 DEF_ARG_SPOT,
916 DEF_ARG_STRIKE,
917 DEF_ARG_TIME_DIVIDEND,
918 DEF_ARG_TIME_EXPIRATION,
919 DEF_ARG_RATE_ANNUALIZED,
920 DEF_ARG_AMOUNT,
921 DEF_ARG_VOLATILITY,
922 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
923 { GNM_FUNC_HELP_END}
926 /* the Barone-Adesi and Whaley (1987) American approximation */
927 static GnmValue *
928 opt_baw_amer (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
930 OptionSide call_put = option_side (value_peek_string (argv[0]));
931 gnm_float s = value_get_as_float (argv[1]);
932 gnm_float x = value_get_as_float (argv[2]);
933 gnm_float t = value_get_as_float (argv[3]);
934 gnm_float r = value_get_as_float (argv[4]);
935 gnm_float v = value_get_as_float (argv[5]);
936 gnm_float b = value_get_as_float (argv[6]);
937 gnm_float gfresult;
939 switch (call_put) {
940 case OS_Call:
941 gfresult = opt_baw_call (s, x, t, r, v, b);
942 break;
943 case OS_Put:
944 gfresult = opt_baw_put (s, x, t, r, v, b);
945 break;
946 default:
947 return value_new_error_NUM (ei->pos);
950 if (gnm_isnan (gfresult))
951 return value_new_error_NUM (ei->pos);
953 return value_new_float (gfresult);
956 static GnmFuncHelp const help_opt_baw_amer[] = {
957 { GNM_FUNC_HELP_NAME, F_("OPT_BAW_AMER:theoretical price of an option according to the Barone Adesie & Whaley approximation")},
958 DEF_ARG_CALL_PUT_FLAG,
959 DEF_ARG_SPOT,
960 DEF_ARG_STRIKE,
961 DEF_ARG_TIME_MATURITY_D,
962 DEF_ARG_RATE_RISKFREE_ANN,
963 DEF_ARG_CC,
964 DEF_ARG_VOLATILITY_SHORT,
965 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
966 { GNM_FUNC_HELP_END}
969 /* American call */
970 static gnm_float
971 opt_baw_call (gnm_float s, gnm_float x, gnm_float t, gnm_float r, gnm_float v, gnm_float b)
973 gnm_float sk, n, k;
974 gnm_float d1, q2, a2;
975 gnm_float gfresult;
976 if (b >= r)
977 gfresult = opt_bs1 (OS_Call, s, x, t, r, v, b);
978 else
980 sk = NRA_c (x, t, r, v, b);
981 n = 2 * b / (v * v);
982 k = 2 * r / ((v * v) * (1.0 - gnm_exp (-r * t)));
983 d1 = (gnm_log (sk / x) + (b + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
984 q2 = (-(n - 1.0) + gnm_sqrt ((n - 1.0) * (n - 1.0) + 4.0 * k)) / 2.0;
985 a2 = (sk / q2) * (1.0 - gnm_exp ((b - r) * t) * ncdf (d1));
986 if (s < sk)
987 gfresult = opt_bs1 (OS_Call, s, x, t, r, v, b) + a2 * gnm_pow (s / sk, q2);
988 else
989 gfresult = s - x;
991 } /*end if statement*/
992 return gfresult;
999 /* Newton Raphson algorithm to solve for the critical commodity price for a Call */
1000 static gnm_float
1001 NRA_c (gnm_float x, gnm_float t, gnm_float r, gnm_float v, gnm_float b)
1003 gnm_float n, m;
1004 gnm_float su, si;
1005 gnm_float h2, k;
1006 gnm_float d1, q2, q2u;
1007 gnm_float LHS, RHS;
1008 gnm_float bi, e;
1010 /* Calculation of seed value, si */
1011 n = 2 * b / (v * v);
1012 m = 2 * r / (v * v);
1013 q2u = (-(n - 1.0) + gnm_sqrt (((n - 1.0) * (n - 1.0)) + 4.0 * m)) / 2.0;
1014 su = x / (1.0 - 1.0/ q2u);
1015 h2 = -(b * t + 2.0 * v * gnm_sqrt (t)) * x / (su - x);
1016 si = x + (su - x) * (1.0 - gnm_exp (h2));
1018 k = 2 * r / ((v * v) * (1.0 - gnm_exp (-r * t)));
1019 d1 = (gnm_log (si / x) + (b + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
1020 q2 = (-(n - 1.0) + gnm_sqrt (((n - 1.0) * (n - 1.0)) + 4.0 * k)) / 2.0;
1021 LHS = si - x;
1022 RHS = opt_bs1 (OS_Call, si, x, t, r, v, b) + (1.0 - gnm_exp ((b - r) * t) * ncdf (d1)) * si / q2;
1023 bi = gnm_exp ((b - r) * t) * ncdf (d1) * (1.0 - 1.0/ q2)
1024 + (1.0 - gnm_exp ((b - r) * t) * ncdf (d1) / (v * gnm_sqrt (t))) / q2;
1025 e = 0.000001;
1027 /* Newton Raphson algorithm for finding critical price si */
1028 while ((gnm_abs (LHS - RHS) / x) > e)
1030 si = (x + RHS - bi * si) / (1.0 - bi);
1031 d1 = (gnm_log (si / x) + (b + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
1032 LHS = si - x;
1033 RHS = opt_bs1 (OS_Call, si, x, t, r, v, b) + (1.0 - gnm_exp ((b - r) * t) * ncdf (d1)) * si / q2;
1034 bi = gnm_exp ((b - r) * t) * ncdf (d1) * (1.0 - 1.0/ q2)
1035 + (1.0 - gnm_exp ((b - r) * t) * npdf (d1) / (v * gnm_sqrt (t))) / q2;
1037 return si;
1040 static gnm_float
1041 opt_baw_put (gnm_float s, gnm_float x, gnm_float t, gnm_float r, gnm_float v, gnm_float b)
1043 gnm_float sk = NRA_p (x, t, r, v, b);
1044 gnm_float n = 2 * b / (v * v);
1045 gnm_float k = 2 * r / ((v * v) * (1.0 - gnm_exp (-r * t)));
1046 gnm_float d1 = (gnm_log (sk / x) + (b + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
1047 gnm_float q1 = (-(n - 1.0) - gnm_sqrt (((n - 1.0) * (n - 1.0)) + 4.0 * k)) / 2.0;
1048 gnm_float a1 = -(sk / q1) * (1.0 - gnm_exp ((b - r) * t) * ncdf (-d1));
1050 if (s > sk)
1051 return opt_bs1 (OS_Put, s, x, t, r, v, b) + a1 * gnm_pow (s/ sk, q1);
1052 else
1053 return x - s;
1056 /* Newton Raphson algorithm to solve for the critical commodity price for a Put*/
1057 static gnm_float
1058 NRA_p (gnm_float x, gnm_float t, gnm_float r, gnm_float v, gnm_float b)
1061 gnm_float n, m;
1062 gnm_float su, si;
1063 gnm_float h1, k;
1064 gnm_float d1, q1u, q1;
1065 gnm_float LHS, RHS;
1066 gnm_float bi, e;
1068 /* Calculation of seed value, si */
1069 n = 2 * b / (v * v);
1070 m = 2 * r / (v * v);
1071 q1u = (-(n - 1.0) - gnm_sqrt (((n - 1.0) * (n - 1.0)) + 4.0 * m)) / 2.0;
1072 su = x / (1.0 - 1.0/ q1u);
1073 h1 = (b * t - 2.0 * v * gnm_sqrt (t)) * x / (x - su);
1074 si = su + (x - su) * gnm_exp (h1);
1076 k = 2 * r / ((v * v) * (1.0 - gnm_exp (-r * t)));
1077 d1 = (gnm_log (si / x) + (b + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
1078 q1 = (-(n - 1.0) - gnm_sqrt (((n - 1.0) * (n - 1.0)) + 4.0 * k)) / 2.0;
1079 LHS = x - si;
1080 RHS = opt_bs1 (OS_Put, si, x, t, r, v, b) - (1.0 - gnm_exp ((b - r) * t) * ncdf (-d1)) * si / q1;
1081 bi = -gnm_exp ((b - r) * t) * ncdf (-d1) * (1.0 - 1.0/ q1)
1082 - (1.0 + gnm_exp ((b - r) * t) * npdf (-d1) / (v * gnm_sqrt (t))) / q1;
1083 e = 0.000001;
1085 /* Newton Raphson algorithm for finding critical price si */
1086 while(gnm_abs (LHS - RHS) / x > e) {
1087 si = (x - RHS + bi * si) / (1.0 + bi);
1088 d1 = (gnm_log (si / x) + (b + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
1089 LHS = x - si;
1090 RHS = opt_bs1 (OS_Put, si, x, t, r, v, b) - (1.0 - gnm_exp ((b - r) * t) * ncdf (-d1)) * si / q1;
1091 bi = -gnm_exp ((b - r) * t) * ncdf (-d1) * (1.0 - 1.0/ q1)
1092 - (1.0 + gnm_exp ((b - r) * t) * ncdf (-d1) / (v * gnm_sqrt (t))) / q1;
1094 return si;
1097 /* the Bjerksund and stensland (1993) American approximation */
1098 static gnm_float
1099 opt_bjer_stens1 (OptionSide side, gnm_float s, gnm_float x, gnm_float t, gnm_float r, gnm_float v, gnm_float b)
1101 switch (side) {
1102 case OS_Call:
1103 return opt_bjer_stens1_c (s, x, t, r, v, b);
1104 case OS_Put:
1105 /* Use the Bjerksund and stensland put-call transformation */
1106 return opt_bjer_stens1_c (x, s, t, r - b, v, -b);
1107 default:
1108 return gnm_nan;
1112 static GnmValue *
1113 opt_bjer_stens (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1115 OptionSide call_put = option_side (value_peek_string (argv[0]));
1116 gnm_float s = value_get_as_float (argv[1]);
1117 gnm_float x = value_get_as_float (argv[2]);
1118 gnm_float t = value_get_as_float (argv[3]);
1119 gnm_float r = value_get_as_float (argv[4]);
1120 gnm_float v = value_get_as_float (argv[5]);
1121 gnm_float b = argv[6] ? value_get_as_float (argv[6]):0;
1122 gnm_float gfresult =
1123 opt_bjer_stens1 (call_put, s, x, t, r, v, b);
1124 return value_new_float (gfresult);
1128 static GnmFuncHelp const help_opt_bjer_stens[] = {
1129 { GNM_FUNC_HELP_NAME, F_("OPT_BJER_STENS:theoretical price of American options according to the Bjerksund & Stensland approximation technique")},
1130 DEF_ARG_CALL_PUT_FLAG,
1131 DEF_ARG_SPOT,
1132 DEF_ARG_STRIKE,
1133 DEF_ARG_TIME_MATURITY_D,
1134 DEF_ARG_RATE_RISKFREE_ANN,
1135 DEF_ARG_VOLATILITY_SHORT,
1136 DEF_ARG_CC_OPT,
1137 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1138 { GNM_FUNC_HELP_END}
1141 static gnm_float
1142 opt_bjer_stens1_c (gnm_float s, gnm_float x, gnm_float t, gnm_float r, gnm_float v, gnm_float b)
1144 if (b >= r) /* Never optimal to exersice before maturity */
1145 return opt_bs1 (OS_Call, s, x, t, r, v, b);
1146 else {
1147 gnm_float Beta =
1148 (1.0/ 2.0 - b / (v * v)) +
1149 gnm_sqrt (gnm_pow (b / (v * v) - 1.0/ 2.0, 2) + 2 * r / (v * v));
1150 gnm_float BInfinity = Beta / (Beta - 1.0) * x;
1151 gnm_float B0 = MAX (x, r / (r - b) * x);
1152 gnm_float ht = -(b * t + 2.0 * v * gnm_sqrt (t)) * B0 / (BInfinity - B0);
1153 gnm_float I = B0 + (BInfinity - B0) * (1.0 - gnm_exp (ht));
1154 if (s >= I)
1155 return s - x;
1156 else {
1157 gnm_float alpha = (I - x) * gnm_pow (I ,-Beta);
1158 return alpha * gnm_pow (s ,Beta) -
1159 alpha * phi (s, t, Beta, I, I, r, v, b) +
1160 phi (s, t, 1.0, I, I, r, v, b) -
1161 phi (s, t, 1.0, x, I, r, v, b) -
1162 x * phi (s, t, 0.0, I, I, r, v, b) +
1163 x * phi (s, t, 0.0, x, I, r, v, b);
1168 static gnm_float
1169 phi (gnm_float s, gnm_float t, gnm_float gamma, gnm_float H, gnm_float I, gnm_float r, gnm_float v, gnm_float b)
1171 gnm_float lambda, kappa;
1172 gnm_float d;
1173 gnm_float gfresult;
1175 lambda = (-r + gamma * b + 0.5 * gamma * (gamma - 1.0) * (v * v)) * t;
1176 d = -(gnm_log (s / H) + (b + (gamma - 0.5) * (v * v)) * t) / (v * gnm_sqrt (t));
1177 kappa = 2 * b / (v * v) + (2.0 * gamma - 1.0);
1178 gfresult = gnm_exp (lambda) * gnm_pow (s, gamma) * (ncdf (d) - gnm_pow (I / s, kappa) * ncdf (d - 2.0 * gnm_log (I / s) / (v * gnm_sqrt (t))));
1180 return gfresult;
1184 /* Executive stock options */
1185 static GnmValue *
1186 opt_exec (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1188 OptionSide call_put = option_side (value_peek_string (argv[0]));
1189 gnm_float s = value_get_as_float (argv[1]);
1190 gnm_float x = value_get_as_float (argv[2]);
1191 gnm_float t = value_get_as_float (argv[3]);
1192 gnm_float r = value_get_as_float (argv[4]);
1193 gnm_float v = value_get_as_float (argv[5]);
1194 gnm_float b = value_get_as_float (argv[6]);
1195 gnm_float lambda = value_get_as_float (argv[7]);
1196 gnm_float gfresult =
1197 gnm_exp (-lambda * t) * opt_bs1 (call_put, s, x, t, r, v, b);
1198 return value_new_float (gfresult);
1202 static GnmFuncHelp const help_opt_exec[] = {
1203 { GNM_FUNC_HELP_NAME, F_("OPT_EXEC:theoretical price of executive stock options")},
1204 DEF_ARG_CALL_PUT_FLAG,
1205 DEF_ARG_SPOT,
1206 DEF_ARG_STRIKE,
1207 DEF_ARG_TIME_MATURITY_D,
1208 DEF_ARG_RATE_RISKFREE_ANN,
1209 DEF_ARG_VOLATILITY_SHORT,
1210 DEF_ARG_CC,
1211 { GNM_FUNC_HELP_ARG, F_("lambda:jump rate for executives")},
1212 { GNM_FUNC_HELP_NOTE, F_("The model assumes executives forfeit their options if they leave the company.")},
1213 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1214 { GNM_FUNC_HELP_END}
1221 /* Forward start options */
1222 static GnmValue *
1223 opt_forward_start(GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1225 OptionSide call_put = option_side (value_peek_string (argv[0]));
1226 gnm_float s = value_get_as_float (argv[1]);
1227 gnm_float alpha = value_get_as_float (argv[2]);
1228 gnm_float t1 = value_get_as_float (argv[3]);
1229 gnm_float t = value_get_as_float (argv[4]);
1230 gnm_float r = value_get_as_float (argv[5]);
1231 gnm_float v = value_get_as_float (argv[6]);
1232 gnm_float b = value_get_as_float (argv[7]);
1233 gnm_float gfresult =
1234 s * gnm_exp ((b - r) * t1) * opt_bs1 (call_put, 1, alpha, t - t1, r, v, b);
1235 return value_new_float (gfresult);
1240 static GnmFuncHelp const help_opt_forward_start[] = {
1241 { GNM_FUNC_HELP_NAME, F_("OPT_FORWARD_START:theoretical price of forward start options")},
1242 DEF_ARG_CALL_PUT_FLAG,
1243 DEF_ARG_SPOT,
1244 { GNM_FUNC_HELP_ARG, F_("alpha:fraction setting the strike price at the future date @{time_start}")},
1245 { GNM_FUNC_HELP_ARG, F_("time_start:time until the option starts in days")},
1246 DEF_ARG_TIME_MATURITY_D,
1247 DEF_ARG_RATE_RISKFREE_ANN,
1248 DEF_ARG_VOLATILITY_SHORT,
1249 DEF_ARG_CC,
1250 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1251 { GNM_FUNC_HELP_END}
1255 /* time switch options (discrete) */
1256 static GnmValue *
1257 opt_time_switch (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1259 OptionSide call_put = option_side (value_peek_string (argv[0]));
1260 gnm_float s = value_get_as_float (argv[1]);
1261 gnm_float x = value_get_as_float (argv[2]);
1262 gnm_float a = value_get_as_float (argv[3]);
1263 gnm_float t = value_get_as_float (argv[4]);
1264 gnm_float m = value_get_as_float (argv[5]);
1265 gnm_float dt = value_get_as_float (argv[6]);
1266 gnm_float r = value_get_as_float (argv[7]);
1267 gnm_float b = value_get_as_float (argv[8]);
1268 gnm_float v = value_get_as_float (argv[9]);
1270 gnm_float gfresult;
1271 gnm_float sum, d;
1272 int i, n, Z = 0;
1274 switch (call_put) {
1275 case OS_Call: Z = +1; break;
1276 case OS_Put: Z = -1; break;
1277 default: return value_new_error_NUM (ei->pos);
1280 sum = 0.0;
1281 n = t / dt;
1282 for (i = 1; i < n; ++i) {
1283 d = (gnm_log (s / x) + (b - (v * v) / 2.0) * i * dt) / (v * gnm_sqrt (i * dt));
1284 sum = sum + ncdf (Z * d) * dt;
1287 gfresult = a * gnm_exp (-r * t) * sum + dt * a * gnm_exp (-r * t) * m;
1288 return value_new_float (gfresult);
1292 static GnmFuncHelp const help_opt_time_switch[] = {
1293 { GNM_FUNC_HELP_NAME, F_("OPT_TIME_SWITCH:theoretical price of time switch options")},
1294 DEF_ARG_CALL_PUT_FLAG,
1295 DEF_ARG_SPOT,
1296 DEF_ARG_STRIKE,
1297 { GNM_FUNC_HELP_ARG, F_("a:amount received for each time period")},
1298 DEF_ARG_TIME_MATURITY_Y,
1299 { GNM_FUNC_HELP_ARG, F_("m:number of time units the option has already met the condition")},
1300 { GNM_FUNC_HELP_ARG, F_("dt:agreed upon discrete time period expressed as "
1301 "a fraction of a year")},
1302 DEF_ARG_RATE_RISKFREE_ANN,
1303 DEF_ARG_CC,
1304 DEF_ARG_VOLATILITY_SHORT,
1305 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_TIME_SWITCH models the theoretical price of time switch options. (Pechtl 1995). "
1306 "The holder receives @{a} * @{dt} for each period that the asset price was "
1307 "greater than @{strike} (for a call) or below it (for a put).")},
1308 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1309 { GNM_FUNC_HELP_END}
1313 /* simple chooser options */
1314 static GnmValue *
1315 opt_simple_chooser (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1317 gnm_float s = value_get_as_float (argv[0]);
1318 gnm_float x = value_get_as_float (argv[1]);
1319 gnm_float t1 = value_get_as_float (argv[2]);
1320 gnm_float t2 = value_get_as_float (argv[3]);
1321 gnm_float r = value_get_as_float (argv[4]);
1322 gnm_float b = value_get_as_float (argv[5]);
1323 gnm_float v = value_get_as_float (argv[6]);
1325 gnm_float d = (gnm_log (s / x) + (b + (v * v) / 2.0) * t2) / (v * gnm_sqrt (t2));
1326 gnm_float y = (gnm_log (s / x) + b * t2 + (v * v) * t1 / 2.0) / (v * gnm_sqrt (t1));
1327 gnm_float gfresult =
1328 s * gnm_exp ((b - r) * t2) * ncdf ( d) - x * gnm_exp (-r * t2) * ncdf ( d - v * gnm_sqrt (t2)) -
1329 s * gnm_exp ((b - r) * t2) * ncdf (-y) + x * gnm_exp (-r * t2) * ncdf (-y + v * gnm_sqrt (t1));
1331 return value_new_float (gfresult);
1334 static GnmFuncHelp const help_opt_simple_chooser[] = {
1335 { GNM_FUNC_HELP_NAME, F_("OPT_SIMPLE_CHOOSER:theoretical price of a simple chooser option")},
1336 DEF_ARG_CALL_PUT_FLAG,
1337 DEF_ARG_SPOT,
1338 DEF_ARG_STRIKE,
1339 { GNM_FUNC_HELP_ARG, F_("time1:time in years until the holder chooses a put or a call option")},
1340 { GNM_FUNC_HELP_ARG, F_("time2:time in years until the chosen option expires")},
1341 DEF_ARG_CC,
1342 DEF_ARG_VOLATILITY_SHORT,
1343 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1344 { GNM_FUNC_HELP_END}
1348 /* Complex chooser options */
1349 static GnmValue *
1350 opt_complex_chooser(GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1352 gnm_float s = value_get_as_float (argv[0]);
1353 gnm_float xc = value_get_as_float (argv[1]);
1354 gnm_float xp = value_get_as_float (argv[2]);
1355 gnm_float t = value_get_as_float (argv[3]);
1356 gnm_float tc = value_get_as_float (argv[4]);
1357 gnm_float tp = value_get_as_float (argv[5]);
1358 gnm_float r = value_get_as_float (argv[6]);
1359 gnm_float b = value_get_as_float (argv[7]);
1360 gnm_float v = value_get_as_float (argv[8]);
1362 gnm_float gfresult;
1364 gnm_float d1, d2, y1, y2;
1365 gnm_float rho1, rho2, I;
1367 I = opt_crit_val_chooser (s, xc, xp, t, tc, tp, r, b, v);
1368 d1 = (gnm_log (s / I) + (b + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
1369 d2 = d1 - v * gnm_sqrt (t);
1370 y1 = (gnm_log (s / xc) + (b + (v * v) / 2.0) * tc) / (v * gnm_sqrt (tc));
1371 y2 = (gnm_log (s / xp) + (b + (v * v) / 2.0) * tp) / (v * gnm_sqrt (tp));
1372 rho1 = gnm_sqrt (t / tc);
1373 rho2 = gnm_sqrt (t / tp);
1375 gfresult = s * gnm_exp ((b - r) * tc) * cum_biv_norm_dist1 (d1, y1, rho1) - xc * gnm_exp (-r * tc)
1376 * cum_biv_norm_dist1 (d2, y1 - v * gnm_sqrt (tc), rho1) - s * gnm_exp ((b - r) * tp)
1377 * cum_biv_norm_dist1 (-d1, -y2, rho2) + xp * gnm_exp (-r * tp) * cum_biv_norm_dist1 (-d2, -y2 + v * gnm_sqrt (tp), rho2);
1379 return value_new_float (gfresult);
1383 static GnmFuncHelp const help_opt_complex_chooser[] = {
1384 { GNM_FUNC_HELP_NAME, F_("OPT_COMPLEX_CHOOSER:theoretical price of a complex chooser option")},
1385 DEF_ARG_SPOT,
1386 { GNM_FUNC_HELP_ARG, F_("strike_call:strike price, if exercised as a call option")},
1387 { GNM_FUNC_HELP_ARG, F_("strike_put:strike price, if exercised as a put option")},
1388 { GNM_FUNC_HELP_ARG, F_("time:time in years until the holder chooses a put or a call option")},
1389 { GNM_FUNC_HELP_ARG, F_("time_call:time in years to maturity of the call option if chosen")},
1390 { GNM_FUNC_HELP_ARG, F_("time_put:time in years to maturity of the put option if chosen")},
1391 DEF_ARG_RATE_RISKFREE_ANN,
1392 DEF_ARG_CC,
1393 DEF_ARG_VOLATILITY,
1394 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1395 { GNM_FUNC_HELP_END}
1401 /* Critical value complex chooser option */
1402 static gnm_float
1403 opt_crit_val_chooser (gnm_float s,gnm_float xc,gnm_float xp,gnm_float t,
1404 gnm_float tc, gnm_float tp, gnm_float r, gnm_float b, gnm_float v)
1406 gnm_float sv, ci, Pi, epsilon;
1407 gnm_float dc, dp, yi, di;
1409 sv = s;
1410 ci = opt_bs1 (OS_Call, sv, xc, tc - t, r, v, b);
1411 Pi = opt_bs1 (OS_Put, sv, xp, tp - t, r, v, b);
1412 dc = opt_bs_delta1 (OS_Call, sv, xc, tc - t, r, v, b);
1413 dp = opt_bs_delta1 (OS_Put, sv, xp, tp - t, r, v, b);
1414 yi = ci - Pi;
1415 di = dc - dp;
1416 epsilon = 0.001;
1417 /* Newton-Raphson */
1418 while (gnm_abs (yi) > epsilon)
1420 sv = sv - (yi) / di;
1421 ci = opt_bs1 (OS_Call, sv, xc, tc - t, r, v, b);
1422 Pi = opt_bs1 (OS_Put, sv, xp, tp - t, r, v, b);
1423 dc = opt_bs_delta1 (OS_Call, sv, xc, tc - t, r, v, b);
1424 dp = opt_bs_delta1 (OS_Put, sv, xp, tp - t, r, v, b);
1425 yi = ci - Pi;
1426 di = dc - dp;
1429 return sv;
1433 /* Options on options */
1434 static GnmValue *
1435 opt_on_options (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1437 char const *type_flag = value_peek_string (argv[0]);
1438 gnm_float s = value_get_as_float (argv[1]);
1439 gnm_float x1 = value_get_as_float (argv[2]);
1440 gnm_float x2 = value_get_as_float (argv[3]);
1441 gnm_float t1 = value_get_as_float (argv[4]);
1442 gnm_float t2 = value_get_as_float (argv[5]);
1443 gnm_float r = value_get_as_float (argv[6]);
1444 gnm_float b = value_get_as_float (argv[7]);
1445 gnm_float v = value_get_as_float (argv[8]);
1447 gnm_float gfresult;
1449 gnm_float y1, y2, z1, z2;
1450 gnm_float I, rho;
1451 OptionSide call_put;
1453 if (!strcmp (type_flag, "cc") || !strcmp (type_flag, "pc"))
1454 call_put = OS_Call;
1455 else
1456 call_put = OS_Put;
1458 I = CriticalValueOptionsOnOptions (call_put, x1, x2, t2 - t1, r, b, v);
1460 rho = gnm_sqrt (t1 / t2);
1461 y1 = (gnm_log (s / I) + (b + (v * v) / 2.0) * t1) / (v * gnm_sqrt (t1));
1462 y2 = y1 - v * gnm_sqrt (t1);
1463 z1 = (gnm_log (s / x1) + (b + (v * v) / 2.0) * t2) / (v * gnm_sqrt (t2));
1464 z2 = z1 - v * gnm_sqrt (t2);
1466 if (!strcmp (type_flag, "cc"))
1467 gfresult = s * gnm_exp ((b - r) * t2) * cum_biv_norm_dist1 (z1, y1, rho) -
1468 x1 * gnm_exp (-r * t2) * cum_biv_norm_dist1 (z2, y2, rho) - x2 * gnm_exp (-r * t1) * ncdf (y2);
1469 else if (!strcmp (type_flag, "pc"))
1470 gfresult = x1 * gnm_exp (-r * t2) * cum_biv_norm_dist1 (z2, -y2, -rho) -
1471 s * gnm_exp ((b - r) * t2) * cum_biv_norm_dist1 (z1, -y1, -rho) + x2 * gnm_exp (-r * t1) * ncdf (-y2);
1472 else if (!strcmp (type_flag, "cp"))
1473 gfresult = x1 * gnm_exp (-r * t2) * cum_biv_norm_dist1 (-z2, -y2, rho) -
1474 s * gnm_exp ((b - r) * t2) * cum_biv_norm_dist1 (-z1, -y1, rho) - x2 * gnm_exp (-r * t1) * ncdf (-y2);
1475 else if (!strcmp (type_flag, "pp"))
1476 gfresult = s * gnm_exp ((b - r) * t2) * cum_biv_norm_dist1 (-z1, y1, -rho) -
1477 x1 * gnm_exp (-r * t2) * cum_biv_norm_dist1 (-z2, y2, -rho) + gnm_exp (-r * t1) * x2 * ncdf (y2);
1478 else
1479 return value_new_error_VALUE (ei->pos);
1481 return value_new_float (gfresult);
1484 static GnmFuncHelp const help_opt_on_options[] = {
1485 { GNM_FUNC_HELP_NAME, F_("OPT_ON_OPTIONS:theoretical price of options on options")},
1486 { GNM_FUNC_HELP_ARG, F_("type_flag:'cc' for calls on calls, 'cp' for calls on puts, and so on for 'pc', and 'pp'")},
1487 DEF_ARG_SPOT,
1488 { GNM_FUNC_HELP_ARG, F_("strike1:strike price at which the option being valued is struck")},
1489 { GNM_FUNC_HELP_ARG, F_("strike2:strike price at which the underlying option is struck")},
1490 { GNM_FUNC_HELP_ARG, F_("time1:time in years to maturity of the option")},
1491 { GNM_FUNC_HELP_ARG, F_("time2:time in years to the maturity of the underlying option")},
1492 DEF_ARG_RATE_RISKFREE_ANN,
1493 { GNM_FUNC_HELP_ARG, F_("cost_of_carry:net cost of holding the underlying asset of the underlying option")},
1494 { GNM_FUNC_HELP_ARG, F_("volatility:annualized volatility in price of the underlying asset of the underlying option")},
1495 { GNM_FUNC_HELP_NOTE, F_("For common stocks, @{cost_of_carry} is the risk free rate less the dividend yield.")},
1496 { GNM_FUNC_HELP_NOTE, F_("@{time2} \xe2\x89\xa5 @{time1}")},
1497 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1498 { GNM_FUNC_HELP_END}
1502 /* Calculation of critical price options on options */
1503 static gnm_float
1504 CriticalValueOptionsOnOptions (OptionSide side, gnm_float x1, gnm_float x2, gnm_float t,
1505 gnm_float r, gnm_float b, gnm_float v)
1507 gnm_float si, ci, di, epsilon;
1509 si = x1;
1510 ci = opt_bs1 (side, si, x1, t, r, v, b);
1511 di = opt_bs_delta1 (side, si, x1, t, r, v, b);
1513 /* Newton-Raphson algorithm */
1514 epsilon = 0.0001;
1515 while (gnm_abs (ci - x2) > epsilon) {
1516 si = si - (ci - x2) / di;
1517 ci = opt_bs1 (side, si, x1, t, r, v, b);
1518 di = opt_bs_delta1 (side, si, x1, t, r, v, b);
1520 return si;
1523 /* Writer extendible options */
1524 static GnmValue *
1525 opt_extendible_writer (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1527 OptionSide call_put = option_side (value_peek_string (argv[0]));
1528 gnm_float s = value_get_as_float (argv[1]);
1529 gnm_float x1 = value_get_as_float (argv[2]);
1530 gnm_float x2 = value_get_as_float (argv[3]);
1531 gnm_float t1 = value_get_as_float (argv[4]);
1532 gnm_float t2 = value_get_as_float (argv[5]);
1533 gnm_float r = value_get_as_float (argv[6]);
1534 gnm_float b = value_get_as_float (argv[7]);
1535 gnm_float v = value_get_as_float (argv[8]);
1537 gnm_float rho = gnm_sqrt (t1 / t2);
1538 gnm_float z1 = (gnm_log (s / x2) + (b + (v * v) / 2.0) * t2) / (v * gnm_sqrt (t2));
1539 gnm_float z2 = (gnm_log (s / x1) + (b + (v * v) / 2.0) * t1) / (v * gnm_sqrt (t1));
1541 gnm_float gfresult;
1543 switch (call_put) {
1544 case OS_Call:
1545 gfresult = opt_bs1 (call_put, s, x1, t1, r, v, b) +
1546 s * gnm_exp ((b - r) * t2) * cum_biv_norm_dist1 (z1, -z2, -rho) -
1547 x2 * gnm_exp (-r * t2) * cum_biv_norm_dist1 (z1 - gnm_sqrt ((v * v) * t2), -z2 + gnm_sqrt ((v * v) * t1), -rho);
1548 break;
1549 case OS_Put:
1550 gfresult = opt_bs1 (call_put, s, x1, t1, r, v, b) +
1551 x2 * gnm_exp (-r * t2) * cum_biv_norm_dist1 (-z1 + gnm_sqrt ((v * v) * t2), z2 - gnm_sqrt ((v * v) * t1), -rho) -
1552 s * gnm_exp ((b - r) * t2) * cum_biv_norm_dist1 (-z1, z2, -rho);
1553 break;
1554 default:
1555 return value_new_error_NUM (ei->pos);
1558 return value_new_float (gfresult);
1561 static GnmFuncHelp const help_opt_extendible_writer[] = {
1562 { GNM_FUNC_HELP_NAME, F_("OPT_EXTENDIBLE_WRITER:theoretical price of extendible writer options")},
1563 DEF_ARG_CALL_PUT_FLAG,
1564 DEF_ARG_SPOT,
1565 { GNM_FUNC_HELP_ARG, F_("strike1:strike price at which the option is struck")},
1566 { GNM_FUNC_HELP_ARG, F_("strike2:strike price at which the option is re-struck if out of the money at @{time1}")},
1567 { GNM_FUNC_HELP_ARG, F_("time1:initial maturity of the option in years")},
1568 { GNM_FUNC_HELP_ARG, F_("time2:extended maturity in years if chosen")},
1569 DEF_ARG_RATE_RISKFREE_ANN,
1570 DEF_ARG_CC,
1571 DEF_ARG_VOLATILITY_SHORT,
1572 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_EXTENDIBLE_WRITER models the theoretical price of extendible "
1573 "writer options. These are options that have their maturity "
1574 "extended to @{time2} if the option is "
1575 "out of the money at @{time1}.")},
1576 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1577 { GNM_FUNC_HELP_END}
1582 /* Two asset correlation options */
1583 static GnmValue *
1584 opt_2_asset_correlation(GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1586 OptionSide call_put = option_side (value_peek_string(argv[0]));
1587 gnm_float s1 = value_get_as_float (argv[1]);
1588 gnm_float s2 = value_get_as_float (argv[2]);
1589 gnm_float x1 = value_get_as_float (argv[3]);
1590 gnm_float x2 = value_get_as_float (argv[4]);
1591 gnm_float t = value_get_as_float (argv[5]);
1592 gnm_float b1 = value_get_as_float (argv[6]);
1593 gnm_float b2 = value_get_as_float (argv[7]);
1594 gnm_float r = value_get_as_float (argv[8]);
1595 gnm_float v1 = value_get_as_float (argv[9]);
1596 gnm_float v2 = value_get_as_float (argv[10]);
1597 gnm_float rho = value_get_as_float (argv[11]);
1599 gnm_float y1 = (gnm_log (s1 / x1) + (b1 - (v1 * v1) / 2.0) * t) / (v1 * gnm_sqrt (t));
1600 gnm_float y2 = (gnm_log (s2 / x2) + (b2 - (v2 * v2) / 2.0) * t) / (v2 * gnm_sqrt (t));
1602 if (call_put == OS_Call) {
1603 return value_new_float (s2 * gnm_exp ((b2 - r) * t)
1604 * cum_biv_norm_dist1 (y2 + v2 * gnm_sqrt (t), y1 + rho * v2 * gnm_sqrt (t), rho)
1605 - x2 * gnm_exp (-r * t) * cum_biv_norm_dist1 (y2, y1, rho));
1606 } else if (call_put == OS_Put) {
1607 return value_new_float (x2 * gnm_exp (-r * t) * cum_biv_norm_dist1 (-y2, -y1, rho)
1608 - s2 * gnm_exp ((b2 - r) * t) * cum_biv_norm_dist1 (-y2 - v2 * gnm_sqrt (t), -y1 - rho * v2 * gnm_sqrt (t), rho));
1609 } else
1610 return value_new_error_NUM (ei->pos);
1613 static GnmFuncHelp const help_opt_2_asset_correlation[] = {
1614 { GNM_FUNC_HELP_NAME, F_("OPT_2_ASSET_CORRELATION:theoretical price of options on 2 assets with correlation @{rho}")},
1615 DEF_ARG_CALL_PUT_FLAG,
1616 { GNM_FUNC_HELP_ARG, F_("spot1:spot price of the underlying asset of the first option")},
1617 { GNM_FUNC_HELP_ARG, F_("spot2:spot price of the underlying asset of the second option")},
1618 { GNM_FUNC_HELP_ARG, F_("strike1:strike prices of the first option")},
1619 { GNM_FUNC_HELP_ARG, F_("strike2:strike prices of the second option")},
1620 DEF_ARG_TIME_MATURITY_Y,
1621 { GNM_FUNC_HELP_ARG, F_("cost_of_carry1:net cost of holding the underlying asset of the first option "
1622 "(for common stocks, the risk free rate less the dividend yield)")},
1623 { GNM_FUNC_HELP_ARG, F_("cost_of_carry2:net cost of holding the underlying asset of the second option "
1624 "(for common stocks, the risk free rate less the dividend yield)")},
1625 DEF_ARG_RATE_RISKFREE_ANN,
1626 { GNM_FUNC_HELP_ARG, F_("volatility1:annualized volatility in price of the underlying asset of the first option")},
1627 { GNM_FUNC_HELP_ARG, F_("volatility2:annualized volatility in price of the underlying asset of the second option")},
1628 { GNM_FUNC_HELP_ARG, F_("rho:correlation between the two underlying assets")},
1629 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_2_ASSET_CORRELATION models the theoretical price of options "
1630 "on 2 assets with correlation @{rho}. The payoff for a call is "
1631 "max(@{spot2} - @{strike2},0) if @{spot1} > @{strike1} or 0 otherwise. "
1632 "The payoff for a put is max (@{strike2} - @{spot2}, 0) if @{spot1} < @{strike1} or 0 otherwise.")},
1633 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1634 { GNM_FUNC_HELP_END}
1638 /* European option to exchange one asset for another */
1639 static GnmValue *
1640 opt_euro_exchange(GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1642 gnm_float s1 = value_get_as_float (argv[0]);
1643 gnm_float s2 = value_get_as_float (argv[1]);
1644 gnm_float q1 = value_get_as_float (argv[2]);
1645 gnm_float q2 = value_get_as_float (argv[3]);
1646 gnm_float t = value_get_as_float (argv[4]);
1647 gnm_float r = value_get_as_float (argv[5]);
1648 gnm_float b1 = value_get_as_float (argv[6]);
1649 gnm_float b2 = value_get_as_float (argv[7]);
1650 gnm_float v1 = value_get_as_float (argv[8]);
1651 gnm_float v2 = value_get_as_float (argv[9]);
1652 gnm_float rho = value_get_as_float (argv[10]);
1653 gnm_float v, d1, d2;
1655 v = gnm_sqrt (v1 * v1 + v2 * v2 - 2 * rho * v1 * v2);
1656 d1 = (gnm_log (q1 * s1 / (q2 * s2)) + (b1 - b2 + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
1657 d2 = d1 - v * gnm_sqrt (t);
1659 return value_new_float (q1 * s1 * gnm_exp ((b1 - r) * t) * ncdf (d1) -
1660 q2 * s2 * gnm_exp ((b2 - r) * t) * ncdf (d2));
1663 static GnmFuncHelp const help_opt_euro_exchange[] = {
1664 { GNM_FUNC_HELP_NAME, F_("OPT_EURO_EXCHANGE:theoretical price of a European option to exchange assets")},
1665 { GNM_FUNC_HELP_ARG, F_("spot1:spot price of asset 1")},
1666 { GNM_FUNC_HELP_ARG, F_("spot2:spot price of asset 2")},
1667 { GNM_FUNC_HELP_ARG, F_("qty1:quantity of asset 1")},
1668 { GNM_FUNC_HELP_ARG, F_("qty2:quantity of asset 2")},
1669 DEF_ARG_TIME_MATURITY_Y,
1670 DEF_ARG_RATE_RISKFREE_ANN,
1671 { GNM_FUNC_HELP_ARG, F_("cost_of_carry1:net cost of holding asset 1 "
1672 "(for common stocks, the risk free rate less the dividend yield)")},
1673 { GNM_FUNC_HELP_ARG, F_("cost_of_carry2:net cost of holding asset 2 "
1674 "(for common stocks, the risk free rate less the dividend yield)")},
1675 { GNM_FUNC_HELP_ARG, F_("volatility1:annualized volatility in price of asset 1")},
1676 { GNM_FUNC_HELP_ARG, F_("volatility2:annualized volatility in price of asset 2")},
1677 { GNM_FUNC_HELP_ARG, F_("rho:correlation between the prices of the two assets")},
1678 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_EURO_EXCHANGE models the theoretical price of a European "
1679 "option to exchange one asset with quantity @{qty2} and spot "
1680 "price @{spot2} for another with quantity @{qty1} and spot price "
1681 "@{spot1}.")},
1682 { GNM_FUNC_HELP_SEEALSO, "OPT_AMER_EXCHANGE,OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1683 { GNM_FUNC_HELP_END}
1687 /* American option to exchange one asset for another */
1688 static GnmValue *
1689 opt_amer_exchange(GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1691 gnm_float s1 = value_get_as_float (argv[0]);
1692 gnm_float s2 = value_get_as_float (argv[1]);
1693 gnm_float q1 = value_get_as_float (argv[2]);
1694 gnm_float q2 = value_get_as_float (argv[3]);
1695 gnm_float t = value_get_as_float (argv[4]);
1696 gnm_float r = value_get_as_float (argv[5]);
1697 gnm_float b1 = value_get_as_float (argv[6]);
1698 gnm_float b2 = value_get_as_float (argv[7]);
1699 gnm_float v1 = value_get_as_float (argv[8]);
1700 gnm_float v2 = value_get_as_float (argv[9]);
1701 gnm_float rho = value_get_as_float (argv[10]);
1702 gnm_float v = gnm_sqrt (v1 * v1 + v2 * v2 - 2 * rho * v1 * v2);
1704 return value_new_float (opt_bjer_stens1 (OS_Call, q1 * s1, q2 * s2, t, r - b2, v,b1 - b2));
1707 static GnmFuncHelp const help_opt_amer_exchange[] = {
1708 { GNM_FUNC_HELP_NAME, F_("OPT_AMER_EXCHANGE:theoretical price of an American option to exchange assets")},
1709 { GNM_FUNC_HELP_ARG, F_("spot1:spot price of asset 1")},
1710 { GNM_FUNC_HELP_ARG, F_("spot2:spot price of asset 2")},
1711 { GNM_FUNC_HELP_ARG, F_("qty1:quantity of asset 1")},
1712 { GNM_FUNC_HELP_ARG, F_("qty2:quantity of asset 2")},
1713 DEF_ARG_TIME_MATURITY_Y,
1714 DEF_ARG_RATE_RISKFREE_ANN,
1715 { GNM_FUNC_HELP_ARG, F_("cost_of_carry1:net cost of holding asset 1 "
1716 "(for common stocks, the risk free rate less the dividend yield)")},
1717 { GNM_FUNC_HELP_ARG, F_("cost_of_carry2:net cost of holding asset 2 "
1718 "(for common stocks, the risk free rate less the dividend yield)")},
1719 { GNM_FUNC_HELP_ARG, F_("volatility1:annualized volatility in price of asset 1")},
1720 { GNM_FUNC_HELP_ARG, F_("volatility2:annualized volatility in price of asset 2")},
1721 { GNM_FUNC_HELP_ARG, F_("rho:correlation between the prices of the two assets")},
1722 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_AMER_EXCHANGE models the theoretical price of an American "
1723 "option to exchange one asset with quantity @{qty2} and spot "
1724 "price @{spot2} for another with quantity @{qty1} and spot price "
1725 "@{spot1}.")},
1726 { GNM_FUNC_HELP_SEEALSO, "OPT_EURO_EXCHANGE,OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1727 { GNM_FUNC_HELP_END}
1731 /* Spread option approximation */
1732 static GnmValue *
1733 opt_spread_approx(GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1735 OptionSide call_put_flag = option_side (value_peek_string(argv[0]));
1736 gnm_float f1 = value_get_as_float (argv[1]);
1737 gnm_float f2 = value_get_as_float (argv[2]);
1738 gnm_float x = value_get_as_float (argv[3]);
1739 gnm_float t = value_get_as_float (argv[4]);
1740 gnm_float r = value_get_as_float (argv[5]);
1741 gnm_float v1 = value_get_as_float (argv[6]);
1742 gnm_float v2 = value_get_as_float (argv[7]);
1743 gnm_float rho = value_get_as_float (argv[8]);
1745 gnm_float v = gnm_sqrt (v1 * v1 + gnm_pow ((v2 * f2 / (f2 + x)), 2) - 2 * rho * v1 * v2 * f2 / (f2 + x));
1746 gnm_float F = f1 / (f2 + x);
1748 return value_new_float (opt_bs1 (call_put_flag, F, 1.0, t, r, v, 0.0) * (f2 + x));
1751 static GnmFuncHelp const help_opt_spread_approx[] = {
1752 { GNM_FUNC_HELP_NAME, F_("OPT_SPREAD_APPROX:theoretical price of a European option on the spread between two futures contracts")},
1753 DEF_ARG_CALL_PUT_FLAG,
1754 { GNM_FUNC_HELP_ARG, F_("fut_price1:price of the first futures contract")},
1755 { GNM_FUNC_HELP_ARG, F_("fut_price2:price of the second futures contract")},
1756 DEF_ARG_STRIKE,
1757 DEF_ARG_TIME_MATURITY_Y,
1758 DEF_ARG_RATE_RISKFREE_ANN,
1759 { GNM_FUNC_HELP_ARG, F_("volatility1:annualized volatility in price of the first underlying futures contract")},
1760 { GNM_FUNC_HELP_ARG, F_("volatility2:annualized volatility in price of the second underlying futures contract")},
1761 { GNM_FUNC_HELP_ARG, F_("rho:correlation between the two futures contracts")},
1762 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1763 { GNM_FUNC_HELP_END}
1767 /* Floating strike lookback options */
1768 static GnmValue *
1769 opt_float_strk_lkbk(GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1771 OptionSide call_put_flag = option_side (value_peek_string(argv[0]));
1772 gnm_float s = value_get_as_float (argv[1]);
1773 gnm_float s_min = value_get_as_float (argv[2]);
1774 gnm_float s_max = value_get_as_float (argv[3]);
1775 gnm_float t = value_get_as_float (argv[4]);
1776 gnm_float r = value_get_as_float (argv[5]);
1777 gnm_float b = value_get_as_float (argv[6]);
1778 gnm_float v = value_get_as_float (argv[7]);
1780 gnm_float a1, a2, m;
1782 if(OS_Call == call_put_flag)
1783 m = s_min;
1784 else if(OS_Put == call_put_flag)
1785 m = s_max;
1786 else
1787 return value_new_error_NUM (ei->pos);
1789 a1 = (gnm_log (s / m) + (b + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
1790 a2 = a1 - v * gnm_sqrt (t);
1792 if(OS_Call == call_put_flag)
1793 return value_new_float (s * gnm_exp ((b - r) * t) * ncdf (a1) -
1794 m * gnm_exp (-r * t) * ncdf (a2) +
1795 gnm_exp (-r * t) * (v * v) / (2 * b) * s * (gnm_pow (s / m, (-2 * b / (v * v))) * ncdf (-a1 + 2 * b / v * gnm_sqrt (t)) -
1796 gnm_exp (b * t) * ncdf (-a1)));
1797 else if(OS_Put == call_put_flag)
1798 return value_new_float (m * gnm_exp (-r * t) * ncdf (-a2) -
1799 s * gnm_exp ((b - r) * t) * ncdf (-a1) +
1800 gnm_exp (-r * t) * (v * v) / (2 * b) * s * (-gnm_pow (s / m, ((-2 * b) / (v * v))) * ncdf (a1 - 2 * b / v * gnm_sqrt (t)) +
1801 gnm_exp (b * t) * ncdf (a1)));
1803 return value_new_error_VALUE (ei->pos);
1806 static GnmFuncHelp const help_opt_float_strk_lkbk[] = {
1807 { GNM_FUNC_HELP_NAME, F_("OPT_FLOAT_STRK_LKBK:theoretical price of floating-strike lookback option")},
1808 DEF_ARG_CALL_PUT_FLAG,
1809 DEF_ARG_SPOT,
1810 { GNM_FUNC_HELP_ARG, F_("spot_min:minimum spot price of the underlying asset so far observed")},
1811 { GNM_FUNC_HELP_ARG, F_("spot_max:maximum spot price of the underlying asset so far observed")},
1812 DEF_ARG_TIME_MATURITY_Y,
1813 DEF_ARG_RATE_RISKFREE_ANN,
1814 DEF_ARG_CC,
1815 DEF_ARG_VOLATILITY_SHORT,
1816 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_FLOAT_STRK_LKBK determines the theoretical price of a "
1817 "floating-strike lookback option where the holder "
1818 "of the option may exercise on expiry at the most favourable price "
1819 "observed during the options life of the underlying asset.")},
1820 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1821 { GNM_FUNC_HELP_END}
1825 /* Fixed strike lookback options */
1827 static GnmValue *
1828 opt_fixed_strk_lkbk(GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1830 OptionSide call_put_flag = option_side (value_peek_string(argv[0]));
1831 gnm_float s = value_get_as_float (argv[1]);
1832 gnm_float s_min = value_get_as_float (argv[2]);
1833 gnm_float s_max = value_get_as_float (argv[3]);
1834 gnm_float x = value_get_as_float (argv[4]);
1835 gnm_float t = value_get_as_float (argv[5]);
1836 gnm_float r = value_get_as_float (argv[6]);
1837 gnm_float b = value_get_as_float (argv[7]);
1838 gnm_float v = value_get_as_float (argv[8]);
1840 gnm_float d1, d2;
1841 gnm_float e1, e2, m;
1843 if (OS_Call == call_put_flag)
1844 m = s_max;
1845 else if (OS_Put == call_put_flag)
1846 m = s_min;
1847 else
1848 return value_new_error_VALUE (ei->pos);
1850 d1 = (gnm_log (s / x) + (b + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
1851 d2 = d1 - v * gnm_sqrt (t);
1852 e1 = (gnm_log (s / m) + (b + (v * v) / 2.0) * t) / (v * gnm_sqrt (t));
1853 e2 = e1 - v * gnm_sqrt (t);
1855 if (OS_Call == call_put_flag && x > m)
1856 return value_new_float (s * gnm_exp ((b - r) * t) * ncdf (d1) - x * gnm_exp (-r * t) * ncdf (d2) + s * gnm_exp (-r * t) * (v * v) / (2 * b) * (-gnm_pow ((s / x), (-2 * b / (v * v))) * ncdf (d1 - 2 * b / v * gnm_sqrt (t)) + gnm_exp (b * t) * ncdf (d1)));
1858 else if (OS_Call == call_put_flag && x <= m)
1859 return value_new_float (gnm_exp (-r * t) * (m - x) + s * gnm_exp ((b - r) * t) * ncdf (e1) - gnm_exp (-r * t) * m * ncdf (e2) + s * gnm_exp (-r * t) * (v * v) / (2 * b) * (-gnm_pow ((s / m), (-2 * b / (v * v))) * ncdf (e1 - 2 * b / v * gnm_sqrt (t)) + gnm_exp (b * t) * ncdf (e1)));
1861 else if (OS_Put == call_put_flag && x < m)
1862 return value_new_float (-s * gnm_exp ((b - r) * t) * ncdf (-d1) + x * gnm_exp (-r * t) * ncdf (-d1 + v * gnm_sqrt (t)) + s * gnm_exp (-r * t) * (v * v) / (2 * b) * (gnm_pow ((s / x), (-2 * b / (v * v))) * ncdf (-d1 + 2 * b / v * gnm_sqrt (t)) - gnm_exp (b * t) * ncdf (-d1)));
1864 else if (OS_Put == call_put_flag && x >= m)
1865 return value_new_float (gnm_exp (-r * t) * (x - m) - s * gnm_exp ((b - r) * t) * ncdf (-e1) + gnm_exp (-r * t) * m * ncdf (-e1 + v * gnm_sqrt (t)) + gnm_exp (-r * t) * (v * v) / (2 * b) * s * (gnm_pow ((s / m), (-2 * b / (v * v))) * ncdf (-e1 + 2 * b / v * gnm_sqrt (t)) - gnm_exp (b * t) * ncdf (-e1)));
1867 return value_new_error_VALUE (ei->pos);
1870 static GnmFuncHelp const help_opt_fixed_strk_lkbk[] = {
1871 { GNM_FUNC_HELP_NAME, F_("OPT_FIXED_STRK_LKBK:theoretical price of a fixed-strike lookback option")},
1872 DEF_ARG_CALL_PUT_FLAG,
1873 DEF_ARG_SPOT,
1874 { GNM_FUNC_HELP_ARG, F_("spot_min:minimum spot price of the underlying asset so far observed")},
1875 { GNM_FUNC_HELP_ARG, F_("spot_max:maximum spot price of the underlying asset so far observed")},
1876 DEF_ARG_STRIKE,
1877 DEF_ARG_TIME_MATURITY_Y,
1878 DEF_ARG_RATE_RISKFREE_ANN,
1879 DEF_ARG_CC,
1880 DEF_ARG_VOLATILITY_SHORT,
1881 { GNM_FUNC_HELP_DESCRIPTION, F_("OPT_FIXED_STRK_LKBK determines the theoretical price of a "
1882 "fixed-strike lookback option where the holder "
1883 "of the option may exercise on expiry at the most favourable price "
1884 "observed during the options life of the underlying asset.")},
1885 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1886 { GNM_FUNC_HELP_END}
1891 /* Binomial Tree valuation */
1892 static GnmValue *
1893 opt_binomial(GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1895 OptionType amer_euro_flag = option_type(value_peek_string(argv[0]));
1896 OptionSide call_put_flag = option_side (value_peek_string(argv[1]));
1897 gnm_float n = gnm_floor (value_get_as_float (argv[2]));
1898 gnm_float s = value_get_as_float (argv[3]);
1899 gnm_float x = value_get_as_float (argv[4]);
1900 gnm_float t = value_get_as_float (argv[5]);
1901 gnm_float r = value_get_as_float (argv[6]);
1902 gnm_float v = value_get_as_float (argv[7]);
1903 gnm_float b = argv[8] ? value_get_as_float (argv[8]) : 0;
1905 gnm_float *value_array;
1906 gnm_float u, d, p, dt, Df, temp1, temp2, gf_result;
1907 gint i, j, z;
1909 if (n < 0 || n > 100000)
1910 return value_new_error_NUM (ei->pos);
1912 if (OS_Call == call_put_flag)
1913 z = 1;
1914 else if (OS_Put == call_put_flag)
1915 z = -1;
1916 else
1917 return value_new_error_NUM (ei->pos);
1919 if (OT_Error == amer_euro_flag)
1920 return value_new_error_NUM (ei->pos);
1922 value_array = (gnm_float *) g_try_malloc ((n + 2)* sizeof(gnm_float));
1923 if (value_array == NULL)
1924 return value_new_error_NUM (ei->pos);
1926 dt = t / n;
1927 u = gnm_exp (v * gnm_sqrt (dt));
1928 d = 1.0 / u;
1929 p = (gnm_exp (b * dt) - d) / (u - d);
1930 Df = gnm_exp (-r * dt);
1932 for (i = 0; i <= n; ++i) {
1933 temp1 = z * (s * gnm_pow (u, i) * gnm_pow (d, (n - i)) - x);
1934 value_array[i] = MAX (temp1, 0.0);
1937 for (j = n - 1; j > -1; --j) {
1938 for (i = 0; i <= j; ++i) {
1939 /*if (0==i)g_printerr("secondloop %d\n",j);*/
1940 if (OT_Euro == amer_euro_flag)
1941 value_array[i] = (p * value_array[i + 1] + (1.0 - p) * value_array[i]) * Df;
1942 else if (OT_Amer == amer_euro_flag) {
1943 temp1 = (z * (s * gnm_pow (u, i) * gnm_pow (d, (gnm_abs (i - j))) - x));
1944 temp2 = (p * value_array[i + 1] + (1.0 - p) * value_array[i]) * Df;
1945 value_array[i] = MAX (temp1, temp2);
1949 gf_result = value_array[0];
1950 g_free (value_array);
1951 return value_new_float (gf_result);
1954 static GnmFuncHelp const help_opt_binomial[] = {
1955 { GNM_FUNC_HELP_NAME, F_("OPT_BINOMIAL:theoretical price of either an American or European style option using a binomial tree")},
1956 { GNM_FUNC_HELP_ARG, F_("amer_euro_flag:'a' for an American style option or 'e' for a European style option")},
1957 DEF_ARG_CALL_PUT_FLAG,
1958 { GNM_FUNC_HELP_ARG, F_("num_time_steps:number of time steps used in the valuation")},
1959 DEF_ARG_SPOT,
1960 DEF_ARG_STRIKE,
1961 DEF_ARG_TIME_MATURITY_Y,
1962 DEF_ARG_RATE_RISKFREE_ANN,
1963 DEF_ARG_VOLATILITY_SHORT,
1964 DEF_ARG_CC,
1965 { GNM_FUNC_HELP_NOTE, F_("A larger @{num_time_steps} yields greater accuracy but OPT_BINOMIAL is slower to calculate.")},
1966 { GNM_FUNC_HELP_SEEALSO, "OPT_BS,OPT_BS_DELTA,OPT_BS_RHO,OPT_BS_THETA,OPT_BS_GAMMA"},
1967 { GNM_FUNC_HELP_END}
1972 GnmFuncDescriptor const derivatives_functions [] = {
1973 { "opt_bs",
1974 "sfffff|f",
1975 help_opt_bs, opt_bs, NULL, NULL, NULL,
1976 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
1978 { "opt_bs_delta",
1979 "sfffff|f",
1980 help_opt_bs_delta, opt_bs_delta, NULL, NULL, NULL,
1981 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
1983 { "opt_bs_rho",
1984 "sfffff|f",
1985 help_opt_bs_rho, opt_bs_rho, NULL, NULL, NULL,
1986 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
1988 { "opt_bs_theta",
1989 "sfffff|f",
1990 help_opt_bs_theta, opt_bs_theta, NULL, NULL, NULL,
1991 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
1993 { "opt_bs_gamma",
1994 "fffff|f",
1995 help_opt_bs_gamma, opt_bs_gamma, NULL, NULL, NULL,
1996 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
1998 { "opt_bs_vega",
1999 "fffff|f",
2000 help_opt_bs_vega, opt_bs_vega, NULL, NULL, NULL,
2001 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2003 { "opt_bs_carrycost",
2004 "sfffff|f",
2005 help_opt_bs_carrycost, opt_bs_carrycost, NULL, NULL, NULL,
2006 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2008 { "cum_biv_norm_dist",
2009 "fff",
2010 help_cum_biv_norm_dist, cum_biv_norm_dist, NULL, NULL, NULL,
2011 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE },
2013 { "opt_garman_kohlhagen",
2014 "sffffff",
2015 help_opt_garman_kohlhagen, opt_garman_kohlhagen, NULL, NULL, NULL,
2016 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2018 { "opt_french",
2019 "sfffffff",
2020 help_opt_french, opt_french, NULL, NULL, NULL,
2021 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2023 { "opt_jump_diff",
2024 "sfffffff",
2025 help_opt_jump_diff, opt_jump_diff, NULL, NULL, NULL,
2026 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2028 { "opt_exec",
2029 "sfffffff",
2030 help_opt_exec, opt_exec, NULL, NULL, NULL,
2031 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2033 { "opt_bjer_stens",
2034 "sffffff",
2035 help_opt_bjer_stens, opt_bjer_stens, NULL, NULL, NULL,
2036 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2038 { "opt_miltersen_schwartz",
2039 "sfffffffffffff",
2040 help_opt_miltersen_schwartz, opt_miltersen_schwartz, NULL, NULL, NULL,
2041 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2043 { "opt_baw_amer",
2044 "sffffff",
2045 help_opt_baw_amer, opt_baw_amer, NULL, NULL, NULL,
2046 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2048 { "opt_rgw",
2049 "fffffff",
2050 help_opt_rgw, opt_rgw, NULL, NULL, NULL,
2051 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2053 { "opt_forward_start",
2054 "sfffffff",
2055 help_opt_forward_start, opt_forward_start, NULL, NULL, NULL,
2056 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2058 { "opt_time_switch",
2059 "sfffffffff",
2060 help_opt_time_switch, opt_time_switch, NULL, NULL, NULL,
2061 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2063 { "opt_simple_chooser",
2064 "fffffff",
2065 help_opt_simple_chooser, opt_simple_chooser, NULL, NULL, NULL,
2066 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2068 { "opt_complex_chooser",
2069 "fffffffff",
2070 help_opt_complex_chooser, opt_complex_chooser, NULL, NULL, NULL,
2071 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2073 { "opt_on_options",
2074 "sffffffff",
2075 help_opt_on_options, opt_on_options, NULL, NULL, NULL,
2076 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2078 { "opt_extendible_writer",
2079 "sffffffff",
2080 help_opt_extendible_writer, opt_extendible_writer, NULL, NULL, NULL,
2081 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2083 { "opt_2_asset_correlation",
2084 "sfffffffffff",
2085 help_opt_2_asset_correlation, opt_2_asset_correlation, NULL, NULL, NULL,
2086 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2088 { "opt_euro_exchange",
2089 "fffffffffff",
2090 help_opt_euro_exchange, opt_euro_exchange, NULL, NULL, NULL,
2091 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2093 { "opt_amer_exchange",
2094 "fffffffffff",
2095 help_opt_amer_exchange, opt_amer_exchange, NULL, NULL, NULL,
2096 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2098 { "opt_spread_approx",
2099 "sffffffff",
2100 help_opt_spread_approx, opt_spread_approx, NULL, NULL, NULL,
2101 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2103 { "opt_float_strk_lkbk",
2104 "sfffffff",
2105 help_opt_float_strk_lkbk, opt_float_strk_lkbk, NULL, NULL, NULL,
2106 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2108 { "opt_fixed_strk_lkbk",
2109 "sffffffff",
2110 help_opt_fixed_strk_lkbk, opt_fixed_strk_lkbk, NULL, NULL, NULL,
2111 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2114 { "opt_binomial",
2115 "ssffffff|f",
2116 help_opt_binomial, opt_binomial, NULL, NULL, NULL,
2117 GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_BASIC },
2119 { NULL}