[Quantlib-dev] LocalvolSurface.cpp

[Michael H. <Mic...@gm...>]

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