QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
perturbativebarrieroptionengine.cpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2008 Lorella Fatone
5 Copyright (C) 2008 Maria Cristina Recchioni
6 Copyright (C) 2008 Francesco Zirilli
7 Copyright (C) 2008 StatPro Italia srl
8
9 This file is part of QuantLib, a free-software/open-source library
10 for financial quantitative analysts and developers - http://quantlib.org/
11
12 QuantLib is free software: you can redistribute it and/or modify it
13 under the terms of the QuantLib license. You should have received a
14 copy of the license along with this program; if not, please email
15 <quantlib-dev@lists.sf.net>. The license is also available online at
16 <http://quantlib.org/license.shtml>.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 FOR A PARTICULAR PURPOSE. See the license for more details.
21*/
22
23#include <ql/errors.hpp>
24#include <ql/exercise.hpp>
25#include <ql/experimental/barrieroption/perturbativebarrieroptionengine.hpp>
26#include <ql/functional.hpp>
27#include <ql/types.hpp>
28#include <algorithm>
29#include <cmath>
30#include <utility>
31
32using namespace std;
33
34namespace {
35
36 inline QuantLib::Real SIGN(const QuantLib::Real& a, const QuantLib::Real& b)
37 {
38 if (b > 0.0)
39 return std::fabs(a);
40 else
41 return -std::fabs(a);
42 }
43
44}
45
46#define PI 3.14159265358979324
47
48namespace QuantLib {
49
50 namespace {
51
52 Real ND2(Real a, Real b, Real rho);
53
54 Real H1, H2, H3, R23, RUA, RUB, AR, RUC;
55 int NUC;
56
57 // standard normal cumulative distribution function
58 Real PHID(Real Z);
59
60 // Functions used to compute the first order approximation
61 Real ff(Real p,Real tt,Real a, Real b, Real gm);
62 Real v(Real p, Real tt,Real a,Real b,Real gm);
63 Real llold(Real p,Real tt, Real a, Real b,
64 Real c, Real gm);
65
66 // Functions used to compute the second order approximation
67 Real derivn3(Real limit[4],Real sigmarho[4], int idx);
68 Real ddvv(Real s, Real p, Real tt, Real a,
69 Real b, Real gm);
70 Real ddff(Real s, Real p,Real tt,Real a,Real b,Real gm);
71 Real dll(Real s,Real p,Real tt,Real a,Real b,
72 Real c,Real gm);
73 Real ddll(Real s,Real p,Real tt, Real ax, Real bx,
74 Real c, Real gm);
75 Real dvv(Real s,Real p,Real tt,Real a,Real b,Real gm);
76 Real dff(Real s, Real p,Real tt,Real a,Real b,Real gm);
77 Real tvtl(int jj, const Real limit[4], const Real sigmarho[4], Real epsi);
78
79 Real BarrierUPD(Real kprice,
80 Real stock,
81 Real hbarr,
82 Real taumin,
83 Real taumax,
84 int iord,
85 int igm,
86 const ext::function<Real(Real, Real)>& integr,
87 const ext::function<Real(Real, Real)>& integalpha,
88 const ext::function<Real(Real, Real)>& integs,
89 const ext::function<Real(Real)>& alpha,
90 const ext::function<Real(Real)>& sigmaq) {
91 Real v0=0.0, v1=0.0, v1p=0.0, v2p=0.0, v2pp=0.0, gm=0.0;
92 int i=0,j=0;
93 Real tmp=0.0, e1=0.0, e2=0.0, e3=0.0, e4=0.0;
94 Real xstar=0.0, s0=0.0;
95 Real sigmat=0.0, disc=0.0, d1=0.0,d2=0.0,d3=0.0,d4=0.0;
96 Real et=0.0,tt=0.0, dt=0.0,p=0.0;
97 int npoint,npoint2;
98 static double pi= 3.14159265358979324;
99 Real dsqpi;
100 Real caux=0.0,ccaux=0.0;
101 Real auxnew=0.0;
102 Real x=0.0,b=0.0,c=0.0;
103
104 if(igm==0) {
105 gm=0.0;
106 } else if(igm==1) {
107 gm=integalpha(taumin,taumax)/(0.5*integs(taumin,taumax));
108 } else {
109 gm=0.0;
110 }
111
112 /*
113 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
114 !! !!
115 !! xstar=min(0,log(kprice/hbarr)) !!
116 !! !!
117 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
118 */
119
120 xstar=log(kprice/hbarr);
121
122 if(xstar>0.0) xstar=0.0;
123 sigmat=integs(taumin,taumax);
124 disc=-integr(taumin,taumax);
125
126 /*
127 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
128 !! Change of variable !!
129 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
130 */
131 s0=stock/hbarr;
132
133 /*
134 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
135 !! !!
136 !! Compute the zero-th order term P_0 !!
137 !! !!
138 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
139 */
140 d1=(xstar-log(s0)+(1.0-gm)*0.5*sigmat)/sqrt(sigmat);
141 d2=(xstar+log(s0)+(1.0-gm)*0.5*sigmat)/sqrt(sigmat);
142 d3=(xstar-log(s0)-(1.0+gm)*0.5*sigmat)/sqrt(sigmat);
143 d4=(xstar+log(s0)-(1.0+gm)*0.5*sigmat)/sqrt(sigmat);
144
145 e1=PHID(d1);
146 e2=PHID(d2);
147 e3=PHID(d3);
148 e4=PHID(d4);
149
150 v0=kprice*e1-kprice*std::pow(s0,(1.0-gm))*e2;
151 v0=v0+exp(gm*0.5*sigmat)*(-hbarr*s0*e3+hbarr*std::pow(s0,-gm)*e4);
152 v0=v0*exp(disc);
153
154 if(iord==0) return v0;
155
156 /*
157 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
158 !! !!
159 !! Compute the first order term P_1 !!
160 !! !!
161 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
162 */
163
164 npoint=1000;
165 npoint2=100;
166
167 dt=(taumax-taumin)/double(npoint);
168
169 tt=0.5*integs(taumin,taumax);
170
171 x=log(s0);
172 et=exp(0.5*(1.0-gm)*x);
173
174 dsqpi=std::pow(pi,0.5);
175
176 v1=0.0;
177 for( i=1;i<=npoint;i++) {
178 v1p=0.0;
179 tmp=taumin+dt*double(2*i-1)*0.5;
180 p=0.5*integs(tmp,taumax);
181 /*
182 !!
183 !! Function E(p,tt,a,b,gm)
184 !!
185 */
186 ccaux=v(p,tt,x,xstar,gm)+v(p,tt,x,-xstar,gm)-v(p,tt,-x,xstar,gm)-v(p,tt,-x,-xstar,gm);
187 auxnew=ccaux*(-kprice*exp(-xstar*0.5*(1.0-gm))+hbarr*exp(xstar*0.5*(1.0+gm)));
188 v1p=v1p+auxnew;
189 /*
190 !!
191 !! Function L(p,tt,a,b,c,gm)
192 !!
193 */
194 b=gm-1.0;
195 c=-xstar;
196 ccaux=llold(p,tt,x,b,c,gm)-llold(p,tt,-x,b,c,gm);
197 auxnew=kprice*(1.0-gm)*ccaux;
198 v1p=v1p+auxnew;
199
200 b=-(gm+1.0);
201 c=xstar;
202 ccaux=llold(p,tt,x,b,c,gm)-llold(p,tt,-x,b,c,gm);
203 auxnew=-exp(gm*p)*hbarr*ccaux;
204 v1p=v1p+auxnew;
205
206 b=(gm+1.0);
207 c=-xstar;
208 ccaux=llold(p,tt,x,b,c,gm)-llold(p,tt,-x,b,c,gm);
209 auxnew=exp(gm*p)*hbarr*gm*ccaux;
210 v1p=v1p+auxnew;
211 /*
212 !!
213 !! Function F(p,tt,a,b,c,gm)
214 !!
215 */
216 b=gm-1.0;
217 auxnew=-kprice*(1.0-gm)*(ff(p,tt,x,b,gm)-ff(p,tt,-x,b,gm));
218 v1p=v1p+auxnew;
219
220 b=gm+1.0;
221 auxnew=-exp(gm*p)*gm*hbarr*(ff(p,tt,x,b,gm)-ff(p,tt,-x,b,gm));
222 v1p=v1p+auxnew;
223
224 v1=v1+(alpha(tmp)-gm*0.5*sigmaq(tmp))*v1p;
225 }
226
227 v1=exp(disc)*et*v1*dt/(dsqpi*2.0);
228
229 if(iord==1) return v0+v1;
230
231 /*
232 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
233 !! !!
234 !! Compute the second order term P_2 !!
235 !! !!
236 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
237 */
238
239 Real v2,dtp, tmp1,s,caux2;
240 v2=0.0;
241
242 for(i=1;i<=npoint;i++) {
243 v2p=0.0;
244 tmp=taumin+dt*(double)(2*i-1)*0.5;
245 p=0.5*integs(tmp,taumax);
246
247 dtp=(taumax-tmp)/(double)(npoint2);
248
249 for(j=1;j<=npoint2; j++) {
250 tmp1=tmp+dtp*(double)(2*j-1)*0.50;
251 s=0.50*integs(tmp1,taumax);
252
253 caux=dll(s,p,tt,-x,-1.0+gm,-xstar,gm)-dll(s,p,tt,x,-1.0+gm,-xstar,gm);
254 v2pp=caux*kprice*(1.0-gm);
255
256 caux=dll(s,p,tt,-x,-1.0-gm,xstar,gm)-dll(s,p,tt,x,-1.0-gm,xstar,gm);
257 v2pp=v2pp-exp(gm*s)*hbarr*caux;
258
259 caux=dll(s,p,tt,-x,1.0+gm,-xstar,gm)-dll(s,p,tt,x,1.0+gm,-xstar,gm);
260 v2pp=v2pp+exp(gm*s)*gm*hbarr*caux;
261
262 caux=+dvv(s,p,tt,-x,xstar,gm)-dvv(s,p,tt,x,xstar,gm);
263 caux=caux+(dvv(s,p,tt,-x,-xstar,gm)-dvv(s,p,tt,x,-xstar,gm));
264 caux2=hbarr*exp(0.5*(1.0+gm)*xstar)-kprice*exp(-0.5*(1.0-gm)*xstar);
265 v2pp=v2pp+caux2*caux;
266
267 caux=dff(s,p,tt,-x,-1.0+gm,gm)-dff(s,p,tt,x,-1.0+gm,gm);
268 v2pp=v2pp-(1.0-gm)*kprice*caux;
269
270 caux=dff(s,p,tt,-x,1.0+gm,gm)-dff(s,p,tt,x,1.0+gm,gm);
271 v2pp=v2pp-exp(gm*s)*gm*hbarr*caux;
272
273 v2pp=v2pp*0.5*(1.0-gm);
274
275 caux=-ddll(s,p,tt,-x,-1.0+gm,-xstar,gm)+ddll(s,p,tt,x,-1.0+gm,-xstar,gm);
276 v2pp=v2pp+caux*kprice*(1.0-gm);
277
278 caux=-ddll(s,p,tt,-x,-1.0-gm,xstar,gm)+ddll(s,p,tt,x,-1.0-gm,xstar,gm);
279 v2pp=v2pp-exp(gm*s)*hbarr*caux;
280
281 caux=-ddll(s,p,tt,-x,1.0+gm,-xstar,gm)+ddll(s,p,tt,x,1.0+gm,-xstar,gm);
282 v2pp=v2pp+exp(gm*s)*gm*hbarr*caux;
283
284 caux=-ddvv(s,p,tt,-x,xstar,gm)+ddvv(s,p,tt,x,xstar,gm);
285 caux=caux+(-dvv(s,p,tt,-x,-xstar,gm)+dvv(s,p,tt,x,-xstar,gm));
286 caux2=hbarr*exp(0.5*(1.0+gm)*xstar)-kprice*exp(-0.5*(1-gm)*xstar);
287
288 v2pp=v2pp+caux2*caux;
289
290 caux=-ddff(s,p,tt,-x,-1+gm,gm)+ddff(s,p,tt,x,-1+gm,gm);
291 v2pp=v2pp-(1.0-gm)*kprice*caux;
292
293 caux=-ddff(s,p,tt,-x,1.0+gm,gm)+ddff(s,p,tt,x,1.0+gm,gm);
294 v2pp=v2pp-exp(gm*s)*gm*hbarr*caux;
295
296 v2p=v2p+(alpha(tmp1)-gm*0.5*sigmaq(tmp1))*v2pp;
297 }
298
299 v2=v2+v2p*(alpha(tmp)-gm*0.5*sigmaq(tmp))*dtp;
300 }
301
302 v2=exp(disc)*et*v2*dt;
303
304 return v0+v1+v2;
305 }
306
307
308 Real PHID(Real Z){
309 /*
310 * Normal distribution probabilities accurate to 1D-15.
311 * Z = number of standard deviations from the mean.
312 *
313 * The software that computes the normal distribution
314 * probabilities has been developed by M.C. Recchioni
315 * based upon algorithm 5666 (Programmer Alan Miller)
316 * for the error function, taken from:
317 * Hart, J.F. et al, 'Computer Approximations', Wiley, 1968
318 *
319 */
320 Real P0, P1, P2, P3, P4, P5, P6;
321 Real Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7;
322 Real P, EXPNTL, CUTOFF, ROOTPI, ZABS;
323
324 P0 = 220.2068679123761;
325 P1 = 221.2135961699311;
326 P2 = 112.0792914978709;
327 P3 = 33.91286607838300;
328 P4 = 6.373962203531650;
329 P5 = 0.7003830644436881;
330 P6 = 0.03526249659989109;
331
332 Q0 = 440.4137358247522;
333 Q1 = 793.8265125199484;
334 Q2 = 637.3336333788311;
335 Q3 = 296.5642487796737;
336 Q4 = 86.78073220294608;
337 Q5 = 16.064177579206950;
338 Q6 = 1.7556671631826420;
339 Q7 = 0.088388347648318440;
340 ROOTPI = 2.506628274631001;
341 CUTOFF = 7.071067811865475;
342
343 ZABS = fabs(Z);
344 /*
345 |Z| > 37
346 */
347 if (ZABS > 37)
348 P = 0;
349 else
350 {
351 /*
352 |Z| <= 37
353 */
354 EXPNTL =exp(-ZABS*ZABS/2);
355 /*
356 |Z| < CUTOFF = 10/SQRT(2)
357 */
358 if ( ZABS < CUTOFF )
359 P = EXPNTL*((((((P6*ZABS + P5)*ZABS + P4)*ZABS + P3)*ZABS+ P2)*ZABS + P1)*ZABS + P0)/(((((((Q7*ZABS + Q6)*ZABS + Q5)*ZABS + Q4)*ZABS + Q3)*ZABS + Q2)*ZABS + Q1)*ZABS + Q0);
360 /*
361 |Z| >= CUTOFF.
362 */
363 else
364 P = EXPNTL/(ZABS + 1/(ZABS + 2/(ZABS + 3/(ZABS + 4/(ZABS + 0.65)))))/ROOTPI;
365
366 }
367 if ( Z > 0 ) P = 1 - P;
368
369 return(P);
370 }
371
372
373 /*
374 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
375 !! !!
376 !! Functions needed to compute the first order term P_1 !!
377 !! !!
378 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
379 */
380
381 /*
382 !!
383 !! Function F(p,tt,a,b,gm)
384 !!
385 */
386 Real ff(Real p,Real tt,Real a, Real b, Real gm) {
387 Real phid;
388 Real aa, caux;
389 Real ppi= 3.14159265358979324;
390
391 aa=-(b*p-b*tt+a)/std::pow(2.0*(tt-p),0.5);
392
393 caux=2.0*std::pow(ppi,0.5)*PHID(aa);
394 aa=b*b-(1.0-gm)*(1.0-gm);
395 aa=aa/4.0;
396 phid=exp(-0.5*a*b)*exp(aa*(tt-p))*caux;
397
398 return phid;
399 }
400
401 /*
402 !!
403 !! Function E(p,tt,a,b,gm)
404 !!
405 */
406 Real v(Real p, Real tt,Real a,Real b,Real gm)
407 {
408 Real result;
409 Real aa,caux;
410
411 aa=-(p*(a-b)+b*tt)/std::pow(2.0*p*tt*(tt-p),0.5);
412 caux=PHID(aa);
413
414 aa=exp(std::pow((a-b),2)/(4.0*tt))*exp(std::pow((1.0-gm),2)*tt/4.0)*std::pow(tt,0.5);
415 result=caux/aa;
416
417 return(result);
418 }
419
420 /*
421 !!
422 !! Fuction L(p,tt,a,b,c,gm)
423 !!
424 */
425 Real llold(Real p,Real tt, Real a, Real b,Real c, Real gm){
426 Real bvnd;
427 Real xx,yy,rho,caux;
428 Real ppi= 3.14159265358979324;
429 Real aa;
430
431 xx=(-a+b*(tt-p))/std::pow(2.0*(tt-p),0.5);
432 yy=(-a+b*tt+c)/std::pow(2.0*tt,0.5);
433 rho=std::pow((tt-p)/tt,0.5);
434 aa=b*b-(1.0-gm)*(1.0-gm);
435 aa=aa/4.0;
436 caux=ND2(-xx,-yy,rho);
437
438 bvnd=2.0*std::pow(ppi,0.5)*exp(-a*b*0.5)*exp(aa*(tt-p))*caux;
439 return(bvnd);
440 }
441
442 /*
443 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
444 !! !!
445 !! Functions needed to compute the second order term P_2 !!
446 !! !!
447 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
448 */
449
450 /*
451 !!
452 !! Function D_E(s,p,tt,a,b,gm)
453 !!
454 */
455 Real dvv(Real s,Real p,Real tt,Real a,Real b,Real gm)
456 {
457 double ppi= 3.14159265358979324;
458 Real result;
459 Real aa,caux,caux1,caux2;
460 Real xx,yy,rho;
461
462 aa=(a*p+b*(tt-p))/std::pow(2.0*p*tt*(tt-p),0.5);
463 caux=PHID(aa);
464
465 aa=exp((a-b)*(a-b)/(4.0*tt))*exp(std::pow((1.0-gm),2)*tt/4.0)*std::pow(tt,0.5);
466 caux=-caux/aa;
467
468 xx=(a*p+b*(tt-p))/std::pow(2.0*tt*p*(tt-p),0.5);
469 yy=(a*s+b*(tt-s))/std::pow(2.0*tt*s*(tt-s),0.5);
470 rho=std::pow((s*(tt-p))/(p*(tt-s)),0.5);
471 caux1=ND2(-xx,-yy,rho);
472 caux1=caux1/aa;
473
474
475 aa=exp((a+b)*(a+b)/(4.0*tt))*exp(std::pow((1.0-gm),2)*tt/4.0)*std::pow(tt,0.5);
476
477 xx=(a*p-b*(tt-p))/std::pow(2.0*tt*p*(tt-p),0.5);
478 yy=(a*s-b*(tt-s))/std::pow(2.0*tt*s*(tt-s),0.5);
479 rho=std::pow((s*(tt-p))/(p*(tt-s)),0.5);
480 caux2=ND2(-xx,-yy,rho);
481 caux2=caux2/aa;
482 result=(caux+caux1+caux2)/(2.0*std::pow(ppi,0.5));
483 return(result);
484 }
485
486 /*
487 !!
488 !! Function D_F(s,p,tt,a,b,gm)
489 !!
490 */
491 Real dff(Real s, Real p,Real tt,Real a,Real b,Real gm)
492 {
493 Real result;
494 Real aa,caux,caux1,caux2;
495 Real xx,yy,rho;
496
497 xx=(a-b*(tt-p))/std::pow(2.0*(tt-p),0.5);
498 caux=-PHID(xx)*exp(-0.5*a*b);
499
500 xx=(a+b*(tt-p))/std::pow(2.0*(tt-p),0.5);
501 yy=(a+b*(tt-s))/std::pow(2.0*(tt-s),0.5);
502 rho=std::pow((tt-p)/(tt-s),0.5);
503 caux1=ND2(-xx,-yy,rho);
504 caux1=exp(0.5*a*b)*caux1;
505
506 xx=(a-b*(tt-p))/std::pow(2.0*(tt-p),0.5);
507 yy=(a-b*(tt-s))/std::pow(2.0*(tt-s),0.5);
508 rho=std::pow((tt-p)/(tt-s),0.5);
509 caux2=ND2(-xx,-yy,rho);
510 caux2=exp(-0.5*a*b)*caux2;
511
512 aa=exp((b*b-(1.0-gm)*(1.0-gm))*(tt-s)/4.0);
513
514 result=(caux+caux1+caux2)*aa;
515 return(result);
516 }
517
518
519 /*
520 !!
521 !! Function D_L(s,p,a,b,c,gm)
522 !!
523 */
524 Real dll(Real s,Real p,Real tt,Real a,Real b,Real c,Real gm)
525 {
526 Real result;
527 Real aa,caux,caux1;
528 Real sigmarho[4],limit[4],epsi;
529
530 epsi=1.e-12;
531 limit[1]=(a+b*(tt-p))/std::pow(2.0*(tt-p),0.5);
532 limit[2]=(a+b*(tt-s))/std::pow(2.0*(tt-s),0.5);
533 limit[3]=(a+b*tt+c)/std::pow(2.0*tt,0.5);
534 sigmarho[1]=std::pow((tt-p)/(tt-s),0.5);
535 sigmarho[2]=std::pow((tt-p)/tt,0.5);
536 sigmarho[3]=std::pow((tt-s)/tt,0.5);
537
538 caux=exp(0.5*a*b)*tvtl(0,limit,sigmarho,epsi);
539
540 limit[1]=(a-b*(tt-p))/std::pow(2.0*(tt-p),0.5);
541 limit[2]=(-a+b*(tt-s))/std::pow(2.0*(tt-s),0.5);
542 limit[3]=(-a+b*tt+c)/std::pow(2.0*tt,0.5);
543 sigmarho[1]=-std::pow((tt-p)/(tt-s),0.5);
544 sigmarho[2]=-std::pow((tt-p)/tt,0.5);
545 sigmarho[3]=std::pow((tt-s)/tt,0.5);
546
547 caux1=-exp(-0.5*a*b)*tvtl(0,limit,sigmarho,epsi);
548
549 aa=exp((b*b-(1.0-gm)*(1.0-gm))*(tt-s)/4.0);
550
551
552 result=(caux+caux1)*aa;
553
554 return(result);
555 }
556
557 /*
558 !!
559 !! Derivative with respect to a of the function D_F(s,p,tt,a,b,gm)
560 !!
561 */
562 Real ddff(Real s, Real p,Real tt,Real a,Real b,Real gm)
563 {
564 Real aa,caux,caux1,caux2,caux3,caux4;
565 Real xx,yy,rho;
566 Real result;
567 double ppi= 3.14159265358979324;
568
569 xx=(a-b*(tt-p))/std::pow(2.0*(tt-p),0.5);
570 caux=PHID(xx)*exp(-0.5*a*b);
571
572 xx=(a+b*(tt-p))/std::pow(2.0*(tt-p),0.5);
573 yy=(a+b*(tt-s))/std::pow(2.0*(tt-s),0.5);
574 rho=std::pow((tt-p)/(tt-s),0.5);
575 caux1=ND2(-xx,-yy,rho);
576 caux1=exp(0.5*a*b)*caux1;
577
578 xx=(a-b*(tt-p))/std::pow(2.0*(tt-p),0.5);
579 yy=(a-b*(tt-s))/std::pow(2.0*(tt-s),0.5);
580 rho=std::pow((tt-p)/(tt-s),0.5);
581 caux2=ND2(-xx,-yy,rho);
582 caux2=-exp(-0.5*a*b)*caux2;
583
584 caux=0.5*b*(caux+caux1+caux2);
585
586 xx=(a+b*(tt-p))/std::pow(2.0*(tt-p),0.5);
587 yy=b*std::pow((p-s),0.5)/std::pow(2.0,0.5);
588 caux1=exp(-0.5*xx*xx)*exp(0.5*a*b)*PHID(yy)/(2.0*std::pow(ppi*(tt-p),0.5));
589
590
591 xx=(a+b*(tt-s))/std::pow(2.0*(tt-s),0.5);
592 yy=a*std::pow((p-s),0.5)/std::pow(2.0*(tt-p)*(tt-s),0.5);
593 caux2=exp(-0.5*xx*xx)*exp(0.5*a*b)*PHID(yy)/(2.0*std::pow(ppi*(tt-s),0.5));
594
595 xx=(a-b*(tt-p))/std::pow(2.0*(tt-p),0.5);
596 yy=b*std::pow((p-s),0.5)/std::pow(2.0,0.5);
597 caux3=-exp(-0.5*xx*xx)*exp(-0.5*a*b)*PHID(yy)/(2.0*std::pow(ppi*(tt-p),0.5));
598
599
600 xx=(a-b*(tt-s))/std::pow(2.0*(tt-s),0.5);
601 yy=a*std::pow((p-s),0.5)/std::pow(2.0*(tt-p)*(tt-s),0.5);
602 caux4=exp(-0.5*xx*xx)*exp(-0.5*a*b)*PHID(yy)/(2.0*std::pow(ppi*(tt-s),0.5));
603
604
605
606 aa=exp((b*b-(1.0-gm)*(1.0-gm))*(tt-p)/4.0);
607
608
609 result=(caux+caux1+caux2+caux3+caux4)*aa;
610 return(result);
611 }
612
613 /*
614 !!
615 !! Derivative with respect to a of the function D_L(s,p,tt,a,b,c,gm)
616 !!
617 */
618 Real ddll(Real s,Real p,Real tt, Real ax, Real bx,Real c, Real gm)
619 {
620 static Real result;
621 static Real aa,caux,caux1;
622 static Real sigmarho[4],limit[4];
623 static int idx;
624 Real epsi;
625
626 epsi=1.e-12;
627 limit[1]=(ax+bx*(tt-p))/std::pow(2.0*(tt-p),0.5);
628 limit[2]=(ax+bx*(tt-s))/std::pow(2.0*(tt-s),0.5);
629 limit[3]=(ax+bx*tt+c)/std::pow(2.0*tt,0.5);
630 sigmarho[1]=std::pow((tt-p)/(tt-s),0.5);
631 sigmarho[2]=std::pow((tt-p)/tt,0.5);
632 sigmarho[3]=std::pow((tt-s)/tt,0.5);
633
634 caux=0.5*bx*tvtl(0,limit,sigmarho,epsi);
635
636
637 idx=1;
638 caux=caux+derivn3(limit,sigmarho,idx)/std::pow(2.0*(tt-p),0.5);
639
640 idx=2;
641 caux=caux+derivn3(limit,sigmarho,idx)/std::pow(2.0*(tt-s),0.5);
642
643 idx=3;
644 caux=caux+derivn3(limit,sigmarho,idx)/std::pow(2.0*tt,0.5);
645
646 caux=exp(0.5*ax*bx)*caux;
647
648 limit[1]=(ax-bx*(tt-p))/std::pow(2.0*(tt-p),0.5);
649 limit[2]=(-ax+bx*(tt-s))/std::pow(2.0*(tt-s),0.5);
650 limit[3]=(-ax+bx*tt+c)/std::pow(2.0*tt,0.5);
651 sigmarho[1]=-std::pow((tt-p)/(tt-s),0.5);
652 sigmarho[2]=-std::pow((tt-p)/tt,0.5);
653 sigmarho[3]=std::pow((tt-s)/tt,0.5);
654
655 caux1=0.5*bx*tvtl(0,limit,sigmarho,epsi);
656
657 idx=1;
658 caux1=caux1-derivn3(limit,sigmarho,idx)/std::pow(2.0*(tt-p),0.5);
659
660 idx=2;
661 caux1=caux1+derivn3(limit,sigmarho,idx)/std::pow(2.0*(tt-s),0.5);
662
663
664 idx=3;
665 caux1=caux1+derivn3(limit,sigmarho,idx)/std::pow(2.0*tt,0.5);
666
667 caux1=exp(-0.5*ax*bx)*caux1;
668
669
670 aa=exp((bx*bx-(1.0-gm)*(1.0-gm))*(tt-s)/4.0);
671
672 result=(caux+caux1)*aa;
673 return(result);
674 }
675
676 /*
677 !!
678 !! Derivative with respect to a of the function D_E(s,p,tt,a,b,gm)
679 !!
680 */
681 Real ddvv(Real s, Real p, Real tt, Real a, Real b, Real gm)
682 {
683 static Real result;
684 static Real aa,caux,caux1,caux2,caux6;
685 static Real caux3,caux4,caux5,aux;
686 static Real xx,yy,rho;
687 static double ppi= 3.14159265358979324;
688
689 aa=(a*p+b*(tt-p))/std::pow(2.0*p*tt*(tt-p),0.5);
690 caux=PHID(aa);
691
692 aa=exp(-(a-b)*(a-b)/(4.0*tt))/tt;
693
694 caux=0.5*aa*caux*(a-b);
695
696 xx=(a*p+b*(tt-p))/std::pow(2.0*tt*p*(tt-p),0.5);
697 yy=(a*s+b*(tt-s))/std::pow(2.0*tt*s*(tt-s),0.5);
698 rho=std::pow((s*(tt-p))/(p*(tt-s)),0.5);
699 caux1=ND2(-xx,-yy,rho);
700 caux1=-0.5*aa*caux1*(a-b);
701
702
703 aa=exp(-(a+b)*(a+b)/(4.0*tt))/tt;
704
705 xx=(a*p-b*(tt-p))/std::pow(2.0*tt*p*(tt-p),0.5);
706 yy=(a*s-b*(tt-s))/std::pow(2.0*tt*s*(tt-s),0.5);
707 rho=std::pow((s*(tt-p))/(p*(tt-s)),0.5);
708 caux2=ND2(-xx,-yy,rho);
709 caux2=-0.5*aa*caux2*(a+b);
710
711 aa=-b*std::pow((p-s)/std::pow(2.0*p*s,0.5),0.5);
712 aux=std::pow(p/(ppi*tt*(tt-p)),0.5)*PHID(aa);
713
714 xx=(a+b)*(a+b)/(4.0*tt);
715 yy=std::pow((a*p-b*(tt-p)),2)/(4.0*p*tt*(tt-p));
716 caux3=aux*exp(-xx)*exp(-yy)/2.0;
717
718
719 xx=(a-b)*(a-b)/(4.0*tt);
720 yy=std::pow((a*p+b*(tt-p)),2)/(4.0*p*tt*(tt-p));
721 caux4=aux*exp(-xx)*exp(-yy)/2.0;
722
723 aa=a*std::pow((p-s)/std::pow(2.0*(tt-p)*(tt-s),0.5),0.5);
724 aux=std::pow(s/(ppi*tt*(tt-s)),0.5)*PHID(aa);
725
726 xx=(a+b)*(a+b)/(4.0*tt);
727 yy=std::pow((a*s-b*(tt-s)),2)/(4.0*s*tt*(tt-s));
728 caux5=aux*exp(-xx)*exp(-yy)/2.0;
729
730 xx=(a-b)*(a-b)/(4.0*tt);
731 yy=std::pow((a*s+b*(tt-s)),2)/(4.0*s*tt*(tt-s));
732 caux6=aux*exp(-xx)*exp(-yy)/2.0;
733
734 aux=exp((1.0-gm)*(1.0-gm)*tt/4.0)*std::pow(tt,0.5);
735
736 result=(caux+caux1+caux2+caux3+caux4+caux5+caux6)/(aux*2.0*std::pow(ppi,0.5));
737 return(result);
738 }
739
740 /*
741 !!
742 !! Derivn3 computes the derivatives of the trivariate cumulative normal
743 !! distribution with respect to one of the integration limits
744 !!
745 */
746 Real derivn3(Real limit[4],Real sigmarho[4], int idx)
747 {
748 static Real aa;
749 static Real xx,yy,rho,sc;
750 static double ppi= 3.14159265358979324;
751 static Real deriv;
752 sc=std::pow(2.0*ppi,0.5);
753
754 if(idx==1)
755 {
756 aa=exp(-0.5*std::pow(limit[1],2));
757 xx=(limit[3]-sigmarho[2]*limit[1])/std::pow((1.0-std::pow(sigmarho[2],2)),0.5);
758 yy=(limit[2]-sigmarho[1]*limit[1])/std::pow((1.0-std::pow(sigmarho[1],2)),0.5);
759 rho=(sigmarho[3]-sigmarho[1]*sigmarho[2])/std::pow((1.0-sigmarho[1]*sigmarho[1])*(1.0-sigmarho[2]*sigmarho[2]),0.5);
760 deriv=aa*ND2(-xx,-yy,rho)/sc;
761 }
762 else
763 {
764 if(idx==2)
765 {
766 aa=exp(-0.5*limit[2]*limit[2]);
767 xx=(limit[1]-sigmarho[1]*limit[2])/std::pow((1.0-std::pow(sigmarho[1],2)),0.5);
768 yy=(limit[3]-sigmarho[3]*limit[2])/std::pow((1.0-std::pow(sigmarho[3],2)),0.5);
769 rho=(sigmarho[2]-sigmarho[1]*sigmarho[3])/ \
770 std::pow((1.0-sigmarho[1]*sigmarho[1])*(1.0-sigmarho[3]*sigmarho[3]),0.5);
771 deriv=aa*ND2(-xx,-yy,rho)/sc;
772 }
773 else
774 {
775 //!!! idx=3
776 aa=exp(-0.5*limit[3]*limit[3]);
777
778 xx=(limit[1]-sigmarho[2]*limit[3])/std::pow((1.0-std::pow(sigmarho[2],2)),0.5);
779 yy=(limit[2]-sigmarho[3]*limit[3])/std::pow((1.0-std::pow(sigmarho[3],2)),0.5);
780 rho=(sigmarho[1]-sigmarho[2]*sigmarho[3])/ \
781 std::pow((1.0-sigmarho[2]*sigmarho[2])*(1.0-sigmarho[3]*sigmarho[3]),0.5);
782 deriv=aa*ND2(-xx,-yy,rho)/sc;
783 }
784
785 }
786 return(deriv);
787 }
788
789
790 Real BVTL(int NU, Real DH, Real DK, Real RRR );
791 Real TVTMFN(Real X, Real H1, Real H2, Real H3,
792 Real R23, Real RUA, Real RUB, Real AR,
793 Real RUC, int NUC);
794 Real ADONET(Real ZRO,Real ONE,Real EPS,
795 Real(*TVTMFN)(Real X, Real H1, Real H2,
796 Real H3, Real R23, Real RUA,
797 Real RUB, Real AR, Real RUC,
798 int NUC));
799
800 Real tvtl(int NU, const Real limit[4], const Real sigmarho[4], Real epsi) {
801 /*
802 A function for computing trivariate normal and t-probabilities.
803
804 This function uses algorithms developed from the ideas
805 described in the papers:
806 R.L. Plackett, Biometrika 41(1954), pp. 351-360.
807 Z. Drezner, Math. Comp. 62(1994), pp. 289-294.
808 and uses adaptive integration.
809
810 The software given here is based on algorithms described in
811 the paper A. Genz: "Numerical Computation of Rectangular
812 Bivariate and Trivariate Normal and t Probabilities",
813 Statistics and Computing 14 (2004) 251-260.
814
815 This software has been developed by M.C. Recchioni based on
816 previous software developed by
817 Alan Genz
818 Department of Mathematics
819 Washington State University
820 Pullman, WA 99164-3113
821 Email : alangenz@wsu.edu
822 The software developed by A. Genz is available free of
823 charge in the website:
824 www.math.wsu.edu/faculty/genz/software/software.html
825
826 The software calculates the probability that
827 X(I) < H(I), for I = 1,2,3
828
829 NU INTEGER degrees of freedom; use NU = 0 for normal cases.
830 LIMIT REAL array of uppoer limits for probability distribution
831 SIGMARHO REAL array of three correlation coefficients, should
832 contain the lower left portion of the correlation matrix.
833 SIGMARHO should contains the values r21, r31, r23 in that order.
834 EPSI REAL required absolute accuracy; maximum accuracy for most
835 computations is approximately 1D-14
836
837 */
838
839 static Real result;
840 static Real ONE=1.0, ZRO=0.0, EPS, TVT;
841 static Real PT, R12, R13;
842 static double ppi= 3.14159265358979324;
843 EPS = max( 1.e-14, epsi );
844 PT=ppi/2.0;
845
846 NUC = NU;
847 H1 = limit[1];
848 H2 = limit[2];
849 H3 = limit[3];
850 R12 = sigmarho[1];
851 R13 = sigmarho[2];
852 R23 = sigmarho[3];
853 /*
854 * Sort R's and check for special cases
855 */
856 if ( fabs(R12) > fabs(R13) ) {
857 H2 = H3;
858 H3 = limit[2];
859 R12 = R13;
860 R13 = sigmarho[1];
861 }
862
863 if ( fabs(R13) > fabs(R23) ) {
864 H1 = H2;
865 H2 = limit[1];
866 R23 = R13;
867 R13 = sigmarho[3];
868 }
869
870 TVT = 0.0;
871 if ( (fabs(H1) + fabs(H2) + fabs(H3)) < EPS ) TVT = ( 1 + ( asin(R12) + asin(R13) + asin(R23) )/PT )/8.0;
872
873 else if ( (NU < 1) && ( (fabs(R12) + fabs(R13)) < EPS) ) TVT = PHID(H1)*BVTL( NU, H2, H3, R23 );
874
875 else if ( (NU < 1) && ((fabs(R13) + fabs(R23))< EPS) ) TVT = PHID(H3)*BVTL( NU, H1, H2, R12 );
876
877 else if( (NU < 1) && ((fabs(R12) + fabs(R23))< EPS) ) TVT = PHID(H2)*BVTL( NU, H1, H3, R13 );
878
879 else if ( (1.0 - R23)< EPS ) TVT = BVTL( NU, H1, min( H2, H3 ), R12 );
880
881 else if ( (R23 + 1.0) <EPS ) {
882 if ( H2 > -H3 ) TVT = BVTL( NU, H1, H2, R12 ) - BVTL( NU, H1, -H3, R12 );}
883 else
884 {
885 /*
886 * Compute singular TVT value
887 */
888 if ( NU < 1 ) TVT = BVTL( NU, H2, H3, R23 )*PHID(H1);
889
890 else if ( R23 > 0 ) TVT = BVTL( NU, H1, min( H2, H3 ), ZRO );
891
892 else if ( H2 > -H3 ) TVT = BVTL( NU, H1, H2, ZRO ) - BVTL( NU, H1, -H3, ZRO );
893
894 /*
895 * Use numerical integration to compute probability
896 *
897 */
898 RUA = asin( R12 );
899 RUB = asin( R13 );
900 AR = asin( R23);
901 RUC = SIGN( PT, AR ) - AR;
902 TVT = TVT + ADONET(ZRO, ONE, EPS, TVTMFN) / (4.0 * PT);
903 }
904 result = max( ZRO, min( TVT, ONE ) );
905
906 return(result);
907 }
908
909 void SINCS(Real v1,Real& v2, Real& v3);
910 Real PNTGND(int , Real ,Real ,Real ,
911 Real ,Real ,Real ,Real );
912
913 Real TVTMFN(Real X, Real H1, Real H2, Real H3, Real R23,
914 Real RUA, Real RUB, Real AR,Real RUC, int NUC ){
915 /*
916 Computes Plackett formula integrands
917 */
918
919 static Real R12=0.0, RR2=0, R13=0.0, RR3=0.0, R=0.0, RR=0.0;
920 const Real ZRO = 0.0;
921 Real result = 0.0;
922
923 SINCS( RUA*X, R12, RR2 );
924 SINCS( RUB*X, R13, RR3 );
925
926 if ( fabs(RUA)> 0 ) result += RUA*PNTGND( NUC, H1,H2,H3, R13,R23,R12,RR2);
927 if( fabs(RUB)>0 ) result += RUB*PNTGND( NUC, H1,H3,H2, R12,R23,R13,RR3 ) ;
928 if ( NUC > 0 )
929 {
930 SINCS( AR + RUC*X, R, RR );
931 result -= RUC*PNTGND( NUC, H2, H3, H1, ZRO, ZRO, R, RR );
932 }
933 return(result);
934 }
935 //
936
937
938 void SINCS(Real X, Real& SX, Real& CS )
939 {
940 /*
941 Computes SIN(X), COS(X)^2, with series approx. for |X| near PI/2
942 */
943 static Real PT, EE;
944 PT = 1.57079632679489661923132169163975;
945 EE = std::pow(( PT - fabs(X) ),2);
946
947 if ( EE < 5e-5 )
948 {
949 SX = SIGN( 1 - EE*( 1 - EE/12 )/2, X );
950 CS = EE*( 1 - EE*( 1 - 2*EE/15 )/3 );
951 }
952 else
953 {
954 SX = sin(X);
955 CS = 1 - SX*SX;
956 }
957 }
958 //
959
960 Real KRNRDT(Real, Real,
961 Real(*TVTMFN)(Real, Real, Real, Real,
962 Real, Real, Real, Real, Real, int),
963 Real& );
964
965 Real ADONET(Real A,Real B, Real TOL,Real(*TVTMFN)(Real X, Real H1, Real H2, Real H3, Real R23, Real RUA, Real RUB, Real AR, Real RUC, int NUC)){
966 //
967 // One Dimensional Globally Adaptive Integration Function
968 //
969 int NL=100, I, IM, IP;
970 static Real EI[101], AI[101], BI[101], FI[101], FIN;
971 static Real result,ERR;
972
973
974 AI[1] = A;
975 BI[1] = B;
976 ERR = 1;
977 IP = 1;
978 IM = 1;
979 while ( ((4*ERR)> TOL) && (IM< NL) )
980 {
981 IM = IM + 1;
982 BI[IM] = BI[IP];
983 AI[IM] = (AI[IP] + BI[IP] )/2.0;
984 BI[IP] = AI[IM];
985 FI[IP] = KRNRDT( AI[IP], BI[IP], *TVTMFN, EI[IP] );
986 FI[IM] = KRNRDT( AI[IM], BI[IM], *TVTMFN, EI[IM] );
987
988 ERR = 0.0;
989 FIN = 0.0;
990 for(I = 1; I<=IM; I++)
991 {
992 if( EI[I] > EI[IP]) IP = I;
993 FIN = FIN + FI[I];
994 ERR = ERR + EI[I]*EI[I];
995 }
996 ERR = std::pow( ERR,0.5 );
997 }
998 result=FIN;
999 // ADONET = FIN
1000 return(result);
1001 }
1002 //
1003
1004 Real KRNRDT(Real A, Real B,Real(*TVTMFN)(Real X, Real H1, Real H2, Real H3, Real R23, Real RUA, Real RUB, Real AR, Real RUC, int NUC),Real& ERR ){
1005
1006 //
1007 // Kronrod Rule
1008 //
1009 static Real T, CEN, FC, WID, RESG, RESK;
1010
1011 static Real result;
1012 //
1013 // The abscissae and weights are given for the interval (-1,1);
1014 // only positive abscissae and corresponding weights are given.
1015 //
1016 // XGK - abscissae of the 2N+1-point Kronrod rule:
1017 // XGK(2), XGK(4), ... N-point Gauss rule abscissae;
1018 // XGK(1), XGK(3), ... optimally added abscissae.
1019 // WGK - weights of the 2N+1-point Kronrod rule.
1020 // WG - weights of the N-point Gauss rule.
1021 //
1022 int J, N=11;
1023
1024 static Real WG[7], WGK[13], XGK[13];
1025
1026 WG[1]= 0.2729250867779007;
1027 WG[2]=0.05566856711617449;
1028 WG[3]=0.1255803694649048;
1029 WG[4]=0.1862902109277352;
1030 WG[5]= 0.2331937645919914;
1031 WG[6]= 0.2628045445102478;
1032 //
1033 XGK[1]= 0.0000000000000000;
1034 XGK[2]= 0.9963696138895427;
1035 XGK[3]= 0.9782286581460570;
1036 XGK[4]= 0.9416771085780681;
1037 XGK[5]= 0.8870625997680953;
1038 XGK[6]= 0.8160574566562211;
1039 XGK[7]= 0.7301520055740492;
1040 XGK[8]= 0.6305995201619651;
1041 XGK[9]= 0.5190961292068118;
1042 XGK[10]= 0.3979441409523776;
1043 XGK[11]= 0.2695431559523450;
1044 XGK[12]= 0.1361130007993617;
1045 //
1046 WGK[1]=0.1365777947111183;
1047 WGK[2]=0.9765441045961290e-02;
1048 WGK[3]=0.2715655468210443e-01;
1049 WGK[4]=0.4582937856442671e-01;
1050 WGK[5]=0.6309742475037484e-01;
1051 WGK[6]=0.7866457193222764e-01;
1052 WGK[7]=0.9295309859690074e-01;
1053 WGK[8]=0.1058720744813894;
1054 WGK[9]=0.1167395024610472;
1055 WGK[10]=0.1251587991003195;
1056 WGK[11]=0.1312806842298057;
1057 WGK[12]=0.1351935727998845;
1058 /*
1059 Major variables
1060
1061 CEN - mid point of the interval
1062 WID - half-length of the interval
1063 RESG - result of the N-point Gauss formula
1064 RESK - result of the 2N+1-point Kronrod formula
1065 Compute the 2N+1-point Kronrod approximation to
1066 the integral, and estimate the absolute error.
1067 */
1068 WID = ( B - A )/2.0;
1069 CEN = ( B + A )/2.0;
1070
1071 FC = TVTMFN(CEN,H1, H2, H3, R23, RUA, RUB, AR, RUC, NUC);
1072
1073 RESG = FC*WG[0+1];
1074 RESK = FC*WGK[0+1];
1075
1076 for (J = 1; J<= N; J++)
1077 {
1078 T = WID*XGK[J+1];
1079 FC = TVTMFN(CEN-T,H1, H2, H3, R23, RUA, RUB, AR, RUC, NUC )+TVTMFN(CEN+T,H1, H2, H3, R23, RUA, RUB, AR, RUC, NUC );
1080 RESK = RESK + WGK[J+1]*FC;
1081 if((J-2*int(J/2)) == 0 ) RESG = RESG + WG[1+J/2]*FC;
1082 }
1083 result = WID*RESK;
1084 ERR = fabs( WID*( RESK - RESG ) );
1085 return(result);
1086 }
1087
1088 //
1089 Real STUDNT(int NU, Real T )
1090 {
1091 /*
1092 Student t Distribution Function
1093 */
1094 static int J;
1095 static Real ZRO=0.0, ONE=1.0;
1096 static Real CSSTHE, SNTHE, POLYN, TT, TS, RN;
1097 static Real result;
1098
1099
1100 if ( NU < 1 ) result= PHID( T );
1101 else if ( NU == 1 ) result = ( 1 + 2.0*atan(T)/PI )/2.0;
1102 else if ( NU == 2 ) result = ( 1 + T/std::pow(( 2.0 + T*T),0.5))/2.0;
1103 else
1104 {
1105 TT = T*T;
1106 CSSTHE = 1/( 1 + TT/double(NU) );
1107 POLYN = 1;
1108 for( J = NU-2; J>= 2; J=J-2)
1109 {
1110 POLYN = 1.0 + ( J - 1.0 )*CSSTHE*POLYN/(double)J;
1111 }
1112 if ((NU-2*int(NU/2) ) == 1 )
1113 {
1114 RN = NU;
1115 TS = T/std::pow(RN,0.5);
1116 result = ( 1.0 + 2.0*( atan(TS) + TS*CSSTHE*POLYN )/PI )/2.0;
1117 }
1118 else
1119 {
1120 SNTHE = T/std::pow(( NU + TT ),0.5);
1121 result = ( 1 + SNTHE*POLYN )/2.0;
1122 }
1123 result = max( ZRO, min( result, ONE ) );
1124 }
1125 return(result);
1126 }
1127
1128 //
1129 Real BVTL(int NU, Real DH, Real DK, Real R )
1130 {
1131 /*
1132 A function for computing bivariate t probabilities.
1133 This function is based on the method described by
1134 Dunnett, C.W. and M. Sobel, (1954),
1135 A bivariate generalization of Student's t-distribution
1136 with tables for certain special cases,
1137 Biometrika 41, pp. 153-169.
1138 The software given here has been developed by M.C. Recchioni based on previous
1139 software developed by
1140 Alan Genz
1141 Department of Mathematics
1142 Washington State University
1143 Pullman, WA 99164-3113
1144 Email : alangenz@wsu.edu
1145 The software developed by A. Genz is available free of charge in
1146 the website: www.math.wsu.edu/faculty/genz/software/software.html
1147 ***
1148
1149 BVTL - calculate the probability that X < DH and Y < DK.
1150
1151 parameters
1152
1153 NU number of degrees of freedom
1154 DH 1st lower integration limit
1155 DK 2nd lower integration limit
1156 R correlation coefficient
1157 */
1158 static int J, HS, KS;
1159 static Real TPI, ORS, HRK, KRH, BVT, SNU;
1160 static Real GMPH, GMPK, XNKH, XNHK, QHRK, HKN, HPK, HKRN;
1161 static Real BTNCKH, BTNCHK, BTPDKH, BTPDHK, ONE, EPS;
1162 static Real result;
1163 ONE = 1;
1164 EPS = 1e-15;
1165 if ( NU <1 ) result = ND2( -DH, -DK, R );
1166
1167 else if ( (1 - R)<= EPS ) result = STUDNT( NU, min( DH, DK ) );
1168
1169 else if( (R + 1)<=EPS )
1170 {
1171 if( DH > -DK ) result = STUDNT( NU, DH ) - STUDNT( NU, -DK );
1172 else
1173 result = 0.0;
1174 }
1175 else
1176 {
1177 TPI = 2.0*PI;
1178 SNU = (double)NU;
1179 SNU = std::pow(SNU,0.5);
1180 ORS = 1.0 - R*R;
1181 HRK = DH - R*DK;
1182 KRH = DK - R*DH;
1183 if((fabs(HRK) + ORS)> 0 )
1184 {
1185 XNHK = HRK*HRK/( HRK*HRK + ORS*( NU + DK*DK ) );
1186 XNKH = KRH*KRH/( KRH*KRH+ ORS*( NU + DH*DH ) );
1187 }
1188 else
1189 {
1190 XNHK = 0.0;
1191 XNKH = 0.0;
1192 }
1193
1194 HS =(int)SIGN( ONE, DH - R*DK );
1195 KS =(int)SIGN( ONE, DK - R*DH );
1196 if((NU-2*(int)(NU/2))==0 )
1197 {
1198 BVT = atan2( std::pow(ORS,0.5), -R )/TPI;
1199 GMPH = DH/std::pow( 16*( NU + DH*DH ),0.5 );
1200 GMPK = DK/std::pow( 16*( NU + DK*DK),0.5);
1201 BTNCKH = 2*atan2( std::pow( XNKH,0.5 ), std::pow(( 1-XNKH),0.5) )/PI;
1202 BTPDKH = 2*std::pow( XNKH*( 1 - XNKH ),0.5 )/PI;
1203 BTNCHK = 2*atan2( std::pow( XNHK,0.5 ), std::pow((1 - XNHK),0.5) )/PI;
1204 BTPDHK = 2*std::pow( XNHK*( 1 - XNHK ),0.5 )/PI;
1205 for( J = 1; J<= NU/2;J++)
1206 {
1207 BVT = BVT + GMPH*( 1 + KS*BTNCKH );
1208 BVT = BVT + GMPK*( 1 + HS*BTNCHK );
1209 BTNCKH = BTNCKH + BTPDKH;
1210 BTPDKH = 2*J*BTPDKH*( 1 - XNKH )/( 2*J + 1 );
1211 BTNCHK = BTNCHK + BTPDHK;
1212 BTPDHK = 2*J*BTPDHK*( 1 - XNHK )/( 2*J + 1 );
1213 GMPH = GMPH*( 2*J - 1 )/( 2*J*( 1 + DH*DH/NU ) );
1214 GMPK = GMPK*( 2*J - 1 )/( 2*J*( 1 + DK*DK/NU ) );
1215 }
1216 }
1217 else
1218 {
1219 QHRK = std::pow((DH*DH + DK*DK - 2*R*DH*DK + NU*ORS),0.5 ) ;
1220 HKRN = DH*DK + R*NU ;
1221 HKN = DH*DK - NU;
1222 HPK = DH + DK;
1223 BVT = atan2( -SNU*( HKN*QHRK + HPK*HKRN ),HKN*HKRN-NU*HPK*QHRK )/TPI;
1224 if ( BVT < -EPS ) BVT = BVT + 1;
1225 GMPH = DH/( TPI*SNU*( 1 + DH*DH/NU ) );
1226 GMPK = DK/( TPI*SNU*( 1 + DK*DK/NU ) );
1227 BTNCKH = std::pow( XNKH,0.5 );
1228 BTPDKH = BTNCKH;
1229 BTNCHK = std::pow( XNHK,0.5 );
1230 BTPDHK = BTNCHK;
1231 for( J = 1;J<= ( NU - 1 )/2; J++)
1232 {
1233 BVT = BVT + GMPH*( 1 + KS*BTNCKH );
1234 BVT = BVT + GMPK*( 1 + HS*BTNCHK );
1235 BTPDKH = ( 2*J - 1 )*BTPDKH*( 1 - XNKH )/( 2*J );
1236 BTNCKH = BTNCKH + BTPDKH;
1237 BTPDHK = ( 2*J - 1 )*BTPDHK*( 1 - XNHK )/( 2*J );
1238 BTNCHK = BTNCHK + BTPDHK;
1239 GMPH = 2*J*GMPH/( ( 2*J + 1 )*( 1 + DH*DH/NU ) );
1240 GMPK = 2*J*GMPK/( ( 2*J + 1 )*( 1 + DK*DK/NU ) );
1241 }
1242 }
1243 result = BVT;
1244 }
1245
1246 return(result);
1247
1248 }
1249
1250
1251
1252
1253 Real PNTGND(int NUC, Real BA, Real BB, Real BC, Real RA, Real RB, Real R, Real RR) {
1254 /*
1255 Computes Plackett formula integrand
1256 */
1257 static Real DT, FT, BT,result;
1258
1259 result = 0.0;
1260 DT = RR*( RR - std::pow(( RA - RB ),2) - 2*RA*RB*( 1 - R ) );
1261 if( DT > 0 ) {
1262 BT = ( BC*RR + BA*( R*RB - RA ) + BB*( R*RA -RB ) )/std::pow(DT,0.5);
1263 FT = std::pow(( BA - R*BB ),0.5)/RR + BB*BB;
1264 if( NUC<1 ) {
1265 if ( (BT > -10) && (FT <100) ) {
1266 result = exp( -FT/2 );
1267 if ( BT <10 ) result= result*PHID(BT);
1268 } else {
1269 FT = std::pow((1 + FT/NUC),0.5);
1270 result = STUDNT( NUC, BT/FT )/std::pow(FT,NUC);
1271 }
1272 }
1273 }
1274 return(result);
1275 }
1276
1277
1278
1279 //***********************************************************
1280 Real ND2(Real a, Real b, Real rho ){
1281 /*
1282 * A function for computing bivariate normal probabilities.
1283 * This function is based on the method described by
1284 * Z. Drezner and G.O. Wesolowsky, (1989),
1285 * On the computation of the bivariate normal integral,
1286 * Journal of Statist. Comput. Simul. 35, pp. 101-107,
1287 * with major modifications for double precision, and for |R| close to 1.
1288 * The software given here has been developed by M.C. Recchioni based on previous
1289 * software developed by:
1290 * Alan Genz
1291 * Department of Mathematics
1292 * Washington State University
1293 * Pullman, WA 99164-3113
1294 * Email : alangenz@wsu.edu
1295 * The software developed by A. Genz is available free of charge in the website:
1296 * www.math.wsu.edu/faculty/genz/software/software.html
1297 *
1298 *
1299 * ND2 calculates the probability that X > DH and Y > DK.
1300 * Note: Prob( X < DH, Y < DK ) = ND2( -DH, -DK, R ).
1301 *
1302 * Parameters
1303 *
1304 * DH DOUBLE PRECISION, integration limit
1305 * DK DOUBLE PRECISION, integration limit
1306 * R DOUBLE PRECISION, correlation coefficient
1307 */
1308 static double TWOPI = 6.283185307179586;
1309 static Real result, DK, DH, R;
1310 static int I, IS, LG, NG;
1311
1312 static Real XL[11][4], WL[11][4], AS, AA, BB, C, D, RS, XS, BVN;
1313 static Real SN, ASR, H, K, BS, HS, HK;
1314 // Gauss Legendre Points and Weights, N = 6
1315 // DATA ( W(I,1), X(I,1), I = 1,3) /
1316 WL[1][1]=0.1713244923791705;
1317 XL[1][1]=-0.9324695142031522;
1318 WL[2][1]= 0.3607615730481384;
1319 XL[2][1]=-0.6612093864662647;
1320 WL[3][1]= 0.4679139345726904;
1321 XL[3][1]=-0.2386191860831970;
1322
1323 // Gauss Legendre Points and Weights, N = 12
1324 // DATA ( W(I,2), X(I,2), I = 1,6) /
1325 WL[1][2]=0.4717533638651177e-01;
1326 XL[1][2]=-0.9815606342467191;
1327 WL[2][2]= 0.1069393259953183;
1328 XL[2][2]=-0.9041172563704750;
1329 WL[3][2]= 0.1600783285433464;
1330 XL[3][2]=-0.7699026741943050;
1331 WL[4][2]= 0.2031674267230659;
1332 XL[4][2]=-0.5873179542866171;
1333 WL[5][2]= 0.2334925365383547;
1334 XL[5][2]=-0.3678314989981802;
1335 WL[6][2] = 0.2491470458134029;
1336 XL[6][2]=-0.1252334085114692;
1337
1338 // Gauss Legendre Points and Weights, N = 20
1339 // DATA ( W(I,3), X(I,3), I = 1, 10 ) /
1340 WL[1][3]=0.1761400713915212e-01;
1341 XL[1][3]=-0.9931285991850949;
1342 WL[2][3]=0.4060142980038694e-01;
1343 XL[2][3]=-0.9639719272779138;
1344 WL[3][3]=0.6267204833410906e-01;
1345 XL[3][3]=-0.9122344282513259;
1346 WL[4][3]=0.8327674157670475e-01;
1347 XL[4][3]=-0.8391169718222188;
1348 WL[5][3]=0.1019301198172404;
1349 XL[5][3]=-0.7463319064601508;
1350 WL[6][3]=0.1181945319615184;
1351 XL[6][3]=-0.6360536807265150;
1352 WL[7][3]=0.1316886384491766;
1353 XL[7][3]=-0.5108670019508271;
1354 WL[8][3]=0.1420961093183821;
1355 XL[8][3]=-0.3737060887154196;
1356 WL[9][3]=0.1491729864726037;
1357 XL[9][3]=-0.2277858511416451;
1358 WL[10][3]=0.1527533871307259;
1359 XL[10][3]=-0.7652652113349733e-01;
1360
1361 R=rho;
1362 DH=a;
1363 DK=b;
1364
1365 if( fabs(R) < 0.3 ) {
1366 NG = 1;
1367 LG = 3; }
1368 else if ( fabs(R) < 0.75 ) {
1369 NG = 2;
1370 LG = 6;}
1371 else{
1372 NG = 3;
1373 LG = 10;
1374 }
1375 H = DH;
1376 K = DK;
1377 HK = H*K;
1378
1379 BVN = 0.0 ;
1380 if( fabs(R) < 0.925 ) {
1381 if( fabs(R) > 0 ) {
1382 HS = ( H*H + K*K )/2;
1383 ASR = asin(R);
1384 for (I = 1;I<= LG; I++){
1385 for( IS = -1; IS<= 1; IS=IS+2){
1386 SN = sin( ASR*( IS*XL[I][NG] + 1 )/2 );
1387 BVN = BVN + WL[I][NG]*exp( ( SN*HK-HS )/( 1.0-SN*SN ) );
1388
1389 }
1390 }
1391 BVN = BVN*ASR/( 2*TWOPI );
1392
1393 }
1394
1395 BVN = BVN + PHID(-H)*PHID(-K);
1396
1397 }
1398 else
1399 {
1400 if ( R < 0 ) {
1401 K = -K;
1402 HK = -HK;
1403 }
1404
1405 if( fabs(R) <1 ) {
1406 AS = ( 1 - R )*( 1 + R );
1407 AA = std::pow(AS,0.5);
1408
1409 BS = std::pow(( H - K ),2);
1410 C = ( 4 - HK )/8 ;
1411 D = ( 12 - HK )/16;
1412 ASR = -( BS/AS + HK )/2;
1413 if( ASR > -100 ) BVN = AA*exp(ASR)*( 1 - C*( BS - AS )*( 1 - D*BS/5 )/3 + C*D*AS*AS/5 );
1414 if( -HK<100 ){
1415 BB = std::pow(BS,0.5);
1416 BVN = BVN - exp( -HK/2 )*std::pow(TWOPI,0.5)*PHID(-BB/AA)*BB*( 1 - C*BS*( 1 - D*BS/5 )/3 );
1417 }
1418 AA = AA/2 ;
1419 for (I = 1; I<= LG;I++){
1420 for( IS = -1; IS<=1; IS=IS+2){
1421 XS =std::pow( ( AA*( IS*XL[I][NG] + 1 ) ),2) ;
1422 RS = std::pow( (1 - XS),2 );
1423 ASR = -( BS/XS + HK )/2;
1424 if ( ASR > -100 ) {
1425
1426 BVN = BVN + AA*WL[I][NG]*exp( ASR )*(exp( -HK*( 1 - RS )/( 2*( 1 + RS ) ) )/RS- ( 1 + C*XS*( 1 + D*XS ) ) );
1427
1428
1429 }
1430 }
1431 }
1432 BVN = -BVN/TWOPI;
1433 }
1434 if ( R > 0 ) {
1435
1436 BVN = BVN + PHID( -max( H, K ) );
1437
1438 }
1439 else
1440 {
1441 BVN = -BVN;
1442 if( K > H ) BVN = BVN + PHID(K) - PHID(H);
1443 }
1444 }
1445
1446
1447 result=BVN;
1448
1449 return(result);
1450
1451 }
1452
1453 struct integr_adapter {
1454 ext::shared_ptr<YieldTermStructure> r;
1455 explicit integr_adapter(const ext::shared_ptr<GeneralizedBlackScholesProcess>& process)
1456 : r(*(process->riskFreeRate())) {}
1457 Real operator()(Real t1,Real t2) const {
1458 return r->forwardRate(t1,t2,Continuous) * (t2-t1);
1459 }
1460 };
1461
1462 struct integalpha_adapter {
1463 ext::shared_ptr<YieldTermStructure> r;
1464 ext::shared_ptr<YieldTermStructure> q;
1465 explicit integalpha_adapter(
1466 const ext::shared_ptr<GeneralizedBlackScholesProcess>& process)
1467 : r(*(process->riskFreeRate())), q(*(process->dividendYield())) {}
1468 Real operator()(Real t1,Real t2) const {
1469 Real alpha = r->forwardRate(t1,t2,Continuous).rate()
1470 - q->forwardRate(t1,t2,Continuous).rate();
1471 return alpha * (t2-t1);
1472 }
1473 };
1474
1475 struct alpha_adapter {
1476 ext::shared_ptr<YieldTermStructure> r;
1477 ext::shared_ptr<YieldTermStructure> q;
1478 explicit alpha_adapter(const ext::shared_ptr<GeneralizedBlackScholesProcess>& process)
1479 : r(*(process->riskFreeRate())), q(*(process->dividendYield())) {}
1480 Real operator()(Real t) const {
1481 return r->forwardRate(t,t,Continuous).rate()
1482 - q->forwardRate(t,t,Continuous).rate();
1483 }
1484 };
1485
1486 struct sigmaq_adapter {
1487 ext::shared_ptr<BlackVolTermStructure> v;
1488 Real s;
1489 explicit sigmaq_adapter(const ext::shared_ptr<GeneralizedBlackScholesProcess>& process)
1490 : v(*(process->blackVolatility())), s(process->x0()) {}
1491 Real operator()(Real t) const {
1492 Real sigma = v->blackForwardVol(t,t,s,true);
1493 return sigma*sigma;
1494 }
1495 };
1496
1497 struct integs_adapter {
1498 ext::shared_ptr<BlackVolTermStructure> v;
1499 Real s;
1500 explicit integs_adapter(const ext::shared_ptr<GeneralizedBlackScholesProcess>& process)
1501 : v(*(process->blackVolatility())), s(process->x0()) {}
1502 Real operator()(Real t1,Real t2) const {
1503 return v->blackForwardVariance(t1,t2,s,true);
1504 }
1505 };
1506
1507 }
1508
1509
1510 PerturbativeBarrierOptionEngine::PerturbativeBarrierOptionEngine(
1511 ext::shared_ptr<GeneralizedBlackScholesProcess> process, Natural order, bool zeroGamma)
1512 : process_(std::move(process)), order_(order), zeroGamma_(zeroGamma) {
1513 registerWith(process_);
1514 }
1515
1516 void PerturbativeBarrierOptionEngine::calculate() const {
1517
1518 QL_REQUIRE(arguments_.barrierType == Barrier::UpOut,
1519 "this engine only manages up-and-out options");
1520
1521 QL_REQUIRE(arguments_.rebate == 0.0,
1522 "this engine does not manage non-null rebates");
1523
1524 ext::shared_ptr<PlainVanillaPayoff> payoff =
1525 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);
1526 QL_REQUIRE(payoff && payoff->optionType() == Option::Put,
1527 "this engine only manages put options");
1528
1529 Real stock = process_->x0();
1530 Real kprice = payoff->strike();
1531 Real hbarr = arguments_.barrier;
1532
1533 Time tauMin = 0.0;
1534 Time tauMax = process_->time(arguments_.exercise->lastDate());
1535
1536 QL_REQUIRE(order_ <= 2, "order must be <= 2");
1537
1538 int igm = zeroGamma_ ? 0 : 1;
1539
1540 results_.value = BarrierUPD(kprice, stock, hbarr,
1541 tauMin, tauMax, order_, igm,
1542 integr_adapter(process_),
1543 integalpha_adapter(process_),
1544 integs_adapter(process_),
1545 alpha_adapter(process_),
1546 sigmaq_adapter(process_));
1547 }
1548
1549}
1550
QuantLib::Observer::registerWith
std::pair< iterator, bool > registerWith(const ext::shared_ptr< Observable > &)
Definition: observable.hpp:228
QuantLib::PerturbativeBarrierOptionEngine::PerturbativeBarrierOptionEngine
PerturbativeBarrierOptionEngine(ext::shared_ptr< GeneralizedBlackScholesProcess >, Natural order=1, bool zeroGamma=false)
Definition: perturbativebarrieroptionengine.cpp:1510
QuantLib::PerturbativeBarrierOptionEngine::process_
ext::shared_ptr< GeneralizedBlackScholesProcess > process_
Definition: perturbativebarrieroptionengine.hpp:48
t
const DefaultType & t
Definition: defaultprobabilitykey.cpp:39
errors.hpp
Classes and functions for error handling.
QL_REQUIRE
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
exercise.hpp
Option exercise classes and payoff function.
b
ext::function< Real(Real)> b
Definition: extendedornsteinuhlenbeckprocess.cpp:30
functional.hpp
Maps function, bind and cref to either the boost or std implementation.
QuantLib::Time
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QuantLib::Real
QL_REAL Real
real number
Definition: types.hpp:50
QuantLib::Natural
unsigned QL_INTEGER Natural
positive integer
Definition: types.hpp:43
v0
Real v0
Definition: hestonrndcalculator.cpp:36
rho
Real rho
Definition: hestonrndcalculator.cpp:36
sigma
Real sigma
Definition: hestonrndcalculator.cpp:36
payoff
ext::shared_ptr< QuantLib::Payoff > payoff
Definition: integralhestonvarianceoptionengine.cpp:350
order_
const Size order_
Definition: lsmbasissystem.cpp:51
QuantLib
Definition: any.hpp:35
std
STL namespace.
q
ext::shared_ptr< YieldTermStructure > q
Definition: perturbativebarrieroptionengine.cpp:1464
r
ext::shared_ptr< YieldTermStructure > r
Definition: perturbativebarrieroptionengine.cpp:1454
v
ext::shared_ptr< BlackVolTermStructure > v
Definition: perturbativebarrieroptionengine.cpp:1487
perturbativebarrieroptionengine.hpp
perturbative barrier-option engine
alpha
Real alpha
Definition: sabr.cpp:200
types.hpp
Custom types.
Generated by Doxygen 1.9.5