QuantLib › quantlib-users

calibration G2++ with defferential evolution

Classic

List

Threaded

11 messages Options

André de Boer

calibration G2++ with defferential evolution

Hello,

Two days ago I posted a question about differential evolution and the shortratemodel G2.

Meanwhile I learned how to minimize the Rosenbrockfunction with differential evolution.

But my goal is to calibrate the G2 shortratemodel.

For the 5 parameters of the G2++ model I want to have boundedconstraints:

For instance:

0 < a < 6

0 < b < 6

0 < sigma < 0.25

0 < eta < 0.25

-1 < rho < 1

To accomplish this I made a class MyConstraint, see code below.

...

MyConstraint mc;

DifferentialEvolution deOptim(conf);

modelG2pp->calibrate(swaptions, deOptim, EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);

...

But still I isn't not working.

Can someone give me a hint of parts of code?

Kind regards,

André

#include <iostream>

#include <vector>

#include <ql/quantlib.hpp>

using namespace QuantLib;

class MyConstraint : public Constraint {

class Impl : public Constraint::Impl {

public:

bool test(const Array& x) const {

Real a1 = x[0], a2 = x[1], a3 = x[2], a4=x[3], a5=x[4];

return (0.0 <= a1 && a1 <= 6.0) &&
(0.0 <= a2 && a2 <= 6.0) &&
(0.0 <= a3 && a3 <= 0.25) &&
(0.0 <= a4 && a4 <= 0.25) &&
(-1.0 <= a5 && a5 <= 1.0);

}

};

public:

MyConstraint()
: Constraint(boost::shared_ptr<Constraint::Impl>(new MyConstraint::Impl)) {}

};

int main(){

...
double a = 0.1;

double b = 0.1;
double sigma = 0.01;
double eta = 0.01;
double rho = -0.75;
boost::shared_ptr<G2> modelG2pp(new G2(yieldCurve, a, b, sigma, eta, rho));
std::cout << "G2 (analytic formulae) calibration" << std::endl;

MyConstraint mc;
DifferentialEvolution::Configuration conf =
DifferentialEvolution::Configuration()
.withStepsizeWeight(0.4)
.withBounds()
.withCrossoverProbability(0.35)

.withPopulationMembers(100)
.withStrategy(DifferentialEvolution::BestMemberWithJitter)
.withCrossoverType(DifferentialEvolution::Normal)
.withAdaptiveCrossover()
.withSeed(0);

MyConstraint mc;
DifferentialEvolution deOptim(conf);
modelG2pp->calibrate(swaptions, deOptim, EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);

}

------------------------------------------------------------------------------
Android apps run on BlackBerry 10
Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
Now with support for Jelly Bean, Bluetooth, Mapview and more.
Get your Android app in front of a whole new audience. Start now.
http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
_______________________________________________
QuantLib-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-users

Peter Caspers-4

Re: calibration G2++ with defferential evolution

Hi André,

can you please post the complete code needed to reproduce the problem
? With a short description of what the problem is (if not obvious) ?

best regards
Peter

On 16 February 2014 17:28, André de Boer <[hidden email]> wrote:

> Hello,
>
> Two days ago I posted a question about differential evolution and the
> shortratemodel G2.
> Meanwhile I learned how to minimize the Rosenbrockfunction with differential
> evolution.
> But my goal is to calibrate the G2 shortratemodel.
>
> For the 5 parameters of the G2++ model I want to have boundedconstraints:
> For instance:
> 0 < a < 6
> 0 0 < sigma < 0.25
> 0 < eta < 0.25
> -1 < rho < 1
>
> To accomplish this I made a class MyConstraint, see code below.
> ...
> MyConstraint mc;
> DifferentialEvolution deOptim(conf);
> modelG2pp->calibrate(swaptions, deOptim,
> EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
> ...
> But still I isn't not working.
> Can someone give me a hint of parts of code?
>
> Kind regards,
> André
>
>
>
> #include <iostream>
> #include <vector>
> #include <ql/quantlib.hpp>
> using namespace QuantLib;
>
> class MyConstraint : public Constraint {
> class Impl : public Constraint::Impl {
> public:
>
> bool test(const Array& x) const {
> Real a1 = x[0], a2 = x[1], a3 = x[2], a4=x[3], a5=x[4];
>
> return (0.0 <= a1 && a1 <= 6.0) &&
> (0.0 <= a2 && a2 <= 6.0) &&
> (0.0 <= a3 && a3 <= 0.25) &&
> (0.0 <= a4 && a4 <= 0.25) &&
> (-1.0 <= a5 && a5 <= 1.0);
>
> }
>
> };
> public:
>
> MyConstraint()
> : Constraint(boost::shared_ptr<Constraint::Impl>(new MyConstraint::Impl)) {}
>
> };
>
> int main(){
>
> ...
> double a = 0.1;
>
> double b = 0.1;
> double sigma = 0.01;
> double eta = 0.01;
> double rho = -0.75;
> boost::shared_ptr<G2> modelG2pp(new G2(yieldCurve, a, b, sigma, eta, rho));
> std::cout << "G2 (analytic formulae) calibration" << std::endl;
>
> MyConstraint mc;
> DifferentialEvolution::Configuration conf =
> DifferentialEvolution::Configuration()
> .withStepsizeWeight(0.4)
> .withBounds()
> .withCrossoverProbability(0.35)
> .withPopulationMembers(100)
> .withStrategy(DifferentialEvolution::BestMemberWithJitter)
> .withCrossoverType(DifferentialEvolution::Normal)
> .withAdaptiveCrossover()
> .withSeed(0);
>
> MyConstraint mc;
> DifferentialEvolution deOptim(conf);
> modelG2pp->calibrate(swaptions, deOptim,
> EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>
> }
>
> ------------------------------------------------------------------------------
> Android apps run on BlackBerry 10
> Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
> Now with support for Jelly Bean, Bluetooth, Mapview and more.
> Get your Android app in front of a whole new audience. Start now.
> http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
> _______________________________________________
> QuantLib-users mailing list
> [hidden email]
> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>

Ralph Schreyer-3

Re: calibration G2++ with defferential evolution

Hi Andre, Peter,

I think the current implementation is buggy. You would do something like

Array lower(5); lower[0] = 0.0; lower[1] = 0.0; lower[2] = 0.0; lower[3] = -1.0; lower[4] = 0.0; Array upper(5); upper[0] = 10.0; upper[1] = 50.0; upper[2] = 10.0; upper[3] = 1.0; upper[4] = 10.0; model->calibrate(options, deOptimizer, EndCriteria(400, 40, 1.0e-8, 1.0e-8, Null<Real>()), NonhomogeneousBoundaryConstraint(lower, upper));

Then, in model.cpp, line 88 ff, the CompositeConstraint is built from the PrivateConstraint of the model and the NonhomogeneousBoundaryConstraint. When you look at the implementation of e.g. upperBound(...) in constraint.hpp, line 156 ff., you find

for (Size iter = 0; iter < c1ub.size(); iter++) { rtrnArray.at(iter) = std::min(c1ub.at(iter), c2ub.at(iter)); }

Because c1ub is the upperBound of the PrivateConstraint and c1ub.size() = 0, the DifferentialEvolution does not get its bounds and therefore e.g. DifferentialEvolution::fillInitialPopulation(...) fails.

Maybe one can change the CompositeConstraint upperBound and lowerBound implementations? Or is there a better place? Or isn't there a bug at all??

Best

Ralph

2014-02-17 13:00 GMT+01:00 Peter Caspers <[hidden email]>:

Hi André,

can you please post the complete code needed to reproduce the problem
? With a short description of what the problem is (if not obvious) ?

best regards
Peter

On 16 February 2014 17:28, André de Boer <[hidden email]> wrote:
> Hello,
>
> Two days ago I posted a question about differential evolution and the
> shortratemodel G2.
> Meanwhile I learned how to minimize the Rosenbrockfunction with differential
> evolution.
> But my goal is to calibrate the G2 shortratemodel.
>
> For the 5 parameters of the G2++ model I want to have boundedconstraints:
> For instance:
> 0 < a < 6
> 0 0 < sigma < 0.25
> 0 < eta < 0.25
> -1 < rho < 1
>
> To accomplish this I made a class MyConstraint, see code below.
> ...
> MyConstraint mc;
> DifferentialEvolution deOptim(conf);
> modelG2pp->calibrate(swaptions, deOptim,
> EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
> ...
> But still I isn't not working.
> Can someone give me a hint of parts of code?
>
> Kind regards,
> André
>
>
>
> #include <iostream>
> #include <vector>
> #include <ql/quantlib.hpp>
> using namespace QuantLib;
>
> class MyConstraint : public Constraint {
> class Impl : public Constraint::Impl {
> public:
>
> bool test(const Array& x) const {
> Real a1 = x[0], a2 = x[1], a3 = x[2], a4=x[3], a5=x[4];
>
> return (0.0 <= a1 && a1 <= 6.0) &&
> (0.0 <= a2 && a2 <= 6.0) &&
> (0.0 <= a3 && a3 <= 0.25) &&
> (0.0 <= a4 && a4 <= 0.25) &&
> (-1.0 <= a5 && a5 <= 1.0);
>
> }
>
> };
> public:
>
> MyConstraint()
> : Constraint(boost::shared_ptr<Constraint::Impl>(new MyConstraint::Impl)) {}
>
> };
>
> int main(){
>
> ...
> double a = 0.1;
>
> double b = 0.1;
> double sigma = 0.01;
> double eta = 0.01;
> double rho = -0.75;
> boost::shared_ptr<G2> modelG2pp(new G2(yieldCurve, a, b, sigma, eta, rho));
> std::cout << "G2 (analytic formulae) calibration" << std::endl;
>
> MyConstraint mc;
> DifferentialEvolution::Configuration conf =
> DifferentialEvolution::Configuration()
> .withStepsizeWeight(0.4)
> .withBounds()
> .withCrossoverProbability(0.35)
> .withPopulationMembers(100)
> .withStrategy(DifferentialEvolution::BestMemberWithJitter)
> .withCrossoverType(DifferentialEvolution::Normal)
> .withAdaptiveCrossover()
> .withSeed(0);
>
> MyConstraint mc;
> DifferentialEvolution deOptim(conf);
> modelG2pp->calibrate(swaptions, deOptim,
> EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>
> }
>
> ------------------------------------------------------------------------------
> Android apps run on BlackBerry 10
> Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
> Now with support for Jelly Bean, Bluetooth, Mapview and more.
> Get your Android app in front of a whole new audience. Start now.
> http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
> _______________________________________________
> QuantLib-users mailing list
> [hidden email]
> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>

------------------------------------------------------------------------------
Android apps run on BlackBerry 10
Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
Now with support for Jelly Bean, Bluetooth, Mapview and more.
Get your Android app in front of a whole new audience. Start now.
http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
_______________________________________________
QuantLib-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-users

Peter Caspers-4

Re: calibration G2++ with defferential evolution

oh, ok. got it. On second thought, instead of "fixing"
CompositeConstraint we should rather

1 add upperBound() and lowerBound() implementations to
CalibratedModel::PrivateConstraint (for this we need access to the
parameters' constraints though to read their bounds, so an additional
inspector might be needed)

2 maybe change the default implementation in Constraint to return an
array with size params.size() and value numeric_limits::max() / min()
? If not overwritten in derived classes, at least the correct
structure is returned indicating no constraint

still does not feel like the perfect solution ...

Peter

On 17 February 2014 14:01, Ralph Schreyer <[hidden email]> wrote:

> Hi Peter,
>
> but the DifferentialEvolution in its minimize(...) implementation of the
> Optimizer interface calls DifferentialEvolution::fillInitialPopulation which
> in turn calls the upperBound() and lowerBound() methods of the problem
> given, no?
>
>
> Best
> Ralph
>
>
>
>
>
> 2014-02-17 13:54 GMT+01:00 Peter Caspers <[hidden email]>:
>
>> Hi Ralph,
>>
>> indeed, this looks like a bug in CompositeConstraint.
>>
>> This should however not cause trouble if André implements his
>> constraints via the test() method directly (as he did ?). The
>> CompositeConstraint then just returns c1_.test(params) &&
>> c2_.test(params) without referring to any upperBound() / lowerBound()
>> implementation ?
>>
>> I am just guessing though.
>>
>> best
>> Peter
>>
>> On 17 February 2014 13:39, Ralph Schreyer <[hidden email]>
>> wrote:
>> > Hi Andre, Peter,
>> >
>> > I think the current implementation is buggy. You would do something like
>> >
>> > Array lower(5); lower[0] = 0.0; lower[1] = 0.0; lower[2] = 0.0; lower[3]
>> > =
>> > -1.0; lower[4] = 0.0; Array upper(5); upper[0] = 10.0; upper[1] = 50.0;
>> > upper[2] = 10.0; upper[3] = 1.0; upper[4] = 10.0;
>> > model->calibrate(options,
>> > deOptimizer, EndCriteria(400, 40, 1.0e-8, 1.0e-8, Null<Real>()),
>> > NonhomogeneousBoundaryConstraint(lower, upper));
>> >
>> > Then, in model.cpp, line 88 ff, the CompositeConstraint is built from
>> > the
>> > PrivateConstraint of the model and the NonhomogeneousBoundaryConstraint.
>> > When you look at the implementation of e.g. upperBound(...) in
>> > constraint.hpp, line 156 ff., you find
>> >
>> > for (Size iter = 0; iter < c1ub.size(); iter++) { rtrnArray.at(iter) =
>> > std::min(c1ub.at(iter), c2ub.at(iter)); }
>> >
>> > Because c1ub is the upperBound of the PrivateConstraint and c1ub.size()
>> > = 0,
>> > the DifferentialEvolution does not get its bounds and therefore e.g.
>> > DifferentialEvolution::fillInitialPopulation(...) fails.
>> >
>> > Maybe one can change the CompositeConstraint upperBound and lowerBound
>> > implementations? Or is there a better place? Or isn't there a bug at
>> > all??
>> >
>> >
>> > Best
>> > Ralph
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> > 2014-02-17 13:00 GMT+01:00 Peter Caspers <[hidden email]>:
>> >
>> >> Hi André,
>> >>
>> >> can you please post the complete code needed to reproduce the problem
>> >> ? With a short description of what the problem is (if not obvious) ?
>> >>
>> >> best regards
>> >> Peter
>> >>
>> >>
>> >> On 16 February 2014 17:28, André de Boer <[hidden email]> wrote:
>> >> > Hello,
>> >> >
>> >> > Two days ago I posted a question about differential evolution and the
>> >> > shortratemodel G2.
>> >> > Meanwhile I learned how to minimize the Rosenbrockfunction with
>> >> > differential
>> >> > evolution.
>> >> > But my goal is to calibrate the G2 shortratemodel.
>> >> >
>> >> > For the 5 parameters of the G2++ model I want to have
>> >> > boundedconstraints:
>> >> > For instance:
>> >> > 0 < a < 6
>> >> > 0 > >> > 0 < sigma < 0.25
>> >> > 0 < eta < 0.25
>> >> > -1 < rho < 1
>> >> >
>> >> > To accomplish this I made a class MyConstraint, see code below.
>> >> > ...
>> >> > MyConstraint mc;
>> >> > DifferentialEvolution deOptim(conf);
>> >> > modelG2pp->calibrate(swaptions, deOptim,
>> >> > EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>> >> > ...
>> >> > But still I isn't not working.
>> >> > Can someone give me a hint of parts of code?
>> >> >
>> >> > Kind regards,
>> >> > André
>> >> >
>> >> >
>> >> >
>> >> > #include <iostream>
>> >> > #include <vector>
>> >> > #include <ql/quantlib.hpp>
>> >> > using namespace QuantLib;
>> >> >
>> >> > class MyConstraint : public Constraint {
>> >> > class Impl : public Constraint::Impl {
>> >> > public:
>> >> >
>> >> > bool test(const Array& x) const {
>> >> > Real a1 = x[0], a2 = x[1], a3 = x[2], a4=x[3], a5=x[4];
>> >> >
>> >> > return (0.0 <= a1 && a1 <= 6.0) &&
>> >> > (0.0 <= a2 && a2 <= 6.0) &&
>> >> > (0.0 <= a3 && a3 <= 0.25) &&
>> >> > (0.0 <= a4 && a4 <= 0.25) &&
>> >> > (-1.0 <= a5 && a5 <= 1.0);
>> >> >
>> >> > }
>> >> >
>> >> > };
>> >> > public:
>> >> >
>> >> > MyConstraint()
>> >> > : Constraint(boost::shared_ptr<Constraint::Impl>(new
>> >> > MyConstraint::Impl)) {}
>> >> >
>> >> > };
>> >> >
>> >> > int main(){
>> >> >
>> >> > ...
>> >> > double a = 0.1;
>> >> >
>> >> > double b = 0.1;
>> >> > double sigma = 0.01;
>> >> > double eta = 0.01;
>> >> > double rho = -0.75;
>> >> > boost::shared_ptr<G2> modelG2pp(new G2(yieldCurve, a, b, sigma, eta,
>> >> > rho));
>> >> > std::cout << "G2 (analytic formulae) calibration" << std::endl;
>> >> >
>> >> > MyConstraint mc;
>> >> > DifferentialEvolution::Configuration conf =
>> >> > DifferentialEvolution::Configuration()
>> >> > .withStepsizeWeight(0.4)
>> >> > .withBounds()
>> >> > .withCrossoverProbability(0.35)
>> >> > .withPopulationMembers(100)
>> >> > .withStrategy(DifferentialEvolution::BestMemberWithJitter)
>> >> > .withCrossoverType(DifferentialEvolution::Normal)
>> >> > .withAdaptiveCrossover()
>> >> > .withSeed(0);
>> >> >
>> >> > MyConstraint mc;
>> >> > DifferentialEvolution deOptim(conf);
>> >> > modelG2pp->calibrate(swaptions, deOptim,
>> >> > EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>> >> >
>> >> > }
>> >> >
>> >> >
>> >> >
>> >> > ------------------------------------------------------------------------------
>> >> > Android apps run on BlackBerry 10
>> >> > Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
>> >> > Now with support for Jelly Bean, Bluetooth, Mapview and more.
>> >> > Get your Android app in front of a whole new audience. Start now.
>> >> >
>> >> >
>> >> > http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
>> >> > _______________________________________________
>> >> > QuantLib-users mailing list
>> >> > [hidden email]
>> >> > https://lists.sourceforge.net/lists/listinfo/quantlib-users
>> >> >
>> >>
>> >>
>> >>
>> >> ------------------------------------------------------------------------------
>> >> Android apps run on BlackBerry 10
>> >> Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
>> >> Now with support for Jelly Bean, Bluetooth, Mapview and more.
>> >> Get your Android app in front of a whole new audience. Start now.
>> >>
>> >>
>> >> http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
>> >> _______________________________________________
>> >> QuantLib-users mailing list
>> >> [hidden email]
>> >> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>> >
>> >
>
>

Ralph Schreyer-3

Re: calibration G2++ with defferential evolution

In reply to this post by Peter Caspers-4

Hi Peter, André,

"fixing" the CompositeConstraint in constraint.hpp like this (this is not a real fix, it is just to get the code running)

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

class CompositeConstraint : public Constraint {

private:

class Impl : public Constraint::Impl {

public:

Impl(const Constraint& c1,

const Constraint& c2)

: c1_(c1), c2_(c2) {}

bool test(const Array& params) const {

return c1_.test(params) && c2_.test(params);

}

Array upperBound(const Array& params) const {

Array c1ub = c1_.upperBound(params);

Array c2ub = c2_.upperBound(params);

Size ubSize = std::max(c1ub.size(), c2ub.size());

Array rtrnArray(ubSize, 0.0);

if (c1ub.size() == c2ub.size()) {

for (Size iter = 0; iter < c1ub.size(); iter++) {

rtrnArray.at(iter) = std::min(c1ub.at(iter), c2ub.at(iter));

}

else if (c1ub.size() > c2ub.size()) {

for (Size iter = 0; iter < c1ub.size(); iter++) {

rtrnArray.at(iter) = c1ub.at(iter);

}

else {

for (Size iter = 0; iter < c2ub.size(); iter++) {

rtrnArray.at(iter) = c2ub.at(iter);

}

return rtrnArray;

}

Array lowerBound(const Array& params) const {

Array c1lb = c1_.lowerBound(params);

Array c2lb = c2_.lowerBound(params);

Size lbSize = std::max(c1lb.size(), c2lb.size());

Array rtrnArray(lbSize, 0.0);

if (c1lb.size() == c2lb.size()) {

for (Size iter = 0; iter < c1lb.size(); iter++) {

rtrnArray.at(iter) = std::min(c1lb.at(iter), c2lb.at(iter));

}

else if (c1lb.size() > c2lb.size()) {

for (Size iter = 0; iter < c1lb.size(); iter++) {

rtrnArray.at(iter) = c1lb.at(iter);

}

else {

for (Size iter = 0; iter < c2lb.size(); iter++) {

rtrnArray.at(iter) = c2lb.at(iter);

}

return rtrnArray;

}

private:

Constraint c1_, c2_;

};

public:

CompositeConstraint(const Constraint& c1, const Constraint& c2)

: Constraint(boost::shared_ptr<Constraint::Impl>(

new CompositeConstraint::Impl(c1,c2))) {}

};

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

and exchanging the LevenbergMarquardt optimizer in method testDAXCalibration() for the DifferentialEvolution like this

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

void HestonModelTest::testDAXCalibration() {

BOOST_TEST_MESSAGE(

"Testing Heston model calibration using DAX volatility data...");

SavedSettings backup;

Date settlementDate(5, July, 2002);

Settings::instance().evaluationDate() = settlementDate;

CalibrationMarketData marketData = getDAXCalibrationMarketData();

const Handle<YieldTermStructure> riskFreeTS = marketData.riskFreeTS;

const Handle<YieldTermStructure> dividendTS = marketData.dividendYield;

const Handle<Quote> s0 = marketData.s0;

const std::vector<boost::shared_ptr<CalibrationHelper> > options

= marketData.options;

const Real v0=0.1;

const Real kappa=1.0;

const Real theta=0.1;

const Real sigma=0.5;

const Real rho=-0.5;

boost::shared_ptr<HestonProcess> process(new HestonProcess(

riskFreeTS, dividendTS, s0, v0, kappa, theta, sigma, rho));

boost::shared_ptr<HestonModel> model(new HestonModel(process));

boost::shared_ptr<PricingEngine> engine(

new AnalyticHestonEngine(model, 64));

for (Size i = 0; i < options.size(); ++i)

options[i]->setPricingEngine(engine);

DifferentialEvolution::Configuration conf =

DifferentialEvolution::Configuration()

.withStepsizeWeight(1.8)

.withBounds()

.withCrossoverProbability(0.9)

.withPopulationMembers(50)

.withStrategy(DifferentialEvolution::Rand1SelfadaptiveWithRotation)

.withCrossoverType(DifferentialEvolution::Normal)

.withAdaptiveCrossover()

.withSeed(3242);

DifferentialEvolution deOptim(conf);

Array lower(5);

lower[0] = 0.001;

lower[1] = 0.001;

lower[2] = 0.001;

lower[3] = -1.0;

lower[4] = 0.001;

Array upper(5);

upper[0] = 10.0;

upper[1] = 50.0;

upper[2] = 10.0;

upper[3] = 1.0;

upper[4] = 10.0;

model->calibrate(options, deOptim, EndCriteria(400, 40, 1.0e-8, 1.0e-8, Null<Real>()),

NonhomogeneousBoundaryConstraint(lower, upper));

Real sse = 0;

for (Size i = 0; i < 13*8; ++i) {

const Real diff = options[i]->calibrationError()*100.0;

sse += diff*diff;

}

Real expected = 177.2; //see article by A. Sepp.

if (std::fabs(sse - expected) > 1.0) {

BOOST_FAIL("Failed to reproduce calibration error"

<< "\n calculated: " << sse

<< "\n expected: " << expected);

}

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

I get the code runnning and the test succeeds.

Best

Ralph

2014-02-17 13:00 GMT+01:00 Peter Caspers <[hidden email]>:

Hi André,

can you please post the complete code needed to reproduce the problem
? With a short description of what the problem is (if not obvious) ?

best regards
Peter

On 16 February 2014 17:28, André de Boer <[hidden email]> wrote:
> Hello,
>
> Two days ago I posted a question about differential evolution and the
> shortratemodel G2.
> Meanwhile I learned how to minimize the Rosenbrockfunction with differential
> evolution.
> But my goal is to calibrate the G2 shortratemodel.
>
> For the 5 parameters of the G2++ model I want to have boundedconstraints:
> For instance:
> 0 < a < 6
> 0 0 < sigma < 0.25
> 0 < eta < 0.25
> -1 < rho < 1
>
> To accomplish this I made a class MyConstraint, see code below.
> ...
> MyConstraint mc;
> DifferentialEvolution deOptim(conf);
> modelG2pp->calibrate(swaptions, deOptim,
> EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
> ...
> But still I isn't not working.
> Can someone give me a hint of parts of code?
>
> Kind regards,
> André
>
>
>
> #include <iostream>
> #include <vector>
> #include <ql/quantlib.hpp>
> using namespace QuantLib;
>
> class MyConstraint : public Constraint {
> class Impl : public Constraint::Impl {
> public:
>
> bool test(const Array& x) const {
> Real a1 = x[0], a2 = x[1], a3 = x[2], a4=x[3], a5=x[4];
>
> return (0.0 <= a1 && a1 <= 6.0) &&
> (0.0 <= a2 && a2 <= 6.0) &&
> (0.0 <= a3 && a3 <= 0.25) &&
> (0.0 <= a4 && a4 <= 0.25) &&
> (-1.0 <= a5 && a5 <= 1.0);
>
> }
>
> };
> public:
>
> MyConstraint()
> : Constraint(boost::shared_ptr<Constraint::Impl>(new MyConstraint::Impl)) {}
>
> };
>
> int main(){
>
> ...
> double a = 0.1;
>
> double b = 0.1;
> double sigma = 0.01;
> double eta = 0.01;
> double rho = -0.75;
> boost::shared_ptr<G2> modelG2pp(new G2(yieldCurve, a, b, sigma, eta, rho));
> std::cout << "G2 (analytic formulae) calibration" << std::endl;
>
> MyConstraint mc;
> DifferentialEvolution::Configuration conf =
> DifferentialEvolution::Configuration()
> .withStepsizeWeight(0.4)
> .withBounds()
> .withCrossoverProbability(0.35)
> .withPopulationMembers(100)
> .withStrategy(DifferentialEvolution::BestMemberWithJitter)
> .withCrossoverType(DifferentialEvolution::Normal)
> .withAdaptiveCrossover()
> .withSeed(0);
>
> MyConstraint mc;
> DifferentialEvolution deOptim(conf);
> modelG2pp->calibrate(swaptions, deOptim,
> EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>
> }
>
> ------------------------------------------------------------------------------
> Android apps run on BlackBerry 10
> Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
> Now with support for Jelly Bean, Bluetooth, Mapview and more.
> Get your Android app in front of a whole new audience. Start now.
> http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
> _______________________________________________
> QuantLib-users mailing list
> [hidden email]
> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>

------------------------------------------------------------------------------
Android apps run on BlackBerry 10
Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
Now with support for Jelly Bean, Bluetooth, Mapview and more.
Get your Android app in front of a whole new audience. Start now.
http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
_______________________________________________
QuantLib-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-users

André de Boer

Re: calibration G2++ with defferential evolution

Hi Ralph, Peter,

Thanks for your comments and code, it looks promising and I will definitely use parts of it.

Meanwhile I also made some progress

Reimplementation of G2 with MyG2 with NoConstraint and defining NonhomogeneousBoundaryConstraint in main.cpp gets the code also running.

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

// file myG2.h

#ifndef MYG2_H

#define MYG2_H

#include <ql/models/shortrate/twofactormodels/g2.hpp>

class MyG2 : public QuantLib::G2

{

public:

MyG2(const QuantLib::Handle<QuantLib::YieldTermStructure>& termStructure,

QuantLib::Real a = 0.1,

QuantLib::Real sigma = 0.01,

QuantLib::Real b = 0.1,

QuantLib::Real eta = 0.01,

QuantLib::Real rho = -0.75);

};

#endif // MYG2_H

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

//file myG2.cpp

#include "myg2.h"

#include <ql/math/optimization/constraint.hpp>

MyG2::MyG2(const QuantLib::Handle<QuantLib::YieldTermStructure> &termStructure,

QuantLib::Real a,

QuantLib::Real sigma,

QuantLib::Real b,

QuantLib::Real eta,

QuantLib::Real rho) : G2(termStructure,a,sigma,b,eta,rho)

{

constraint_=boost::shared_ptr<QuantLib::NoConstraint>(new QuantLib::NoConstraint);

}

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

// main.cpp

....

boost::shared_ptr<MyG2> modelG2pp(new MyG2(yieldCurve, a, sigma, b, eta, rho));
DifferentialEvolution::Configuration conf =
DifferentialEvolution::Configuration()
.withStepsizeWeight(0.4)
.withBounds()
.withCrossoverProbability(0.35)
.withPopulationMembers(50)
.withStrategy(DifferentialEvolution::BestMemberWithJitter)
.withCrossoverType(DifferentialEvolution::Normal)
.withAdaptiveCrossover()
.withSeed(0);

DifferentialEvolution deOptim(conf);

QuantLib::Array lowerBound(5);
QuantLib::Array upperBound(5);
lowerBound[0]=0.0;
lowerBound[1]=0.0;
lowerBound[2]=0.0;
lowerBound[3]=0.0;
lowerBound[4]=-1.0;

upperBound[0]=20.0;
upperBound[1]=0.1;
upperBound[2]=20.0;
upperBound[3]=0.1;
upperBound[4]=1.0;

QuantLib::NonhomogeneousBoundaryConstraint constraint(lowerBound,upperBound);

modelG2pp->calibrate(swaptions, deOptim, EndCriteria(400,40,1e-10,1e-8,Null<Real>()),constraint);

...

2014-02-17 17:01 GMT+01:00 Ralph Schreyer <[hidden email]>:

Hi Peter, André,

"fixing" the CompositeConstraint in constraint.hpp like this (this is not a real fix, it is just to get the code running)

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

class CompositeConstraint : public Constraint {
private:
class Impl : public Constraint::Impl {
public:
Impl(const Constraint& c1,

const Constraint& c2)
: c1_(c1), c2_(c2) {}
bool test(const Array& params) const {
return c1_.test(params) && c2_.test(params);

}
Array upperBound(const Array& params) const {
Array c1ub = c1_.upperBound(params);
Array c2ub = c2_.upperBound(params);

Size ubSize = std::max(c1ub.size(), c2ub.size());
Array rtrnArray(ubSize, 0.0);
if (c1ub.size() == c2ub.size()) {

for (Size iter = 0; iter < c1ub.size(); iter++) {
rtrnArray.at(iter) = std::min(c1ub.at(iter), c2ub.at(iter));

}
}
else if (c1ub.size() > c2ub.size()) {

for (Size iter = 0; iter < c1ub.size(); iter++) {
rtrnArray.at(iter) = c1ub.at(iter);

}
}
else {
for (Size iter = 0; iter < c2ub.size(); iter++) {

rtrnArray.at(iter) = c2ub.at(iter);
}
}

return rtrnArray;
}
Array lowerBound(const Array& params) const {
Array c1lb = c1_.lowerBound(params);
Array c2lb = c2_.lowerBound(params);

Size lbSize = std::max(c1lb.size(), c2lb.size());
Array rtrnArray(lbSize, 0.0);
if (c1lb.size() == c2lb.size()) {

for (Size iter = 0; iter < c1lb.size(); iter++) {
rtrnArray.at(iter) = std::min(c1lb.at(iter), c2lb.at(iter));

}
}
else if (c1lb.size() > c2lb.size()) {

for (Size iter = 0; iter < c1lb.size(); iter++) {
rtrnArray.at(iter) = c1lb.at(iter);

}
}
else {
for (Size iter = 0; iter < c2lb.size(); iter++) {

rtrnArray.at(iter) = c2lb.at(iter);
}
}

return rtrnArray;
}
private:
Constraint c1_, c2_;
};
public:
CompositeConstraint(const Constraint& c1, const Constraint& c2)

: Constraint(boost::shared_ptr<Constraint::Impl>(
new CompositeConstraint::Impl(c1,c2))) {}
};

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

and exchanging the LevenbergMarquardt optimizer in method testDAXCalibration() for the DifferentialEvolution like this

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

void HestonModelTest::testDAXCalibration() {

BOOST_TEST_MESSAGE(
"Testing Heston model calibration using DAX volatility data...");

SavedSettings backup;

Date settlementDate(5, July, 2002);
Settings::instance().evaluationDate() = settlementDate;

CalibrationMarketData marketData = getDAXCalibrationMarketData();

const Handle<YieldTermStructure> riskFreeTS = marketData.riskFreeTS;
const Handle<YieldTermStructure> dividendTS = marketData.dividendYield;
const Handle<Quote> s0 = marketData.s0;

const std::vector<boost::shared_ptr<CalibrationHelper> > options
= marketData.options;

const Real v0=0.1;

const Real kappa=1.0;
const Real theta=0.1;
const Real sigma=0.5;
const Real rho=-0.5;

boost::shared_ptr<HestonProcess> process(new HestonProcess(

riskFreeTS, dividendTS, s0, v0, kappa, theta, sigma, rho));

boost::shared_ptr<HestonModel> model(new HestonModel(process));

boost::shared_ptr<PricingEngine> engine(

new AnalyticHestonEngine(model, 64));

for (Size i = 0; i < options.size(); ++i)
options[i]->setPricingEngine(engine);

DifferentialEvolution::Configuration conf =
DifferentialEvolution::Configuration()
.withStepsizeWeight(1.8)
.withBounds()
.withCrossoverProbability(0.9)

.withPopulationMembers(50)
.withStrategy(DifferentialEvolution::Rand1SelfadaptiveWithRotation)
.withCrossoverType(DifferentialEvolution::Normal)
.withAdaptiveCrossover()

.withSeed(3242);
DifferentialEvolution deOptim(conf);

Array lower(5);
lower[0] = 0.001;

lower[1] = 0.001;
lower[2] = 0.001;
lower[3] = -1.0;

lower[4] = 0.001;
Array upper(5);
upper[0] = 10.0;

upper[1] = 50.0;
upper[2] = 10.0;
upper[3] = 1.0;

upper[4] = 10.0;
model->calibrate(options, deOptim, EndCriteria(400, 40, 1.0e-8, 1.0e-8, Null<Real>()),

NonhomogeneousBoundaryConstraint(lower, upper));

Real sse = 0;
for (Size i = 0; i < 13*8; ++i) {
const Real diff = options[i]->calibrationError()*100.0;

sse += diff*diff;
}
Real expected = 177.2; //see article by A. Sepp.
if (std::fabs(sse - expected) > 1.0) {
BOOST_FAIL("Failed to reproduce calibration error"

<< "\n calculated: " << sse
<< "\n expected: " << expected);
}
}

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<:

Hi André,

can you please post the complete code needed to reproduce the problem
? With a short description of what the problem is (if not obvious) ?

best regards
Peter

On 16 February 2014 17:28, André de Boer <[hidden email]> wrote:
> Hello,
>
> Two days ago I posted a question about differential evolution and the
> shortratemodel G2.
> Meanwhile I learned how to minimize the Rosenbrockfunction with differential
> evolution.
> But my goal is to calibrate the G2 shortratemodel.
>
> For the 5 parameters of the G2++ model I want to have boundedconstraints:
> For instance:
> 0 < a < 6
> 0 0 < sigma < 0.25
> 0 < eta < 0.25
> -1 < rho < 1
>
> To accomplish this I made a class MyConstraint, see code below.
> ...
> MyConstraint mc;
> DifferentialEvolution deOptim(conf);
> modelG2pp->calibrate(swaptions, deOptim,
> EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
> ...
> But still I isn't not working.
> Can someone give me a hint of parts of code?
>
> Kind regards,
> André
>
>
>
> #include <iostream>
> #include <vector>
> #include <ql/quantlib.hpp>
> using namespace QuantLib;
>
> class MyConstraint : public Constraint {
> class Impl : public Constraint::Impl {
> public:
>
> bool test(const Array& x) const {
> Real a1 = x[0], a2 = x[1], a3 = x[2], a4=x[3], a5=x[4];
>
> return (0.0 <= a1 && a1 <= 6.0) &&
> (0.0 <= a2 && a2 <= 6.0) &&
> (0.0 <= a3 && a3 <= 0.25) &&
> (0.0 <= a4 && a4 <= 0.25) &&
> (-1.0 <= a5 && a5 <= 1.0);
>
> }
>
> };
> public:
>
> MyConstraint()
> : Constraint(boost::shared_ptr<Constraint::Impl>(new MyConstraint::Impl)) {}
>
> };
>
> int main(){
>
> ...
> double a = 0.1;
>
> double b = 0.1;
> double sigma = 0.01;
> double eta = 0.01;
> double rho = -0.75;
> boost::shared_ptr<G2> modelG2pp(new G2(yieldCurve, a, b, sigma, eta, rho));
> std::cout << "G2 (analytic formulae) calibration" << std::endl;
>
> MyConstraint mc;
> DifferentialEvolution::Configuration conf =
> DifferentialEvolution::Configuration()
> .withStepsizeWeight(0.4)
> .withBounds()
> .withCrossoverProbability(0.35)
> .withPopulationMembers(100)
> .withStrategy(DifferentialEvolution::BestMemberWithJitter)
> .withCrossoverType(DifferentialEvolution::Normal)
> .withAdaptiveCrossover()
> .withSeed(0);
>
> MyConstraint mc;
> DifferentialEvolution deOptim(conf);
> modelG2pp->calibrate(swaptions, deOptim,
> EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>
> }
>
> ------------------------------------------------------------------------------
> Android apps run on BlackBerry 10
> Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
> Now with support for Jelly Bean, Bluetooth, Mapview and more.
> Get your Android app in front of a whole new audience. Start now.
> http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
> _______________________________________________
> QuantLib-users mailing list
> [hidden email]
> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>

------------------------------------------------------------------------------
Android apps run on BlackBerry 10
Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
Now with support for Jelly Bean, Bluetooth, Mapview and more.
Get your Android app in front of a whole new audience. Start now.
http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
_______________________________________________
QuantLib-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-users

------------------------------------------------------------------------------
Managing the Performance of Cloud-Based Applications
Take advantage of what the Cloud has to offer - Avoid Common Pitfalls.
Read the Whitepaper.
http://pubads.g.doubleclick.net/gampad/clk?id=121054471&iu=/4140/ostg.clktrk
_______________________________________________
QuantLib-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-users

André de Boer

Re: calibration G2++ with defferential evolution

Hi Peter, Ralph,

It still puzzles me that the CompositeConstraint has to be "fixed" as you states in your answer:

"fixing" the CompositeConstraint in constraint.hpp like this...".

Why is it not possible to use the class "NonhomogeneousBoundaryConstraint" in constraint.hpp directly?

It looks like a bug to me, but I am not an expert.

2014-02-17 17:51 GMT+01:00 André de Boer <[hidden email]>:

Hi Ralph, Peter,

Thanks for your comments and code, it looks promising and I will definitely use parts of it.

Meanwhile I also made some progress

Reimplementation of G2 with MyG2 with NoConstraint and defining NonhomogeneousBoundaryConstraint in main.cpp gets the code also running.

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

// file myG2.h

#ifndef MYG2_H

#define MYG2_H

#include <ql/models/shortrate/twofactormodels/g2.hpp>

class MyG2 : public QuantLib::G2

{

public:

MyG2(const QuantLib::Handle<QuantLib::YieldTermStructure>& termStructure,

QuantLib::Real a = 0.1,

QuantLib::Real sigma = 0.01,

QuantLib::Real b = 0.1,

QuantLib::Real eta = 0.01,

QuantLib::Real rho = -0.75);

};

#endif // MYG2_H

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

//file myG2.cpp

#include "myg2.h"

#include <ql/math/optimization/constraint.hpp>

MyG2::MyG2(const QuantLib::Handle<QuantLib::YieldTermStructure> &termStructure,

QuantLib::Real a,

QuantLib::Real sigma,

QuantLib::Real b,

QuantLib::Real eta,

QuantLib::Real rho) : G2(termStructure,a,sigma,b,eta,rho)

{

constraint_=boost::shared_ptr<QuantLib::NoConstraint>(new QuantLib::NoConstraint);

}

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

// main.cpp

....
boost::shared_ptr<MyG2> modelG2pp(new MyG2(yieldCurve, a, sigma, b, eta, rho));

DifferentialEvolution::Configuration conf =
DifferentialEvolution::Configuration()
.withStepsizeWeight(0.4)
.withBounds()
.withCrossoverProbability(0.35)
.withPopulationMembers(50)

.withStrategy(DifferentialEvolution::BestMemberWithJitter)
.withCrossoverType(DifferentialEvolution::Normal)
.withAdaptiveCrossover()
.withSeed(0);

DifferentialEvolution deOptim(conf);

QuantLib::Array lowerBound(5);
QuantLib::Array upperBound(5);
lowerBound[0]=0.0;
lowerBound[1]=0.0;
lowerBound[2]=0.0;
lowerBound[3]=0.0;
lowerBound[4]=-1.0;

upperBound[0]=20.0;
upperBound[1]=0.1;
upperBound[2]=20.0;
upperBound[3]=0.1;
upperBound[4]=1.0;

QuantLib::NonhomogeneousBoundaryConstraint constraint(lowerBound,upperBound);

modelG2pp->calibrate(swaptions, deOptim, EndCriteria(400,40,1e-10,1e-8,Null<Real>()),constraint);

...

2014-02-17 17:01 GMT+01:00 Ralph Schreyer <[hidden email]>:

Hi Peter, André,

"fixing" the CompositeConstraint in constraint.hpp like this (this is not a real fix, it is just to get the code running)

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

class CompositeConstraint : public Constraint {

private:

class Impl : public Constraint::Impl {

public:

Impl(const Constraint& c1,

const Constraint& c2)

: c1_(c1), c2_(c2) {}

bool test(const Array& params) const {

return c1_.test(params) && c2_.test(params);

}

Array upperBound(const Array& params) const {

Array c1ub = c1_.upperBound(params);

Array c2ub = c2_.upperBound(params);

Size ubSize = std::max(c1ub.size(), c2ub.size());

Array rtrnArray(ubSize, 0.0);

if (c1ub.size() == c2ub.size()) {

for (Size iter = 0; iter < c1ub.size(); iter++) {

rtrnArray.at(iter) = std::min(c1ub.at(iter), c2ub.at(iter));

}

}

else if (c1ub.size() > c2ub.size()) {

for (Size iter = 0; iter < c1ub.size(); iter++) {

rtrnArray.at(iter) = c1ub.at(iter);

}

}

else {

for (Size iter = 0; iter < c2ub.size(); iter++) {

rtrnArray.at(iter) = c2ub.at(iter);

}

}

return rtrnArray;

}

Array lowerBound(const Array& params) const {

Array c1lb = c1_.lowerBound(params);

Array c2lb = c2_.lowerBound(params);

Size lbSize = std::max(c1lb.size(), c2lb.size());

Array rtrnArray(lbSize, 0.0);

if (c1lb.size() == c2lb.size()) {

for (Size iter = 0; iter < c1lb.size(); iter++) {

rtrnArray.at(iter) = std::min(c1lb.at(iter), c2lb.at(iter));

}

}

else if (c1lb.size() > c2lb.size()) {

for (Size iter = 0; iter < c1lb.size(); iter++) {

rtrnArray.at(iter) = c1lb.at(iter);

}

}

else {

for (Size iter = 0; iter < c2lb.size(); iter++) {

rtrnArray.at(iter) = c2lb.at(iter);

}

}

return rtrnArray;

}

private:

Constraint c1_, c2_;

};

public:

CompositeConstraint(const Constraint& c1, const Constraint& c2)

: Constraint(boost::shared_ptr<Constraint::Impl>(

new CompositeConstraint::Impl(c1,c2))) {}

};

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

and exchanging the LevenbergMarquardt optimizer in method testDAXCalibration() for the DifferentialEvolution like this

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

void HestonModelTest::testDAXCalibration() {

BOOST_TEST_MESSAGE(

"Testing Heston model calibration using DAX volatility data...");

SavedSettings backup;

Date settlementDate(5, July, 2002);

Settings::instance().evaluationDate() = settlementDate;

CalibrationMarketData marketData = getDAXCalibrationMarketData();

const Handle<YieldTermStructure> riskFreeTS = marketData.riskFreeTS;

const Handle<YieldTermStructure> dividendTS = marketData.dividendYield;

const Handle<Quote> s0 = marketData.s0;

const std::vector<boost::shared_ptr<CalibrationHelper> > options

= marketData.options;

const Real v0=0.1;

const Real kappa=1.0;

const Real theta=0.1;

const Real sigma=0.5;

const Real rho=-0.5;

boost::shared_ptr<HestonProcess> process(new HestonProcess(

riskFreeTS, dividendTS, s0, v0, kappa, theta, sigma, rho));

boost::shared_ptr<HestonModel> model(new HestonModel(process));

boost::shared_ptr<PricingEngine> engine(

new AnalyticHestonEngine(model, 64));

for (Size i = 0; i < options.size(); ++i)

options[i]->setPricingEngine(engine);

DifferentialEvolution::Configuration conf =

DifferentialEvolution::Configuration()

.withStepsizeWeight(1.8)

.withBounds()

.withCrossoverProbability(0.9)

.withPopulationMembers(50)

.withStrategy(DifferentialEvolution::Rand1SelfadaptiveWithRotation)

.withCrossoverType(DifferentialEvolution::Normal)

.withAdaptiveCrossover()

.withSeed(3242);

DifferentialEvolution deOptim(conf);

Array lower(5);

lower[0] = 0.001;

lower[1] = 0.001;

lower[2] = 0.001;

lower[3] = -1.0;

lower[4] = 0.001;

Array upper(5);

upper[0] = 10.0;

upper[1] = 50.0;

upper[2] = 10.0;

upper[3] = 1.0;

upper[4] = 10.0;

model->calibrate(options, deOptim, EndCriteria(400, 40, 1.0e-8, 1.0e-8, Null<Real>()),

NonhomogeneousBoundaryConstraint(lower, upper));

Real sse = 0;

for (Size i = 0; i < 13*8; ++i) {

const Real diff = options[i]->calibrationError()*100.0;

sse += diff*diff;

}

Real expected = 177.2; //see article by A. Sepp.

if (std::fabs(sse - expected) > 1.0) {

BOOST_FAIL("Failed to reproduce calibration error"

<< "\n calculated: " << sse

<< "\n expected: " << expected);

}

}

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<:

Hi André,

can you please post the complete code needed to reproduce the problem
? With a short description of what the problem is (if not obvious) ?

best regards
Peter

On 16 February 2014 17:28, André de Boer <[hidden email]> wrote:
> Hello,
>
> Two days ago I posted a question about differential evolution and the
> shortratemodel G2.
> Meanwhile I learned how to minimize the Rosenbrockfunction with differential
> evolution.
> But my goal is to calibrate the G2 shortratemodel.
>
> For the 5 parameters of the G2++ model I want to have boundedconstraints:
> For instance:
> 0 < a < 6
> 0 0 < sigma < 0.25
> 0 < eta < 0.25
> -1 < rho < 1
>
> To accomplish this I made a class MyConstraint, see code below.
> ...
> MyConstraint mc;
> DifferentialEvolution deOptim(conf);
> modelG2pp->calibrate(swaptions, deOptim,
> EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
> ...
> But still I isn't not working.
> Can someone give me a hint of parts of code?
>
> Kind regards,
> André
>
>
>
> #include <iostream>
> #include <vector>
> #include <ql/quantlib.hpp>
> using namespace QuantLib;
>
> class MyConstraint : public Constraint {
> class Impl : public Constraint::Impl {
> public:
>
> bool test(const Array& x) const {
> Real a1 = x[0], a2 = x[1], a3 = x[2], a4=x[3], a5=x[4];
>
> return (0.0 <= a1 && a1 <= 6.0) &&
> (0.0 <= a2 && a2 <= 6.0) &&
> (0.0 <= a3 && a3 <= 0.25) &&
> (0.0 <= a4 && a4 <= 0.25) &&
> (-1.0 <= a5 && a5 <= 1.0);
>
> }
>
> };
> public:
>
> MyConstraint()

> : Constraint(boost::shared_ptr<Constraint::Impl>(new MyConstraint::Impl)) {}
>
> };
>
> int main(){
>
> ...
> double a = 0.1;
>
> double b = 0.1;
> double sigma = 0.01;
> double eta = 0.01;
> double rho = -0.75;
> boost::shared_ptr<G2> modelG2pp(new G2(yieldCurve, a, b, sigma, eta, rho));
> std::cout << "G2 (analytic formulae) calibration" << std::endl;
>
> MyConstraint mc;
> DifferentialEvolution::Configuration conf =
> DifferentialEvolution::Configuration()
> .withStepsizeWeight(0.4)
> .withBounds()
> .withCrossoverProbability(0.35)
> .withPopulationMembers(100)
> .withStrategy(DifferentialEvolution::BestMemberWithJitter)
> .withCrossoverType(DifferentialEvolution::Normal)
> .withAdaptiveCrossover()
> .withSeed(0);
>
> MyConstraint mc;
> DifferentialEvolution deOptim(conf);
> modelG2pp->calibrate(swaptions, deOptim,
> EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>
> }
>
> ------------------------------------------------------------------------------
> Android apps run on BlackBerry 10
> Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
> Now with support for Jelly Bean, Bluetooth, Mapview and more.
> Get your Android app in front of a whole new audience. Start now.
> http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
> _______________________________________________
> QuantLib-users mailing list
> [hidden email]
> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>

------------------------------------------------------------------------------
Android apps run on BlackBerry 10
Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
Now with support for Jelly Bean, Bluetooth, Mapview and more.
Get your Android app in front of a whole new audience. Start now.
http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
_______________________________________________
QuantLib-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-users

Peter Caspers-4

Re: calibration G2++ with defferential evolution

I tried to write a fix, which you can find here

https://github.com/lballabio/quantlib/pull/80

It compiles and the test suite does not break, but I did not test
anything else. Maybe you can try to apply this patch and see if it
solves the
original problem ?

best regards
Peter

On 18 February 2014 09:15, André de Boer <[hidden email]> wrote:

> Hi Peter, Ralph,
>
> It still puzzles me that the CompositeConstraint has to be "fixed" as you
> states in your answer:
> "fixing" the CompositeConstraint in constraint.hpp like this...".
> Why is it not possible to use the class "NonhomogeneousBoundaryConstraint"
> in constraint.hpp directly?
> It looks like a bug to me, but I am not an expert.
>
> 2014-02-17 17:51 GMT+01:00 André de Boer <[hidden email]>:
>>
>> Hi Ralph, Peter,
>>
>> Thanks for your comments and code, it looks promising and I will
>> definitely use parts of it.
>> Meanwhile I also made some progress
>>
>> Reimplementation of G2 with MyG2 with NoConstraint and defining
>> NonhomogeneousBoundaryConstraint in main.cpp gets the code also running.
>>
>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>> // file myG2.h
>> #ifndef MYG2_H
>> #define MYG2_H
>>
>> #include <ql/models/shortrate/twofactormodels/g2.hpp>
>>
>> class MyG2 : public QuantLib::G2
>> {
>> public:
>> MyG2(const QuantLib::Handle<QuantLib::YieldTermStructure>&
>> termStructure,
>> QuantLib::Real a = 0.1,
>> QuantLib::Real sigma = 0.01,
>> QuantLib::Real b = 0.1,
>> QuantLib::Real eta = 0.01,
>> QuantLib::Real rho = -0.75);
>> };
>>
>> #endif // MYG2_H
>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>
>> //file myG2.cpp
>> #include "myg2.h"
>> #include <ql/math/optimization/constraint.hpp>
>>
>> MyG2::MyG2(const QuantLib::Handle<QuantLib::YieldTermStructure>
>> &termStructure,
>> QuantLib::Real a,
>> QuantLib::Real sigma,
>> QuantLib::Real b,
>> QuantLib::Real eta,
>> QuantLib::Real rho) : G2(termStructure,a,sigma,b,eta,rho)
>> {
>> constraint_=boost::shared_ptr<QuantLib::NoConstraint>(new
>> QuantLib::NoConstraint);
>> }
>>
>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>
>> // main.cpp
>>
>> ....
>> boost::shared_ptr<MyG2> modelG2pp(new MyG2(yieldCurve, a, sigma, b, eta,
>> rho));
>>
>> DifferentialEvolution::Configuration conf =
>> DifferentialEvolution::Configuration()
>> .withStepsizeWeight(0.4)
>> .withBounds()
>> .withCrossoverProbability(0.35)
>> .withPopulationMembers(50)
>>
>> .withStrategy(DifferentialEvolution::BestMemberWithJitter)
>> .withCrossoverType(DifferentialEvolution::Normal)
>> .withAdaptiveCrossover()
>> .withSeed(0);
>>
>> DifferentialEvolution deOptim(conf);
>>
>> QuantLib::Array lowerBound(5);
>> QuantLib::Array upperBound(5);
>> lowerBound[0]=0.0;
>> lowerBound[1]=0.0;
>> lowerBound[2]=0.0;
>> lowerBound[3]=0.0;
>> lowerBound[4]=-1.0;
>>
>> upperBound[0]=20.0;
>> upperBound[1]=0.1;
>> upperBound[2]=20.0;
>> upperBound[3]=0.1;
>> upperBound[4]=1.0;
>>
>> QuantLib::NonhomogeneousBoundaryConstraint
>> constraint(lowerBound,upperBound);
>> modelG2pp->calibrate(swaptions, deOptim,
>> EndCriteria(400,40,1e-10,1e-8,Null<Real>()),constraint);
>> ...
>>
>>
>> 2014-02-17 17:01 GMT+01:00 Ralph Schreyer <[hidden email]>:
>>
>>> Hi Peter, André,
>>>
>>> "fixing" the CompositeConstraint in constraint.hpp like this (this is not
>>> a real fix, it is just to get the code running)
>>>
>>> >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
>>> class CompositeConstraint : public Constraint {
>>> private:
>>> class Impl : public Constraint::Impl {
>>> public:
>>> Impl(const Constraint& c1,
>>> const Constraint& c2)
>>> : c1_(c1), c2_(c2) {}
>>> bool test(const Array& params) const {
>>> return c1_.test(params) && c2_.test(params);
>>> }
>>> Array upperBound(const Array& params) const {
>>> Array c1ub = c1_.upperBound(params);
>>> Array c2ub = c2_.upperBound(params);
>>> Size ubSize = std::max(c1ub.size(), c2ub.size());
>>> Array rtrnArray(ubSize, 0.0);
>>> if (c1ub.size() == c2ub.size()) {
>>> for (Size iter = 0; iter < c1ub.size(); iter++) {
>>> rtrnArray.at(iter) = std::min(c1ub.at(iter), c2ub.at(iter));
>>> }
>>> }
>>> else if (c1ub.size() > c2ub.size()) {
>>> for (Size iter = 0; iter < c1ub.size(); iter++) {
>>> rtrnArray.at(iter) = c1ub.at(iter);
>>> }
>>> }
>>> else {
>>> for (Size iter = 0; iter < c2ub.size(); iter++) {
>>> rtrnArray.at(iter) = c2ub.at(iter);
>>> }
>>> }
>>> return rtrnArray;
>>> }
>>> Array lowerBound(const Array& params) const {
>>> Array c1lb = c1_.lowerBound(params);
>>> Array c2lb = c2_.lowerBound(params);
>>> Size lbSize = std::max(c1lb.size(), c2lb.size());
>>> Array rtrnArray(lbSize, 0.0);
>>> if (c1lb.size() == c2lb.size()) {
>>> for (Size iter = 0; iter < c1lb.size(); iter++) {
>>> rtrnArray.at(iter) = std::min(c1lb.at(iter), c2lb.at(iter));
>>> }
>>> }
>>> else if (c1lb.size() > c2lb.size()) {
>>> for (Size iter = 0; iter < c1lb.size(); iter++) {
>>> rtrnArray.at(iter) = c1lb.at(iter);
>>> }
>>> }
>>> else {
>>> for (Size iter = 0; iter < c2lb.size(); iter++) {
>>> rtrnArray.at(iter) = c2lb.at(iter);
>>> }
>>> }
>>> return rtrnArray;
>>> }
>>> private:
>>> Constraint c1_, c2_;
>>> };
>>> public:
>>> CompositeConstraint(const Constraint& c1, const Constraint& c2)
>>> : Constraint(boost::shared_ptr<Constraint::Impl>(
>>> new
>>> CompositeConstraint::Impl(c1,c2))) {}
>>> };
>>>
>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>>
>>> and exchanging the LevenbergMarquardt optimizer in method
>>> testDAXCalibration() for the DifferentialEvolution like this
>>>
>>> >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
>>>
>>> void HestonModelTest::testDAXCalibration() {
>>>
>>> BOOST_TEST_MESSAGE(
>>> "Testing Heston model calibration using DAX volatility
>>> data...");
>>>
>>> SavedSettings backup;
>>>
>>> Date settlementDate(5, July, 2002);
>>> Settings::instance().evaluationDate() = settlementDate;
>>>
>>> CalibrationMarketData marketData = getDAXCalibrationMarketData();
>>>
>>> const Handle<YieldTermStructure> riskFreeTS = marketData.riskFreeTS;
>>> const Handle<YieldTermStructure> dividendTS =
>>> marketData.dividendYield;
>>> const Handle<Quote> s0 = marketData.s0;
>>>
>>> const std::vector<boost::shared_ptr<CalibrationHelper> > options
>>> = marketData.options;
>>>
>>> const Real v0=0.1;
>>> const Real kappa=1.0;
>>> const Real theta=0.1;
>>> const Real sigma=0.5;
>>> const Real rho=-0.5;
>>>
>>> boost::shared_ptr<HestonProcess> process(new HestonProcess(
>>> riskFreeTS, dividendTS, s0, v0, kappa, theta, sigma, rho));
>>>
>>> boost::shared_ptr<HestonModel> model(new HestonModel(process));
>>>
>>> boost::shared_ptr<PricingEngine> engine(
>>> new AnalyticHestonEngine(model,
>>> 64));
>>>
>>> for (Size i = 0; i < options.size(); ++i)
>>> options[i]->setPricingEngine(engine);
>>>
>>> DifferentialEvolution::Configuration conf =
>>> DifferentialEvolution::Configuration()
>>> .withStepsizeWeight(1.8)
>>> .withBounds()
>>> .withCrossoverProbability(0.9)
>>> .withPopulationMembers(50)
>>>
>>> .withStrategy(DifferentialEvolution::Rand1SelfadaptiveWithRotation)
>>> .withCrossoverType(DifferentialEvolution::Normal)
>>> .withAdaptiveCrossover()
>>> .withSeed(3242);
>>> DifferentialEvolution deOptim(conf);
>>>
>>> Array lower(5);
>>> lower[0] = 0.001;
>>> lower[1] = 0.001;
>>> lower[2] = 0.001;
>>> lower[3] = -1.0;
>>> lower[4] = 0.001;
>>> Array upper(5);
>>> upper[0] = 10.0;
>>> upper[1] = 50.0;
>>> upper[2] = 10.0;
>>> upper[3] = 1.0;
>>> upper[4] = 10.0;
>>> model->calibrate(options, deOptim, EndCriteria(400, 40, 1.0e-8, 1.0e-8,
>>> Null<Real>()),
>>> NonhomogeneousBoundaryConstraint(lower, upper));
>>>
>>> Real sse = 0;
>>> for (Size i = 0; i < 13*8; ++i) {
>>> const Real diff = options[i]->calibrationError()*100.0;
>>> sse += diff*diff;
>>> }
>>> Real expected = 177.2; //see article by A. Sepp.
>>> if (std::fabs(sse - expected) > 1.0) {
>>> BOOST_FAIL("Failed to reproduce calibration error"
>>> << "\n calculated: " << sse
>>> << "\n expected: " << expected);
>>> }
>>> }
>>>
>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>>
>>> I get the code runnning and the test succeeds.
>>>
>>> Best
>>> Ralph
>>>
>>>
>>> 2014-02-17 13:00 GMT+01:00 Peter Caspers <[hidden email]>:
>>>>
>>>> Hi André,
>>>>
>>>>
>>>> can you please post the complete code needed to reproduce the problem
>>>> ? With a short description of what the problem is (if not obvious) ?
>>>>
>>>> best regards
>>>> Peter
>>>>
>>>>
>>>> On 16 February 2014 17:28, André de Boer <[hidden email]> wrote:
>>>> > Hello,
>>>> >
>>>> > Two days ago I posted a question about differential evolution and the
>>>> > shortratemodel G2.
>>>> > Meanwhile I learned how to minimize the Rosenbrockfunction with
>>>> > differential
>>>> > evolution.
>>>> > But my goal is to calibrate the G2 shortratemodel.
>>>> >
>>>> > For the 5 parameters of the G2++ model I want to have
>>>> > boundedconstraints:
>>>> > For instance:
>>>> > 0 < a < 6
>>>> > 0 >>> > 0 < sigma < 0.25
>>>> > 0 < eta < 0.25
>>>> > -1 < rho < 1
>>>> >
>>>> > To accomplish this I made a class MyConstraint, see code below.
>>>> > ...
>>>> > MyConstraint mc;
>>>> > DifferentialEvolution deOptim(conf);
>>>> > modelG2pp->calibrate(swaptions, deOptim,
>>>> > EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>>>> > ...
>>>> > But still I isn't not working.
>>>> > Can someone give me a hint of parts of code?
>>>> >
>>>> > Kind regards,
>>>> > André
>>>> >
>>>> >
>>>> >
>>>> > #include <iostream>
>>>> > #include <vector>
>>>> > #include <ql/quantlib.hpp>
>>>> > using namespace QuantLib;
>>>> >
>>>> > class MyConstraint : public Constraint {
>>>> > class Impl : public Constraint::Impl {
>>>> > public:
>>>> >
>>>> > bool test(const Array& x) const {
>>>> > Real a1 = x[0], a2 = x[1], a3 = x[2], a4=x[3], a5=x[4];
>>>> >
>>>> > return (0.0 <= a1 && a1 <= 6.0) &&
>>>> > (0.0 <= a2 && a2 <= 6.0) &&
>>>> > (0.0 <= a3 && a3 <= 0.25) &&
>>>> > (0.0 <= a4 && a4 <= 0.25) &&
>>>> > (-1.0 <= a5 && a5 <= 1.0);
>>>> >
>>>> > }
>>>> >
>>>> > };
>>>> > public:
>>>> >
>>>> > MyConstraint()
>>>> > : Constraint(boost::shared_ptr<Constraint::Impl>(new
>>>> > MyConstraint::Impl)) {}
>>>> >
>>>> > };
>>>> >
>>>> > int main(){
>>>> >
>>>> > ...
>>>> > double a = 0.1;
>>>> >
>>>> > double b = 0.1;
>>>> > double sigma = 0.01;
>>>> > double eta = 0.01;
>>>> > double rho = -0.75;
>>>> > boost::shared_ptr<G2> modelG2pp(new G2(yieldCurve, a, b, sigma, eta,
>>>> > rho));
>>>> > std::cout << "G2 (analytic formulae) calibration" << std::endl;
>>>> >
>>>> > MyConstraint mc;
>>>> > DifferentialEvolution::Configuration conf =
>>>> > DifferentialEvolution::Configuration()
>>>> > .withStepsizeWeight(0.4)
>>>> > .withBounds()
>>>> > .withCrossoverProbability(0.35)
>>>> > .withPopulationMembers(100)
>>>> > .withStrategy(DifferentialEvolution::BestMemberWithJitter)
>>>> > .withCrossoverType(DifferentialEvolution::Normal)
>>>> > .withAdaptiveCrossover()
>>>> > .withSeed(0);
>>>> >
>>>> > MyConstraint mc;
>>>> > DifferentialEvolution deOptim(conf);
>>>> > modelG2pp->calibrate(swaptions, deOptim,
>>>> > EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>>>> >
>>>> > }
>>>> >
>>>> >
>>>> > ------------------------------------------------------------------------------
>>>> > Android apps run on BlackBerry 10
>>>> > Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
>>>> > Now with support for Jelly Bean, Bluetooth, Mapview and more.
>>>> > Get your Android app in front of a whole new audience. Start now.
>>>> >
>>>> > http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
>>>> > _______________________________________________
>>>> > QuantLib-users mailing list
>>>> > [hidden email]
>>>> > https://lists.sourceforge.net/lists/listinfo/quantlib-users
>>>> >
>>>>
>>>>
>>>> ------------------------------------------------------------------------------
>>>> Android apps run on BlackBerry 10
>>>> Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
>>>> Now with support for Jelly Bean, Bluetooth, Mapview and more.
>>>> Get your Android app in front of a whole new audience. Start now.
>>>>
>>>> http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
>>>> _______________________________________________
>>>> QuantLib-users mailing list
>>>> [hidden email]
>>>> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>>>
>>>
>>
>

André de Boer

Re: calibration G2++ with defferential evolution

Hoi Peter,

Thanks for your contribution to solve this issue.

I will try this patch and let you know.

Kind regards,

André

2014-02-22 7:59 GMT+01:00 Peter Caspers <[hidden email]>:

I tried to write a fix, which you can find here

https://github.com/lballabio/quantlib/pull/80

It compiles and the test suite does not break, but I did not test
anything else. Maybe you can try to apply this patch and see if it
solves the
original problem ?

best regards
Peter

On 18 February 2014 09:15, André de Boer <[hidden email]> wrote:
> Hi Peter, Ralph,
>
> It still puzzles me that the CompositeConstraint has to be "fixed" as you
> states in your answer:
> "fixing" the CompositeConstraint in constraint.hpp like this...".
> Why is it not possible to use the class "NonhomogeneousBoundaryConstraint"
> in constraint.hpp directly?
> It looks like a bug to me, but I am not an expert.
>
> 2014-02-17 17:51 GMT+01:00 André de Boer <[hidden email]>:
>>
>> Hi Ralph, Peter,
>>
>> Thanks for your comments and code, it looks promising and I will
>> definitely use parts of it.
>> Meanwhile I also made some progress
>>
>> Reimplementation of G2 with MyG2 with NoConstraint and defining
>> NonhomogeneousBoundaryConstraint in main.cpp gets the code also running.
>>
>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>> // file myG2.h
>> #ifndef MYG2_H
>> #define MYG2_H
>>
>> #include <ql/models/shortrate/twofactormodels/g2.hpp>
>>
>> class MyG2 : public QuantLib::G2
>> {
>> public:
>> MyG2(const QuantLib::Handle<QuantLib::YieldTermStructure>&
>> termStructure,
>> QuantLib::Real a = 0.1,
>> QuantLib::Real sigma = 0.01,
>> QuantLib::Real b = 0.1,
>> QuantLib::Real eta = 0.01,
>> QuantLib::Real rho = -0.75);
>> };
>>
>> #endif // MYG2_H
>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>
>> //file myG2.cpp
>> #include "myg2.h"
>> #include <ql/math/optimization/constraint.hpp>
>>
>> MyG2::MyG2(const QuantLib::Handle<QuantLib::YieldTermStructure>
>> &termStructure,
>> QuantLib::Real a,
>> QuantLib::Real sigma,
>> QuantLib::Real b,
>> QuantLib::Real eta,
>> QuantLib::Real rho) : G2(termStructure,a,sigma,b,eta,rho)
>> {
>> constraint_=boost::shared_ptr<QuantLib::NoConstraint>(new
>> QuantLib::NoConstraint);
>> }
>>
>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>
>> // main.cpp
>>
>> ....
>> boost::shared_ptr<MyG2> modelG2pp(new MyG2(yieldCurve, a, sigma, b, eta,
>> rho));
>>
>> DifferentialEvolution::Configuration conf =
>> DifferentialEvolution::Configuration()
>> .withStepsizeWeight(0.4)
>> .withBounds()
>> .withCrossoverProbability(0.35)
>> .withPopulationMembers(50)
>>
>> .withStrategy(DifferentialEvolution::BestMemberWithJitter)
>> .withCrossoverType(DifferentialEvolution::Normal)
>> .withAdaptiveCrossover()
>> .withSeed(0);
>>
>> DifferentialEvolution deOptim(conf);
>>
>> QuantLib::Array lowerBound(5);
>> QuantLib::Array upperBound(5);
>> lowerBound[0]=0.0;
>> lowerBound[1]=0.0;
>> lowerBound[2]=0.0;
>> lowerBound[3]=0.0;
>> lowerBound[4]=-1.0;
>>
>> upperBound[0]=20.0;
>> upperBound[1]=0.1;
>> upperBound[2]=20.0;
>> upperBound[3]=0.1;
>> upperBound[4]=1.0;
>>
>> QuantLib::NonhomogeneousBoundaryConstraint
>> constraint(lowerBound,upperBound);
>> modelG2pp->calibrate(swaptions, deOptim,
>> EndCriteria(400,40,1e-10,1e-8,Null<Real>()),constraint);
>> ...
>>
>>
>> 2014-02-17 17:01 GMT+01:00 Ralph Schreyer <[hidden email]>:
>>
>>> Hi Peter, André,
>>>
>>> "fixing" the CompositeConstraint in constraint.hpp like this (this is not
>>> a real fix, it is just to get the code running)
>>>
>>> >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
>>> class CompositeConstraint : public Constraint {
>>> private:
>>> class Impl : public Constraint::Impl {
>>> public:
>>> Impl(const Constraint& c1,
>>> const Constraint& c2)
>>> : c1_(c1), c2_(c2) {}
>>> bool test(const Array& params) const {
>>> return c1_.test(params) && c2_.test(params);
>>> }
>>> Array upperBound(const Array& params) const {
>>> Array c1ub = c1_.upperBound(params);
>>> Array c2ub = c2_.upperBound(params);
>>> Size ubSize = std::max(c1ub.size(), c2ub.size());
>>> Array rtrnArray(ubSize, 0.0);
>>> if (c1ub.size() == c2ub.size()) {
>>> for (Size iter = 0; iter < c1ub.size(); iter++) {
>>> rtrnArray.at(iter) = std::min(c1ub.at(iter), c2ub.at(iter));
>>> }
>>> }
>>> else if (c1ub.size() > c2ub.size()) {
>>> for (Size iter = 0; iter < c1ub.size(); iter++) {
>>> rtrnArray.at(iter) = c1ub.at(iter);
>>> }
>>> }
>>> else {
>>> for (Size iter = 0; iter < c2ub.size(); iter++) {
>>> rtrnArray.at(iter) = c2ub.at(iter);
>>> }
>>> }
>>> return rtrnArray;
>>> }
>>> Array lowerBound(const Array& params) const {
>>> Array c1lb = c1_.lowerBound(params);
>>> Array c2lb = c2_.lowerBound(params);
>>> Size lbSize = std::max(c1lb.size(), c2lb.size());
>>> Array rtrnArray(lbSize, 0.0);
>>> if (c1lb.size() == c2lb.size()) {
>>> for (Size iter = 0; iter < c1lb.size(); iter++) {
>>> rtrnArray.at(iter) = std::min(c1lb.at(iter), c2lb.at(iter));
>>> }
>>> }
>>> else if (c1lb.size() > c2lb.size()) {
>>> for (Size iter = 0; iter < c1lb.size(); iter++) {
>>> rtrnArray.at(iter) = c1lb.at(iter);
>>> }
>>> }
>>> else {
>>> for (Size iter = 0; iter < c2lb.size(); iter++) {
>>> rtrnArray.at(iter) = c2lb.at(iter);
>>> }
>>> }
>>> return rtrnArray;
>>> }
>>> private:
>>> Constraint c1_, c2_;
>>> };
>>> public:
>>> CompositeConstraint(const Constraint& c1, const Constraint& c2)
>>> : Constraint(boost::shared_ptr<Constraint::Impl>(
>>> new
>>> CompositeConstraint::Impl(c1,c2))) {}
>>> };
>>>
>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>>
>>> and exchanging the LevenbergMarquardt optimizer in method
>>> testDAXCalibration() for the DifferentialEvolution like this
>>>
>>> >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
>>>
>>> void HestonModelTest::testDAXCalibration() {
>>>
>>> BOOST_TEST_MESSAGE(
>>> "Testing Heston model calibration using DAX volatility
>>> data...");
>>>
>>> SavedSettings backup;
>>>
>>> Date settlementDate(5, July, 2002);
>>> Settings::instance().evaluationDate() = settlementDate;
>>>
>>> CalibrationMarketData marketData = getDAXCalibrationMarketData();
>>>
>>> const Handle<YieldTermStructure> riskFreeTS = marketData.riskFreeTS;
>>> const Handle<YieldTermStructure> dividendTS =
>>> marketData.dividendYield;
>>> const Handle<Quote> s0 = marketData.s0;
>>>
>>> const std::vector<boost::shared_ptr<CalibrationHelper> > options
>>> = marketData.options;
>>>
>>> const Real v0=0.1;
>>> const Real kappa=1.0;
>>> const Real theta=0.1;
>>> const Real sigma=0.5;
>>> const Real rho=-0.5;
>>>
>>> boost::shared_ptr<HestonProcess> process(new HestonProcess(
>>> riskFreeTS, dividendTS, s0, v0, kappa, theta, sigma, rho));
>>>
>>> boost::shared_ptr<HestonModel> model(new HestonModel(process));
>>>
>>> boost::shared_ptr<PricingEngine> engine(
>>> new AnalyticHestonEngine(model,
>>> 64));
>>>
>>> for (Size i = 0; i < options.size(); ++i)
>>> options[i]->setPricingEngine(engine);
>>>
>>> DifferentialEvolution::Configuration conf =
>>> DifferentialEvolution::Configuration()
>>> .withStepsizeWeight(1.8)
>>> .withBounds()
>>> .withCrossoverProbability(0.9)
>>> .withPopulationMembers(50)
>>>
>>> .withStrategy(DifferentialEvolution::Rand1SelfadaptiveWithRotation)
>>> .withCrossoverType(DifferentialEvolution::Normal)
>>> .withAdaptiveCrossover()
>>> .withSeed(3242);
>>> DifferentialEvolution deOptim(conf);
>>>
>>> Array lower(5);
>>> lower[0] = 0.001;
>>> lower[1] = 0.001;
>>> lower[2] = 0.001;
>>> lower[3] = -1.0;
>>> lower[4] = 0.001;
>>> Array upper(5);
>>> upper[0] = 10.0;
>>> upper[1] = 50.0;
>>> upper[2] = 10.0;
>>> upper[3] = 1.0;
>>> upper[4] = 10.0;
>>> model->calibrate(options, deOptim, EndCriteria(400, 40, 1.0e-8, 1.0e-8,
>>> Null<Real>()),
>>> NonhomogeneousBoundaryConstraint(lower, upper));
>>>
>>> Real sse = 0;
>>> for (Size i = 0; i < 13*8; ++i) {
>>> const Real diff = options[i]->calibrationError()*100.0;
>>> sse += diff*diff;
>>> }
>>> Real expected = 177.2; //see article by A. Sepp.
>>> if (std::fabs(sse - expected) > 1.0) {
>>> BOOST_FAIL("Failed to reproduce calibration error"
>>> << "\n calculated: " << sse
>>> << "\n expected: " << expected);
>>> }
>>> }
>>>
>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>>
>>> I get the code runnning and the test succeeds.
>>>
>>> Best
>>> Ralph
>>>
>>>
>>> 2014-02-17 13:00 GMT+01:00 Peter Caspers <[hidden email]>:
>>>>
>>>> Hi André,
>>>>
>>>>
>>>> can you please post the complete code needed to reproduce the problem
>>>> ? With a short description of what the problem is (if not obvious) ?
>>>>
>>>> best regards
>>>> Peter
>>>>
>>>>
>>>> On 16 February 2014 17:28, André de Boer <[hidden email]> wrote:
>>>> > Hello,
>>>> >
>>>> > Two days ago I posted a question about differential evolution and the
>>>> > shortratemodel G2.
>>>> > Meanwhile I learned how to minimize the Rosenbrockfunction with
>>>> > differential
>>>> > evolution.
>>>> > But my goal is to calibrate the G2 shortratemodel.
>>>> >
>>>> > For the 5 parameters of the G2++ model I want to have
>>>> > boundedconstraints:
>>>> > For instance:
>>>> > 0 < a < 6
>>>> > 0 >>> > 0 < sigma < 0.25
>>>> > 0 < eta < 0.25
>>>> > -1 < rho < 1
>>>> >
>>>> > To accomplish this I made a class MyConstraint, see code below.
>>>> > ...
>>>> > MyConstraint mc;
>>>> > DifferentialEvolution deOptim(conf);
>>>> > modelG2pp->calibrate(swaptions, deOptim,
>>>> > EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>>>> > ...
>>>> > But still I isn't not working.
>>>> > Can someone give me a hint of parts of code?
>>>> >
>>>> > Kind regards,
>>>> > André
>>>> >
>>>> >
>>>> >
>>>> > #include <iostream>
>>>> > #include <vector>
>>>> > #include <ql/quantlib.hpp>
>>>> > using namespace QuantLib;
>>>> >
>>>> > class MyConstraint : public Constraint {
>>>> > class Impl : public Constraint::Impl {
>>>> > public:
>>>> >
>>>> > bool test(const Array& x) const {
>>>> > Real a1 = x[0], a2 = x[1], a3 = x[2], a4=x[3], a5=x[4];
>>>> >
>>>> > return (0.0 <= a1 && a1 <= 6.0) &&
>>>> > (0.0 <= a2 && a2 <= 6.0) &&
>>>> > (0.0 <= a3 && a3 <= 0.25) &&
>>>> > (0.0 <= a4 && a4 <= 0.25) &&
>>>> > (-1.0 <= a5 && a5 <= 1.0);
>>>> >
>>>> > }
>>>> >
>>>> > };
>>>> > public:
>>>> >
>>>> > MyConstraint()
>>>> > : Constraint(boost::shared_ptr<Constraint::Impl>(new
>>>> > MyConstraint::Impl)) {}
>>>> >
>>>> > };
>>>> >
>>>> > int main(){
>>>> >
>>>> > ...
>>>> > double a = 0.1;
>>>> >
>>>> > double b = 0.1;
>>>> > double sigma = 0.01;
>>>> > double eta = 0.01;
>>>> > double rho = -0.75;
>>>> > boost::shared_ptr<G2> modelG2pp(new G2(yieldCurve, a, b, sigma, eta,
>>>> > rho));
>>>> > std::cout << "G2 (analytic formulae) calibration" << std::endl;
>>>> >
>>>> > MyConstraint mc;
>>>> > DifferentialEvolution::Configuration conf =
>>>> > DifferentialEvolution::Configuration()
>>>> > .withStepsizeWeight(0.4)
>>>> > .withBounds()
>>>> > .withCrossoverProbability(0.35)
>>>> > .withPopulationMembers(100)
>>>> > .withStrategy(DifferentialEvolution::BestMemberWithJitter)
>>>> > .withCrossoverType(DifferentialEvolution::Normal)
>>>> > .withAdaptiveCrossover()
>>>> > .withSeed(0);
>>>> >
>>>> > MyConstraint mc;
>>>> > DifferentialEvolution deOptim(conf);
>>>> > modelG2pp->calibrate(swaptions, deOptim,
>>>> > EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>>>> >
>>>> > }
>>>> >
>>>> >
>>>> > ------------------------------------------------------------------------------
>>>> > Android apps run on BlackBerry 10
>>>> > Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
>>>> > Now with support for Jelly Bean, Bluetooth, Mapview and more.
>>>> > Get your Android app in front of a whole new audience. Start now.
>>>> >
>>>> > http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
>>>> > _______________________________________________
>>>> > QuantLib-users mailing list
>>>> > [hidden email]
>>>> > https://lists.sourceforge.net/lists/listinfo/quantlib-users
>>>> >
>>>>
>>>>
>>>> ------------------------------------------------------------------------------
>>>> Android apps run on BlackBerry 10
>>>> Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
>>>> Now with support for Jelly Bean, Bluetooth, Mapview and more.
>>>> Get your Android app in front of a whole new audience. Start now.
>>>>
>>>> http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
>>>> _______________________________________________
>>>> QuantLib-users mailing list
>>>> [hidden email]
>>>> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>>>
>>>
>>
>

André de Boer

Re: calibration G2++ with defferential evolution

Hi Peter,

I tried the fix but I get the message:
Error: parameter array too big!

This is my code:

#include <iostream>
#include <vector>
#include <ql/quantlib.hpp>
using namespace QuantLib;
int main()
{
Date today(13, Sep, 2013);
Natural settlementDays = 0;
Calendar calendar = UnitedStates(UnitedStates::Settlement);
DayCounter dayCounter = Actual360();
DayCounter dayCounterFixed = Thirty360();
Date settlementDate = calendar.advance(today, settlementDays, Days);
// SWAP Curve
std::vector<Date> dates; std::vector<Rate> rates;
dates.push_back(settlementDate); rates.push_back(0.445/100.0);
dates.push_back(calendar.advance(settlementDate,1,Years));
rates.push_back(0.445/100.0);
dates.push_back(calendar.advance(settlementDate,2,Years));
rates.push_back(0.6426/100.0);
dates.push_back(calendar.advance(settlementDate,3,Years));
rates.push_back(0.883/100.0);
dates.push_back(calendar.advance(settlementDate,4,Years));
rates.push_back(1.148/100.0);
dates.push_back(calendar.advance(settlementDate,5,Years));
rates.push_back(1.3976/100.0);
dates.push_back(calendar.advance(settlementDate,6,Years));
rates.push_back(1.6065/100.0);
dates.push_back(calendar.advance(settlementDate,7,Years));
rates.push_back(1.7906/100.0);
dates.push_back(calendar.advance(settlementDate,8,Years));
rates.push_back(1.9501/100.0);
dates.push_back(calendar.advance(settlementDate,9,Years));
rates.push_back(2.091/100.0);
dates.push_back(calendar.advance(settlementDate,10,Years));
rates.push_back(2.22/100.0);
dates.push_back(calendar.advance(settlementDate,12,Years));
rates.push_back(2.425/100.0);
dates.push_back(calendar.advance(settlementDate,15,Years));
rates.push_back(2.6227/100.0);
dates.push_back(calendar.advance(settlementDate,20,Years));
rates.push_back(2.7328/100.0);
dates.push_back(calendar.advance(settlementDate,25,Years));
rates.push_back(2.7428/100.0);
dates.push_back(calendar.advance(settlementDate,30,Years));
rates.push_back(2.73/100.0);
dates.push_back(calendar.advance(settlementDate,40,Years));
rates.push_back(2.7525/100.0);
dates.push_back(calendar.advance(settlementDate,50,Years));
rates.push_back(2.7945/100.0);
dates.push_back(calendar.advance(settlementDate,60,Years));
rates.push_back(2.7945/100.0);
dates.push_back(calendar.advance(settlementDate,61,Years));
rates.push_back(2.7945/100.0);
Handle<YieldTermStructure>
yieldCurve(boost::shared_ptr<YieldTermStructure>(new ZeroCurve(dates,
rates, dayCounter) ));
boost::shared_ptr<IborIndex> index(new USDLibor(Period(3, Months), yieldCurve));
Integer swapLengths[] = {1,2,3,4,5,6,7,8,9,10,15,20,25,30};
Integer swaptionLengths[] = {1,2,3,4,5,6,7,10,15,20,25,30};
Volatility swaptionVols[] = {
83.28, 65.22, 57.38, 51.62, 46.60, 41.99, 39.29, 36.28, 34.10, 32.67,
28.11, 26.00, 25.67, 25.30,
66.81, 52.88, 46.91, 42.32, 39.09, 36.35, 34.21, 32.34, 30.82, 29.85,
26.83, 25.46, 25.35, 25.11,
51.74, 43.57, 39.70, 36.58, 34.09, 32.31, 30.88, 29.53, 28.52, 27.87,
25.62, 24.67, 24.71, 24.50,
43.72, 37.47, 34.40, 32.28, 30.46, 29.25, 28.32, 27.44, 26.79, 26.25,
24.50, 24.17, 24.23, 23.90,
36.08, 32.42, 30.64, 29.11, 28.07, 26.95, 26.38, 25.91, 25.53, 25.26,
23.64, 23.68, 23.73, 23.36,
31.96, 29.30, 27.97, 26.77, 25.77, 25.28, 24.91, 24.66, 24.50, 24.40,
23.33, 23.01, 23.01, 22.57,
29.19, 26.94, 25.86, 25.02, 24.34, 24.06, 23.79, 23.60, 23.55, 23.52,
22.89, 22.51, 22.30, 21.92,
23.50, 22.61, 22.39, 22.25, 22.18, 22.31, 22.11, 22.08, 22.14, 22.24,
21.11, 21.20, 20.61, 19.90,
21.70, 22.10, 22.31, 22.51, 22.80, 22.87, 22.18, 22.31, 22.23, 22.23,
21.17, 19.40, 18.30, 17.40,
22.51, 23.07, 23.25, 23.42, 23.57, 23.45, 22.67, 22.45, 22.09, 21.84,
19.59, 17.35, 16.24, 15.38,
23.30, 23.47, 23.39, 23.25, 23.04, 22.54, 22.02, 21.48, 20.95, 20.39,
17.23, 15.13, 14.14, 13.43,
20.85, 20.61, 20.22, 19.88, 19.54, 19.06, 18.61, 18.20, 17.81, 17.41,
14.94, 13.20, 12.53, 11.84};

try {
// construct SwaptionHelpers
std::vector<boost::shared_ptr<CalibrationHelper> > swaptions;
std::list<Time> times;
for (Size i=0; i<12; i++) { // Swaption
for (Size j=0; j<14; j++) { // Swap
boost::shared_ptr<Quote> vol(new SimpleQuote(swaptionVols[i*14+j]/100.0));
swaptions.push_back(boost::shared_ptr<CalibrationHelper> (new
SwaptionHelper(Period(swaptionLengths[i], Years),
Period(swapLengths[j], Years), Handle<Quote>(vol),
index, // floating index
Period(6, Months), // fixed leg tenor
dayCounterFixed, // fixed leg dayCounter
dayCounter, // floating leg dayCounter
yieldCurve)));
swaptions.back()->addTimesTo(times);
}
}
double a = 0.1;
double sigma = 0.01;
double b = 0.1;
double eta = 0.01;
double rho = -0.75;
boost::shared_ptr<G2> modelG2pp(new G2(yieldCurve));
std::cout << "G2 (analytic formulae) calibration" << std::endl;
for (Size i=0; i<swaptions.size(); i++)
swaptions[i]->setPricingEngine(boost::shared_ptr<PricingEngine>(new
G2SwaptionEngine(modelG2pp, 6.0, 16)));
std::cout << "\nCalibration with DE" << std::endl;
DifferentialEvolution::Configuration conf =
DifferentialEvolution::Configuration()
.withStepsizeWeight(0.5)
.withBounds()
.withCrossoverProbability(0.4)
.withPopulationMembers(100) //greater than 209 will erase an error
.withStrategy(DifferentialEvolution::BestMemberWithJitter)
.withCrossoverType(DifferentialEvolution::Normal)
.withAdaptiveCrossover()
.withSeed(314);
DifferentialEvolution deOptim(conf);
QuantLib::Array lowerBound(5);
QuantLib::Array upperBound(5);
lowerBound[0]=0.0;
lowerBound[1]=0.0;
lowerBound[2]=0.0;
lowerBound[3]=0.0;
lowerBound[4]=-1.0;
upperBound[0]=5.0;
upperBound[1]=0.25;
upperBound[2]=5.0;
upperBound[3]=0.25;
upperBound[4]=1.0;
QuantLib::NonhomogeneousBoundaryConstraint constraints(lowerBound,upperBound);
modelG2pp->calibrate(swaptions, deOptim,
EndCriteria(50,20,1.0e-5,1.0e-5,Null<Real>()), constraints);
std::cout << "a = " << modelG2pp->params()[0]
<< "\nsigma = " << modelG2pp->params()[1]
<< "\nb = " << modelG2pp->params()[2]
<< "\neta = " << modelG2pp->params()[3]
<< "\nrho = " << modelG2pp->params()[4]
<< std::endl;
std::cout << "\nCriterium: " << modelG2pp->endCriteria() << std::endl;
Real sse = 0;
for (Size i = 0; i < 14*12; ++i) {
const Real diff = swaptions[i]->calibrationError()*100.0;
sse += diff*diff;
}
std::cout << "\ncalculated calibration error: " << sse << std::endl;
}
catch(std::exception &e)
{
std::cout << "Error: " << e.what() << std::endl;
}
return 0;
}

2014-02-23 22:42 GMT+01:00 André de Boer <[hidden email]>:

> Hoi Peter,
>
> Thanks for your contribution to solve this issue.
> I will try this patch and let you know.
>
> Kind regards,
> André
>
>
> 2014-02-22 7:59 GMT+01:00 Peter Caspers <[hidden email]>:
>
>> I tried to write a fix, which you can find here
>>
>> https://github.com/lballabio/quantlib/pull/80
>>
>> It compiles and the test suite does not break, but I did not test
>> anything else. Maybe you can try to apply this patch and see if it
>> solves the
>> original problem ?
>>
>> best regards
>> Peter
>>
>>
>> On 18 February 2014 09:15, André de Boer <[hidden email]> wrote:
>> > Hi Peter, Ralph,
>> >
>> > It still puzzles me that the CompositeConstraint has to be "fixed" as
>> > you
>> > states in your answer:
>> > "fixing" the CompositeConstraint in constraint.hpp like this...".
>> > Why is it not possible to use the class
>> > "NonhomogeneousBoundaryConstraint"
>> > in constraint.hpp directly?
>> > It looks like a bug to me, but I am not an expert.
>> >
>> > 2014-02-17 17:51 GMT+01:00 André de Boer <[hidden email]>:
>> >>
>> >> Hi Ralph, Peter,
>> >>
>> >> Thanks for your comments and code, it looks promising and I will
>> >> definitely use parts of it.
>> >> Meanwhile I also made some progress
>> >>
>> >> Reimplementation of G2 with MyG2 with NoConstraint and defining
>> >> NonhomogeneousBoundaryConstraint in main.cpp gets the code also
>> >> running.
>> >>
>> >> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>> >> // file myG2.h
>> >> #ifndef MYG2_H
>> >> #define MYG2_H
>> >>
>> >> #include <ql/models/shortrate/twofactormodels/g2.hpp>
>> >>
>> >> class MyG2 : public QuantLib::G2
>> >> {
>> >> public:
>> >> MyG2(const QuantLib::Handle<QuantLib::YieldTermStructure>&
>> >> termStructure,
>> >> QuantLib::Real a = 0.1,
>> >> QuantLib::Real sigma = 0.01,
>> >> QuantLib::Real b = 0.1,
>> >> QuantLib::Real eta = 0.01,
>> >> QuantLib::Real rho = -0.75);
>> >> };
>> >>
>> >> #endif // MYG2_H
>> >> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>> >>
>> >> //file myG2.cpp
>> >> #include "myg2.h"
>> >> #include <ql/math/optimization/constraint.hpp>
>> >>
>> >> MyG2::MyG2(const QuantLib::Handle<QuantLib::YieldTermStructure>
>> >> &termStructure,
>> >> QuantLib::Real a,
>> >> QuantLib::Real sigma,
>> >> QuantLib::Real b,
>> >> QuantLib::Real eta,
>> >> QuantLib::Real rho) : G2(termStructure,a,sigma,b,eta,rho)
>> >> {
>> >> constraint_=boost::shared_ptr<QuantLib::NoConstraint>(new
>> >> QuantLib::NoConstraint);
>> >> }
>> >>
>> >> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>> >>
>> >> // main.cpp
>> >>
>> >> ....
>> >> boost::shared_ptr<MyG2> modelG2pp(new MyG2(yieldCurve, a, sigma, b,
>> >> eta,
>> >> rho));
>> >>
>> >> DifferentialEvolution::Configuration conf =
>> >> DifferentialEvolution::Configuration()
>> >> .withStepsizeWeight(0.4)
>> >> .withBounds()
>> >> .withCrossoverProbability(0.35)
>> >> .withPopulationMembers(50)
>> >>
>> >> .withStrategy(DifferentialEvolution::BestMemberWithJitter)
>> >> .withCrossoverType(DifferentialEvolution::Normal)
>> >> .withAdaptiveCrossover()
>> >> .withSeed(0);
>> >>
>> >> DifferentialEvolution deOptim(conf);
>> >>
>> >> QuantLib::Array lowerBound(5);
>> >> QuantLib::Array upperBound(5);
>> >> lowerBound[0]=0.0;
>> >> lowerBound[1]=0.0;
>> >> lowerBound[2]=0.0;
>> >> lowerBound[3]=0.0;
>> >> lowerBound[4]=-1.0;
>> >>
>> >> upperBound[0]=20.0;
>> >> upperBound[1]=0.1;
>> >> upperBound[2]=20.0;
>> >> upperBound[3]=0.1;
>> >> upperBound[4]=1.0;
>> >>
>> >> QuantLib::NonhomogeneousBoundaryConstraint
>> >> constraint(lowerBound,upperBound);
>> >> modelG2pp->calibrate(swaptions, deOptim,
>> >> EndCriteria(400,40,1e-10,1e-8,Null<Real>()),constraint);
>> >> ...
>> >>
>> >>
>> >> 2014-02-17 17:01 GMT+01:00 Ralph Schreyer <[hidden email]>:
>> >>
>> >>> Hi Peter, André,
>> >>>
>> >>> "fixing" the CompositeConstraint in constraint.hpp like this (this is
>> >>> not
>> >>> a real fix, it is just to get the code running)
>> >>>
>> >>> >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
>> >>> class CompositeConstraint : public Constraint {
>> >>> private:
>> >>> class Impl : public Constraint::Impl {
>> >>> public:
>> >>> Impl(const Constraint& c1,
>> >>> const Constraint& c2)
>> >>> : c1_(c1), c2_(c2) {}
>> >>> bool test(const Array& params) const {
>> >>> return c1_.test(params) && c2_.test(params);
>> >>> }
>> >>> Array upperBound(const Array& params) const {
>> >>> Array c1ub = c1_.upperBound(params);
>> >>> Array c2ub = c2_.upperBound(params);
>> >>> Size ubSize = std::max(c1ub.size(), c2ub.size());
>> >>> Array rtrnArray(ubSize, 0.0);
>> >>> if (c1ub.size() == c2ub.size()) {
>> >>> for (Size iter = 0; iter < c1ub.size(); iter++) {
>> >>> rtrnArray.at(iter) = std::min(c1ub.at(iter), c2ub.at(iter));
>> >>> }
>> >>> }
>> >>> else if (c1ub.size() > c2ub.size()) {
>> >>> for (Size iter = 0; iter < c1ub.size(); iter++) {
>> >>> rtrnArray.at(iter) = c1ub.at(iter);
>> >>> }
>> >>> }
>> >>> else {
>> >>> for (Size iter = 0; iter < c2ub.size(); iter++) {
>> >>> rtrnArray.at(iter) = c2ub.at(iter);
>> >>> }
>> >>> }
>> >>> return rtrnArray;
>> >>> }
>> >>> Array lowerBound(const Array& params) const {
>> >>> Array c1lb = c1_.lowerBound(params);
>> >>> Array c2lb = c2_.lowerBound(params);
>> >>> Size lbSize = std::max(c1lb.size(), c2lb.size());
>> >>> Array rtrnArray(lbSize, 0.0);
>> >>> if (c1lb.size() == c2lb.size()) {
>> >>> for (Size iter = 0; iter < c1lb.size(); iter++) {
>> >>> rtrnArray.at(iter) = std::min(c1lb.at(iter), c2lb.at(iter));
>> >>> }
>> >>> }
>> >>> else if (c1lb.size() > c2lb.size()) {
>> >>> for (Size iter = 0; iter < c1lb.size(); iter++) {
>> >>> rtrnArray.at(iter) = c1lb.at(iter);
>> >>> }
>> >>> }
>> >>> else {
>> >>> for (Size iter = 0; iter < c2lb.size(); iter++) {
>> >>> rtrnArray.at(iter) = c2lb.at(iter);
>> >>> }
>> >>> }
>> >>> return rtrnArray;
>> >>> }
>> >>> private:
>> >>> Constraint c1_, c2_;
>> >>> };
>> >>> public:
>> >>> CompositeConstraint(const Constraint& c1, const Constraint&
>> >>> c2)
>> >>> : Constraint(boost::shared_ptr<Constraint::Impl>(
>> >>> new
>> >>> CompositeConstraint::Impl(c1,c2))) {}
>> >>> };
>> >>>
>> >>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>> >>>
>> >>> and exchanging the LevenbergMarquardt optimizer in method
>> >>> testDAXCalibration() for the DifferentialEvolution like this
>> >>>
>> >>> >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
>> >>>
>> >>> void HestonModelTest::testDAXCalibration() {
>> >>>
>> >>> BOOST_TEST_MESSAGE(
>> >>> "Testing Heston model calibration using DAX volatility
>> >>> data...");
>> >>>
>> >>> SavedSettings backup;
>> >>>
>> >>> Date settlementDate(5, July, 2002);
>> >>> Settings::instance().evaluationDate() = settlementDate;
>> >>>
>> >>> CalibrationMarketData marketData = getDAXCalibrationMarketData();
>> >>>
>> >>> const Handle<YieldTermStructure> riskFreeTS =
>> >>> marketData.riskFreeTS;
>> >>> const Handle<YieldTermStructure> dividendTS =
>> >>> marketData.dividendYield;
>> >>> const Handle<Quote> s0 = marketData.s0;
>> >>>
>> >>> const std::vector<boost::shared_ptr<CalibrationHelper> > options
>> >>> =
>> >>> marketData.options;
>> >>>
>> >>> const Real v0=0.1;
>> >>> const Real kappa=1.0;
>> >>> const Real theta=0.1;
>> >>> const Real sigma=0.5;
>> >>> const Real rho=-0.5;
>> >>>
>> >>> boost::shared_ptr<HestonProcess> process(new HestonProcess(
>> >>> riskFreeTS, dividendTS, s0, v0, kappa, theta, sigma,
>> >>> rho));
>> >>>
>> >>> boost::shared_ptr<HestonModel> model(new HestonModel(process));
>> >>>
>> >>> boost::shared_ptr<PricingEngine> engine(
>> >>> new
>> >>> AnalyticHestonEngine(model,
>> >>> 64));
>> >>>
>> >>> for (Size i = 0; i < options.size(); ++i)
>> >>> options[i]->setPricingEngine(engine);
>> >>>
>> >>> DifferentialEvolution::Configuration conf =
>> >>> DifferentialEvolution::Configuration()
>> >>> .withStepsizeWeight(1.8)
>> >>> .withBounds()
>> >>> .withCrossoverProbability(0.9)
>> >>> .withPopulationMembers(50)
>> >>>
>> >>> .withStrategy(DifferentialEvolution::Rand1SelfadaptiveWithRotation)
>> >>> .withCrossoverType(DifferentialEvolution::Normal)
>> >>> .withAdaptiveCrossover()
>> >>> .withSeed(3242);
>> >>> DifferentialEvolution deOptim(conf);
>> >>>
>> >>> Array lower(5);
>> >>> lower[0] = 0.001;
>> >>> lower[1] = 0.001;
>> >>> lower[2] = 0.001;
>> >>> lower[3] = -1.0;
>> >>> lower[4] = 0.001;
>> >>> Array upper(5);
>> >>> upper[0] = 10.0;
>> >>> upper[1] = 50.0;
>> >>> upper[2] = 10.0;
>> >>> upper[3] = 1.0;
>> >>> upper[4] = 10.0;
>> >>> model->calibrate(options, deOptim, EndCriteria(400, 40, 1.0e-8,
>> >>> 1.0e-8,
>> >>> Null<Real>()),
>> >>> NonhomogeneousBoundaryConstraint(lower, upper));
>> >>>
>> >>> Real sse = 0;
>> >>> for (Size i = 0; i < 13*8; ++i) {
>> >>> const Real diff = options[i]->calibrationError()*100.0;
>> >>> sse += diff*diff;
>> >>> }
>> >>> Real expected = 177.2; //see article by A. Sepp.
>> >>> if (std::fabs(sse - expected) > 1.0) {
>> >>> BOOST_FAIL("Failed to reproduce calibration error"
>> >>> << "\n calculated: " << sse
>> >>> << "\n expected: " << expected);
>> >>> }
>> >>> }
>> >>>
>> >>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>> >>>
>> >>> I get the code runnning and the test succeeds.
>> >>>
>> >>> Best
>> >>> Ralph
>> >>>
>> >>>
>> >>> 2014-02-17 13:00 GMT+01:00 Peter Caspers <[hidden email]>:
>> >>>>
>> >>>> Hi André,
>> >>>>
>> >>>>
>> >>>> can you please post the complete code needed to reproduce the problem
>> >>>> ? With a short description of what the problem is (if not obvious) ?
>> >>>>
>> >>>> best regards
>> >>>> Peter
>> >>>>
>> >>>>
>> >>>> On 16 February 2014 17:28, André de Boer <[hidden email]> wrote:
>> >>>> > Hello,
>> >>>> >
>> >>>> > Two days ago I posted a question about differential evolution and
>> >>>> > the
>> >>>> > shortratemodel G2.
>> >>>> > Meanwhile I learned how to minimize the Rosenbrockfunction with
>> >>>> > differential
>> >>>> > evolution.
>> >>>> > But my goal is to calibrate the G2 shortratemodel.
>> >>>> >
>> >>>> > For the 5 parameters of the G2++ model I want to have
>> >>>> > boundedconstraints:
>> >>>> > For instance:
>> >>>> > 0 < a < 6
>> >>>> > 0 > >>>> > 0 < sigma < 0.25
>> >>>> > 0 < eta < 0.25
>> >>>> > -1 < rho < 1
>> >>>> >
>> >>>> > To accomplish this I made a class MyConstraint, see code below.
>> >>>> > ...
>> >>>> > MyConstraint mc;
>> >>>> > DifferentialEvolution deOptim(conf);
>> >>>> > modelG2pp->calibrate(swaptions, deOptim,
>> >>>> > EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>> >>>> > ...
>> >>>> > But still I isn't not working.
>> >>>> > Can someone give me a hint of parts of code?
>> >>>> >
>> >>>> > Kind regards,
>> >>>> > André
>> >>>> >
>> >>>> >
>> >>>> >
>> >>>> > #include <iostream>
>> >>>> > #include <vector>
>> >>>> > #include <ql/quantlib.hpp>
>> >>>> > using namespace QuantLib;
>> >>>> >
>> >>>> > class MyConstraint : public Constraint {
>> >>>> > class Impl : public Constraint::Impl {
>> >>>> > public:
>> >>>> >
>> >>>> > bool test(const Array& x) const {
>> >>>> > Real a1 = x[0], a2 = x[1], a3 = x[2], a4=x[3], a5=x[4];
>> >>>> >
>> >>>> > return (0.0 <= a1 && a1 <= 6.0) &&
>> >>>> > (0.0 <= a2 && a2 <= 6.0) &&
>> >>>> > (0.0 <= a3 && a3 <= 0.25) &&
>> >>>> > (0.0 <= a4 && a4 <= 0.25) &&
>> >>>> > (-1.0 <= a5 && a5 <= 1.0);
>> >>>> >
>> >>>> > }
>> >>>> >
>> >>>> > };
>> >>>> > public:
>> >>>> >
>> >>>> > MyConstraint()
>> >>>> > : Constraint(boost::shared_ptr<Constraint::Impl>(new
>> >>>> > MyConstraint::Impl)) {}
>> >>>> >
>> >>>> > };
>> >>>> >
>> >>>> > int main(){
>> >>>> >
>> >>>> > ...
>> >>>> > double a = 0.1;
>> >>>> >
>> >>>> > double b = 0.1;
>> >>>> > double sigma = 0.01;
>> >>>> > double eta = 0.01;
>> >>>> > double rho = -0.75;
>> >>>> > boost::shared_ptr<G2> modelG2pp(new G2(yieldCurve, a, b, sigma,
>> >>>> > eta,
>> >>>> > rho));
>> >>>> > std::cout << "G2 (analytic formulae) calibration" << std::endl;
>> >>>> >
>> >>>> > MyConstraint mc;
>> >>>> > DifferentialEvolution::Configuration conf =
>> >>>> > DifferentialEvolution::Configuration()
>> >>>> > .withStepsizeWeight(0.4)
>> >>>> > .withBounds()
>> >>>> > .withCrossoverProbability(0.35)
>> >>>> > .withPopulationMembers(100)
>> >>>> > .withStrategy(DifferentialEvolution::BestMemberWithJitter)
>> >>>> > .withCrossoverType(DifferentialEvolution::Normal)
>> >>>> > .withAdaptiveCrossover()
>> >>>> > .withSeed(0);
>> >>>> >
>> >>>> > MyConstraint mc;
>> >>>> > DifferentialEvolution deOptim(conf);
>> >>>> > modelG2pp->calibrate(swaptions, deOptim,
>> >>>> > EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>> >>>> >
>> >>>> > }
>> >>>> >
>> >>>> >
>> >>>> >
>> >>>> > ------------------------------------------------------------------------------
>> >>>> > Android apps run on BlackBerry 10
>> >>>> > Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
>> >>>> > Now with support for Jelly Bean, Bluetooth, Mapview and more.
>> >>>> > Get your Android app in front of a whole new audience. Start now.
>> >>>> >
>> >>>> >
>> >>>> > http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
>> >>>> > _______________________________________________
>> >>>> > QuantLib-users mailing list
>> >>>> > [hidden email]
>> >>>> > https://lists.sourceforge.net/lists/listinfo/quantlib-users
>> >>>> >
>> >>>>
>> >>>>
>> >>>>
>> >>>> ------------------------------------------------------------------------------
>> >>>> Android apps run on BlackBerry 10
>> >>>> Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
>> >>>> Now with support for Jelly Bean, Bluetooth, Mapview and more.
>> >>>> Get your Android app in front of a whole new audience. Start now.
>> >>>>
>> >>>>
>> >>>> http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
>> >>>> _______________________________________________
>> >>>> QuantLib-users mailing list
>> >>>> [hidden email]
>> >>>> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>> >>>
>> >>>
>> >>
>> >
>
>

Peter Caspers-4

Re: calibration G2++ with defferential evolution

Hi Andre,

I forgot to amend the projected constraint, you can find the missing part here

https://github.com/lballabio/quantlib/pull/82

With that your code runs (forever ... :-) ) on my machine and returns
some calibrated model parameters finally.

Can you try on your side ?

Thanks
Peter

On 23 February 2014 23:49, André de Boer <[hidden email]> wrote:

> Hi Peter,
>
> I tried the fix but I get the message:
> Error: parameter array too big!
>
> This is my code:
>
> #include <iostream>
> #include <vector>
> #include <ql/quantlib.hpp>
> using namespace QuantLib;
> int main()
> {
> Date today(13, Sep, 2013);
> Natural settlementDays = 0;
> Calendar calendar = UnitedStates(UnitedStates::Settlement);
> DayCounter dayCounter = Actual360();
> DayCounter dayCounterFixed = Thirty360();
> Date settlementDate = calendar.advance(today, settlementDays, Days);
> // SWAP Curve
> std::vector<Date> dates; std::vector<Rate> rates;
> dates.push_back(settlementDate); rates.push_back(0.445/100.0);
> dates.push_back(calendar.advance(settlementDate,1,Years));
> rates.push_back(0.445/100.0);
> dates.push_back(calendar.advance(settlementDate,2,Years));
> rates.push_back(0.6426/100.0);
> dates.push_back(calendar.advance(settlementDate,3,Years));
> rates.push_back(0.883/100.0);
> dates.push_back(calendar.advance(settlementDate,4,Years));
> rates.push_back(1.148/100.0);
> dates.push_back(calendar.advance(settlementDate,5,Years));
> rates.push_back(1.3976/100.0);
> dates.push_back(calendar.advance(settlementDate,6,Years));
> rates.push_back(1.6065/100.0);
> dates.push_back(calendar.advance(settlementDate,7,Years));
> rates.push_back(1.7906/100.0);
> dates.push_back(calendar.advance(settlementDate,8,Years));
> rates.push_back(1.9501/100.0);
> dates.push_back(calendar.advance(settlementDate,9,Years));
> rates.push_back(2.091/100.0);
> dates.push_back(calendar.advance(settlementDate,10,Years));
> rates.push_back(2.22/100.0);
> dates.push_back(calendar.advance(settlementDate,12,Years));
> rates.push_back(2.425/100.0);
> dates.push_back(calendar.advance(settlementDate,15,Years));
> rates.push_back(2.6227/100.0);
> dates.push_back(calendar.advance(settlementDate,20,Years));
> rates.push_back(2.7328/100.0);
> dates.push_back(calendar.advance(settlementDate,25,Years));
> rates.push_back(2.7428/100.0);
> dates.push_back(calendar.advance(settlementDate,30,Years));
> rates.push_back(2.73/100.0);
> dates.push_back(calendar.advance(settlementDate,40,Years));
> rates.push_back(2.7525/100.0);
> dates.push_back(calendar.advance(settlementDate,50,Years));
> rates.push_back(2.7945/100.0);
> dates.push_back(calendar.advance(settlementDate,60,Years));
> rates.push_back(2.7945/100.0);
> dates.push_back(calendar.advance(settlementDate,61,Years));
> rates.push_back(2.7945/100.0);
> Handle<YieldTermStructure>
> yieldCurve(boost::shared_ptr<YieldTermStructure>(new ZeroCurve(dates,
> rates, dayCounter) ));
> boost::shared_ptr<IborIndex> index(new USDLibor(Period(3, Months), yieldCurve));
> Integer swapLengths[] = {1,2,3,4,5,6,7,8,9,10,15,20,25,30};
> Integer swaptionLengths[] = {1,2,3,4,5,6,7,10,15,20,25,30};
> Volatility swaptionVols[] = {
> 83.28, 65.22, 57.38, 51.62, 46.60, 41.99, 39.29, 36.28, 34.10, 32.67,
> 28.11, 26.00, 25.67, 25.30,
> 66.81, 52.88, 46.91, 42.32, 39.09, 36.35, 34.21, 32.34, 30.82, 29.85,
> 26.83, 25.46, 25.35, 25.11,
> 51.74, 43.57, 39.70, 36.58, 34.09, 32.31, 30.88, 29.53, 28.52, 27.87,
> 25.62, 24.67, 24.71, 24.50,
> 43.72, 37.47, 34.40, 32.28, 30.46, 29.25, 28.32, 27.44, 26.79, 26.25,
> 24.50, 24.17, 24.23, 23.90,
> 36.08, 32.42, 30.64, 29.11, 28.07, 26.95, 26.38, 25.91, 25.53, 25.26,
> 23.64, 23.68, 23.73, 23.36,
> 31.96, 29.30, 27.97, 26.77, 25.77, 25.28, 24.91, 24.66, 24.50, 24.40,
> 23.33, 23.01, 23.01, 22.57,
> 29.19, 26.94, 25.86, 25.02, 24.34, 24.06, 23.79, 23.60, 23.55, 23.52,
> 22.89, 22.51, 22.30, 21.92,
> 23.50, 22.61, 22.39, 22.25, 22.18, 22.31, 22.11, 22.08, 22.14, 22.24,
> 21.11, 21.20, 20.61, 19.90,
> 21.70, 22.10, 22.31, 22.51, 22.80, 22.87, 22.18, 22.31, 22.23, 22.23,
> 21.17, 19.40, 18.30, 17.40,
> 22.51, 23.07, 23.25, 23.42, 23.57, 23.45, 22.67, 22.45, 22.09, 21.84,
> 19.59, 17.35, 16.24, 15.38,
> 23.30, 23.47, 23.39, 23.25, 23.04, 22.54, 22.02, 21.48, 20.95, 20.39,
> 17.23, 15.13, 14.14, 13.43,
> 20.85, 20.61, 20.22, 19.88, 19.54, 19.06, 18.61, 18.20, 17.81, 17.41,
> 14.94, 13.20, 12.53, 11.84};
>
> try {
> // construct SwaptionHelpers
> std::vector<boost::shared_ptr<CalibrationHelper> > swaptions;
> std::list<Time> times;
> for (Size i=0; i<12; i++) { // Swaption
> for (Size j=0; j<14; j++) { // Swap
> boost::shared_ptr<Quote> vol(new SimpleQuote(swaptionVols[i*14+j]/100.0));
> swaptions.push_back(boost::shared_ptr<CalibrationHelper> (new
> SwaptionHelper(Period(swaptionLengths[i], Years),
> Period(swapLengths[j], Years), Handle<Quote>(vol),
> index, // floating index
> Period(6, Months), // fixed leg tenor
> dayCounterFixed, // fixed leg dayCounter
> dayCounter, // floating leg dayCounter
> yieldCurve)));
> swaptions.back()->addTimesTo(times);
> }
> }
> double a = 0.1;
> double sigma = 0.01;
> double b = 0.1;
> double eta = 0.01;
> double rho = -0.75;
> boost::shared_ptr<G2> modelG2pp(new G2(yieldCurve));
> std::cout << "G2 (analytic formulae) calibration" << std::endl;
> for (Size i=0; i<swaptions.size(); i++)
> swaptions[i]->setPricingEngine(boost::shared_ptr<PricingEngine>(new
> G2SwaptionEngine(modelG2pp, 6.0, 16)));
> std::cout << "\nCalibration with DE" << std::endl;
> DifferentialEvolution::Configuration conf =
> DifferentialEvolution::Configuration()
> .withStepsizeWeight(0.5)
> .withBounds()
> .withCrossoverProbability(0.4)
> .withPopulationMembers(100) //greater than 209 will erase an error
> .withStrategy(DifferentialEvolution::BestMemberWithJitter)
> .withCrossoverType(DifferentialEvolution::Normal)
> .withAdaptiveCrossover()
> .withSeed(314);
> DifferentialEvolution deOptim(conf);
> QuantLib::Array lowerBound(5);
> QuantLib::Array upperBound(5);
> lowerBound[0]=0.0;
> lowerBound[1]=0.0;
> lowerBound[2]=0.0;
> lowerBound[3]=0.0;
> lowerBound[4]=-1.0;
> upperBound[0]=5.0;
> upperBound[1]=0.25;
> upperBound[2]=5.0;
> upperBound[3]=0.25;
> upperBound[4]=1.0;
> QuantLib::NonhomogeneousBoundaryConstraint constraints(lowerBound,upperBound);
> modelG2pp->calibrate(swaptions, deOptim,
> EndCriteria(50,20,1.0e-5,1.0e-5,Null<Real>()), constraints);
> std::cout << "a = " << modelG2pp->params()[0]
> << "\nsigma = " << modelG2pp->params()[1]
> << "\nb = " << modelG2pp->params()[2]
> << "\neta = " << modelG2pp->params()[3]
> << "\nrho = " << modelG2pp->params()[4]
> << std::endl;
> std::cout << "\nCriterium: " << modelG2pp->endCriteria() << std::endl;
> Real sse = 0;
> for (Size i = 0; i < 14*12; ++i) {
> const Real diff = swaptions[i]->calibrationError()*100.0;
> sse += diff*diff;
> }
> std::cout << "\ncalculated calibration error: " << sse << std::endl;
> }
> catch(std::exception &e)
> {
> std::cout << "Error: " << e.what() << std::endl;
> }
> return 0;
> }
>
> 2014-02-23 22:42 GMT+01:00 André de Boer <[hidden email]>:
>> Hoi Peter,
>>
>> Thanks for your contribution to solve this issue.
>> I will try this patch and let you know.
>>
>> Kind regards,
>> André
>>
>>
>> 2014-02-22 7:59 GMT+01:00 Peter Caspers <[hidden email]>:
>>
>>> I tried to write a fix, which you can find here
>>>
>>> https://github.com/lballabio/quantlib/pull/80
>>>
>>> It compiles and the test suite does not break, but I did not test
>>> anything else. Maybe you can try to apply this patch and see if it
>>> solves the
>>> original problem ?
>>>
>>> best regards
>>> Peter
>>>
>>>
>>> On 18 February 2014 09:15, André de Boer <[hidden email]> wrote:
>>> > Hi Peter, Ralph,
>>> >
>>> > It still puzzles me that the CompositeConstraint has to be "fixed" as
>>> > you
>>> > states in your answer:
>>> > "fixing" the CompositeConstraint in constraint.hpp like this...".
>>> > Why is it not possible to use the class
>>> > "NonhomogeneousBoundaryConstraint"
>>> > in constraint.hpp directly?
>>> > It looks like a bug to me, but I am not an expert.
>>> >
>>> > 2014-02-17 17:51 GMT+01:00 André de Boer <[hidden email]>:
>>> >>
>>> >> Hi Ralph, Peter,
>>> >>
>>> >> Thanks for your comments and code, it looks promising and I will
>>> >> definitely use parts of it.
>>> >> Meanwhile I also made some progress
>>> >>
>>> >> Reimplementation of G2 with MyG2 with NoConstraint and defining
>>> >> NonhomogeneousBoundaryConstraint in main.cpp gets the code also
>>> >> running.
>>> >>
>>> >> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>> >> // file myG2.h
>>> >> #ifndef MYG2_H
>>> >> #define MYG2_H
>>> >>
>>> >> #include <ql/models/shortrate/twofactormodels/g2.hpp>
>>> >>
>>> >> class MyG2 : public QuantLib::G2
>>> >> {
>>> >> public:
>>> >> MyG2(const QuantLib::Handle<QuantLib::YieldTermStructure>&
>>> >> termStructure,
>>> >> QuantLib::Real a = 0.1,
>>> >> QuantLib::Real sigma = 0.01,
>>> >> QuantLib::Real b = 0.1,
>>> >> QuantLib::Real eta = 0.01,
>>> >> QuantLib::Real rho = -0.75);
>>> >> };
>>> >>
>>> >> #endif // MYG2_H
>>> >> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>> >>
>>> >> //file myG2.cpp
>>> >> #include "myg2.h"
>>> >> #include <ql/math/optimization/constraint.hpp>
>>> >>
>>> >> MyG2::MyG2(const QuantLib::Handle<QuantLib::YieldTermStructure>
>>> >> &termStructure,
>>> >> QuantLib::Real a,
>>> >> QuantLib::Real sigma,
>>> >> QuantLib::Real b,
>>> >> QuantLib::Real eta,
>>> >> QuantLib::Real rho) : G2(termStructure,a,sigma,b,eta,rho)
>>> >> {
>>> >> constraint_=boost::shared_ptr<QuantLib::NoConstraint>(new
>>> >> QuantLib::NoConstraint);
>>> >> }
>>> >>
>>> >> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>> >>
>>> >> // main.cpp
>>> >>
>>> >> ....
>>> >> boost::shared_ptr<MyG2> modelG2pp(new MyG2(yieldCurve, a, sigma, b,
>>> >> eta,
>>> >> rho));
>>> >>
>>> >> DifferentialEvolution::Configuration conf =
>>> >> DifferentialEvolution::Configuration()
>>> >> .withStepsizeWeight(0.4)
>>> >> .withBounds()
>>> >> .withCrossoverProbability(0.35)
>>> >> .withPopulationMembers(50)
>>> >>
>>> >> .withStrategy(DifferentialEvolution::BestMemberWithJitter)
>>> >> .withCrossoverType(DifferentialEvolution::Normal)
>>> >> .withAdaptiveCrossover()
>>> >> .withSeed(0);
>>> >>
>>> >> DifferentialEvolution deOptim(conf);
>>> >>
>>> >> QuantLib::Array lowerBound(5);
>>> >> QuantLib::Array upperBound(5);
>>> >> lowerBound[0]=0.0;
>>> >> lowerBound[1]=0.0;
>>> >> lowerBound[2]=0.0;
>>> >> lowerBound[3]=0.0;
>>> >> lowerBound[4]=-1.0;
>>> >>
>>> >> upperBound[0]=20.0;
>>> >> upperBound[1]=0.1;
>>> >> upperBound[2]=20.0;
>>> >> upperBound[3]=0.1;
>>> >> upperBound[4]=1.0;
>>> >>
>>> >> QuantLib::NonhomogeneousBoundaryConstraint
>>> >> constraint(lowerBound,upperBound);
>>> >> modelG2pp->calibrate(swaptions, deOptim,
>>> >> EndCriteria(400,40,1e-10,1e-8,Null<Real>()),constraint);
>>> >> ...
>>> >>
>>> >>
>>> >> 2014-02-17 17:01 GMT+01:00 Ralph Schreyer <[hidden email]>:
>>> >>
>>> >>> Hi Peter, André,
>>> >>>
>>> >>> "fixing" the CompositeConstraint in constraint.hpp like this (this is
>>> >>> not
>>> >>> a real fix, it is just to get the code running)
>>> >>>
>>> >>> >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
>>> >>> class CompositeConstraint : public Constraint {
>>> >>> private:
>>> >>> class Impl : public Constraint::Impl {
>>> >>> public:
>>> >>> Impl(const Constraint& c1,
>>> >>> const Constraint& c2)
>>> >>> : c1_(c1), c2_(c2) {}
>>> >>> bool test(const Array& params) const {
>>> >>> return c1_.test(params) && c2_.test(params);
>>> >>> }
>>> >>> Array upperBound(const Array& params) const {
>>> >>> Array c1ub = c1_.upperBound(params);
>>> >>> Array c2ub = c2_.upperBound(params);
>>> >>> Size ubSize = std::max(c1ub.size(), c2ub.size());
>>> >>> Array rtrnArray(ubSize, 0.0);
>>> >>> if (c1ub.size() == c2ub.size()) {
>>> >>> for (Size iter = 0; iter < c1ub.size(); iter++) {
>>> >>> rtrnArray.at(iter) = std::min(c1ub.at(iter), c2ub.at(iter));
>>> >>> }
>>> >>> }
>>> >>> else if (c1ub.size() > c2ub.size()) {
>>> >>> for (Size iter = 0; iter < c1ub.size(); iter++) {
>>> >>> rtrnArray.at(iter) = c1ub.at(iter);
>>> >>> }
>>> >>> }
>>> >>> else {
>>> >>> for (Size iter = 0; iter < c2ub.size(); iter++) {
>>> >>> rtrnArray.at(iter) = c2ub.at(iter);
>>> >>> }
>>> >>> }
>>> >>> return rtrnArray;
>>> >>> }
>>> >>> Array lowerBound(const Array& params) const {
>>> >>> Array c1lb = c1_.lowerBound(params);
>>> >>> Array c2lb = c2_.lowerBound(params);
>>> >>> Size lbSize = std::max(c1lb.size(), c2lb.size());
>>> >>> Array rtrnArray(lbSize, 0.0);
>>> >>> if (c1lb.size() == c2lb.size()) {
>>> >>> for (Size iter = 0; iter < c1lb.size(); iter++) {
>>> >>> rtrnArray.at(iter) = std::min(c1lb.at(iter), c2lb.at(iter));
>>> >>> }
>>> >>> }
>>> >>> else if (c1lb.size() > c2lb.size()) {
>>> >>> for (Size iter = 0; iter < c1lb.size(); iter++) {
>>> >>> rtrnArray.at(iter) = c1lb.at(iter);
>>> >>> }
>>> >>> }
>>> >>> else {
>>> >>> for (Size iter = 0; iter < c2lb.size(); iter++) {
>>> >>> rtrnArray.at(iter) = c2lb.at(iter);
>>> >>> }
>>> >>> }
>>> >>> return rtrnArray;
>>> >>> }
>>> >>> private:
>>> >>> Constraint c1_, c2_;
>>> >>> };
>>> >>> public:
>>> >>> CompositeConstraint(const Constraint& c1, const Constraint&
>>> >>> c2)
>>> >>> : Constraint(boost::shared_ptr<Constraint::Impl>(
>>> >>> new
>>> >>> CompositeConstraint::Impl(c1,c2))) {}
>>> >>> };
>>> >>>
>>> >>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>> >>>
>>> >>> and exchanging the LevenbergMarquardt optimizer in method
>>> >>> testDAXCalibration() for the DifferentialEvolution like this
>>> >>>
>>> >>> >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
>>> >>>
>>> >>> void HestonModelTest::testDAXCalibration() {
>>> >>>
>>> >>> BOOST_TEST_MESSAGE(
>>> >>> "Testing Heston model calibration using DAX volatility
>>> >>> data...");
>>> >>>
>>> >>> SavedSettings backup;
>>> >>>
>>> >>> Date settlementDate(5, July, 2002);
>>> >>> Settings::instance().evaluationDate() = settlementDate;
>>> >>>
>>> >>> CalibrationMarketData marketData = getDAXCalibrationMarketData();
>>> >>>
>>> >>> const Handle<YieldTermStructure> riskFreeTS =
>>> >>> marketData.riskFreeTS;
>>> >>> const Handle<YieldTermStructure> dividendTS =
>>> >>> marketData.dividendYield;
>>> >>> const Handle<Quote> s0 = marketData.s0;
>>> >>>
>>> >>> const std::vector<boost::shared_ptr<CalibrationHelper> > options
>>> >>> =
>>> >>> marketData.options;
>>> >>>
>>> >>> const Real v0=0.1;
>>> >>> const Real kappa=1.0;
>>> >>> const Real theta=0.1;
>>> >>> const Real sigma=0.5;
>>> >>> const Real rho=-0.5;
>>> >>>
>>> >>> boost::shared_ptr<HestonProcess> process(new HestonProcess(
>>> >>> riskFreeTS, dividendTS, s0, v0, kappa, theta, sigma,
>>> >>> rho));
>>> >>>
>>> >>> boost::shared_ptr<HestonModel> model(new HestonModel(process));
>>> >>>
>>> >>> boost::shared_ptr<PricingEngine> engine(
>>> >>> new
>>> >>> AnalyticHestonEngine(model,
>>> >>> 64));
>>> >>>
>>> >>> for (Size i = 0; i < options.size(); ++i)
>>> >>> options[i]->setPricingEngine(engine);
>>> >>>
>>> >>> DifferentialEvolution::Configuration conf =
>>> >>> DifferentialEvolution::Configuration()
>>> >>> .withStepsizeWeight(1.8)
>>> >>> .withBounds()
>>> >>> .withCrossoverProbability(0.9)
>>> >>> .withPopulationMembers(50)
>>> >>>
>>> >>> .withStrategy(DifferentialEvolution::Rand1SelfadaptiveWithRotation)
>>> >>> .withCrossoverType(DifferentialEvolution::Normal)
>>> >>> .withAdaptiveCrossover()
>>> >>> .withSeed(3242);
>>> >>> DifferentialEvolution deOptim(conf);
>>> >>>
>>> >>> Array lower(5);
>>> >>> lower[0] = 0.001;
>>> >>> lower[1] = 0.001;
>>> >>> lower[2] = 0.001;
>>> >>> lower[3] = -1.0;
>>> >>> lower[4] = 0.001;
>>> >>> Array upper(5);
>>> >>> upper[0] = 10.0;
>>> >>> upper[1] = 50.0;
>>> >>> upper[2] = 10.0;
>>> >>> upper[3] = 1.0;
>>> >>> upper[4] = 10.0;
>>> >>> model->calibrate(options, deOptim, EndCriteria(400, 40, 1.0e-8,
>>> >>> 1.0e-8,
>>> >>> Null<Real>()),
>>> >>> NonhomogeneousBoundaryConstraint(lower, upper));
>>> >>>
>>> >>> Real sse = 0;
>>> >>> for (Size i = 0; i < 13*8; ++i) {
>>> >>> const Real diff = options[i]->calibrationError()*100.0;
>>> >>> sse += diff*diff;
>>> >>> }
>>> >>> Real expected = 177.2; //see article by A. Sepp.
>>> >>> if (std::fabs(sse - expected) > 1.0) {
>>> >>> BOOST_FAIL("Failed to reproduce calibration error"
>>> >>> << "\n calculated: " << sse
>>> >>> << "\n expected: " << expected);
>>> >>> }
>>> >>> }
>>> >>>
>>> >>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>> >>>
>>> >>> I get the code runnning and the test succeeds.
>>> >>>
>>> >>> Best
>>> >>> Ralph
>>> >>>
>>> >>>
>>> >>> 2014-02-17 13:00 GMT+01:00 Peter Caspers <[hidden email]>:
>>> >>>>
>>> >>>> Hi André,
>>> >>>>
>>> >>>>
>>> >>>> can you please post the complete code needed to reproduce the problem
>>> >>>> ? With a short description of what the problem is (if not obvious) ?
>>> >>>>
>>> >>>> best regards
>>> >>>> Peter
>>> >>>>
>>> >>>>
>>> >>>> On 16 February 2014 17:28, André de Boer <[hidden email]> wrote:
>>> >>>> > Hello,
>>> >>>> >
>>> >>>> > Two days ago I posted a question about differential evolution and
>>> >>>> > the
>>> >>>> > shortratemodel G2.
>>> >>>> > Meanwhile I learned how to minimize the Rosenbrockfunction with
>>> >>>> > differential
>>> >>>> > evolution.
>>> >>>> > But my goal is to calibrate the G2 shortratemodel.
>>> >>>> >
>>> >>>> > For the 5 parameters of the G2++ model I want to have
>>> >>>> > boundedconstraints:
>>> >>>> > For instance:
>>> >>>> > 0 < a < 6
>>> >>>> > 0 >> >>>> > 0 < sigma < 0.25
>>> >>>> > 0 < eta < 0.25
>>> >>>> > -1 < rho < 1
>>> >>>> >
>>> >>>> > To accomplish this I made a class MyConstraint, see code below.
>>> >>>> > ...
>>> >>>> > MyConstraint mc;
>>> >>>> > DifferentialEvolution deOptim(conf);
>>> >>>> > modelG2pp->calibrate(swaptions, deOptim,
>>> >>>> > EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>>> >>>> > ...
>>> >>>> > But still I isn't not working.
>>> >>>> > Can someone give me a hint of parts of code?
>>> >>>> >
>>> >>>> > Kind regards,
>>> >>>> > André
>>> >>>> >
>>> >>>> >
>>> >>>> >
>>> >>>> > #include <iostream>
>>> >>>> > #include <vector>
>>> >>>> > #include <ql/quantlib.hpp>
>>> >>>> > using namespace QuantLib;
>>> >>>> >
>>> >>>> > class MyConstraint : public Constraint {
>>> >>>> > class Impl : public Constraint::Impl {
>>> >>>> > public:
>>> >>>> >
>>> >>>> > bool test(const Array& x) const {
>>> >>>> > Real a1 = x[0], a2 = x[1], a3 = x[2], a4=x[3], a5=x[4];
>>> >>>> >
>>> >>>> > return (0.0 <= a1 && a1 <= 6.0) &&
>>> >>>> > (0.0 <= a2 && a2 <= 6.0) &&
>>> >>>> > (0.0 <= a3 && a3 <= 0.25) &&
>>> >>>> > (0.0 <= a4 && a4 <= 0.25) &&
>>> >>>> > (-1.0 <= a5 && a5 <= 1.0);
>>> >>>> >
>>> >>>> > }
>>> >>>> >
>>> >>>> > };
>>> >>>> > public:
>>> >>>> >
>>> >>>> > MyConstraint()
>>> >>>> > : Constraint(boost::shared_ptr<Constraint::Impl>(new
>>> >>>> > MyConstraint::Impl)) {}
>>> >>>> >
>>> >>>> > };
>>> >>>> >
>>> >>>> > int main(){
>>> >>>> >
>>> >>>> > ...
>>> >>>> > double a = 0.1;
>>> >>>> >
>>> >>>> > double b = 0.1;
>>> >>>> > double sigma = 0.01;
>>> >>>> > double eta = 0.01;
>>> >>>> > double rho = -0.75;
>>> >>>> > boost::shared_ptr<G2> modelG2pp(new G2(yieldCurve, a, b, sigma,
>>> >>>> > eta,
>>> >>>> > rho));
>>> >>>> > std::cout << "G2 (analytic formulae) calibration" << std::endl;
>>> >>>> >
>>> >>>> > MyConstraint mc;
>>> >>>> > DifferentialEvolution::Configuration conf =
>>> >>>> > DifferentialEvolution::Configuration()
>>> >>>> > .withStepsizeWeight(0.4)
>>> >>>> > .withBounds()
>>> >>>> > .withCrossoverProbability(0.35)
>>> >>>> > .withPopulationMembers(100)
>>> >>>> > .withStrategy(DifferentialEvolution::BestMemberWithJitter)
>>> >>>> > .withCrossoverType(DifferentialEvolution::Normal)
>>> >>>> > .withAdaptiveCrossover()
>>> >>>> > .withSeed(0);
>>> >>>> >
>>> >>>> > MyConstraint mc;
>>> >>>> > DifferentialEvolution deOptim(conf);
>>> >>>> > modelG2pp->calibrate(swaptions, deOptim,
>>> >>>> > EndCriteria(100,10,1e-10,1e-8,Null<Real>()), mc);
>>> >>>> >
>>> >>>> > }
>>> >>>> >
>>> >>>> >
>>> >>>> >
>>> >>>> > ------------------------------------------------------------------------------
>>> >>>> > Android apps run on BlackBerry 10
>>> >>>> > Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
>>> >>>> > Now with support for Jelly Bean, Bluetooth, Mapview and more.
>>> >>>> > Get your Android app in front of a whole new audience. Start now.
>>> >>>> >
>>> >>>> >
>>> >>>> > http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
>>> >>>> > _______________________________________________
>>> >>>> > QuantLib-users mailing list
>>> >>>> > [hidden email]
>>> >>>> > https://lists.sourceforge.net/lists/listinfo/quantlib-users
>>> >>>> >
>>> >>>>
>>> >>>>
>>> >>>>
>>> >>>> ------------------------------------------------------------------------------
>>> >>>> Android apps run on BlackBerry 10
>>> >>>> Introducing the new BlackBerry 10.2.1 Runtime for Android apps.
>>> >>>> Now with support for Jelly Bean, Bluetooth, Mapview and more.
>>> >>>> Get your Android app in front of a whole new audience. Start now.
>>> >>>>
>>> >>>>
>>> >>>> http://pubads.g.doubleclick.net/gampad/clk?id=124407151&iu=/4140/ostg.clktrk
>>> >>>> _______________________________________________
>>> >>>> QuantLib-users mailing list
>>> >>>> [hidden email]
>>> >>>> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>>> >>>
>>> >>>
>>> >>
>>> >
>>
>>

------------------------------------------------------------------------------
Flow-based real-time traffic analytics software. Cisco certified tool.
Monitor traffic, SLAs, QoS, Medianet, WAAS etc. with NetFlow Analyzer
Customize your own dashboards, set traffic alerts and generate reports.
Network behavioral analysis & security monitoring. All-in-one tool.
http://pubads.g.doubleclick.net/gampad/clk?id=126839071&iu=/4140/ostg.clktrk
_______________________________________________
QuantLib-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-users