gaussianorthogonalpolynomial.c
Locked
gaussianorthogonalpolynomial.c
 
|
I compiled QuantLib with ICC for x64 target. Everything went well
except a glitch with Null class (I reported it in my previous e-mail) and some issues in gaussianorthogonalpolynomial.cpp. Current implementations of GaussJacobiPolynomial::alpha and GaussJacobiPolynomial::beta have several issues. 1) On x64 targets Size is 64bit and conversion to double may generate a warning. For some strange reason ICC refuses to compile this line Real denom = (2.0*i+alpha_+beta_)*(2.0*i+alpha_+beta_+2); ICC bails out with an internal error. I don't know why it happens (ICC bug?), but casting i from Size into int fixes the problem. 2) Current implementation checks if num and denom variables are zero and then tries to use l'Hospital rule. But denom==0 can be true only if alpha_=0 and beta_=0 and/or i=0. In that case num=0 and l'Hospital rule is never used. In anyway, direct conversion of floating point numbers into bool is a fishy business at best. We should consider using fuzzy comparison instead of direct comparison, something like this: //----------------------------------------------------------------------------------- Real GaussJacobiPolynomial::alpha(Size iIn) const { int i = static_cast<int>(iIn); Real num = beta_*beta_ - alpha_*alpha_; Real denom = (2.0*i+alpha_+beta_)*(2.0*i+alpha_+beta_+2); if (std::fabs(denom)<QL_EPSILON) { if (std::fabs(num)<QL_EPSILON) { QL_FAIL("can't compute a_k for jacobi integration\n"); } else { // l'Hospital num = 2*beta_; denom= 2*(2.0*i+alpha_+beta_+1); QL_ASSERT(std::fabs(num)>QL_EPSILON, "can't compute a_k for jacobi integration\n"); } } return num / denom; } Real GaussJacobiPolynomial::beta(Size iIn) const { int i = static_cast<int>(iIn); Real num = 4.0*i*(i+alpha_)*(i+beta_)*(i+alpha_+beta_); Real denom = (2.0*i+alpha_+beta_)*(2.0*i+alpha_+beta_) * ((2.0*i+alpha_+beta_)*(2.0*i+alpha_+beta_)-1); if (std::fabs(denom)<QL_EPSILON) { if (std::fabs(num)<QL_EPSILON) { QL_FAIL("can't compute b_k for jacobi integration\n"); } else { // l'Hospital num = 4.0*i*(i+beta_)* (2.0*i+2*alpha_+beta_); denom= 2.0*(2.0*i+alpha_+beta_); denom*=denom-1; QL_ASSERT(std::fabs(num)>QL_EPSILON, "can't compute b_k for jacobi integration\n"); } } return num / denom; } //----------------------------------------------------------------------------------- QL_EPSILON can be replaced with something more suitable, but at least it should be OK in most of the cases. Regards, Kakhkhor Abdijalilov. ------------------------------------------------------------------------------ This SF.net email is sponsored by Make an app they can't live without Enter the BlackBerry Developer Challenge http://p.sf.net/sfu/RIM-dev2dev _______________________________________________ QuantLib-dev mailing list [hidden email] https://lists.sourceforge.net/lists/listinfo/quantlib-dev |
Re: gaussianorthogonalpolynomial.c
|
|
Hi Kakhkhor, Kakhkhor Abdijalilov schrieb: > I compiled QuantLib with ICC for x64 target. Everything went well > except a glitch with Null class (I reported it in my previous e-mail) > and some issues in gaussianorthogonalpolynomial.cpp. many polynomials functions are already implemented in boost math toolkit. Can you try using these functions? Kim > Current > implementations of GaussJacobiPolynomial::alpha and > GaussJacobiPolynomial::beta have several issues. > > 1) On x64 targets Size is 64bit and conversion to double may generate > a warning. For some strange reason ICC refuses to compile this line > Real denom = (2.0*i+alpha_+beta_)*(2.0*i+alpha_+beta_+2); > ICC bails out with an internal error. I don't know why it happens (ICC > bug?), but casting i from Size into int fixes the problem. > 2) Current implementation checks if num and denom variables are zero > and then tries to use l'Hospital rule. But denom==0 can be true only > if alpha_=0 and beta_=0 and/or i=0. In that case num=0 and l'Hospital > rule is never used. In anyway, direct conversion of floating point > numbers into bool is a fishy business at best. We should consider > using fuzzy comparison instead of direct comparison, something like > this: > > //----------------------------------------------------------------------------------- > Real GaussJacobiPolynomial::alpha(Size iIn) const { > int i = static_cast<int>(iIn); > Real num = beta_*beta_ - alpha_*alpha_; > Real denom = (2.0*i+alpha_+beta_)*(2.0*i+alpha_+beta_+2); > > if (std::fabs(denom)<QL_EPSILON) { > if (std::fabs(num)<QL_EPSILON) { > QL_FAIL("can't compute a_k for jacobi integration\n"); > } > else { > // l'Hospital > num = 2*beta_; > denom= 2*(2.0*i+alpha_+beta_+1); > QL_ASSERT(std::fabs(num)>QL_EPSILON, "can't compute > a_k for jacobi integration\n"); > } > } > > return num / denom; > } > > Real GaussJacobiPolynomial::beta(Size iIn) const { > int i = static_cast<int>(iIn); > Real num = 4.0*i*(i+alpha_)*(i+beta_)*(i+alpha_+beta_); > Real denom = (2.0*i+alpha_+beta_)*(2.0*i+alpha_+beta_) > * ((2.0*i+alpha_+beta_)*(2.0*i+alpha_+beta_)-1); > > if (std::fabs(denom)<QL_EPSILON) { > if (std::fabs(num)<QL_EPSILON) { > QL_FAIL("can't compute b_k for jacobi integration\n"); > } else { > // l'Hospital > num = 4.0*i*(i+beta_)* (2.0*i+2*alpha_+beta_); > denom= 2.0*(2.0*i+alpha_+beta_); > denom*=denom-1; > QL_ASSERT(std::fabs(num)>QL_EPSILON, "can't compute > b_k for jacobi integration\n"); > } > } > return num / denom; > } > //----------------------------------------------------------------------------------- > > QL_EPSILON can be replaced with something more suitable, but at least > it should be OK in most of the cases. > > > Regards, > Kakhkhor Abdijalilov. > > ------------------------------------------------------------------------------ > This SF.net email is sponsored by > > Make an app they can't live without > Enter the BlackBerry Developer Challenge > http://p.sf.net/sfu/RIM-dev2dev > _______________________________________________ > QuantLib-dev mailing list > [hidden email] > https://lists.sourceforge.net/lists/listinfo/quantlib-dev > > ------------------------------------------------------------------------------ This SF.net email is sponsored by Make an app they can't live without Enter the BlackBerry Developer Challenge http://p.sf.net/sfu/RIM-dev2dev _______________________________________________ QuantLib-dev mailing list [hidden email] https://lists.sourceforge.net/lists/listinfo/quantlib-dev |
Re: gaussianorthogonalpolynomial.c
|
|
They are hardwired into Longstaff-Schwartz pricing engine. I am coding
my own LS method and using QL for sanity check. Btw, can anyone explain why LS path pricer does calibration and pricing separately? Current implementation simply discards all paths used for calibration. It is not only wasteful, but also doesn't agree with Longstaff and Schwartz recipe. The proper way to do LS is to use all path both for calibration and pricing. EquityOption.cpp example uses only 4096 paths for calibration and tries to price with tolerance=0.02. I think 4096 samples are not adequate to achieve such a small tolerance. Any thoughts? ------------------------------------------------------------------------------ This SF.net email is sponsored by Make an app they can't live without Enter the BlackBerry Developer Challenge http://p.sf.net/sfu/RIM-dev2dev _______________________________________________ QuantLib-dev mailing list [hidden email] https://lists.sourceforge.net/lists/listinfo/quantlib-dev |
Re: gaussianorthogonalpolynomial.c
|
|
In reply to this post by Kakhkhor Abdijalilov
On Tue, 2010-08-10 at 09:49 -0500, Kakhkhor Abdijalilov wrote:
> I compiled QuantLib with ICC for x64 target. Everything went well > except a glitch with Null class (I reported it in my previous e-mail) > and some issues in gaussianorthogonalpolynomial.cpp. Current > implementations of GaussJacobiPolynomial::alpha and > GaussJacobiPolynomial::beta have several issues. Kakhkhor, may you resend your fix as a patch? It would make it easier for me to apply it correctly. Thanks, Luigi -- Better to remain silent and be thought a fool than to speak out and remove all doubt. -- Abraham Lincoln ------------------------------------------------------------------------------ This SF.net email is sponsored by Make an app they can't live without Enter the BlackBerry Developer Challenge http://p.sf.net/sfu/RIM-dev2dev _______________________________________________ QuantLib-dev mailing list [hidden email] https://lists.sourceforge.net/lists/listinfo/quantlib-dev |
Re: gaussianorthogonalpolynomial.c
|
|
In reply to this post by Kakhkhor Abdijalilov
Hi
On Thursday 12 August 2010 10:33:16 Kakhkhor Abdijalilov wrote: > Btw, can anyone explain why LS path pricer does calibration and > pricing separately? Current implementation simply discards all paths > used for calibration. It is not only wasteful... If you are using the same paths for both calibration and valuation the resulting estimator includes a "foresight bias". A standard approach to get rid of the foresight bias is to use simple discard all paths used for calibration. Please find more details and other algorithms to remove the foresight bias here http://www.christian-fries.de/finmath/foresightbias/ > EquityOption.cpp example > uses only 4096 paths for calibration and tries to price with > tolerance=0.02 Values are taken (more or less) from Glasserman, Monte-Carlo-Methods in Financial Engineering. Within the Monte-Carlo error (0.02) the Longstaff Schwartz price (4.481675) is consistent with the other american pricer (especially with the finite different pricer, 4.486118). Increasing the number of calibration paths towards e.g. 65535 does not change the price significantly (4.464887). IMO 4096 calibration paths are enough for this example. best regards Klaus ------------------------------------------------------------------------------ This SF.net email is sponsored by Make an app they can't live without Enter the BlackBerry Developer Challenge http://p.sf.net/sfu/RIM-dev2dev _______________________________________________ QuantLib-dev mailing list [hidden email] https://lists.sourceforge.net/lists/listinfo/quantlib-dev |
Re: gaussianorthogonalpolynomial.c
|
|
>may you resend your fix as a patch?
Done. >A standard approach to get of the foresight bias is to use simple discard all paths used for >calibration That makes sense. Thanks for explanation. ------------------------------------------------------------------------------ This SF.net email is sponsored by Make an app they can't live without Enter the BlackBerry Developer Challenge http://p.sf.net/sfu/RIM-dev2dev _______________________________________________ QuantLib-dev mailing list [hidden email] https://lists.sourceforge.net/lists/listinfo/quantlib-dev |
«
Return to quantlib-dev
|
1 view|%1 views
| Free forum by Nabble | Edit this page |