QuantLib: Black / BSM processes and pricing via volatility surface. Different results?

Question

I start this question with a couple of C++ functions that will be useful to show some results. So start your Visual Studio C++ Express or Ceemple or whatever you want and copy & paste this:

#include <ql/quantlib.hpp>
#include <boost/timer.hpp>
#include <iostream>
#include <iomanip>

using namespace QuantLib;

#if defined(QL_ENABLE_SESSIONS)
namespace QuantLib {

    Integer sessionId() { return 0; }

}
#endif

After standard introduction, the first function acts like a small wrapper: taken a shared_ptr of BlackVolTermStructure template and some data, a zero-flat risk free rate curve is built and a BlackProcess is built; hence the option NPV is returned.

double EurVanillaSurfacePricerBlack(boost::shared_ptr<BlackVolTermStructure> forwardVolSurface, Option::Type type, Real underlying, Real strike, Date maturity)
{
  Rate riskFreeRate = 0.00;
  DayCounter dayCounter = Actual365Fixed();
  Calendar calendar = TARGET();
  Natural settlementDays = 3;

  // exercise
  boost::shared_ptr<Exercise> europeanExercise(
      new EuropeanExercise(maturity));

  // underlying
  Handle<Quote> underlyingH(boost::shared_ptr<Quote>(
      new SimpleQuote(underlying)));

  // bootstrap the yield curve
  Handle<YieldTermStructure> flatTermStructure(
      boost::shared_ptr<YieldTermStructure>(
          new FlatForward(settlementDays, calendar, riskFreeRate, dayCounter)));

  // payoff
  boost::shared_ptr<StrikedTypePayoff> payoff(
      new PlainVanillaPayoff(type, strike));

  // process
  boost::shared_ptr<BlackProcess> blackProcess(
      new BlackProcess(underlyingH, flatTermStructure, Handle<BlackVolTermStructure>(forwardVolSurface)));

  // options
  VanillaOption europeanOption(payoff, europeanExercise);
  europeanOption.setPricingEngine(boost::shared_ptr<PricingEngine>(
                                       new AnalyticEuropeanEngine(blackProcess)));
  double optionValue = europeanOption.NPV();

  return(optionValue);
}

The same function can be written using BlackScholesMertonProcess instead of BlackProcess and of course specifying a dividend yield term structure due to we're not using anymore underlying forward price... but here we set everything to zero, both risk free rate and dividend yield - thus giving the function flat and meaningless term structures:

double EurVanillaSurfacePricerBSM(boost::shared_ptr<BlackVolTermStructure> forwardVolSurface, Option::Type type, Real underlying, Real strike, Date maturity)
{
  Spread dividendYield = 0.00;
  Rate riskFreeRate = 0.00;
  DayCounter dayCounter = Actual365Fixed();
  Calendar calendar = TARGET();
  Natural settlementDays = 3;

  // exercise
  boost::shared_ptr<Exercise> europeanExercise(
      new EuropeanExercise(maturity));

  // underlying
  Handle<Quote> underlyingH(boost::shared_ptr<Quote>(
      new SimpleQuote(underlying)));

  // bootstrap the yield curve and the dividend curve
  Handle<YieldTermStructure> flatTermStructure(
      boost::shared_ptr<YieldTermStructure>(
          new FlatForward(settlementDays, calendar, riskFreeRate, dayCounter)));
  Handle<YieldTermStructure> flatDividendTS(
      boost::shared_ptr<YieldTermStructure>(
          new FlatForward(settlementDays, calendar, dividendYield, dayCounter)));


  // payoff
  boost::shared_ptr<StrikedTypePayoff> payoff(
      new PlainVanillaPayoff(type, strike));

  // process
  boost::shared_ptr<BlackScholesMertonProcess> bsmProcess(
      new BlackScholesMertonProcess(underlyingH, flatDividendTS, flatTermStructure, Handle<BlackVolTermStructure>(forwardVolSurface)));

  // options
  VanillaOption europeanOption(payoff, europeanExercise);
  europeanOption.setPricingEngine(boost::shared_ptr<PricingEngine>(
                                       new AnalyticEuropeanEngine(bsmProcess)));
  double optionValue = europeanOption.NPV();

  return(optionValue);
}

According to the little of option theory I know, there shouldn't be any difference between the two functions: $$F(t)=S(0)e^{[r(t)-q(t)]t},$$ where $q$ and $r$ are zero for every maturity thus $F=S(0)$ for every maturity.

The third function is a little variation on the same theme: instead of using standard constructors which require term structures, we plug a constant volatility value into a flat surface:

double EurVanillaPricer(Volatility volatility, Option::Type type, Real underlying, Real strike, Date maturity)
{
  Spread dividendYield = 0.00;
  Rate riskFreeRate = 0.00;
  DayCounter dayCounter = Actual365Fixed();
  Calendar calendar = TARGET();
  // This was "Natural settlementDays = 3;" before Luigi Ballabio's correction
  Natural settlementDays = 0;

  // exercise
  boost::shared_ptr<Exercise> europeanExercise(
      new EuropeanExercise(maturity));

  // underlying
  Handle<Quote> underlyingH(boost::shared_ptr<Quote>(
      new SimpleQuote(underlying)));

  // bootstrap the yield/dividend/vol curves
  Handle<YieldTermStructure> flatTermStructure(
      boost::shared_ptr<YieldTermStructure>(
          new FlatForward(settlementDays, calendar, riskFreeRate, dayCounter)));
  Handle<YieldTermStructure> flatDividendTS(
      boost::shared_ptr<YieldTermStructure>(
          new FlatForward(settlementDays, calendar, dividendYield, dayCounter)));
  Handle<BlackVolTermStructure> flatVolTS(
      boost::shared_ptr<BlackVolTermStructure>(
          new BlackConstantVol(settlementDays, calendar, volatility, dayCounter)));

  // payoff
  boost::shared_ptr<StrikedTypePayoff> payoff(
      new PlainVanillaPayoff(type, strike));

  // process
  boost::shared_ptr<BlackScholesMertonProcess> bsmProcess(
      new BlackScholesMertonProcess(underlyingH, flatDividendTS, flatTermStructure, flatVolTS));

  // options
  VanillaOption europeanOption(payoff, europeanExercise);
  europeanOption.setPricingEngine(boost::shared_ptr<PricingEngine>(
                                       new AnalyticEuropeanEngine(bsmProcess)));
  double optionValue = europeanOption.NPV();

  return(optionValue);
}

Last but not least, let me introduce a forward volatility surface wrapper:

boost::shared_ptr<BlackVolTermStructure> ForwardImpliedVolSurface(Date todaysDate, Date forwardDate, Calendar calendar, std::vector<Date> maturityArray, std::vector<Real> strikeArray, Matrix volatilityMatrix)
{
  // Handle to boost::shared_ptr
  DayCounter dayCounter = Actual365Fixed();
  boost::shared_ptr<BlackVarianceSurface> volatilitySurface(new BlackVarianceSurface(todaysDate, calendar, maturityArray, strikeArray, volatilityMatrix, dayCounter));
  Handle<BlackVolTermStructure> volatilitySurfaceH(volatilitySurface);

  // Volatility surface interpolation
  volatilitySurface->enableExtrapolation(true);

  // Change interpolator to bicubic splines
  volatilitySurface->setInterpolation<Bicubic>(Bicubic());

  // Forward implied volatility surface
  boost::shared_ptr<BlackVolTermStructure> forwardVolSurface(new ImpliedVolTermStructure(volatilitySurfaceH, forwardDate));

  return(forwardVolSurface);
}

What are the output of such functions? Let'see:

