LocalvolSurface.cpp

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
QuantLib-dev@...
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Hi Michael

> 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.

Last week Klaus fixed time derivative (now performed at constant
moneyness instead of constant strike) and added tests.

> 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”. [...]
> the problem lies in the second derivative of the
> black variance with respect to the strike.

you are right. I will fix it, probably moving time/strike (numerical)
derivatives in the BlackVolTermStructure base interface, allowing for
overloading in derived classes which might provide exact derivatives.

BTW I have a related question. Black ATM variance must be increasing
in time; I've always taken for granted that this is also true for
every (not just ATM) constant moneyness section of the Black surface.
Is this a non-arbitrage result or shaky common sense?

ciao -- Nando

------------------------------------------------------------------------------
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
QuantLib-dev@...
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Hi Michael,

you wrote
> 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.

To me the root of the problem was the line

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

For ve y small y this code leads to unrealistic small dy and to numerical
problems during the calculation of the difference quotient. Therfore I've
changed it into

        dy = ((std::fabs(y) > 0.001) ? y*0.0001 : 0.000001);

and at least for my tests the numerical problems with the difference quotient
disappeared. (pls see the latest version of localvolsurface.cpp in the SVN
repository). Could you test this fix using your test cases?

regards
 Klaus


------------------------------------------------------------------------------
Crystal Reports - New Free Runtime and 30 Day Trial
Check out the new simplified licensign option that enables unlimited
royalty-free distribution of the report engine for externally facing
server and web deployment.
http://p.sf.net/sfu/businessobjects
_______________________________________________
QuantLib-dev mailing list
QuantLib-dev@...
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Hallo Klaus,

I checked out your new version of the localvolsurface.cpp from the SVN and
run some extensive tests.

First of all I found some minor errors in your code which are:

1. (Line 116): Real strikept = strike*dr*dq/(drpt*dqpt);
This should be:
        Real strikept = strike*dr*dqpt/(drpt*dq);

2. (same in line 130 and 131):
        Real strikept = strike*dr*dq/(drpt*dqpt);
        Real strikemt = strike*dr*dq/(drmt*dqmt);
This should be:
        Real strikept = strike*dr*dqpt/(drpt*dq);
        Real strikemt = strike*dr*dqmt/(drmt*dq);


Just minor things, but think it over. You divide through the risk free
discount and multiply with the dividend discount. And if you do this from t
to t+dt then this is what you get.


Now a few words what I figured out by extensive testing. I ran tests with MC
with 100 thousand of paths. And your Bug-Fix brought some slight
improvements but the problem with the second derivative still remains.

It happens very often that the program crashes and tells me "negative local
vol^2 at ..."

The black vol surface I am using is very smooth. So this is not the reason
for the problems. Much more is the problem again the numerical instability
when taking the second derivative of the implied variance surface with
respect to the log-moneyness.

Let me introduce 2 more variables:

Real z1,z2;

Where

z1=wp-wm;

z2=wp-2.0*w+wm;


Then the first derivative of w with respect to log-moneyness is:

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

And the second is:

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


It happens at different strikes and moneyness under some conditions that the
d2wdy2 blows up and becomes something like -4231,12

Like a really big negative number. This obviously makes the denominator
which consists of  den1+den2+den3  negative and the whole program crashes.

The numerator is always positive since we make sure that the implied
variance is monotone increasing.

So the only thing we have to take care about is that our denominator is not
getting negative.

I figured out (by extensive testing and only empirically) that

0.4 < den1+den2<1.3

So at least for the surfaces I was testing, this was always given. So
basically we have to make sure that den3 is not getting smaller then -0.4 in
order for the program not to crash.

But my tests showed that den3 is either in the green zone or it blows up
dramatically which gives me the suppression that we encounter some numerical
problems under certain conditions.

The only way how I managed for the program to run stable is by controlling
z1,z2 by setting:

if ((std::abs(z1) < std::abs(dy)/1000 && z1!=0.0) || std::abs(z1) >
2*std::abs(dy))
                        z1=0;

and


if ((std::abs(z2) < dy*dy/1000 && z2!=0.0) || std::abs(z2) > dy*dy)
                        z2=0;


So what I basically do is setting the first derivative to zero when it is
getting so small that it doesn't really affect our result anymore anyways.
Same with the second derivative.

But the more critical part is setting the derivative to zero when it blows
up. I just cut off all critical values.

This makes (obviously) the program run stable.

On the off-side I cant really say for sure (at least till now) how big the
impacts of this manipulation are for the precision of our results. I figured
that we basically never get problems with z1 which means with the first
derivative. But we set the second derivative which means z2 quit often to
zero when it falls out of the good range.

For some reason I don't get rid of the feeling that we would be better of to
work completely without the second derivative since it seems to be
impossible to get this numerically under control. But would that probably
make the dupire formula useless to us?

Do you have any Ideas? Or can I help you with this issue any further?

Greetings,
Michael


PS: is it possible to write a paper which examines under which conditions
(what range for what variables and what interdependencies of the variables
and what restrictions to the shape of the black vol surface) the formula is
theoretically/mathematically possible (i.e. doesn't get negative) and then
use this new gained knowledge to make our code stable with the smallest
impact possible to the precision of the result?



-----Original Message-----
From: Klaus Spanderen [mailto:klaus@...]
Sent: Freitag, 24. April 2009 00:07
To: quantlib-dev@...
Cc: Michael Heckl; Ferdinando Ametrano
Subject: Re: [Quantlib-dev] LocalvolSurface.cpp

Hi Michael,

you wrote
> 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.

To me the root of the problem was the line

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

For ve y small y this code leads to unrealistic small dy and to numerical
problems during the calculation of the difference quotient. Therfore I've
changed it into

        dy = ((std::fabs(y) > 0.001) ? y*0.0001 : 0.000001);

and at least for my tests the numerical problems with the difference
quotient
disappeared. (pls see the latest version of localvolsurface.cpp in the SVN
repository). Could you test this fix using your test cases?

regards
 Klaus


------------------------------------------------------------------------------
Crystal Reports - New Free Runtime and 30 Day Trial
Check out the new simplified licensign option that enables unlimited
royalty-free distribution of the report engine for externally facing
server and web deployment.
http://p.sf.net/sfu/businessobjects
_______________________________________________
QuantLib-dev mailing list
QuantLib-dev@...
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Hi Micheal

> First of all I found some minor errors in your code which are:

yes, you are obviously right. I've changed the code in the SVN repository
accordingly. Thanks for the hint!

>
> It happens at different strikes and moneyness under some conditions that
> the d2wdy2 blows up and becomes something like -4231,12
>

Which interpolation scheme do you use for the volatility surface. The standard
interpolation is linear in the variance which very often leads to problems
with the second derivative. Therefore I've used BicubicSplineInterpolation

    volTS->setInterpolation<Bicubic>();

 
Does this happen only for deep ITM or OTM paths? (Is this a extrapolation
problem for very large (very small) spots?) At least for "extrem"
extrapolations" I found simillar problems and "solved" it be setting negative
variances to zero.

Can I get access to the parameters of your test-surface?

> For some reason I don't get rid of the feeling that we would be better of
> to work completely without the second derivative since it seems to be
> impossible to get this numerically under control.

No, IMO the second derivative is needed.

> But would that probably
> make the dupire formula useless to us?

People your using e.g. splines together with an optimization technique like in

Reconstructing The Unknown Local Volatility Function, Thomas F. Coleman,
Yuying Li, ARUN VERMA,
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.41.6202

to solve the instability. But a lot of code is needed to implemented this;-(.
For the time being I think we should give Nando's approach a try and use the
first and second derivatives taken from the interpolation method


regards
 Klaus

------------------------------------------------------------------------------
Crystal Reports - New Free Runtime and 30 Day Trial
Check out the new simplified licensign option that enables unlimited
royalty-free distribution of the report engine for externally facing
server and web deployment.
http://p.sf.net/sfu/businessobjects
_______________________________________________
QuantLib-dev mailing list
QuantLib-dev@...
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Hallo Klaus,

I ran my tests with both. BiLinear and BiCubic intrapolation. I encountered
the instability issues with both.

You are right that the problems appear far ITM and OTM. But I cant really
figure out why, since I set up a constant extrapolation for both, Strike and
Maturity in my BlackVarianceSurface.

I attached the cpp file that I use for my tests as well as a small jpg of
the surface that is included in my test cases.

I hope this gives you all the information you need and you can give me some
feedback on my test cases.

Thank you also for the link. To me it seem like if we want to use the dupire
formula efficient, stable and productive, we wont get around implementing
some fancy optimization-splines and smoothing algorithm.

I am looking forward to hear back from you.

Greetings,
Michael



-----Original Message-----
From: Klaus Spanderen [mailto:klaus@...]
Sent: Sonntag, 26. April 2009 22:01
To: quantlib-dev@...
Cc: Michael Heckl
Subject: Re: [Quantlib-dev] LocalvolSurface.cpp

>
> It happens at different strikes and moneyness under some conditions that
> the d2wdy2 blows up and becomes something like -4231,12
>

Which interpolation scheme do you use for the volatility surface. The
standard
interpolation is linear in the variance which very often leads to problems
with the second derivative. Therefore I've used BicubicSplineInterpolation

    volTS->setInterpolation<Bicubic>();

 
Does this happen only for deep ITM or OTM paths? (Is this a extrapolation
problem for very large (very small) spots?) At least for "extrem"
extrapolations" I found simillar problems and "solved" it be setting
negative
variances to zero.

Can I get access to the parameters of your test-surface?

> But would that probably
> make the dupire formula useless to us?

People your using e.g. splines together with an optimization technique like
in

Reconstructing The Unknown Local Volatility Function, Thomas F. Coleman,
Yuying Li, ARUN VERMA,
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.41.6202

to solve the instability. But a lot of code is needed to implemented
this;-(.
For the time being I think we should give Nando's approach a try and use the

first and second derivatives taken from the interpolation method


regards
 Klaus

// localSurfaceIssue.cpp : Defines the entry point for the console application.
//


#include <ql/quantlib.hpp>


#include <boost/timer.hpp>
#include <iostream>
#include <iomanip>
#include <string>

using namespace std;
using namespace QuantLib;


/*
        This Class is programmed to track down numerical problems
        with the local Vol Surface (localvolsurface.cpp)
*/



// ####################################################################
// method deklaration -- prototyp (find definitions below; this is just to avoid a header file)



/*
        Use the analytic formula to compare prices of the localvolsurface formula
        if our localvolsurface impl works fine we should get the same price
*/
Real calculatePlainVanillaPriceWithBSAnalytic(Time Laufzeit,
                                                          Real Strike,
                                                          Option::Type type,
                                                          Date SpotDate,
                                                          const boost::shared_ptr<BlackScholesMertonProcess>& process);


/*
        this is the monte Carlo Engine that simulates thousand of pathes
        it crashes if our localvolsurface encounters numerical issues
*/
Real calculatePlainVanillaPriceWithMC(Time Laufzeit,
                                                          Real Strike,
                                                          Option::Type type,
                                                          const boost::shared_ptr<BlackScholesMertonProcess>& process);





/*
        this was the critical method in our localvolsurface class before we started to work on it and improve it
*/
Volatility OldDupireLocalVolSurfaceImpl(Time t,
                                                          Real underlyingLevel,
                                                          const Handle<BlackVolTermStructure>& blackTS,
                              const Handle<YieldTermStructure>& riskFreeTS,
                              const Handle<YieldTermStructure>& dividendTS,
                                                          Real underlying);


/*
        this is the critical method in our localvolsurface as you can find it in svn currently
        it has some improvements but still crashes at far ITM or OTM
*/
Volatility NewDupireLocalVolSurfaceImpl(Time t,
                                                          Real underlyingLevel,
                                                          const Handle<BlackVolTermStructure>& blackTS,
                              const Handle<YieldTermStructure>& riskFreeTS,
                              const Handle<YieldTermStructure>& dividendTS,
                                                          Real underlying);


/*
        this is the critical method modified experimentally to see  if numerical issues still
        occur after some modifications. Here I try out all sorts of modifications.
*/
Volatility ExperimentalDupireLocalVolSurfaceImpl(Time t,
                                                          Real underlyingLevel,
                                                          const Handle<BlackVolTermStructure>& blackTS,
                              const Handle<YieldTermStructure>& riskFreeTS,
                              const Handle<YieldTermStructure>& dividendTS,
                                                          Real underlying);




// ####################################################################
// if debugMode is set to true you get more output on the console
bool debugMode = false;


// ####################################################################
// Main Methode (including the whole test)
int main(int, char* [])
{
        try {
        /*
                        Just a timer to control the runtime
                */
                boost::timer timer;
               

               
                Size widths[] = { 45, 20 };
                Size i;

                DayCounter dayCounter = Actual360();
                Calendar calendar = NullCalendar();

               
                // Setting up our Option Parameters
                Date SpotDate = Date::todaysDate();
                Time Laufzeit = 10;
                Real SpotPreis = 100;
                Handle<Quote> x0(boost::shared_ptr<Quote>(new SimpleQuote(SpotPreis)));
                boost::shared_ptr<SimpleQuote> riskFreeRate(new SimpleQuote(0.03));
                Handle<YieldTermStructure> r(boost::shared_ptr<YieldTermStructure>(new FlatForward(0, NullCalendar(), Handle<Quote>(riskFreeRate), Actual360())));
                boost::shared_ptr<SimpleQuote> dividentYield(new SimpleQuote(0.0));
                Handle<YieldTermStructure> q(boost::shared_ptr<YieldTermStructure>(new FlatForward(0, NullCalendar(), Handle<Quote>(dividentYield), Actual360())));
               
                Real Strike = x0->value()*1.3;

                // build up vector with maturities in month
                Integer t[] = { 12, 24, 36, 48, 60, 72, 84, 96, 108, 120 };
                std::vector<Date> dates;
                for (i = 0; i < 10; ++i) {
                                Date pushinDatum = calendar.advance(SpotDate, t[i], Months);
                                if (debugMode) {
                                        std::cout << "We push into the dates vector: "    << pushinDatum << std::endl;
                                        std::cout << "YearFraction is: "    << dayCounter.yearFraction(SpotDate,pushinDatum) << std::endl;
                                        std::cout << std::endl;
                                }
                                dates.push_back(pushinDatum);
                }
               

               
               
               

                //#######################################
                /*
                        Printing out our Setup:
                */
                std::cout << "Option Setup:"<< std::endl;
                std::cout << std::setw(widths[0]) << std::left << "SpotDate: " << std::setw(widths[1]) << std::left << SpotDate  << std::endl;
                std::cout << std::setw(widths[0]) << std::left << "Options Maturity in years: " << std::setw(widths[1]) << std::left << Laufzeit  << std::endl;
                std::cout << std::setw(widths[0]) << std::left << "Exercise Date: " << std::setw(widths[1]) << std::left << NullCalendar().advance(SpotDate, Laufzeit, Years)  << std::endl;
                std::cout << std::setw(widths[0]) << std::left << "Strike Price: " << std::setw(widths[1]) << std::left << Strike  << std::endl;
                std::cout << std::endl;
                std::cout << std::setw(widths[0]) << std::left << "Used DayCount Convention: " << std::setw(widths[1]) << std::left << dayCounter.name()  << std::endl;
                std::cout << std::setw(widths[0]) << std::left << "Used Calender: " << std::setw(widths[1]) << std::left << calendar.name()  << std::endl;
                std::cout << std::endl;
                std::cout << std::setw(widths[0]) << std::left << "constant risk free rate: " << std::setw(widths[1]) << std::left << riskFreeRate->value()  << std::endl;
                std::cout << std::setw(widths[0]) << std::left << "constant dividend yield: " << std::setw(widths[1]) << std::left << dividentYield->value()  << std::endl;
                std::cout << "-----------------------------" << std::endl;
                std::cout << "-----------------------------" << std::endl;
                std::cout << std::endl;
                std::cout << std::endl;
                std::cout << std::endl;


               

               
               
                //#######################################
                // Local Vol Surface  Tests with BlackVarianceSurface
                       
                        // this is the surface, it is quite smooth
                        // it was actually generated by heston with the paramters
                        /*const Real v0 = 0.16;
                        const Real kappa = 1.0;
                        const Real theta = 0.04;
                        const Real sigma = 0.8;
                        const Real rho = -0.4;*/
                        // feel free to try any surface, just copy here

                         Volatility volas[] =
                          { 0.410908,0.346389,0.309471,0.287156,0.272528,0.262294,0.254745,0.24899,0.244452,0.240788,
                                0.362241,0.306316,0.27681,0.259683,0.248766,0.241307,0.235923,0.2319,0.228788,0.22632,
                                0.333088,0.283417,0.258434,0.244316,0.235512,0.22962,0.225452,0.222398,0.220082,0.21828,
                                0.319807,0.273194,0.250247,0.237459,0.229587,0.224386,0.220755,0.21813,0.216166,0.21466,
                                0.307689,0.263893,0.242759,0.231155,0.224117,0.219538,0.216393,0.214158,0.212517,0.211282,
                                0.302161,0.259622,0.239289,0.228214,0.221552,0.217256,0.214334,0.212279,0.210787,0.209678,
                                0.297044,0.255625,0.236009,0.225416,0.2191,0.215068,0.212355,0.210469,0.209118,0.208129,
                                0.292382,0.251916,0.232927,0.222762,0.216762,0.212973,0.210454,0.208727,0.207509,0.206634,
                                0.288214,0.248511,0.230046,0.220256,0.214538,0.21097,0.20863,0.207052,0.205959,0.20519,
                                0.28148,0.242663,0.224915,0.215694,0.210434,0.207239,0.205212,0.203896,0.203028,0.202454,
                                0.276956,0.238149,0.22065,0.211744,0.20679,0.203874,0.202093,0.200994,0.200316,0.199911,
                                0.273893,0.233039,0.214729,0.205686,0.200879,0.198219,0.196728,0.195919,0.195518,0.19537,
                                0.278921,0.23313,0.211625,0.200898,0.195283,0.192299,0.190748,0.19002,0.189772,0.189813 };

                                       
                       
                        Real s[] = { 0.6, 0.75, 0.85, 0.90, 0.95, 0.975, 1.0, 1.025, 1.05, 1.10, 1.15, 1.25, 1.40 };
                        std::vector<Real> strikes;
                       
                        // build up vektor with strikes
                        for (i = 0; i < 13; ++i) {
                                Real tempStrike = SpotPreis*s[i];
                                if (debugMode) {
                                        std::cout << "We push into the strikes vector: "    << tempStrike << std::endl;
                                        std::cout << std::endl;
                                }
                                strikes.push_back(tempStrike);
                        }

                        // figure out correct dimensions of matrix for implied volatilities
                        Size anzahlStrikes = strikes.size();
                        Size anzahlMaturities = sizeof(t)/sizeof(t[1]);
                        if (debugMode) {
                                std::cout << "number of strikes we have: "    << anzahlStrikes << std::endl;
                                std::cout << "number of maturities we have: "    << anzahlMaturities << std::endl;
                        }

                       

                        // build up matrix with implied volatilities
                        Matrix volSurface(anzahlStrikes,anzahlMaturities);
                        for (Size s = 0; s < anzahlStrikes; ++s) {
                                for (Size m = 0; m < anzahlMaturities; ++m) {
                                        volSurface[s][m]=volas[s*anzahlMaturities+m];
                                }
                        }



                        // We use constant extrapolation at strikes and maturities. I get problems with non-increasing variances far ITM/OTM if i do not extrapolate constant.
                        boost::shared_ptr<BlackVolTermStructure> impliedSurface(new BlackVarianceSurface(SpotDate,calendar,dates,strikes,volSurface,dayCounter,BlackVarianceSurface::ConstantExtrapolation,BlackVarianceSurface::ConstantExtrapolation));
                        Handle<BlackVolTermStructure> impliedVolSurface(impliedSurface);


                        // here we can change interpolation from Bilinear to Bicubic. I tryed it both, but still get problems.
                        (boost::dynamic_pointer_cast<BlackVarianceSurface>(impliedSurface))->setInterpolation<Bicubic>();


                        // here i build the localvolsurface class explicitly to access the currently implemented methods
                        boost::shared_ptr<LocalVolTermStructure> localSurface(new LocalVolSurface(impliedVolSurface,r,q,x0));

                        // just increase the precision of the output
                        std::cout << std::setprecision(17) << std::endl;



                        // Here you can call the single implementations and track down issues at a certain strike and time
                        // merely comment/uncoment the impl you want to run tests on
                        // the only impl that runs stable is my experimental one, but it cuts of problems ... (try it out!)
                        Real strike = 60.0033;
                        Real time = 1.64722;

                        /*std::cout << std::endl;
                        std::cout << "Local Vol at strike " << strike << " and time " << time << " (old Impl): "    << OldDupireLocalVolSurfaceImpl(time,strike,impliedVolSurface,r,q,x0->value()) << std::endl;
                        std::cout << std::endl;
                        std::cout << std::endl;
                        std::cout << "Local Vol at strike " << strike << " and time " << time << " (new Impl): "    << NewDupireLocalVolSurfaceImpl(time,strike,impliedVolSurface,r,q,x0->value()) << std::endl;
                        std::cout << std::endl;
                        std::cout << std::endl;
                        std::cout << "Local Vol at strike " << strike << " and time " << time << " (experimental Impl): "    << ExperimentalDupireLocalVolSurfaceImpl(time,strike,impliedVolSurface,r,q,x0->value()) << std::endl;
                        std::cout << std::endl;
                        std::cout << std::endl;
                        std::cout << "Local Vol at strike " << strike << " and time " << time << " (current Impl): "    << localSurface->localVol(time,strike) << std::endl;
                        std::cout << std::endl;
                        std::cout << std::endl;*/

                       



                        // build up a BlackScholes process with BlackVarianceSurface as BlackVolTermStructure
                        // in the function localVolatility the process automaticall builds up a localvolsurface
                        // so the diffusion is calculated instantly in the localvolsurface class
                        boost::shared_ptr<BlackScholesMertonProcess> bsmProcess3(
                                         new BlackScholesMertonProcess(x0, q, r, impliedVolSurface));





                        /*
                                Calculating analytic BS Preis for a Call Option and a Put Option with Local Vol Surface
                                this is simply to compare the results and see how much impact different modifications
                                of the localvolsurface class have on the precision of our result.
                        */
                        std::cout << "Calculate price for a Plain Vanilla Option with BS Analytic and MC with Local Vol Surface!:"<< std::endl;
                        std::cout << std::endl;
                        std::cout << "analytical prices:"<< std::endl;
                        std::cout << "analytic BS price for a Call with Local Vol Surface: "    << calculatePlainVanillaPriceWithBSAnalytic(Laufzeit,Strike,Option::Call,SpotDate,bsmProcess3) << std::endl;
                        std::cout << "analytic BS price for a Put with Local Vol Surface: "    << calculatePlainVanillaPriceWithBSAnalytic(Laufzeit,Strike,Option::Put,SpotDate,bsmProcess3) << std::endl;
                        std::cout << "-----------------------------" << std::endl;
                        std::cout << std::endl;
                        std::cout << std::endl;

                        /*
                                Calculating MC simulated prices for a Call Option and a Put Option with Local Vol Surface
                        */
                        std::cout << "MC simulated prices:"<< std::endl;
                        std::cout << "MC price for a Call with Local Vol Surface: "    << calculatePlainVanillaPriceWithMC(Laufzeit,Strike,Option::Call,bsmProcess3)*r->discount(Laufzeit) << std::endl;
                        std::cout << "MC price for a Put with Local Vol Surface: "    << calculatePlainVanillaPriceWithMC(Laufzeit,Strike,Option::Put,bsmProcess3)*r->discount(Laufzeit) << std::endl;
                        std::cout << "-----------------------------" << std::endl;
                        std::cout << "-----------------------------" << std::endl;
                        std::cout << std::endl;
                        std::cout << std::endl;
                        std::cout << std::endl;
               





               
               

               
                //#######################################
                // End test
        Real seconds = timer.elapsed();
        Integer hours = int(seconds/3600);
        seconds -= hours * 3600;
        Integer minutes = int(seconds/60);
        seconds -= minutes * 60;
        std::cout << " \nRun completed in ";
        if (hours > 0)
            std::cout << hours << " h ";
        if (hours > 0 || minutes > 0)
            std::cout << minutes << " m ";
        std::cout << std::fixed << std::setprecision(0)
                  << seconds << " s\n" << std::endl;
        return 0;

    } catch (std::exception& e) {
        std::cout << e.what() << std::endl;
        return 1;
    } catch (...) {
        std::cout << "unknown error" << std::endl;
        return 1;
    }
}


// ####################################################################
// calculate price analytically
Real calculatePlainVanillaPriceWithBSAnalytic(Time Laufzeit,
                                                          Real Strike,
                                                          Option::Type type,
                                                          Date SpotDate,
                                                          const boost::shared_ptr<BlackScholesMertonProcess>& process) {

                /*
                        Calculating analytic MC Preis for a Call Option and a Put Option
                */
                Date ExerciseDate = NullCalendar().advance(SpotDate, Laufzeit, Years);
                boost::shared_ptr<StrikedTypePayoff> payoff(new PlainVanillaPayoff(type, Strike));
                boost::shared_ptr<Exercise> europeanExercise(new EuropeanExercise(ExerciseDate));
        VanillaOption europeanOption(payoff, europeanExercise);

                if (debugMode) {
                        std::cout << std::endl;
                        std::cout << "Variance used in the analytic BlackScholes Calculator: "    << process->blackVolatility()->blackVariance(europeanExercise->lastDate(),payoff->strike()) << std::endl;
                        std::cout << "corresponding const implied vol: "    << std::sqrt(process->blackVolatility()->blackVariance(europeanExercise->lastDate(),payoff->strike())/Actual360().yearFraction(SpotDate,ExerciseDate)) << std::endl;
                        std::cout << "last date: "    << europeanExercise->lastDate() << std::endl;
                        std::cout << "strike: "    << payoff->strike() << std::endl;
                        std::cout << "time between dates: "    << Actual360().yearFraction(SpotDate,ExerciseDate) << std::endl;
                }

                europeanOption.setPricingEngine(boost::shared_ptr<PricingEngine>(
                                     new AnalyticEuropeanEngine(process)));
                return europeanOption.NPV();
               
}



// ####################################################################
// Berechne BS Preis mit MC
Real calculatePlainVanillaPriceWithMC(Time Laufzeit,
                                                          Real Strike,
                                                          Option::Type type,
                                                          const boost::shared_ptr<BlackScholesMertonProcess>& process) {

               
               
                //#######################################
                /*
                        Setting up our MC machine
                */

                typedef PseudoRandom::rsg_type rsg_type;
        typedef PathGenerator<rsg_type>::sample_type sample_type;

        // Seed;
                BigNatural seed = 43;
                // how many discretization points do we want? It is a good Question.
                // Since we need a fine discretization in order to get a good result for local vol surface
                // i suggest about one point per day, which means a lot of calculation time
                Size discrPointsPerYear = 360;
        Size absNumOfDiscrPoints = discrPointsPerYear*Laufzeit;
                // build randomSequenceGenerator (mit MersenneTwisterUniformRng und InverseCumulativeNormal)
        rsg_type rsg = PseudoRandom::make_sequence_generator(absNumOfDiscrPoints, seed);
                // baue Pathgenerator für 1-dimensionalen Underlying Process
        PathGenerator<rsg_type> generator(process, Laufzeit, absNumOfDiscrPoints,
                                          rsg, false);
       
                // Number of pathes that we want to generate
                const Size nrTrails = 50000;

                // container for calculating the mean of MC simulated values
                GeneralStatistics stat;
               
               

                //#######################################
                /*
                        Printing out our settings to the console
                */
                if (debugMode) {
                        std::cout << "Settings of our MC Simulation: " << std::endl;
                        // write column headings
                        Size widths[] = { 45, 20 };
                        std::cout << std::setw(widths[0]) << std::left << "Time of the Simulation in years: " << std::setw(widths[1]) << std::left << Laufzeit << std::endl;
                        std::cout << std::setw(widths[0]) << std::left << "discretization points per year: " << std::setw(widths[1]) << std::left << discrPointsPerYear  << std::endl;
                        std::cout << std::setw(widths[0]) << std::left << "overall dicretization points: " << std::setw(widths[1]) << std::left << (absNumOfDiscrPoints + 1)  << std::endl;
                        std::cout << std::setw(widths[0]) << std::left << "Seed Number: " << std::setw(widths[1]) << std::left << seed  << std::endl;
                        std::cout << std::setw(widths[0]) << std::left << "number of pathes we are going to simulate: " << std::setw(widths[1]) << std::left << nrTrails  << std::endl;
                        std::cout << std::setw(widths[0]) << std::left << "Option Type: " << std::setw(widths[1]) << std::left << type  << std::endl;
                        std::cout << std::endl;
                }

               

                //#######################################
                /*
                        starting MC Simulation and return results
                */
                Size i;

                if(type==Option::Call) {
                        for (i=0; i<nrTrails; i++) {
                                sample_type path = (i%2) ? generator.antithetic()
                                                                                                 : generator.next();

                                if (i % 1000 == 0)
                                        std::cout << "simulated path : "  << i << "(of " << nrTrails << ")" << std::endl;

                                stat.add(max(Real(path.value.back()-Strike),Real(0)));
                        }
                } else {
                        for (i=0; i<nrTrails; i++) {
                                sample_type path = (i%2) ? generator.antithetic()
                                                                                                 : generator.next();

                                if (i % 1000 == 0)
                                        std::cout << "simulated path : " << i << "(of " << nrTrails << ")" << std::endl;
                               
                                stat.add(max(Real(Strike-path.value.back()),Real(0)));
                        }
                }

                return stat.mean();

}








// ####################################################################
// Old implementation before we started working on it. That is from QuantLib-0.9.7
Volatility OldDupireLocalVolSurfaceImpl(Time t,
                                                          Real underlyingLevel,
                                                          const Handle<BlackVolTermStructure>& blackTS,
                              const Handle<YieldTermStructure>& riskFreeTS,
                              const Handle<YieldTermStructure>& dividendTS,
                                                          Real underlying)  {


                 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;

                        if (result<0.0) {
                        std::cout << "strike: "    << strike << std::endl;
                        std::cout << "forwardValue: "    << forwardValue << std::endl;
                        std::cout << "y: "    << y << std::endl;
                        std::cout << "dy: "    << dy << std::endl;
                        std::cout << "strikep: "    << strikep << std::endl;
                        std::cout << "strikem: "    << strikem << std::endl;
                        std::cout << "w: "    << w << std::endl;
                        std::cout << "wp: "    << wp << std::endl;
                        std::cout << "wm: "    << wm << std::endl;
                        std::cout << "(wp-wm): "    << (wp-wm) << std::endl;
                        std::cout << "dwdy: "    << dwdy << std::endl;
                        std::cout << "(wp-2.0*w+wm): "    << (wp-2.0*w+wm) << std::endl;
                        std::cout << "d2wdy2: "    << d2wdy2 << std::endl;
                        std::cout << std::endl;
                        std::cout << "dt: "    << dt << std::endl;
                        std::cout << "dwdt: "    << dwdt << std::endl;
                        std::cout << "den: "    << den << std::endl;
                        std::cout << "den1: "    << den1 << std::endl;
                        std::cout << "den2: "    << den2 << std::endl;
                        std::cout << "den3: "    << den3 << std::endl;
                        std::cout << "result is: "    << result << std::endl;
                        }

            QL_ENSURE(result>=0.0,
                      "(old implementation!) 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));
        }
    }




// ####################################################################
// New implementation. That is the latest impl from the svn (from Klaus Spanderen)
Volatility NewDupireLocalVolSurfaceImpl(Time t,
                                                          Real underlyingLevel,
                                                          const Handle<BlackVolTermStructure>& blackTS_,
                              const Handle<YieldTermStructure>& riskFreeTS_,
                              const Handle<YieldTermStructure>& dividendTS_,
                                                          Real underlying)  {


                DiscountFactor dr = riskFreeTS_->discount(t, true);
        DiscountFactor dq = dividendTS_->discount(t, true);
        Real forwardValue = underlying*dq/dr;
       
        // strike derivatives
        Real strike, y, dy, strikep, strikem;
        Real w, wp, wm, dwdy, d2wdy2;
        strike = underlyingLevel;
        y = std::log(strike/forwardValue);
        dy = ((std::fabs(y) > 0.001) ? y*0.0001 : 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;
            DiscountFactor drpt = riskFreeTS_->discount(t+dt, true);
            DiscountFactor dqpt = dividendTS_->discount(t+dt, true);          
            Real strikept = strike*dr*dqpt/(drpt*dq);
       
            wpt = blackTS_->blackVariance(t+dt, strikept, 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);
            DiscountFactor drpt = riskFreeTS_->discount(t+dt, true);
            DiscountFactor drmt = riskFreeTS_->discount(t-dt, true);
            DiscountFactor dqpt = dividendTS_->discount(t+dt, true);
            DiscountFactor dqmt = dividendTS_->discount(t-dt, true);
           
            Real strikept = strike*dr*dqpt/(drpt*dq);
            Real strikemt = strike*dr*dqmt/(drmt*dq);
           
            wpt = blackTS_->blackVariance(t+dt, strikept, true);
            wmt = blackTS_->blackVariance(t-dt, strikemt, 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;

                        if (result<0.0) {
                        std::cout << "strike: "    << strike << std::endl;
                        std::cout << "forwardValue: "    << forwardValue << std::endl;
                        std::cout << "y: "    << y << std::endl;
                        std::cout << "dy: "    << dy << std::endl;
                        std::cout << "strikep: "    << strikep << std::endl;
                        std::cout << "strikem: "    << strikem << std::endl;
                        std::cout << "w: "    << w << std::endl;
                        std::cout << "wp: "    << wp << std::endl;
                        std::cout << "wm: "    << wm << std::endl;
                        std::cout << "(wp-wm): "    << (wp-wm) << std::endl;
                        std::cout << "dwdy: "    << dwdy << std::endl;
                        std::cout << "(wp-2.0*w+wm): "    << (wp-2.0*w+wm) << std::endl;
                        std::cout << "d2wdy2: "    << d2wdy2 << std::endl;
                        std::cout << std::endl;
                        std::cout << "dt: "    << dt << std::endl;
                        std::cout << "dwdt: "    << dwdt << std::endl;
                        std::cout << "den: "    << den << std::endl;
                        std::cout << "den1: "    << den1 << std::endl;
                        std::cout << "den2: "    << den2 << std::endl;
                        std::cout << "den3: "    << den3 << std::endl;
                        std::cout << "result is: "    << result << std::endl;
                        }

            QL_ENSURE(result>=0.0,
                                "(new implementation!) negative local vol^2 at strike " << strike
                      << " and time " << t
                      << "; the black vol surface is not smooth enough");

            return std::sqrt(result);
        }
    }



// ####################################################################
// Experimental Implementation. here i try out all sorts of modifications
// and see what impact they have on the stability and results
Volatility ExperimentalDupireLocalVolSurfaceImpl(Time t,
                                                          Real underlyingLevel,
                                                          const Handle<BlackVolTermStructure>& blackTS_,
                              const Handle<YieldTermStructure>& riskFreeTS_,
                              const Handle<YieldTermStructure>& dividendTS_,
                                                          Real underlying)  {


                 DiscountFactor dr = riskFreeTS_->discount(t, true);
        DiscountFactor dq = dividendTS_->discount(t, true);
        Real forwardValue = underlying*dq/dr;

        // strike derivatives
        Real strike, y, dy, strikep, strikem;
        Real w, wp, wm, dwdy, d2wdy2;
                Real z1,z2;

        strike = underlyingLevel;
        y = std::log(strike/forwardValue);
        dy = ((std::fabs(y) > 0.01) ? 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);
                z1=wp-wm;
                if ((std::abs(z1) < std::abs(dy)/1000 && z1!=0.0) || std::abs(z1) > 2*std::abs(dy)) {
                        std::cout << "Setting z1 to Zero! z1 was: " << z1 << std::endl;
                        z1=0;
                }
        dwdy = (z1)/(2.0*dy);
                z2=wp-2.0*w+wm;
                if ((std::abs(z2) < dy*dy/1000 && z2!=0.0) || std::abs(z2) > dy*dy) {
                        std::cout << "Setting z2 to Zero! z2 was: " << z2 << std::endl;
                        z2=0;
                }
        d2wdy2 = (z2)/(dy*dy);

        // time derivative
        Real dt, wpt, wmt, dwdt;
        if (t==0.0) {
            dt = 0.0001;
            DiscountFactor drpt = riskFreeTS_->discount(t+dt, true);
            DiscountFactor dqpt = dividendTS_->discount(t+dt, true);          
            Real strikept = strike*dr*dqpt/(drpt*dq);

            wpt = blackTS_->blackVariance(t+dt, strikept, 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);
            DiscountFactor drpt = riskFreeTS_->discount(t+dt, true);  //--> risk discount von t bis t+dt = drpt/dr
            DiscountFactor drmt = riskFreeTS_->discount(t-dt, true);  //--> risk discount von t-dt bis t = dr/drmt
            DiscountFactor dqpt = dividendTS_->discount(t+dt, true);  //--> dividend discount von t bis t+dt = dqpt/dq
            DiscountFactor dqmt = dividendTS_->discount(t-dt, true);  //--> dividend discount von t-dt bis t = dq/dqmt
           
            Real strikept = strike*dr*dqpt/(drpt*dq);
            Real strikemt = strike*dr*dqmt/(drmt*dq);
           
            wpt = blackTS_->blackVariance(t+dt, strikept, true);
            wmt = blackTS_->blackVariance(t-dt, strikemt, 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;
                       
                        if (result<0.0) {
                        std::cout << "strike: "    << strike << std::endl;
                        std::cout << "forwardValue: "    << forwardValue << std::endl;
                        std::cout << "y: "    << y << std::endl;
                        std::cout << "dy: "    << dy << std::endl;
                        std::cout << "strikep: "    << strikep << std::endl;
                        std::cout << "strikem: "    << strikem << std::endl;
                        std::cout << "w: "    << w << std::endl;
                        std::cout << "wp: "    << wp << std::endl;
                        std::cout << "wm: "    << wm << std::endl;
                        std::cout << "(wp-wm): "    << (wp-wm) << std::endl;
                        std::cout << "dwdy: "    << dwdy << std::endl;
                        std::cout << "(wp-2.0*w+wm): "    << (wp-2.0*w+wm) << std::endl;
                        std::cout << "d2wdy2: "    << d2wdy2 << std::endl;
                        std::cout << std::endl;
                        std::cout << "dt: "    << dt << std::endl;
                        std::cout << "dwdt: "    << dwdt << std::endl;
                        std::cout << "den: "    << den << std::endl;
                        std::cout << "den1: "    << den1 << std::endl;
                        std::cout << "den2: "    << den2 << std::endl;
                        std::cout << "den3: "    << den3 << std::endl;
                        std::cout << "result is: "    << result << std::endl;
                        }

            QL_ENSURE(result>=0.0,
                      "(experimental implementation!) 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));
        }
    }

------------------------------------------------------------------------------
Register Now & Save for Velocity, the Web Performance & Operations
Conference from O'Reilly Media. Velocity features a full day of
expert-led, hands-on workshops and two days of sessions from industry
leaders in dedicated Performance & Operations tracks. Use code vel09scf
and Save an extra 15% before 5/3. http://p.sf.net/sfu/velocityconf
_______________________________________________
QuantLib-dev mailing list
QuantLib-dev@...
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

On Mon, Apr 27, 2009 at 11:20 AM, Michael Heckl <Michael.Heckl@...> wrote:
> I set up a constant extrapolation for both, Strike and
> Maturity in my BlackVarianceSurface.

I don't work on equities but constant variance extrapolation in strike
and maturity seems plain wrong to me. In time it implies zero forward
volatility, in strike it violates concavity smile requirement

ciao -- Nando

------------------------------------------------------------------------------
Register Now & Save for Velocity, the Web Performance & Operations
Conference from O'Reilly Media. Velocity features a full day of
expert-led, hands-on workshops and two days of sessions from industry
leaders in dedicated Performance & Operations tracks. Use code vel09scf
and Save an extra 15% before 5/3. http://p.sf.net/sfu/velocityconf
_______________________________________________
QuantLib-dev mailing list
QuantLib-dev@...
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Hallo Nando,



-----Original Message-----
On Mon, Apr 27, 2009 at 11:20 AM, Michael Heckl <Michael.Heckl@...>
wrote:
> I set up a constant extrapolation for both, Strike and
> Maturity in my BlackVarianceSurface.

I don't work on equities but constant variance extrapolation in strike
and maturity seems plain wrong to me. In time it implies zero forward
volatility, in strike it violates concavity smile requirement

ciao -- Nando

--------------------------


What is done with the time extrapolation is the following:

if (t<=times_.back())
            return varianceSurface_(t, strike, true);
else // t>times_.back() || extrapolate
           return varianceSurface_(times_.back(), strike, true) *
                t/times_.back();


I.e. the implied variance surface is not extrapolated constant in time, but
the implied volatility surface.

What exactly do you mean by zero forward volatility. IMHO the forward
volatility would be in that case (between two dates past at extrapolation)
exactly the constant extrapolated volatility (analog as with interest rates
for instance).


And at Strikes past the maximum/minimum Strike the implied variance surface
is indeed extrapolated flat, i.e.


// enforce constant extrapolation when required
        if (strike < strikes_.front()
            && lowerExtrapolation_ == ConstantExtrapolation)
            strike = strikes_.front();
        if (strike > strikes_.back()
            && upperExtrapolation_ == ConstantExtrapolation)
            strike = strikes_.back();



Well, I don't know what you mean with your concavity smile requirement. But
I cant see what the problem should be with this extrapolation?

If you look at the No-Arbitrage requirements for the implied volatility
surface they only give limites to the slope.

But constant extrapolation means zero slope. So I cant see any restrictions
to that.

You can find no-arbitrage restrictions for the implied volatility for
instance at Roger W. Lee ("Implied Volatility: Statics, Dynamics, and
Probabalistic Interpretation")

The paper is online here:

http://www.math.uchicago.edu/~rl/



Can you send me a paper which works out the arguments you stated?


Greetings,
Michael


------------------------------------------------------------------------------
Register Now & Save for Velocity, the Web Performance & Operations
Conference from O'Reilly Media. Velocity features a full day of
expert-led, hands-on workshops and two days of sessions from industry
leaders in dedicated Performance & Operations tracks. Use code vel09scf
and Save an extra 15% before 5/3. http://p.sf.net/sfu/velocityconf
_______________________________________________
QuantLib-dev mailing list
QuantLib-dev@...
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Hi

> I ran my tests with both. BiLinear and BiCubic intrapolation. I encountered
> the instability issues with both.

BiLinear interpolation doesn't work due to the jumps in the first derivative.
BiCubic should be much better.

> You are right that the problems appear far ITM and OTM. But I cant really
> figure out why, since I set up a constant extrapolation for both, Strike
> and Maturity in my BlackVarianceSurface.

Constant extrapolation could introduce arbitrage violations far ITM or far OTM
and these arbitrage violations can lead to negative variances. Is this part
of the problem in your tests?

> Thank you also for the link. To me it seem like if we want to use the
> dupire formula efficient, stable and productive, we wont get around
> implementing some fancy optimization-splines and smoothing algorithm.

yes.;-)

regards
 Klaus

------------------------------------------------------------------------------
Register Now & Save for Velocity, the Web Performance & Operations
Conference from O'Reilly Media. Velocity features a full day of
expert-led, hands-on workshops and two days of sessions from industry
leaders in dedicated Performance & Operations tracks. Use code vel09scf
and Save an extra 15% before 5/3. http://p.sf.net/sfu/velocityconf
_______________________________________________
QuantLib-dev mailing list
QuantLib-dev@...
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Just a remark regarding interpolation. I've implemented kernel
------------------------------------------------------------------------------
Register Now & Save for Velocity, the Web Performance & Operations
Conference from O'Reilly Media. Velocity features a full day of
expert-led, hands-on workshops and two days of sessions from industry
leaders in dedicated Performance & Operations tracks. Use code vel09scf
and Save an extra 15% before 5/3. http://p.sf.net/sfu/velocityconf
_______________________________________________
QuantLib-dev mailing list
QuantLib-dev@...
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Hallo Klaus,


> BiLinear interpolation doesn't work due to the jumps in the first
derivative.
> BiCubic should be much better.

I totally agree that BiLinear interpolation is not the smoothest
interpolation out there. And you are right that in theorie you would get
points of discontinuity exactly everywhere where the linear slope changes.
Moreover the first derivative would look like a step function. That would in
theory lead to a sum of dirac delta functions in the second derivative. But
this is only in theory. Since we do discrete approximations in quantlib, I
doubt that we have this effect ...

Anyways, I do also have problems with BiCubic Interpolation. My hope is now
that the newly implemented kernel interpolation (thanks to Dima) gives
better test results.



> Constant extrapolation could introduce arbitrage violations far ITM or far
OTM
> and these arbitrage violations can lead to negative variances. Is this
part
> of the problem in your tests?


I would love to not use the constant extrapolation but instead the
InterpolatorDefaultExtrapolation. But this doesn't work since I run into
another problem by doing this. That is very far ITM/OTM the monoton variance
criteria is violated at some point and the program crashes due to this
issue. I have now Idea how to encounter this problem?? So I cant really say
if that would make it better ...


Anyways, can you send me a link or a paper with more information about the
arbitrage violations due to constant extrapolation? I still cant see why
constant extrapolation is violating the arbitrage criteria. At least it
doesn't violate any of the arbitrage criterias that I know (see Roger Lee or
Gatheral or Musiela/Rutkowski) ...


Greetings from Munich
Michael


------------------------------------------------------------------------------
Register Now & Save for Velocity, the Web Performance & Operations
Conference from O'Reilly Media. Velocity features a full day of
expert-led, hands-on workshops and two days of sessions from industry
leaders in dedicated Performance & Operations tracks. Use code vel09scf
and Save an extra 15% before 5/3. http://p.sf.net/sfu/velocityconf
_______________________________________________
QuantLib-dev mailing list
QuantLib-dev@...
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Hallo Dima,

 

I am very curious about your new kernel interpolation. Sounds like this is something that can help me a lot right now. I would like to try around a bit with different interpolation methods to overcome my numerical problems. If I can make my surface smoother with kernel interpolation I will give it a go and would like to do some testing on it.

 

How do I get your new kernel interpolation running?  Can I check it out from the SVN?  I checked the trunk and I found a file called kernelinterpolation.hpp.   Does it work for 2-dimensions or just for one because for the surface I need it for two dimensions. Can you please also send me a link or a paper with more informations about it so I can build up some theoretical knowledge before I start testing.

 

Greetings,

Michael

 

---

From: Dima [mailto:dimathematician@...]
Sent: Mittwoch, 29. April 2009 11:55
To: Klaus Spanderen
Cc: quantlib-dev@...
Subject: Re: [Quantlib-dev] LocalvolSurface.cpp

 

 

Just a remark regarding interpolation. I've implemented kernel

interpolation (its in the trunk), which can be made sufficiently smooth 

by choosing a proper standard deviation of the gaussian kernel. I've heard 

that this smoothness property makes it a good choice for local vol calibrations.

I don't have any personal experience with that though.

 

 

 

 

------------------------------------------------------------------------------
Register Now & Save for Velocity, the Web Performance & Operations
Conference from O'Reilly Media. Velocity features a full day of
expert-led, hands-on workshops and two days of sessions from industry
leaders in dedicated Performance & Operations tracks. Use code vel09scf
and Save an extra 15% before 5/3. http://p.sf.net/sfu/velocityconf
_______________________________________________
QuantLib-dev mailing list
QuantLib-dev@...
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Hi Michael,

> Anyways, can you send me a link or a paper with more information about the
> arbitrage violations due to constant extrapolation? I still cant see why
> constant extrapolation is violating the arbitrage criteria.

Please find attached a small program, where the constant extrapolation as
implemented in BlackVarianceSurface generates an arbitrage violation  -
negative call spread price when the maturity becomes large enough. To get it
running you have to enable extrapolation in analyticeuopeanengine.hpp at line
45. (Hope I got everything right with the example;-)
 
regards
 Klaus

[attachment removed]
------------------------------------------------------------------------------
Register Now & Save for Velocity, the Web Performance & Operations
Conference from O'Reilly Media. Velocity features a full day of
expert-led, hands-on workshops and two days of sessions from industry
leaders in dedicated Performance & Operations tracks. Use code vel09scf
and Save an extra 15% before 5/3. http://p.sf.net/sfu/velocityconf
_______________________________________________
QuantLib-dev mailing list
QuantLib-dev@...
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Hi.
Was it any feedback from Dima?
I'm also curious about this issue with respect to extrapolation problem during local vol calculation.

By the way, regarding main topic - you should keep in mind, that Duperi formula is used in continuous case calculation, so dt in lattice method must be really short (for instance, 3000 per year). I use LocalVolSurface class in QL PDE framework and it produces good results (see also http://www.nabble.com/LocalVolSurface-class-to20896055.html#a20896055)

MH_quant wrote:

Hallo Dima,

I am very curious about your new kernel interpolation. Sounds like this is
something that can help me a lot right now. I would like to try around a bit
with different interpolation methods to overcome my numerical problems. If I
can make my surface smoother with kernel interpolation I will give it a go
and would like to do some testing on it.

 
How do I get your new kernel interpolation running?  Can I check it out from
the SVN?  I checked the trunk and I found a file called
kernelinterpolation.hpp.   Does it work for 2-dimensions or just for one
because for the surface I need it for two dimensions. Can you please also
send me a link or a paper with more informations about it so I can build up
some theoretical knowledge before I start testing.