Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added moment matching test case for the Kluge model #728

Merged
merged 3 commits into from Dec 23, 2019
Merged
Show file tree
Hide file tree
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
Expand Up @@ -43,9 +43,11 @@ namespace QuantLib {
const ext::shared_ptr<ExtOUWithJumpsProcess>& process,
const ext::shared_ptr<YieldTermStructure>& rTS,
Size tGrid, Size xGrid, Size yGrid,
const ext::shared_ptr<Shape>& shape,
const FdmSchemeDesc& schemeDesc)
: process_(process),
rTS_(rTS),
shape_(shape),
tGrid_(tGrid),
xGrid_(xGrid),
yGrid_(yGrid),
Expand All @@ -58,9 +60,9 @@ namespace QuantLib {
= rTS_->dayCounter().yearFraction(rTS_->referenceDate(),
arguments_.exercise->lastDate());
const ext::shared_ptr<StochasticProcess1D> ouProcess(
process_->getExtendedOrnsteinUhlenbeckProcess());
process_->getExtendedOrnsteinUhlenbeckProcess());
const ext::shared_ptr<Fdm1dMesher> xMesher(
new FdmSimpleProcess1dMesher(xGrid_, ouProcess,maturity));
new FdmSimpleProcess1dMesher(xGrid_, ouProcess,maturity));

const ext::shared_ptr<Fdm1dMesher> yMesher(
new ExponentialJump1dMesher(yGrid_,
Expand All @@ -73,7 +75,7 @@ namespace QuantLib {

// 2. Calculator
const ext::shared_ptr<FdmInnerValueCalculator> calculator(
new FdmExtOUJumpModelInnerValue(arguments_.payoff, mesher));
new FdmExtOUJumpModelInnerValue(arguments_.payoff, mesher, shape_));

// 3. Step conditions
const ext::shared_ptr<FdmStepConditionComposite> conditions =
Expand Down
14 changes: 9 additions & 5 deletions ql/experimental/finitedifferences/fdextoujumpvanillaengine.hpp
Expand Up @@ -28,6 +28,7 @@
#include <ql/pricingengine.hpp>
#include <ql/instruments/vanillaoption.hpp>
#include <ql/methods/finitedifferences/solvers/fdmbackwardsolver.hpp>
#include <ql/experimental/finitedifferences/fdmextoujumpmodelinnervalue.hpp>

namespace QuantLib {

Expand All @@ -38,17 +39,20 @@ namespace QuantLib {
: public GenericEngine<VanillaOption::arguments,
VanillaOption::results> {
public:
FdExtOUJumpVanillaEngine(
const ext::shared_ptr<ExtOUWithJumpsProcess>& p,
const ext::shared_ptr<YieldTermStructure>& rTS,
Size tGrid = 50, Size xGrid = 200, Size yGrid = 50,
const FdmSchemeDesc& schemeDesc=FdmSchemeDesc::Hundsdorfer());
typedef FdmExtOUJumpModelInnerValue::Shape Shape;
FdExtOUJumpVanillaEngine(
const ext::shared_ptr<ExtOUWithJumpsProcess>& p,
const ext::shared_ptr<YieldTermStructure>& rTS,
Size tGrid = 50, Size xGrid = 200, Size yGrid = 50,
const ext::shared_ptr<Shape>& shape = ext::shared_ptr<Shape>(),
const FdmSchemeDesc& schemeDesc=FdmSchemeDesc::Hundsdorfer());

void calculate() const;

private:
const ext::shared_ptr<ExtOUWithJumpsProcess> process_;
const ext::shared_ptr<YieldTermStructure> rTS_;
const ext::shared_ptr<Shape> shape_;
const Size tGrid_, xGrid_, yGrid_;
const FdmSchemeDesc schemeDesc_;
};
Expand Down
10 changes: 7 additions & 3 deletions ql/methods/finitedifferences/meshers/exponentialjump1dmesher.cpp
Expand Up @@ -28,6 +28,8 @@
#include <ql/methods/finitedifferences/meshers/exponentialjump1dmesher.hpp>
#include <ql/functional.hpp>

#include <iostream>
klausspanderen marked this conversation as resolved.
Show resolved Hide resolved

namespace QuantLib {
ExponentialJump1dMesher::ExponentialJump1dMesher(
Size steps, Real beta, Real jumpIntensity, Real eta, Real eps)
Expand All @@ -38,13 +40,15 @@ namespace QuantLib {
QL_REQUIRE(steps > 1, "minimum number of steps is two");

const Real start = 0.0;
const Real end = 1.0-eps;
const Real end = 1.0-eps;
const Real dx = (end-start)/(steps-1);

const Real scale = 1/(1-std::exp(-beta/jumpIntensity));

for (Size i=0; i < steps; ++i) {
const Real p = start + i*dx;
locations_[i] = -1.0/eta*std::log(1.0-p);
locations_[i] = scale*(-1.0/eta*std::log(1.0-p));
}

for (Size i=0; i < steps-1; ++i) {
dminus_[i+1] = dplus_[i] = locations_[i+1]-locations_[i];
}
Expand Down
165 changes: 163 additions & 2 deletions test-suite/swingoption.cpp
Expand Up @@ -20,7 +20,10 @@
#include "swingoption.hpp"
#include "utilities.hpp"

#include <ql/math/factorial.hpp>
#include <ql/math/functional.hpp>
#include <ql/math/richardsonextrapolation.hpp>
#include <ql/math/distributions/normaldistribution.hpp>
#include <ql/quotes/simplequote.hpp>
#include <ql/time/daycounters/actualactual.hpp>
#include <ql/instruments/vanillaoption.hpp>
Expand All @@ -34,6 +37,7 @@
#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
#include <ql/experimental/processes/extouwithjumpsprocess.hpp>
#include <ql/experimental/processes/extendedornsteinuhlenbeckprocess.hpp>
#include <ql/pricingengines/blackformula.hpp>
#include <ql/pricingengines/vanilla/fdsimplebsswingengine.hpp>
#include <ql/experimental/finitedifferences/fdextoujumpvanillaengine.hpp>
#include <ql/experimental/finitedifferences/fdsimpleextoujumpswingengine.hpp>
Expand All @@ -60,8 +64,7 @@ namespace {
new ExtendedOrnsteinUhlenbeckProcess(speed, volatility, x0[0],
constant<Real, Real>(x0[0])));
return ext::make_shared<ExtOUWithJumpsProcess>(
ouProcess, x0[1], beta,
jumpIntensity, eta);
ouProcess, x0[1], beta, jumpIntensity, eta);
}
}

Expand Down Expand Up @@ -434,6 +437,162 @@ void SwingOptionTest::testExtOUJumpSwingOption() {
}
}

namespace {
class SwingPdePricing {
public:
typedef FdSimpleExtOUJumpSwingEngine::Shape Shape;

SwingPdePricing(const ext::shared_ptr<ExtOUWithJumpsProcess>& process,
const ext::shared_ptr<VanillaOption>& option,
const ext::shared_ptr<Shape>& shape)
: process_(process), option_(option), shape_(shape) {}

Real operator()(Real x) const {
const ext::shared_ptr<YieldTermStructure> rTS(
flatRate(0.0, Actual365Fixed()));

const Size gridX = 200;
const Size gridY = 100;
const Size gridT = 100;

option_->setPricingEngine(
ext::make_shared<FdExtOUJumpVanillaEngine>(
process_, rTS,
Size(gridT)/x, Size(gridX/x), Size(gridY/x), shape_));

return option_->NPV();
}

private:
const ext::shared_ptr<ExtOUWithJumpsProcess> process_;
const ext::shared_ptr<VanillaOption> option_;
const ext::shared_ptr<Shape> shape_;
};
}

void SwingOptionTest::testKlugeChFVanillaPricing() {
BOOST_TEST_MESSAGE("Testing Kluge PDE Vanilla Pricing in"
" comparison to moment matching...");

SavedSettings backup;

Date settlementDate = Date(22, November, 2019);
Settings::instance().evaluationDate() = settlementDate;
DayCounter dayCounter = Actual365Fixed();
Date maturityDate = settlementDate + Period(6, Months);
const Time t = dayCounter.yearFraction(settlementDate, maturityDate);

const Real f0 = 30;

const Real x0 = 0.0;
const Real y0 = 0.0;

const Real beta = 5.0;
const Real eta = 5.0;
const Real lambda = 4.0;
const Real alpha = 4.0;
const Real sig = 1.0;

const ext::shared_ptr<ExtOUWithJumpsProcess> klugeProcess =
ext::make_shared<ExtOUWithJumpsProcess>(
ext::make_shared<ExtendedOrnsteinUhlenbeckProcess>(
alpha, sig, x0, constant<Real, Real>(0.0)),
y0, beta, lambda, eta);

const Real strike = f0;

const ext::shared_ptr<VanillaOption> option =
ext::make_shared<VanillaOption>(
ext::make_shared<PlainVanillaPayoff>(Option::Call, strike),
ext::make_shared<EuropeanExercise>(maturityDate));

typedef FdSimpleExtOUJumpSwingEngine::Shape Shape;
const ext::shared_ptr<Shape> shape(ext::make_shared<Shape>());

const Real ps = std::log(f0)
- sig*sig/(4*alpha)*(1-std::exp(-2*alpha*t))
- lambda/beta*std::log((eta-std::exp(-beta*t))/(eta-1.0));

shape->push_back(Shape::value_type(t, ps));

const Real expected =
RichardsonExtrapolation(
SwingPdePricing(klugeProcess, option, shape), 4.0)(2.0, 1.5);

const Real stdDev = std::sqrt((((2 - 2*std::exp(-2*beta*t))*lambda)
/(beta*eta*eta) + ((1 - std::exp(-2*alpha*t))*sig*sig)/alpha)/2.);

const Real bsNPV = blackFormula(Option::Call, strike, f0, stdDev);

const Real g1 = ((2 - 2*std::exp(-3*beta*t))*lambda)/(beta*eta*eta*eta)
/ (stdDev*stdDev*stdDev);

const Real g2 = 3*(std::exp((alpha + beta)*t)
* square<Real>()(2*alpha*std::exp(2*alpha*t)*(-1 + std::exp(2*beta*t))
*lambda + beta*std::exp(2*beta*t)*(-1 + std::exp(2*alpha*t))
*eta*eta*sig*sig)
+ 16*alpha*alpha*beta*std::exp((5*alpha + 3*beta)*t)*lambda
*std::sinh(2*beta*t))
/ (4.*alpha*alpha*beta*beta
*std::exp(5*(alpha + beta)*t)*eta*eta*eta*eta)
/ (stdDev*stdDev*stdDev*stdDev) - 3.0;

const Real d = (std::log(f0/strike) + 0.5*stdDev*stdDev)/stdDev;

// Jurczenko E., Maillet B. and Negrea B.,
// Multi-Moment Approximate Option Pricing Models:
// A General Comparison (Part 1)
// https://papers.ssrn.com/sol3/papers.cfm?abstract_id=300922
const NormalDistribution n;
const Real q3 = 1/Factorial::get(3)*f0*stdDev*(2*stdDev - d)*n(d);
const Real q4 = 1/Factorial::get(4)*f0*stdDev*(d*d - 3*d*stdDev - 1)*n(d);
const Real q5 = 10/Factorial::get(6)*f0*stdDev*(
d*d*d*d - 5*d*d*d*stdDev - 6*d*d + 15*d*stdDev + 3)*n(d);

// Corrado C. and T. Su, (1996-b),
// “Skewness and Kurtosis in S&P 500 IndexReturns Implied by Option Prices”,
// Journal of Financial Research 19 (2), 175-192.
const Real ccs3 = bsNPV + g1*q3;
const Real ccs4 = ccs3 + g2*q4;

// Rubinstein M., (1998), “Edgeworth Binomial Trees”,
// Journal of Derivatives 5 (3), 20-27.
const Real cr = ccs4 + g1*g1*q5;

const Volatility expectedImplVol = blackFormulaImpliedStdDevLiRS(
Option::Call, strike, f0, expected, 1.0)/std::sqrt(t);

const Volatility bsImplVol = blackFormulaImpliedStdDevLiRS(
Option::Call, strike, f0, bsNPV, 1.0)/std::sqrt(t);

const Volatility ccs3ImplVol = blackFormulaImpliedStdDevLiRS(
Option::Call, strike, f0, ccs3, 1.0)/std::sqrt(t);

const Volatility ccs4ImplVol = blackFormulaImpliedStdDevLiRS(
Option::Call, strike, f0, ccs4, 1.0)/std::sqrt(t);

const Volatility crImplVol = blackFormulaImpliedStdDevLiRS(
Option::Call, strike, f0, cr, 1.0)/std::sqrt(t);

const Real tol[] = {0.01, 0.0075, 0.005, 0.004};
const std::string methods[] = {
"Second Order", "Third Order", "Fourth Order", "Rubinstein"};

const Real calculated[] = {bsImplVol, ccs3ImplVol, ccs4ImplVol, crImplVol};

for (Size i=0; i < 4; ++i) {
const Real diff = std::fabs(calculated[i] - expectedImplVol);
if (diff > tol[i]) {
BOOST_ERROR("failed to reproduce vanilla option implied volatility "
"with moment matching"
<< "\n calculated: " << calculated[i]
<< "\n expected: " << expectedImplVol
<< "\n difference: " << diff
<< "\n tolerance: " << tol[i]
<< "\n method: " << methods[i]);
}
}
}

test_suite* SwingOptionTest::suite(SpeedLevel speed) {
test_suite* suite = BOOST_TEST_SUITE("Swing-Option Test");
Expand All @@ -443,6 +602,8 @@ test_suite* SwingOptionTest::suite(SpeedLevel speed) {
suite->add(QUANTLIB_TEST_CASE(&SwingOptionTest::testFdBSSwingOption));
suite->add(QUANTLIB_TEST_CASE(
&SwingOptionTest::testFdmExponentialJump1dMesher));
suite->add(QUANTLIB_TEST_CASE(
&SwingOptionTest::testKlugeChFVanillaPricing));

if (speed <= Fast) {
suite->add(QUANTLIB_TEST_CASE(
Expand Down
1 change: 1 addition & 0 deletions test-suite/swingoption.hpp
Expand Up @@ -33,6 +33,7 @@ class SwingOptionTest {
static void testExtOUJumpSwingOption();
static void testFdmExponentialJump1dMesher();
static void testExtOUJumpVanillaEngine();
static void testKlugeChFVanillaPricing();
static boost::unit_test_framework::test_suite* suite(SpeedLevel);
};

Expand Down