LocalvolSurface.cpp

Posted by MH_quant on
URL: http://quantlib.414.s1.nabble.com/LocalvolSurface-cpp-tp12806.html

Hello everybody,

 

i recently did a lot of testing with the localvolSurface class (which contains Gatherals Dupire Formula) and i am not very happy with the results. Ok, the class is (regarding to the documentation) untested, so I guess I cant expect it to work properly. But I figured that numerical problems cause this class to return the error message “negative local vol …  the black vol surface is not smooth enough”. This also happens for very smooth black vol surfaces. The Problem lies in this code:

 

 

 

 

 

Real forwardValue = underlying *

            (dividendTS->discount(t, true)/

             riskFreeTS->discount(t, true));

 

        // strike derivatives

        Real strike, y, dy, strikep, strikem;

        Real w, wp, wm, dwdy, d2wdy2;

        strike = underlyingLevel;

        y = std::log(strike/forwardValue);

        dy = ((y!=0.0) ? y*0.000001 : 0.000001);

        strikep=strike*std::exp(dy);

        strikem=strike/std::exp(dy);

        w  = blackTS->blackVariance(t, strike,  true);

        wp = blackTS->blackVariance(t, strikep, true);

        wm = blackTS->blackVariance(t, strikem, true);

        dwdy = (wp-wm)/(2.0*dy);

        d2wdy2 = (wp-2.0*w+wm)/(dy*dy);

 

        // time derivative

        Real dt, wpt, wmt, dwdt;

        if (t==0.0) {

            dt = 0.0001;

            wpt = blackTS->blackVariance(t+dt, strike, true);

            QL_ENSURE(wpt>=w,

                      "decreasing variance at strike " << strike

                      << " between time " << t << " and time " << t+dt);

            dwdt = (wpt-w)/dt;

        } else {

            dt = std::min<Time>(0.0001, t/2.0);

            wpt = blackTS->blackVariance(t+dt, strike, true);

            wmt = blackTS->blackVariance(t-dt, strike, true);

            QL_ENSURE(wpt>=w,

                      "decreasing variance at strike " << strike

                      << " between time " << t << " and time " << t+dt);

            QL_ENSURE(w>=wmt,

                      "decreasing variance at strike " << strike

                      << " between time " << t-dt << " and time " << t);

            dwdt = (wpt-wmt)/(2.0*dt);

        }

 

        if (dwdy==0.0 && d2wdy2==0.0) { // avoid /w where w might be 0.0

            return std::sqrt(dwdt);

        } else {

            Real den1 = 1.0 - y/w*dwdy;

            Real den2 = 0.25*(-0.25 - 1.0/w + y*y/w/w)*dwdy*dwdy;

            Real den3 = 0.5*d2wdy2;

            Real den = den1+den2+den3;

            Real result = dwdt / den;

            QL_ENSURE(result>=0.0,

                      "negative local vol^2 at strike " << strike

                      << " and time " << t

                      << "; the black vol surface is not smooth enough");

            return std::sqrt(result);

            // return std::sqrt(dwdt / (1.0 - y/w*dwdy +

            //    0.25*(-0.25 - 1.0/w + y*y/w/w)*dwdy*dwdy + 0.5*d2wdy2));

        }

 

 

 

 

To be a bit more precise  the problem lies in the second derivative of the black variance with respect to the strike. Even if I take surfaces without Smile/Skew, i.e. flat ones, I still get the problem there. This is because (wp-2.0*w+wm) is not exactly zero but very very small. And this gets devided by 0. 000000000001. So if it is not exactly zero but very very little negative, it can blow up the whole calculations and we end up with this error message. To avoid these numerical problems I now decreased the accuracy in the derivatives a little bit and also introduced a check that sets this term to zero if it is very very very small (smaller than 1.0e-12). Furthermore I changed the derivative with respect to the time. This was not necessary but I wanted to make it equal to the localvolcurve.cpp code which works just fine. The modified code now looks like:

 

 

 

 

 

Real forwardValue = underlying *

          (dividendTS->discount(t, true)/

           riskFreeTS->discount(t, true));

 

 

        // strike derivatives

        Real strike, y, dy, strikep, strikem;

        Real w, wp, wm, dwdy, d2wdy2;

            Real z1,z2;

       

            // strike ist gegeben

            strike = underlyingLevel;

       

            // log(strike/forwardValue)

            y = std::log(strike/forwardValue);

            std::cout << "y: "    << y << std::endl;

       

            // wir leiten blackScholesVariance nach y ab, bilde daher diskrete kleine unterteilung

            dy = ((y!=0.0) ? y*0.0001 : 0.0001);

            std::cout << "dy: "    << dy << std::endl;

 

 

        strikep=strike*std::exp(dy);

        strikem=strike/std::exp(dy);

        w  = blackTS->blackVariance(t, strike,  true);

        wp = blackTS->blackVariance(t, strikep, true);

        wm = blackTS->blackVariance(t, strikem, true);

        z1=wp-wm;

            if (std::abs(z1) < 1.0e-12)

                  z1=0;

        dwdy = (z1)/(2.0*dy);

        z2=wp-2.0*w+wm;

        if (std::abs(z2) < 1.0e-12)

                  z2=0;

        d2wdy2 = (z2)/(dy*dy);

           

           

 

 

        // time derivative

        Real dt, wpt, wmt, dwdt;

 

            dt = (1.0/365.0);

        wpt = blackTS->blackVariance(t+dt, strike, true);

            QL_ENSURE(wpt>=w,

                      "decreasing variance at strike " << strike

                      << " between time " << t << " and time " << t+dt);

        dwdt = (wpt-w)/dt;

       

 

        if (dwdy==0.0 && d2wdy2==0.0) { // avoid /w where w might be 0.0

            return std::sqrt(dwdt);

        } else {

            Real den1 = 1.0 - y/w*dwdy;

            Real den2 = 0.25*(-0.25 - 1.0/w + y*y/w/w)*dwdy*dwdy;

            Real den3 = 0.5*d2wdy2;

            Real den = den1+den2+den3;

            Real result = dwdt / den;

            QL_ENSURE(result>=0.0,

                      "negative local vol^2 at strike " << strike

                      << " and time " << t

                      << "; the black vol surface is not smooth enough");

            return std::sqrt(result);

            // return std::sqrt(dwdt / (1.0 - y/w*dwdy +

            //    0.25*(-0.25 - 1.0/w + y*y/w/w)*dwdy*dwdy + 0.5*d2wdy2));

        }

 

 

 

 

So far all the testing I did I received much better and more stable results. But since I am not an expert in computational and numerical methods regarding precision, can anybody who is a bit more experienced with c++ numerical issues give me some advice and probably check this out?

 

 

Greetings

Michael


------------------------------------------------------------------------------
Stay on top of everything new and different, both inside and
around Java (TM) technology - register by April 22, and save
$200 on the JavaOne (SM) conference, June 2-5, 2009, San Francisco.
300 plus technical and hands-on sessions. Register today.
Use priority code J9JMT32. http://p.sf.net/sfu/p
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev