QuantLib › quantlib-dev

Binomial Tree Fever

Classic

List

Threaded

4 messages Options

John Maiden

Binomial Tree Fever

Think I found two errors in the Binomial Tree class. The problems are in
LeisenReimer:: LeisenReimer and Joshi4::Joshi4. In both constructors the
variance should be

Real variance = process->variance(0.0, x0_, dt_);

Instead of

Real variance = process->variance(0.0, x0_, end);

-------------------------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc.
Still grepping through log files to find problems? Stop.
Now Search log events and configuration files using AJAX and a browser.
Download your FREE copy of Splunk now >> http://get.splunk.com/
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Toyin Akin

Re: Binomial Tree Fever

Hi,

Perhaps not.

The variance looks like it is divided by the number of steps before it is being used.

Whereas in the other cases, the variance for a single dt step is taken and used directly.

Just an observation...

Toy out.

> To: [hidden email]
> From: [hidden email]
> Date: Mon, 27 Aug 2007 13:40:57 +0000
> Subject: [Quantlib-dev] Binomial Tree Fever
>
> Think I found two errors in the Binomial Tree class. The problems are in
> LeisenReimer:: LeisenReimer and Joshi4::Joshi4. In both constructors the
> variance should be
>
> Real variance = process->variance(0.0, x0_, dt_);
>
> Instead of
>
> Real variance = process->variance(0.0, x0_, end);
>
>
> -------------------------------------------------------------------------
> This SF.net email is sponsored by: Splunk Inc.
> Still grepping through log files to find problems? Stop.
> Now Search log events and configuration files using AJAX and a browser.
> Download your FREE copy of Splunk now >> http://get.splunk.com/
> _______________________________________________
> QuantLib-dev mailing list
> [hidden email]
> https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Play Movie Mash-up and win BIG prizes!
-------------------------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc.
Still grepping through log files to find problems? Stop.
Now Search log events and configuration files using AJAX and a browser.
Download your FREE copy of Splunk now >> http://get.splunk.com/
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

John Maiden

Re: Binomial Tree Fever

Yeah, I think you might be right. Thanks for noticing and replying so quickly.

-------------------------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc.
Still grepping through log files to find problems? Stop.
Now Search log events and configuration files using AJAX and a browser.
Download your FREE copy of Splunk now >> http://get.splunk.com/
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

John Maiden

Re: Binomial Tree Fever

Working on making extensions to the binomial tree classes.
The main change is that the tree now can use a time-dependent
stochastic process, i.e. a process that has risk-free rate and
volatility curves. I know that there are recombination issues
with a volatility curve, but as far as I know there isn't
anything wrong with incorporating a yield curve.

I've tested the trees with flat vol and different yield curves
(flat, treasury zero, and treasury strip), and half of the trees
(Jarrow-Rudd, Additive EQP, and Tian) produce normal results for
the flat and zero curve, but get very noisy for the strip curve.
The other trees are well-behaved for all the curves. Hoping
someone might be able to point out what's going wrong. I have
a spreadsheet of the results if anyone is interested in seeing
them.

Here's the code. I've renamed all the changed classes with the
extension "Extended":

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
Copyright (C) 2003 Ferdinando Ametrano
Copyright (C) 2001, 2002, 2003 Sadruddin Rejeb
Copyright (C) 2005 StatPro Italia srl

This file is part of QuantLib, a free-software/open-source library
for financial quantitative analysts and developers - http://quantlib.org/

QuantLib is free software: you can redistribute it and/or modify it
under the terms of the QuantLib license. You should have received a
copy of the license along with this program; if not, please email
<[hidden email]>. The license is also available online at
<http://quantlib.org/license.shtml>.

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

/*! \file binomialtree.hpp
\brief Binomial tree class
*/

#ifndef extended_binomial_tree_hpp
#define extended_binomial_tree_hpp

#include <ql/stochasticprocess.hpp>
#include <ql/methods/lattices/tree.hpp>
#include <ql/instruments/dividendschedule.hpp>
#include <ql/processes/blackscholesprocess.hpp>

#include <iostream>

namespace QuantLib {

//! Binomial tree base class
/*! \ingroup lattices */
template <class T>
class ExtendedBinomialTree : public Tree<T> {
public:
enum Branches { branches = 2 };
ExtendedBinomialTree(const boost::shared_ptr<StochasticProcess1D>& process,
Time end,
Size steps)
: Tree<T>(steps+1), treeProcess_(process) {
x0_ = process->x0();
dt_ = end/steps;
driftPerStep_ = process->drift(0.0, x0_) * dt_;
}
Size size(Size i) const {
return i+1;
}
Size descendant(Size, Size index, Size branch) const {
return index + branch;
}

const boost::shared_ptr<StochasticProcess1D> treeProcess() const {return
treeProcess_;}

//the amount the stock price moves up in one step
virtual Real upMove(Time stepTime) const = 0;
//the amount the stock price moves down in one step
virtual Real downMove(Time stepTime) const = 0;

protected:
//time dependent drift per step
Real driftStep(Time driftTime) const {
return this->treeProcess_->drift(driftTime, x0_) * dt_;
}

Real x0_, driftPerStep_;
Time dt_;

protected:
boost::shared_ptr<StochasticProcess1D> treeProcess_;
};

//! Base class for equal probabilities binomial tree
/*! \ingroup lattices */
template <class T>
class ExtendedEqualProbabilitiesBinomialTree : public ExtendedBinomialTree<T> {
public:
ExtendedEqualProbabilitiesBinomialTree(
const boost::shared_ptr<StochasticProcess1D>& process,
Time end,
Size steps)
: ExtendedBinomialTree<T>(process, end, steps) {}
//Real underlying(Size i, Size index) const {
// BigInteger j = 2*BigInteger(index) - BigInteger(i);
// // exploiting the forward value tree centering
// return this->x0_*std::exp(i*this->driftPerStep_ + j*this->up_);
//}
Real underlying(Size i, Size index, Time stepTime) const {
BigInteger j = 2*BigInteger(index) - BigInteger(i);
// exploiting the forward value tree centering
return this->x0_*std::exp(i*this->driftStep(stepTime) +
j*this->upStep(stepTime));
}

Real probability(Size, Size, Size, Time) const { return 0.5; }

Real upMove(Time stepTime) const {
return std::exp(this->driftStep(stepTime) + this->upStep(stepTime));
}

Real downMove(Time stepTime) const {
return std::exp(this->driftStep(stepTime) - this->upStep(stepTime));
}

//the tree dependent up move term at time stepTime
virtual Real upStep(Time stepTime) const = 0;

protected:
Real up_;
};

//! Base class for equal jumps binomial tree
/*! \ingroup lattices */
template <class T>
class ExtendedEqualJumpsBinomialTree : public ExtendedBinomialTree<T> {
public:
ExtendedEqualJumpsBinomialTree(
const boost::shared_ptr<StochasticProcess1D>& process,
Time end,
Size steps)
: ExtendedBinomialTree<T>(process, end, steps) {}
//Real underlying(Size i, Size index) const {
// BigInteger j = 2*BigInteger(index) - BigInteger(i);
// // exploiting equal jump and the x0_ tree centering
// return this->x0_*std::exp(j*this->dx_);
//}
Real underlying(Size i, Size index, Time stepTime) const {
BigInteger j = 2*BigInteger(index) - BigInteger(i);
// exploiting equal jump and the x0_ tree centering
return this->x0_*std::exp(j*this->dxStep(stepTime));
}

/* Real probability(Size, Size, Size branch) const {
return (branch == 1 ? pu_ : pd_);
}*/

Real probability(Size, Size, Size branch, Time stepTime) const {
Real upProb = this->probUp(stepTime);
Real downProb = 1 - upProb;
return (branch == 1 ? upProb : downProb);
}

Real upMove(Time stepTime) const {
return std::exp(this->dxStep(stepTime));
}

Real downMove(Time stepTime) const {
return (1/upMove(stepTime));
}

protected:
//probability of a up move
virtual Real probUp(Time stepTime) const = 0;
//time dependent term dx_
virtual Real dxStep(Time stepTime) const = 0;

Real dx_, pu_, pd_;
};

//! Jarrow-Rudd (multiplicative) equal probabilities binomial tree
/*! \ingroup lattices */
class ExtendedJarrowRudd : public
ExtendedEqualProbabilitiesBinomialTree<ExtendedJarrowRudd> {
public:
ExtendedJarrowRudd(const boost::shared_ptr<StochasticProcess1D>&,
Time end,
Size steps,
Real strike);
protected:
Real upStep(Time stepTime) const;
};

//! Cox-Ross-Rubinstein (multiplicative) equal jumps binomial tree
/*! \ingroup lattices */
class ExtendedCoxRossRubinstein
: public ExtendedEqualJumpsBinomialTree<ExtendedCoxRossRubinstein> {
public:
ExtendedCoxRossRubinstein(const boost::shared_ptr<StochasticProcess1D>&,
Time end,
Size steps,
Real strike);

protected:
Real dxStep(Time stepTime) const;
Real probUp(Time stepTime) const;
};

//! Additive equal probabilities binomial tree
/*! \ingroup lattices */
class ExtendedAdditiveEQPBinomialTree
: public
ExtendedEqualProbabilitiesBinomialTree<ExtendedAdditiveEQPBinomialTree> {
public:
ExtendedAdditiveEQPBinomialTree(
const boost::shared_ptr<StochasticProcess1D>&,
Time end,
Size steps,
Real strike);

protected:
Real upStep(Time stepTime) const;
};

//! %Trigeorgis (additive equal jumps) binomial tree
/*! \ingroup lattices */
class ExtendedTrigeorgis : public
ExtendedEqualJumpsBinomialTree<ExtendedTrigeorgis> {
public:
ExtendedTrigeorgis(const boost::shared_ptr<StochasticProcess1D>&,
Time end,
Size steps,
Real strike);
protected:
Real dxStep(Time stepTime) const;
Real probUp(Time stepTime) const;
};

//! %Tian tree: third moment matching, multiplicative approach
/*! \ingroup lattices */
class ExtendedTian : public ExtendedBinomialTree<ExtendedTian> {
public:
ExtendedTian(const boost::shared_ptr<StochasticProcess1D>&,
Time end,
Size steps,
Real strike);
/*Real underlying(Size i, Size index) const {
return x0_ * std::pow(down_, Real(BigInteger(i)-BigInteger(index)))
* std::pow(up_, Real(index));
};
Real probability(Size, Size, Size branch) const {
return (branch == 1 ? pu_ : pd_);
}*/

Real underlying(Size i, Size index, Time stepTime) const;
Real probability(Size, Size, Size branch, Time stepTime) const;

Real upMove(Time stepTime) const;
Real downMove(Time stepTime) const;

protected:
Real up_, down_, pu_, pd_;
};

//! Leisen & Reimer tree: multiplicative approach
/*! \ingroup lattices */
class ExtendedLeisenReimer : public ExtendedBinomialTree<ExtendedLeisenReimer> {
public:
ExtendedLeisenReimer(const boost::shared_ptr<StochasticProcess1D>&,
Time end,
Size steps,
Real strike);
/*Real underlying(Size i, Size index) const {
return x0_ * std::pow(down_, Real(BigInteger(i)-BigInteger(index)))
* std::pow(up_, Real(index));
}
Real probability(Size, Size, Size branch) const {
return (branch == 1 ? pu_ : pd_);
}*/

Real underlying(Size i, Size index, Time stepTime) const;
Real probability(Size, Size, Size branch, Time stepTime) const;

Real upMove(Time stepTime) const;
Real downMove(Time stepTime) const;

protected:
Real up_, down_, pu_, pd_, OddSteps_, Strike_, End_;
};

class ExtendedJoshi4 : public ExtendedBinomialTree<ExtendedJoshi4> {
public:
ExtendedJoshi4(const boost::shared_ptr<StochasticProcess1D>&,
Time end,
Size steps,
Real strike);
/*Real underlying(Size i, Size index) const {
return x0_ * std::pow(down_, Real(BigInteger(i)-BigInteger(index)))
* std::pow(up_, Real(index));
}
Real probability(Size, Size, Size branch) const {
return (branch == 1 ? pu_ : pd_);
}*/

Real underlying(Size i, Size index, Time stepTime) const;
Real probability(Size, Size, Size branch, Time stepTime) const;

Real upMove(Time stepTime) const;
Real downMove(Time stepTime) const;

protected:
Real computeUpProb(Real k, Real dj) const;
Real up_, down_, pu_, pd_, OddSteps_, Strike_, End_;
};

}

#endif

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
Copyright (C) 2003 Ferdinando Ametrano
Copyright (C) 2001, 2002, 2003 Sadruddin Rejeb
Copyright (C) 2005 StatPro Italia srl

This file is part of QuantLib, a free-software/open-source library
for financial quantitative analysts and developers - http://quantlib.org/

QuantLib is free software: you can redistribute it and/or modify it
under the terms of the QuantLib license. You should have received a
copy of the license along with this program; if not, please email
<[hidden email]>. The license is also available online at
<http://quantlib.org/license.shtml>.

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

#include "extendedbinomialtree.hpp"
#include <ql/math/distributions/binomialdistribution.hpp>

namespace QuantLib {

ExtendedJarrowRudd::ExtendedJarrowRudd(
const boost::shared_ptr<StochasticProcess1D>& process,
Time end, Size steps, Real)
: ExtendedEqualProbabilitiesBinomialTree<ExtendedJarrowRudd>(process, end,
steps) {
// drift removed
up_ = process->stdDeviation(0.0, x0_, dt_);
}

Real ExtendedJarrowRudd::upStep(Time stepTime) const {
return treeProcess_->stdDeviation(stepTime, x0_, dt_);
}

ExtendedCoxRossRubinstein::ExtendedCoxRossRubinstein(
const boost::shared_ptr<StochasticProcess1D>& process,
Time end, Size steps, Real)
: ExtendedEqualJumpsBinomialTree<ExtendedCoxRossRubinstein>(process, end,
steps) {

dx_ = process->stdDeviation(0.0, x0_, dt_);
pu_ = 0.5 + 0.5*this->driftStep(0.0)/dx_;
/*pu_ = 0.5 + 0.5*driftPerStep_/dx_;*/
pd_ = 1.0 - pu_;

QL_REQUIRE(pu_<=1.0, "negative probability");
QL_REQUIRE(pu_>=0.0, "negative probability");
}

Real ExtendedCoxRossRubinstein::dxStep(Time stepTime) const {
return this->treeProcess_->stdDeviation(stepTime, x0_, dt_);
}

Real ExtendedCoxRossRubinstein::probUp(Time stepTime) const {
return 0.5 + 0.5*this->driftStep(stepTime)/dxStep(stepTime);
}

ExtendedAdditiveEQPBinomialTree::ExtendedAdditiveEQPBinomialTree(
const boost::shared_ptr<StochasticProcess1D>& process,
Time end, Size steps, Real)
:
ExtendedEqualProbabilitiesBinomialTree<ExtendedAdditiveEQPBinomialTree>(process,
end, steps) {
/*up_ = - 0.5 * driftPerStep_ + 0.5 *
std::sqrt(4.0*process->variance(0.0, x0_, dt_)-
3.0*driftPerStep_*driftPerStep_);*/
up_ = - 0.5 * this->driftStep(0.0) + 0.5 *
std::sqrt(4.0*process->variance(0.0, x0_, dt_)-
3.0*this->driftStep(0.0)*this->driftStep(0.0));
}

Real ExtendedAdditiveEQPBinomialTree::upStep(Time stepTime) const {
return (- 0.5 * this->driftStep(stepTime) + 0.5 *
std::sqrt(4.0*this->treeProcess_->variance(stepTime, x0_, dt_)-
3.0*this->driftStep(stepTime)*this->driftStep(stepTime)));
}

ExtendedTrigeorgis::ExtendedTrigeorgis(
const boost::shared_ptr<StochasticProcess1D>& process,
Time end, Size steps, Real)
: ExtendedEqualJumpsBinomialTree<ExtendedTrigeorgis>(process, end, steps) {

/*dx_ = std::sqrt(process->variance(0.0, x0_, dt_)+
driftPerStep_*driftPerStep_);
pu_ = 0.5 + 0.5*driftPerStep_/dx_;*/
dx_ = std::sqrt(process->variance(0.0, x0_, dt_)+
this->driftStep(0.0)*this->driftStep(0.0));
pu_ = 0.5 + 0.5*this->driftStep(0.0)/this->dxStep(0.0);
pd_ = 1.0 - pu_;

QL_REQUIRE(pu_<=1.0, "negative probability");
QL_REQUIRE(pu_>=0.0, "negative probability");
}

Real ExtendedTrigeorgis::dxStep(Time stepTime) const {
return std::sqrt(this->treeProcess_->variance(stepTime, x0_, dt_)+
this->driftStep(stepTime)*this->driftStep(stepTime));
}

Real ExtendedTrigeorgis::probUp(Time stepTime) const {
return 0.5 + 0.5*this->driftStep(stepTime)/dxStep(stepTime);
}

ExtendedTian::ExtendedTian(const boost::shared_ptr<StochasticProcess1D>&
process,
Time end, Size steps, Real)
: ExtendedBinomialTree<ExtendedTian>(process, end, steps) {

Real q = std::exp(process->variance(0.0, x0_, dt_));
/*Real r = std::exp(driftPerStep_)*std::sqrt(q);*/

Real r = std::exp(this->driftStep(0.0))*std::sqrt(q);

up_ = 0.5 * r * q * (q + 1 + std::sqrt(q * q + 2 * q - 3));
down_ = 0.5 * r * q * (q + 1 - std::sqrt(q * q + 2 * q - 3));

pu_ = (r - down_) / (up_ - down_);
pd_ = 1.0 - pu_;

// doesn't work
// treeCentering_ = (up_+down_)/2.0;
// up_ = up_-treeCentering_;

QL_REQUIRE(pu_<=1.0, "negative probability");
QL_REQUIRE(pu_>=0.0, "negative probability");
}

Real ExtendedTian::underlying(Size i, Size index, Time stepTime) const {
Real q = std::exp(this->treeProcess_->variance(stepTime, x0_, dt_));
Real r = std::exp(this->driftStep(stepTime))*std::sqrt(q);

Real up = 0.5 * r * q * (q + 1 + std::sqrt(q * q + 2 * q - 3));
Real down = 0.5 * r * q * (q + 1 - std::sqrt(q * q + 2 * q - 3));

return x0_ * std::pow(down, Real(BigInteger(i)-BigInteger(index)))
* std::pow(up, Real(index));
}

Real ExtendedTian::probability(Size, Size, Size branch, Time stepTime) const {
Real q = std::exp(this->treeProcess_->variance(stepTime, x0_, dt_));
Real r = std::exp(this->driftStep(stepTime))*std::sqrt(q);

Real up = 0.5 * r * q * (q + 1 + std::sqrt(q * q + 2 * q - 3));
Real down = 0.5 * r * q * (q + 1 - std::sqrt(q * q + 2 * q - 3));

Real pu = (r - down) / (up - down);
Real pd = 1.0 - pu;

return (branch == 1 ? pu : pd);
}

Real ExtendedTian::upMove(Time stepTime) const {
Real q = std::exp(this->treeProcess_->variance(stepTime, x0_, dt_));
Real r = std::exp(this->driftStep(stepTime))*std::sqrt(q);

return 0.5 * r * q * (q + 1 + std::sqrt(q * q + 2 * q - 3));
}

Real ExtendedTian::downMove(Time stepTime) const {
Real q = std::exp(this->treeProcess_->variance(stepTime, x0_, dt_));
Real r = std::exp(this->driftStep(stepTime))*std::sqrt(q);

return 0.5 * r * q * (q + 1 - std::sqrt(q * q + 2 * q - 3));
}

ExtendedLeisenReimer::ExtendedLeisenReimer(
const boost::shared_ptr<StochasticProcess1D>& process,
Time end, Size steps, Real strike)
: ExtendedBinomialTree<ExtendedLeisenReimer>(process, end, (steps%2 ? steps
: steps+1)),
Strike_(strike), End_(end) {

QL_REQUIRE(strike>0.0, "strike must be positive");
Size oddSteps = (steps%2 ? steps : steps+1);
OddSteps_ = oddSteps;
Real variance = process->variance(0.0, x0_, end);
/*Real ermqdt = std::exp(driftPerStep_ + 0.5*variance/oddSteps);
Real d2 = (std::log(x0_/strike) + driftPerStep_*oddSteps ) /
std::sqrt(variance);*/
Real ermqdt = std::exp(this->driftStep(0.0) + 0.5*variance/oddSteps);
Real d2 = (std::log(x0_/strike) + this->driftStep(0.0)*oddSteps ) /
std::sqrt(variance);

pu_ = PeizerPrattMethod2Inversion(d2, oddSteps);
pd_ = 1.0 - pu_;
Real pdash = PeizerPrattMethod2Inversion(d2+std::sqrt(variance),
oddSteps);
up_ = ermqdt * pdash / pu_;
down_ = (ermqdt - pu_ * up_) / (1.0 - pu_);

}

Real ExtendedLeisenReimer::underlying(Size i, Size index, Time stepTime) const {
Real variance = this->treeProcess_->variance(stepTime, x0_, End_);
Real ermqdt = std::exp(this->driftStep(stepTime) + 0.5*variance/OddSteps_);
Real d2 = (std::log(x0_/Strike_) + this->driftStep(stepTime)*OddSteps_ ) /
std::sqrt(variance);

Real pu = PeizerPrattMethod2Inversion(d2, OddSteps_);
Real pd = 1.0 - pu;
Real pdash = PeizerPrattMethod2Inversion(d2+std::sqrt(variance),
OddSteps_);
Real up = ermqdt * pdash / pu;
Real down = (ermqdt - pu * up) / (1.0 - pu);

return x0_ * std::pow(down, Real(BigInteger(i)-BigInteger(index)))
* std::pow(up, Real(index));
}
Real ExtendedLeisenReimer::probability(Size, Size, Size branch, Time stepTime)
const {
Real variance = this->treeProcess_->variance(stepTime, x0_, End_);
Real ermqdt = std::exp(this->driftStep(stepTime) + 0.5*variance/OddSteps_);
Real d2 = (std::log(x0_/Strike_) + this->driftStep(stepTime)*OddSteps_ ) /
std::sqrt(variance);

Real pu = PeizerPrattMethod2Inversion(d2, OddSteps_);
Real pd = 1.0 - pu;

return (branch == 1 ? pu : pd);
}

Real ExtendedLeisenReimer::upMove(Time stepTime) const {
Real variance = this->treeProcess_->variance(stepTime, x0_, End_);
Real ermqdt = std::exp(this->driftStep(stepTime) + 0.5*variance/OddSteps_);
Real d2 = (std::log(x0_/Strike_) + this->driftStep(stepTime)*OddSteps_ ) /
std::sqrt(variance);

Real pu = PeizerPrattMethod2Inversion(d2, OddSteps_);
Real pd = 1.0 - pu;
Real pdash = PeizerPrattMethod2Inversion(d2+std::sqrt(variance),
OddSteps_);
return ermqdt * pdash / pu;
}

Real ExtendedLeisenReimer::downMove(Time stepTime) const {
Real variance = this->treeProcess_->variance(stepTime, x0_, End_);
Real ermqdt = std::exp(this->driftStep(stepTime) + 0.5*variance/OddSteps_);
Real d2 = (std::log(x0_/Strike_) + this->driftStep(stepTime)*OddSteps_ ) /
std::sqrt(variance);

Real pu = PeizerPrattMethod2Inversion(d2, OddSteps_);
Real pd = 1.0 - pu;
Real pdash = PeizerPrattMethod2Inversion(d2+std::sqrt(variance),
OddSteps_);
Real up = ermqdt * pdash / pu;
return (ermqdt - pu * up) / (1.0 - pu);
}

Real ExtendedJoshi4::computeUpProb(Real k, Real dj) const {
Real alpha = dj/(sqrt(8.0));
Real alpha2 = alpha*alpha;
Real alpha3 = alpha*alpha2;
Real alpha5 = alpha3*alpha2;
Real alpha7 = alpha5*alpha2;
Real beta = -0.375*alpha-alpha3;
Real gamma = (5.0/6.0)*alpha5 + (13.0/12.0)*alpha3
+(25.0/128.0)*alpha;
Real delta = -0.1025 *alpha- 0.9285 *alpha3
-1.43 *alpha5 -0.5 *alpha7;
Real p =0.5;
Real rootk= sqrt(k);
p+= alpha/rootk;
p+= beta /(k*rootk);
p+= gamma/(k*k*rootk);
// delete next line to get results for j three tree
p+= delta/(k*k*k*rootk);
return p;
}

ExtendedJoshi4::ExtendedJoshi4(const boost::shared_ptr<StochasticProcess1D>&
process,
Time end, Size steps, Real strike)
: ExtendedBinomialTree<ExtendedJoshi4>(process, end, (steps%2 ? steps :
steps+1)),
Strike_(strike), End_(end) {

QL_REQUIRE(strike>0.0, "strike must be positive");
Size oddSteps = (steps%2 ? steps : steps+1);
OddSteps_ = oddSteps;
Real variance = process->variance(0.0, x0_, end);
/*Real ermqdt = std::exp(driftPerStep_ + 0.5*variance/oddSteps);
Real d2 = (std::log(x0_/strike) + driftPerStep_*oddSteps ) /
std::sqrt(variance);*/
Real ermqdt = std::exp(this->driftStep(0.0) + 0.5*variance/oddSteps);
Real d2 = (std::log(x0_/strike) + this->driftStep(0.0)*oddSteps ) /
std::sqrt(variance);

pu_ = computeUpProb((oddSteps-1.0)/2.0,d2 );
pd_ = 1.0 - pu_;
Real pdash = computeUpProb((oddSteps-1.0)/2.0,d2+std::sqrt(variance));
up_ = ermqdt * pdash / pu_;
down_ = (ermqdt - pu_ * up_) / (1.0 - pu_);
}

Real ExtendedJoshi4::underlying(Size i, Size index, Time stepTime) const {
Real variance = this->treeProcess_->variance(stepTime, x0_, End_);
Real ermqdt = std::exp(this->driftStep(stepTime) + 0.5*variance/OddSteps_);
Real d2 = (std::log(x0_/Strike_) + this->driftStep(stepTime)*OddSteps_ ) /
std::sqrt(variance);

Real pu = computeUpProb((OddSteps_-1.0)/2.0,d2 );
Real pd = 1.0 - pu;
Real pdash = computeUpProb((OddSteps_-1.0)/2.0,d2+std::sqrt(variance));
Real up = ermqdt * pdash / pu;
Real down = (ermqdt - pu * up) / (1.0 - pu);

return x0_ * std::pow(down, Real(BigInteger(i)-BigInteger(index)))
* std::pow(up, Real(index));
}

Real ExtendedJoshi4::probability(Size, Size, Size branch, Time stepTime) const {
Real variance = this->treeProcess_->variance(stepTime, x0_, End_);
Real ermqdt = std::exp(this->driftStep(stepTime) + 0.5*variance/OddSteps_);
Real d2 = (std::log(x0_/Strike_) + this->driftStep(stepTime)*OddSteps_ ) /
std::sqrt(variance);

Real pu = computeUpProb((OddSteps_-1.0)/2.0,d2 );
Real pd = 1.0 - pu;

return (branch == 1 ? pu : pd);
}

Real ExtendedJoshi4::upMove(Time stepTime) const {
Real variance = this->treeProcess_->variance(stepTime, x0_, End_);
Real ermqdt = std::exp(this->driftStep(stepTime) + 0.5*variance/OddSteps_);
Real d2 = (std::log(x0_/Strike_) + this->driftStep(stepTime)*OddSteps_ ) /
std::sqrt(variance);

Real pu = computeUpProb((OddSteps_-1.0)/2.0,d2 );
Real pd = 1.0 - pu;
Real pdash = computeUpProb((OddSteps_-1.0)/2.0,d2+std::sqrt(variance));
return ermqdt * pdash / pu;
}

Real ExtendedJoshi4::downMove(Time stepTime) const {
Real variance = this->treeProcess_->variance(stepTime, x0_, End_);
Real ermqdt = std::exp(this->driftStep(stepTime) + 0.5*variance/OddSteps_);
Real d2 = (std::log(x0_/Strike_) + this->driftStep(stepTime)*OddSteps_ ) /
std::sqrt(variance);

Real pu = computeUpProb((OddSteps_-1.0)/2.0,d2 );
Real pd = 1.0 - pu;
Real pdash = computeUpProb((OddSteps_-1.0)/2.0,d2+std::sqrt(variance));
Real up = ermqdt * pdash / pu;
return (ermqdt - pu * up) / (1.0 - pu);
}

}

-------------------------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc.
Still grepping through log files to find problems? Stop.
Now Search log events and configuration files using AJAX and a browser.
Download your FREE copy of Splunk now >> http://get.splunk.com/
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev