« Return to Thread: LocalvolSurface.cpp
Re: LocalvolSurface.cpp
by MH_quant
::
Rate this Message:
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
qlSurface.jpg (1M) « Return to Thread: LocalvolSurface.cpp
| Free embeddable forum powered by Nabble | Forum Help |
