Pricing options using the Heston model using QuantLib

./HestonEx --help
Example of pricing a simple Heston model using QuantLib

Option duration is fixed at six months, other parameters can be specified as listed below
Allowed options:
  --help                Produce this help message
  --S arg (=1)          Spot price
  --K arg (=1.05)       Strike price
  --v0 arg (=0.1)       Spot variance
  --kappa arg (=3.16)   Damping of mean reversion of the variance
  --theta arg (=0.09)   Mean reversion level of the variance
  --sigma arg (=0.4)    Volatility of the volatility
  --rho arg (=-0.2)     Correlation between spot and volatility processes

./HestonEx --theta 0.5
Option value: 0.130476

/**
   Copyright (C) 2009 Bojan Nikolic <bojan@bnikolic.co.uk>

   This program is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
   FOR A PARTICULAR PURPOSE.

   A bare-bones example of using the Heston Model in QuantLib. See
   also test-suite/hestonmodel.cpp in QuantLib

   Comments welcome at bojan@bnikolic.co.uk
*/
#include <iostream>

#include <boost/program_options.hpp>

#include <ql/instruments/dividendvanillaoption.hpp>
#include <ql/processes/hestonprocess.hpp>
#include <ql/models/equity/hestonmodel.hpp>
#include <ql/time/daycounters/actualactual.hpp>
#include <ql/time/period.hpp>
#include <ql/termstructures/yield/flatforward.hpp>
#include <ql/termstructures/yield/zerocurve.hpp>
#include <ql/pricingengines/vanilla/analytichestonengine.hpp>
#include <ql/time/calendars/nullcalendar.hpp>
#include <ql/quotes/simplequote.hpp>

using namespace QuantLib;

namespace po = boost::program_options;

boost::shared_ptr<YieldTermStructure>
flatRate(const Date& today,
         const boost::shared_ptr<Quote>& forward,
         const DayCounter &dc)
{
  return boost::shared_ptr<YieldTermStructure>
    (new FlatForward(today,
                     Handle<Quote>(forward),
                     dc));
}

boost::shared_ptr<YieldTermStructure>
flatRate(const boost::shared_ptr<Quote>& forward,
         const DayCounter &dc)
{
  return boost::shared_ptr<YieldTermStructure>
    (new FlatForward(0,
                     NullCalendar(),
                     Handle<Quote>(forward),
                     dc));
}

boost::shared_ptr<YieldTermStructure>
flatRate(Rate forward,
         const DayCounter &dc)
{
  return flatRate(boost::shared_ptr<Quote>(new SimpleQuote(forward)),
                  dc);
}

void Price(const po::variables_map & vm)
{

  Date settlementDate(1, January, 2009);
  Date exerciseDate(1, July, 2009);
  Settings::instance().evaluationDate() = settlementDate;
  DayCounter dayCounter = ActualActual();


  boost::shared_ptr<StrikedTypePayoff>
    payoff(new PlainVanillaPayoff(Option::Call,
                                  vm["K"].as<double>()));
  boost::shared_ptr<Exercise>
    exercise(new EuropeanExercise(exerciseDate));

  Handle<YieldTermStructure> riskFreeTS(flatRate(0.0225,
                                                 dayCounter));
  Handle<YieldTermStructure> dividendTS(flatRate(0.02,
                                                 dayCounter));

  Handle<Quote>
    s0(boost::shared_ptr<Quote>(new SimpleQuote(vm["S"].as<double>())));

  boost::shared_ptr<HestonProcess>
    process(new HestonProcess(riskFreeTS,
                              dividendTS,
                              s0,
                              vm["v0"].as<double>(),
                              vm["kappa"].as<double>(),
                              vm["theta"].as<double>(),
                              vm["sigma"].as<double>(),
                              vm["rho"].as<double>()));

  VanillaOption option(payoff, exercise);

  boost::shared_ptr<AnalyticHestonEngine>
    engine(new AnalyticHestonEngine(boost::shared_ptr<HestonModel>(new HestonModel(process)),
                                    64));

  option.setPricingEngine(engine);

  double val = option.NPV();

  std::cout<<"Option value: "<<val
           <<std::endl;

}

int main(int ac, char* av[])
{

  po::options_description desc("Allowed options");

  desc.add_options()
    ("help", "Produce this help message")
    ("S",
     po::value<double>()->default_value(1.0),
     "Spot price")
    ("K",
     po::value<double>()->default_value(1.05),
     "Strike price")
    ("v0",
     po::value<double>()->default_value(0.1),
     "Spot variance")
    ("kappa",
     po::value<double>()->default_value(3.16),
     "Damping of mean reversion of the variance"
     )
    ("theta",
     po::value<double>()->default_value(0.09),
     "Mean reversion level of the variance"
     )
    ("sigma",
     po::value<double>()->default_value(0.4),
     "Volatility of the volatility")
    ("rho",
     po::value<double>()->default_value(-0.2),
     "Correlation between spot and volatility processes");


  po::variables_map vm;
  po::store(po::parse_command_line(ac, av, desc),
            vm);
  po::notify(vm);

  if (vm.count("help"))
  {
    std::cout<<"Example of pricing a simple Heston model using QuantLib"
             <<std::endl
             <<std::endl
             <<"Option duration is fixed at six months, other parameters can be specified as listed below"
             <<std::endl
             <<desc;
  }
  else
  {
    // Real work
    Price(vm);
  }

  return 0;
}