int main() {

try {

    boost::timer timer;
    std::cout << std::endl;

    /* +---------------------------------------------------------------------------------------------------
     * | Date and calendars parameters
     * +---------------------------------------------------------------------------------------------------
     * */

    // set up dates
    Calendar calendar = TARGET();
    Date todaysDate(03, Jul, 2014);
    Date settlementDate = calendar.advance(todaysDate, 3, Days);
    Settings::instance().evaluationDate() = todaysDate;

    // Maturity dates array
    Date expiry1(15, Aug, 2014);
    Date expiry2(19, Sep, 2014);
    Date expiry3(19, Dec, 2014);
    Date expiry4(20, Mar, 2015);
    Date expiry5(19, Jun, 2015);

    std::vector<Date> maturityArray;
    maturityArray.push_back(expiry1);
    maturityArray.push_back(expiry2);
    maturityArray.push_back(expiry3);
    maturityArray.push_back(expiry4);
    maturityArray.push_back(expiry5);

    // Strikes array
    std::vector<Real> strikeArray;
    for(int i = 2975; i < 2975 + (26 * 25); i = i + 25)
    {
      strikeArray.push_back(i);
    }

    // Implied volatility matrix
        Matrix volatilityMatrix(26, 5);

        volatilityMatrix[0][0]  = 0.198989  ; volatilityMatrix[0][1]  = 0.182889 ; volatilityMatrix[0][2]  = 0.182256 ; volatilityMatrix[0][3]  = 0.183319 ; volatilityMatrix[0][4]  = 0.202197 ;
        volatilityMatrix[1][0]  = 0.192338  ; volatilityMatrix[1][1]  = 0.178463 ; volatilityMatrix[1][2]  = 0.17982  ; volatilityMatrix[1][3]  = 0.181494 ; volatilityMatrix[1][4]  = 0.201261 ;
        volatilityMatrix[2][0]  = 0.185184  ; volatilityMatrix[2][1]  = 0.174239 ; volatilityMatrix[2][2]  = 0.177315 ; volatilityMatrix[2][3]  = 0.179669 ; volatilityMatrix[2][4]  = 0.200291 ;
        volatilityMatrix[3][0]  = 0.178718  ; volatilityMatrix[3][1]  = 0.170046 ; volatilityMatrix[3][2]  = 0.175143 ; volatilityMatrix[3][3]  = 0.177845 ; volatilityMatrix[3][4]  = 0.19928  ;
        volatilityMatrix[4][0]  = 0.172647  ; volatilityMatrix[4][1]  = 0.166123 ; volatilityMatrix[4][2]  = 0.172826 ; volatilityMatrix[4][3]  = 0.176046 ; volatilityMatrix[4][4]  = 0.198271 ;
        volatilityMatrix[5][0]  = 0.166556  ; volatilityMatrix[5][1]  = 0.162275 ; volatilityMatrix[5][2]  = 0.170328 ; volatilityMatrix[5][3]  = 0.174391 ; volatilityMatrix[5][4]  = 0.19764  ;
        volatilityMatrix[6][0]  = 0.160933  ; volatilityMatrix[6][1]  = 0.158344 ; volatilityMatrix[6][2]  = 0.16825  ; volatilityMatrix[6][3]  = 0.172892 ; volatilityMatrix[6][4]  = 0.197454 ;
        volatilityMatrix[7][0]  = 0.155747  ; volatilityMatrix[7][1]  = 0.154688 ; volatilityMatrix[7][2]  = 0.166199 ; volatilityMatrix[7][3]  = 0.17105  ; volatilityMatrix[7][4]  = 0.196211 ;
        volatilityMatrix[8][0]  = 0.150464  ; volatilityMatrix[8][1]  = 0.151097 ; volatilityMatrix[8][2]  = 0.164325 ; volatilityMatrix[8][3]  = 0.16875  ; volatilityMatrix[8][4]  = 0.193533 ;
        volatilityMatrix[9][0]  = 0.145234  ; volatilityMatrix[9][1]  = 0.147602 ; volatilityMatrix[9][2]  = 0.16217  ; volatilityMatrix[9][3]  = 0.16793  ; volatilityMatrix[9][4]  = 0.195104 ;
        volatilityMatrix[10][0] = 0.140751  ; volatilityMatrix[10][1] = 0.144357 ; volatilityMatrix[10][2] = 0.160261 ; volatilityMatrix[10][3] = 0.169107 ; volatilityMatrix[10][4] = 0.202441 ;
        volatilityMatrix[11][0] = 0.136502  ; volatilityMatrix[11][1] = 0.141208 ; volatilityMatrix[11][2] = 0.158546 ; volatilityMatrix[11][3] = 0.165058 ; volatilityMatrix[11][4] = 0.194346 ;
        volatilityMatrix[12][0] = 0.13342   ; volatilityMatrix[12][1] = 0.138357 ; volatilityMatrix[12][2] = 0.156949 ; volatilityMatrix[12][3] = 0.15057  ; volatilityMatrix[12][4] = 0.155503 ;
        volatilityMatrix[13][0] = 0.104896  ; volatilityMatrix[13][1] = 0.119273 ; volatilityMatrix[13][2] = 0.128517 ; volatilityMatrix[13][3] = 0.136208 ; volatilityMatrix[13][4] = 0.116855 ;
        volatilityMatrix[14][0] = 0.10099   ; volatilityMatrix[14][1] = 0.115047 ; volatilityMatrix[14][2] = 0.125638 ; volatilityMatrix[14][3] = 0.132476 ; volatilityMatrix[14][4] = 0.109273 ;
        volatilityMatrix[15][0] = 0.100313  ; volatilityMatrix[15][1] = 0.114395 ; volatilityMatrix[15][2] = 0.125642 ; volatilityMatrix[15][3] = 0.133834 ; volatilityMatrix[15][4] = 0.117099 ;
        volatilityMatrix[16][0] = 0.0981065 ; volatilityMatrix[16][1] = 0.112273 ; volatilityMatrix[16][2] = 0.124137 ; volatilityMatrix[16][3] = 0.132863 ; volatilityMatrix[16][4] = 0.118885 ;
        volatilityMatrix[17][0] = 0.0962976 ; volatilityMatrix[17][1] = 0.109955 ; volatilityMatrix[17][2] = 0.122498 ; volatilityMatrix[17][3] = 0.130647 ; volatilityMatrix[17][4] = 0.116549 ;
        volatilityMatrix[18][0] = 0.0950343 ; volatilityMatrix[18][1] = 0.107924 ; volatilityMatrix[18][2] = 0.121311 ; volatilityMatrix[18][3] = 0.129627 ; volatilityMatrix[18][4] = 0.116142 ;
        volatilityMatrix[19][0] = 0.094729  ; volatilityMatrix[19][1] = 0.106211 ; volatilityMatrix[19][2] = 0.119952 ; volatilityMatrix[19][3] = 0.12918  ; volatilityMatrix[19][4] = 0.116873 ;
        volatilityMatrix[20][0] = 0.0952533 ; volatilityMatrix[20][1] = 0.104712 ; volatilityMatrix[20][2] = 0.118585 ; volatilityMatrix[20][3] = 0.128231 ; volatilityMatrix[20][4] = 0.116804 ;
        volatilityMatrix[21][0] = 0.0977423 ; volatilityMatrix[21][1] = 0.103553 ; volatilityMatrix[21][2] = 0.117229 ; volatilityMatrix[21][3] = 0.126978 ; volatilityMatrix[21][4] = 0.116249 ;
        volatilityMatrix[22][0] = 0.0992171 ; volatilityMatrix[22][1] = 0.102743 ; volatilityMatrix[22][2] = 0.115987 ; volatilityMatrix[22][3] = 0.125834 ; volatilityMatrix[22][4] = 0.115905 ;
        volatilityMatrix[23][0] = 0.102137  ; volatilityMatrix[23][1] = 0.1025   ; volatilityMatrix[23][2] = 0.114716 ; volatilityMatrix[23][3] = 0.124794 ; volatilityMatrix[23][4] = 0.115759 ;
        volatilityMatrix[24][0] = 0.108426  ; volatilityMatrix[24][1] = 0.102351 ; volatilityMatrix[24][2] = 0.113496 ; volatilityMatrix[24][3] = 0.123768 ; volatilityMatrix[24][4] = 0.115648 ;
        volatilityMatrix[25][0] = 0.111779  ; volatilityMatrix[25][1] = 0.102869 ; volatilityMatrix[25][2] = 0.112514 ; volatilityMatrix[25][3] = 0.12274  ; volatilityMatrix[25][4] = 0.11554  ;

/* +---------------------------------------------------------------------------------------------------
 * | Forward volatility (ref. pag. 154-157 of "Dynamic Hedging - Managing Vanilla and Exotic Options")
 * +---------------------------------------------------------------------------------------------------
 * */

    // As instance, go 15 days forward
    Date forwardDate = calendar.advance(todaysDate, 15, Days);

    boost::shared_ptr<BlackVolTermStructure> forwardVolSurface = ForwardImpliedVolSurface(todaysDate, forwardDate, calendar, maturityArray, strikeArray, volatilityMatrix);
    Option::Type typeCall(Option::Call);
    Option::Type typePut(Option::Put);

    Real underlying = 3289.75;

    double myOption4;
    myOption4 = EurVanillaSurfacePricerBlack(forwardVolSurface, typeCall, underlying, 3300, expiry3);
    //disp(myOption4);
    double myOption5;
    myOption5 = EurVanillaSurfacePricerBSM(forwardVolSurface, typeCall, underlying, 3300, expiry3);
    //disp(myOption5);
    double myOption6;
    myOption6 = EurVanillaPricer(forwardVolSurface->blackVol(expiry3, 3300), typeCall, underlying, 3300, expiry3);
    //disp(myOption6);

    // Amend evaluation date...
    Settings::instance().evaluationDate() = forwardDate;

    double myOption1;
    myOption1 = EurVanillaSurfacePricerBlack(forwardVolSurface, typeCall, underlying, 3300, expiry3);
    //disp(myOption1);
    double myOption2;
    myOption2 = EurVanillaSurfacePricerBSM(forwardVolSurface, typeCall, underlying, 3300, expiry3);
    //disp(myOption2);
    double myOption3;
    myOption3 = EurVanillaPricer(forwardVolSurface->blackVol(expiry3, 3300), typeCall, underlying, 3300, expiry3);
    //disp(myOption3);

        return 0;

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

What are the output? Well, of course Black model and BSM model with flat and null term structures return the same values, that are $105.743$ for both the options. Despite of this, however, selecting "manually" the implied volatility to be used from the implied volatility surface via forwardVolSurface->blackVol(expiry3, 3300) returns a very different value, that is, $103.858$.

How would you explain this difference?

What I am afraid of is... automatic shift of volatility surface when evaluation date is amended. I tried to write code to avoid such behavior retaining relinkable handles into functions, so that they're destroyed once function ends.

But I am not sure it works as intended.

Answer 1

It's because of the settlement days you passed when you initialized the flat volatility curve. You're creating the spot, forward and flat volatilities as:

boost::shared_ptr<BlackVarianceSurface> volatilitySurface(
    new BlackVarianceSurface(todaysDate, calendar,
                             maturityArray, strikeArray,
                             volatilityMatrix, dayCounter));

boost::shared_ptr<BlackVolTermStructure> forwardVolSurface(
    new ImpliedVolTermStructure(volatilitySurfaceH, forwardDate));

boost::shared_ptr<BlackVolTermStructure> flatVolCurve(
      new BlackConstantVol(settlementDays, calendar,
                           forwardVolSurface->blackVol(expiry3, 3300),
                           dayCounter)));

(the first two are verbatim from your code; in the last one, I joined your call to retrieve the volatility in main to the constructor call in EurVanillaPricer). In the last one, you've set settlementDays to 3.

Now, the forward and flat curve have the same volatility for the exercise and strike you specified; if you execute

std::cout << forwardVolSurface->blackVol(expiry3, 3300) << std::endl;
std::cout << flatVolTS->blackVol(expiry3, 3300) << std::endl;

you'll get back

0.132405
0.132405

that is, the same. But that's not the whole story, unfortunately. The pricing engine for the European option doesn't ask them for the volatility , but directly for the variance ($\sigma^2 T$). If you do that:

std::cout << forwardVolSurface->blackVariance(expiry3, 3300) << std::endl;
std::cout << flatVolTS->blackVariance(expiry3, 3300) << std::endl;

you'll get back different figures:

0.00710851
0.00686836

Why? Because the two curves have different reference dates:

std::cout << forwardVolSurface->referenceDate() << std::endl;
std::cout << flatVolTS->referenceDate() << std::endl;

gives

July 24th, 2014
July 29th, 2014

In the constructor of forwardVolSurface, you're passing the reference date directly: it's forwardDate, which you also set as the evaluation date. In the constructor of flatVolTS, you're passing settlementDays and calendar, which equal 3 and TARGET(), respectively. This means that the reference date is three business days after the evaluation date, which turns out to be 5 days after when you factor in the weekend.

Thus, when the two curves calculate the variance $\sigma^2 T$, the volatility $\sigma$ is the same; but for the first curve, $T$ is the time between July 24th and the expiry, while for the second curve $T$ is the time between July 29th and the expiry and is 5 days shorter.

To solve your problem, just pass 0 as settlementDays when you create the flat curve. It will have the same reference date as the forward curve, and the option prices will be the same. (This doesn't need to be the same as the settlement days you pass to the risk-free and dividend curves.)

Note 1: In general, it might not be a good idea to pass a non-null number of settlement days to a volatility term structure (I hope that's not too many negatives) (darn, I did it again). The possibility to have settlement days is there because it was inherited from the base TermStructure class and we didn't think to restrict it, but it was added mostly for interest-rate curves where one might want to work spot (i.e., two days from today, which is the settlement date of most quoted interest-rate instruments). For equities, you probably want to start from the evaluation date anyway.

Note 2: the various ways to initialize term structures and how they work are explained in more detail at http://www.implementingquantlib.com/2013/09/chapter-3-part-1-of-n-term-structures.html.