quantlib-dev — Discussion of QuantLib development

Re: [Quantlib-dev] Projecting Forward Volatility Surfaces

[Luigi B. <lui...@gm...>]

Hi David,
	I'm afraid I can't help with the problem. However...

On Wed, 2007-08-29 at 22:03 -0400, David Brown wrote:
>    So I've raised this question here, because I think this is the only
> place where there is a possibility of this problem being solved given
> the expertise of the parties involved in similar projects.

...you might try the Wilmott forums as well.

Later,
	Luigi


-- 

Green's Law of Debate: 
Anything is possible if you don't know what you're talking about. 

Re: [Quantlib-dev] Separating notification channels ?

[DU V. DE V. F. G. P. <fra...@ca...>]

Hi Luigi,



>> In our implementation market data and evaluation date changes are
>> propagated through the same channel (observers update methods).=20
>>=20
>> As we are properly trying to handle properly the evaluation date
>> changes this coupling is becoming more and more an hindrance (eg: the
>> Capstripper can no longer notifying lazily its observers).

>Can you elaborate? What is the problem?

The problem is mainly due to the fact that the CapsStripper structure is
creating some huge reconvergencies in the notification chain. (It
observes hundreds of Caps which all observe the same YC for example). As
a result I suggested to make the LazyObject propagate notifications
lazily.
The issue now comes from the fact that the reference date is managed in
the TermStructure from which the CapsStripper inherits. This class does
not notifiy lazily so we are back to the original problem: a change in
the market data triggers thousands of notification events (even millions
sometimes!).=20

Francois

Re: [Quantlib-dev] Minor bug in date output

[Ferdinando A. <na...@am...>]

On 9/3/07, Robert Lopez <rz...@gm...> wrote:
> I recently installed QuantLib-0.8.1 and think I've run into a minor
> bug in the ostream operator << for short dates. The output adds an
> extra "/" between the day and year fields. From date.cpp, lines
> 407-409:

Thank you Robert, it has been fixed in the repository

ciao -- Nando

[Quantlib-dev] Minor bug in date output

[Robert L. <rz...@gm...>]

I recently installed QuantLib-0.8.1 and think I've run into a minor
bug in the ostream operator << for short dates. The output adds an
extra "/" between the day and year fields. From date.cpp, lines
407-409:

>                out << std::setw(2) << std::setfill('0') << mm << "/";
>                out << std::setw(2) << std::setfill('0') << dd << "/";
>                out << "/" << yyyy;

I think that this code should read

>                out << std::setw(2) << std::setfill('0') << mm << "/";
>                out << std::setw(2) << std::setfill('0') << dd << "/";
>                out << yyyy;

Hope this is helpful.

-- 

 - Rob

Re: [Quantlib-dev] about calendars

[Richard G. <rgo...@ya...>]

Luigi Ballabio wrote:
> On Mon, 2007-09-03 at 12:57 +0100, Toyin Akin wrote:
>> If you take a look at the US or UK calendar classes, there are
>> examples of implementing more than one holiday schedule within the
>> same financial location.
> 
> Right. In your case, there would be Brazil::BACEN and Brazil::BACESPA
> calendars.
>  
>> You can, if you wish, use this as an example for adding a new holiday 
>> schedule.
> 
> And of course, we'll be happy to add your patch to the library if you
> were to contribute it.
> 
> Later,
> 	Luigi
> 
> 

I'll be proud to contribute.

Thanks

Richard

Re: [Quantlib-dev] about calendars

[Toyin A. <toy...@ho...>]

Hi Richard,
=20
If you take a look at the US or UK calendar classes, there are examples of =
implementing more than one holiday schedule within the same financial locat=
ion.
=20
In the case of UK (financial location), I believe there are settlement as w=
ell as Metals holiday rules.
=20
You can, if you wish, use this as an example for adding a new holiday sched=
ule.
=20
Toy out.
> To: qua...@li...> From: rgo...@ya...> Dat=
e: Sun, 2 Sep 2007 21:13:03 +0100> Subject: [Quantlib-dev] about calendars>=
 > Hi guys,> > I implemented an Option valuation procedure written in Java =
and I needed to> calculate the working days. As I'm brazilian, I started wi=
th Brazilian> calendar.> > In QuantLib, it is defined here> http://quantlib=
.org/reference/class_quant_lib_1_1_brazil.html> > It is correct, except tha=
t BOVESPA (the brazilian exchange house) will be> closed on these dates as =
well:> > 25-January : a municipal holiday> 09-July : a state holiday> > My =
question is: Shouldn't we consider all dates when the exchange house will> =
be closed?> > Thanks> > -- > Richard Gomes> > > ---------------------------=
----------------------------------------------> This SF.net email is sponso=
red by: Splunk Inc.> Still grepping through log files to find problems? Sto=
p.> 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> QuantLi=
b-...@li...> https://lists.sourceforge.net/lists/listinfo/q=
uantlib-dev
_________________________________________________________________
100=92s of Music vouchers to be won with MSN Music
https://www.musicmashup.co.uk=

Re: [Quantlib-dev] about calendars

[Luigi B. <lui...@gm...>]

On Mon, 2007-09-03 at 12:57 +0100, Toyin Akin wrote:
> If you take a look at the US or UK calendar classes, there are
> examples of implementing more than one holiday schedule within the
> same financial location.

Right. In your case, there would be Brazil::BACEN and Brazil::BACESPA
calendars.
 
> You can, if you wish, use this as an example for adding a new holiday 
> schedule.

And of course, we'll be happy to add your patch to the library if you
were to contribute it.

Later,
	Luigi


-- 

Do the right thing. It will gratify some people and astonish the rest. 
-- Mark Twain 

Re: [Quantlib-dev] about calendars

[Toyin A. <toy...@ho...>]

Hi Richard,
If you take a look at the US or UK calendar classes, there are examples of =
implementing more than one holiday schedule within the same financial locat=
ion. In the case of UK (financial location), I believe there are settlement=
 as well as Metals holiday rules.You can, if you wish, use this as an examp=
le for adding a new holiday schedule.Toy out.> To: qua...@li...=
eforge.net> From: rgo...@ya...> Date: Sun, 2 Sep 2007 21:13:03 +0=
100> Subject: [Quantlib-dev] about calendars> > Hi guys,> > I implemented a=
n Option valuation procedure written in Java and I needed to> calculate the=
 working days. As I'm brazilian, I started with Brazilian> calendar.> > In =
QuantLib, it is defined here> http://quantlib.org/reference/class_quant_lib=
_1_1_brazil.html> > It is correct, except that BOVESPA (the brazilian excha=
nge house) will be> closed on these dates as well:> > 25-January : a munici=
pal holiday> 09-July : a state holiday> > My question is: Shouldn't we cons=
ider all dates when the exchange house will> be closed?> > Thanks> > -- > R=
ichard Gomes> > > ---------------------------------------------------------=
----------------> This SF.net email is sponsored by: Splunk Inc.> Still gre=
pping through log files to find problems? Stop.> Now Search log events and =
configuration files using AJAX and a browser.> Download your FREE copy of S=
plunk now >> http://get.splunk.com/> ______________________________________=
_________> QuantLib-dev mailing list> Qua...@li...> h=
ttps://lists.sourceforge.net/lists/listinfo/quantlib-dev
_________________________________________________________________
Get free emoticon packs and customisation from Windows Live.=20
http://www.pimpmylive.co.uk=

Re: [Quantlib-dev] Projecting Forward Volatility Surfaces

[Joseph W. <jo...@gn...>]

One way of doing this is the Carr-Mandan technique.  If you have a set of 
liquid vanilla european option prices at a particular time, you can use a 
damped Fast Fourier Transform to get the volatility surface for a particular 
time and a given set of parameter values.  You then minimize least squares 
over the parameter values.

The problem is that this only works if you have a set of liquid European 
option.  The major piece of infrastructure that isn't in QuantLib is the 
ability to do fast fourier transforms plus maybe some C++ data structures to 
handle the data.

-- 
-------------------------------------------------------------------------------
Joseph Wang Ph.D. - jo...@gn...  
China Derivatives Researcher and Software Developer
http://en.wikiversity.org/wiki/User:Roadrunner

[Quantlib-dev] about calendars

[Piter D. <pit...@ma...>]

Richard,

I made this calendar based on Banco Central´s one in order to work with some
fixed income instruments
You could add São Paulo one (people already do it for another countries). So
we would have Bacen and Bovespa ones. That is the easy part.

However, fixed income always follows Banco Central rules. So the interest
should follow Banco Central one.
I another hand you will have stock price following Bovespa one.

I never thought about that once BM&F (derivatives Exchange) prices always
follow Banco Central for calculation. I never played to much with Brazilian
stock options.

Let me know how do you plan to deal with this two calendars in your model.

Regards,

Piter Dias
pit...@ca...

[Quantlib-dev] about calendars

[Piter D. <pit...@ca...>]

Richard,

I made this calendar based on Banco Central´s one in order to work with some
fixed income instruments
You could add São Paulo one (people already do it for another countries). So
we would have Bacen and Bovespa ones. That is the easy part.

However, fixed income always follows Banco Central rules. So the interest
should follow Banco Central one.
I another hand you will have stock price following Bovespa one.

I never thought about that once BM&F (derivatives Exchange) prices always
follow Banco Central for calculation. I never played to much with Brazilian
stock options.

Let me know how do you plan to deal with this two calendars in your model.

Regards,

Piter Dias
pit...@ca...

[Quantlib-dev] swig and java bindings

[Richard G. <rgo...@ya...>]

Hi Chaps,

I have some code already written in Java, using third party math libs and
I'm evaluating the possibility of adopting QuantLib as far as possible.

I've noticed that swig did an incomplete job, for instance:

//
// This is a Java code trying to use Matrix from QuantLib
//

// Create 2 matrices
Matrix a = new Matrix(4, 4);
Matrix b = new Matrix(4, 4);

// SWIG did not converted C++: operator[] to something like
// Java: Array get(int)
a[0][3] = 5; // it does not compile in Java

// SWIG did not converted C++: operator* to something like
// Java: Matrix multiply(Matrix)
Matrix c = a * b; // it does not compile in Java

How it can be done?

Thanks in advance.

-- 
Richard Gomes

[Quantlib-dev] about calendars

[Richard G. <rgo...@ya...>]

Hi guys,

I implemented an Option valuation procedure written in Java and I needed to
calculate the working days. As I'm brazilian, I started with Brazilian
calendar.

In QuantLib, it is defined here
http://quantlib.org/reference/class_quant_lib_1_1_brazil.html

It is correct, except that BOVESPA (the brazilian exchange house) will be
closed on these dates as well:

25-January : a municipal holiday
09-July    : a state holiday

My question is: Shouldn't we consider all dates when the exchange house will
be closed?

Thanks

-- 
Richard Gomes

Re: [Quantlib-dev] Help! newbie question : autogen

[Richard G. <rgo...@ya...>]

Richard Gomes wrote:

> Dirk Eddelbuettel wrote:
> 
>> 
>> On 2 September 2007 at 15:40, Richard Gomes wrote:
>> | You are right... libtool was not installed. :(
>> 
>> You can use Debian's meta-info for that, i.e.
>> 
>> $ apt-get build-dep quantlib
>> 
>> gets you all packages needed to build the main quantlib packages.
>> And
>> 
>> $ apt-get source quantlib
>> 
>> will get the source for you --- if you point to Debian unstable you'd
>> always get the latest.
>> 
>> Hth, Dirk
>> 
> 
> Dirk, Thank you very much.
> 
> In fact, I'm currently using Kubuntu 7.04 and Debian's unstable packages
> do not work well on it. Yesterday I had to use killall to kill thousands
> of apt-cache processes slowing down my box.
> 
> I'm planning to switch to Debian4 in the future, but I'll need to have
> enough spare time and courage to do it.
> 
> Kind Regards
> 

Hi Chaps,

I managed to get everything running.
This is a quick procedure, specially intended to AMD64 addicted ones.

-----------------------------------------------------------------------

echo Setting JAVA_HOME to your JDK flavored for AMD64.
# To be on the safe side, I tested against Java5.
export JAVA_HOME=/opt/JavaIDE/jdk1.5.0_12-amd64

cd $HOME
mkdir svn_working-copy
cd svn_working-copy

echo Download QuantLib
svn checkout
https://quantlib.svn.sourceforge.net/svnroot/quantlib/trunk/QuantLib

echo Download QuantLib-SWIG
svn checkout
https://quantlib.svn.sourceforge.net/svnroot/quantlib/trunk/QuantLib-SWIG

echo Install required packages
sudo apt-get install binutils gcc make autoconf automake libtool swig
libboost-dev

cd $HOME/svn_working-copy
cd QuantLib
echo Configure QuantLib
sh autogen.sh
./configure
make
sudo make install

cd $HOME/svn_working-copy
cd QuantLib-SWIG
echo Configure QuantLib-SWIG for Java environment
sh autogen.sh
echo "================================"
echo "Setting CXXFLAGS=-fPIC for AMD64"
echo "================================"
export CXXFLAGS=-fPIC
echo
./configure \
   --with-jdk-include=${JAVA_HOME}/include \
   --with-jdk-system-include=${JAVA_HOME}/include/linux
make -C Java
sudo make -C Java install

echo Verify installed files
ls -ald /usr/local/lib/*QuantLib*

echo Verify file contents
file /usr/local/lib/*QuantLib*



You should see something similar to...

/usr/local/lib/libQuantLib-0.8.1.so: ELF 64-bit LSB shared object, x86-64,
version 1 (SYSV), not stripped
/usr/local/lib/libQuantLibJNI.so:    ELF 64-bit LSB shared object, x86-64,
version 1 (SYSV), not stripped
/usr/local/lib/libQuantLib.so:       symbolic link to `libQuantLib-0.8.1.so'
/usr/local/lib/libQuantLib-0.8.1.a:  current ar archive
/usr/local/lib/libQuantLib.a:        symbolic link to `libQuantLib-0.8.1.a'
/usr/local/lib/QuantLib.jar:         Zip archive data, at least v2.0 to
extract


To be honest, I've build from released packages and not from the trunk SVN.
If you follow this procedure, probably some file names can be slightly
different.


Now, include /usr/local/lib/QuantLib.jar in the build path of your Java
projects.
Use this test program:


package org.quantlib.examples;

import org.quantlib.QuantLib;
import org.quantlib.Weekday;

public class Hello {
    static {
        try {
            System.loadLibrary("QuantLibJNI");
        } catch (RuntimeException e) {
            e.printStackTrace();
        }
    }
    public static void main(String[] args) {
        QuantLib ql = new QuantLib();
        double nd = QuantLib.nullDouble();

        System.out.println("ql nullDouble is "+nd);
        System.out.println("and sunday is "+Weekday.Sunday.toString());
    }
}

Go to the command line and execute...

$ export LD_LIBRARY_PATH=/usr/local/lib
$ java -cp .:/usr/local/lib/QuantLib.jar org.quantlib.examples.Hello

-----------------------------------------------------------------------

Good Luck! :)

-- 
Richard Gomes

Re: [Quantlib-dev] Help! newbie question : autogen

[Richard G. <rgo...@ya...>]

Dirk Eddelbuettel wrote:

> 
> On 2 September 2007 at 15:40, Richard Gomes wrote:
> | You are right... libtool was not installed. :(
> 
> You can use Debian's meta-info for that, i.e.
> 
> $ apt-get build-dep quantlib
> 
> gets you all packages needed to build the main quantlib packages.
> And
> 
> $ apt-get source quantlib
> 
> will get the source for you --- if you point to Debian unstable you'd
> always get the latest.
> 
> Hth, Dirk
> 

Dirk, Thank you very much.

In fact, I'm currently using Kubuntu 7.04 and Debian's unstable packages do
not work well on it. Yesterday I had to use killall to kill thousands of
apt-cache processes slowing down my box.

I'm planning to switch to Debian4 in the future, but I'll need to have
enough spare time and courage to do it.

Kind Regards

-- 
Richard Gomes

Re: [Quantlib-dev] Help! newbie question : autogen

[Dirk E. <ed...@de...>]

On 2 September 2007 at 15:40, Richard Gomes wrote:
| You are right... libtool was not installed. :(

You can use Debian's meta-info for that, i.e.

	$ apt-get build-dep quantlib

gets you all packages needed to build the main quantlib packages.
And

	$ apt-get source quantlib

will get the source for you --- if you point to Debian unstable you'd always
get the latest.

Hth, Dirk

-- 
Three out of two people have difficulties with fractions.

Re: [Quantlib-dev] [Quantlib-users] multidimensional integration

[Luigi B. <lui...@gm...>]

On Sep 1, 2007, at 9:03 PM, Warren Chou wrote:
> Does Quantlib have any multidimensional integration routines?

No, it doesn't.

> If not, do the 1-D integration routines lend themselves to 
> multidimensional integration?

Possibly. For 2-D integration of a function f(x,y), you can use an 
intermediate function object to return the 1-D integral of f(x,y) along 
y for a fixed x---let's call this function F(x). Then you can integrate 
F(x). For 3-D integration, use the above to get the 2-D integral G(z) 
along x,y for a given z and integrate G(z). Recurse as many times as 
needed.

Later,
	Luigi

Re: [Quantlib-dev] Help! newbie question : autogen

[Richard G. <rgo...@ya...>]

Luigi Ballabio wrote:

> 
> On Sep 2, 2007, at 3:52 PM, Richard Gomes wrote:
>>> I compiled/installed the latest boost release.
>>> I just checked out QuantLib from SVN.
>>> I installed autoconf, automake and m4 packages.
>>>
>>> and then...
>>>
>>> $ sh autogen.sh
>>> /usr/share/aclocal/libmcrypt.m4:17: warning: underquoted definition of
>>> AM_PATH_LIBMCRYPT
>>> /usr/share/aclocal/libmcrypt.m4:17:   run info '(automake)Extending
>>> aclocal'
>>> /usr/share/aclocal/libmcrypt.m4:17:   or see
>>> http://sources.redhat.com/automake/automake.html#Extending-aclocal
>>> configure.ac:63: error: possibly undefined macro: AC_PROG_LIBTOOL
>>>       If this token and others are legitimate, please use
>>>       m4_pattern_allow. See the Autoconf documentation.
>>> autoreconf: /usr/bin/autoconf failed with exit status: 1
> 
> Did you install libtool as well as autoconf and automake? it seems that
> autoconf does not find an m4 macro file provided with libtool.
> 
> 
>> This is a followup to the group.
>>
>> Version 0.8.0 has the c++ sources and swig sources in download area at
>> sourceforge.net.
> 
> Also, it's probably not very clear from the release notes, but
> QuantLib-SWIG 0.8.0 will work with QuantLib 0.8.1.
> 
> Later,
> Luigi
> 
> 
> 
> -------------------------------------------------------------------------
> 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/

Hi Luigi,

I managed to get v0.8.0 compiled.
Yes... I guessed that swig-0.8.0 could fit v0.8.1 as well.
But I decided to be on the safe side.
I'm compiling v0.8.1 now.

You are right... libtool was not installed. :(
After I installed it, I managed to get autogen.sh running fine, in spite of
some warnings. At least "configure" was created and runs OK.
I'm also compiling sources from trunk.

Soon I will need a 2 head amd64x2 :)

Kind Regards

-- 
Richard Gomes

Re: [Quantlib-dev] Help! newbie question : autogen

[Luigi B. <lui...@gm...>]

On Sep 2, 2007, at 3:52 PM, Richard Gomes wrote:
>> I compiled/installed the latest boost release.
>> I just checked out QuantLib from SVN.
>> I installed autoconf, automake and m4 packages.
>>
>> and then...
>>
>> $ sh autogen.sh
>> /usr/share/aclocal/libmcrypt.m4:17: warning: underquoted definition of
>> AM_PATH_LIBMCRYPT
>> /usr/share/aclocal/libmcrypt.m4:17:   run info '(automake)Extending
>> aclocal'
>> /usr/share/aclocal/libmcrypt.m4:17:   or see
>> http://sources.redhat.com/automake/automake.html#Extending-aclocal
>> configure.ac:63: error: possibly undefined macro: AC_PROG_LIBTOOL
>>       If this token and others are legitimate, please use
>>       m4_pattern_allow. See the Autoconf documentation.
>> autoreconf: /usr/bin/autoconf failed with exit status: 1

Did you install libtool as well as autoconf and automake? it seems that 
autoconf does not find an m4 macro file provided with libtool.


> This is a followup to the group.
>
> Version 0.8.0 has the c++ sources and swig sources in download area at
> sourceforge.net.

Also, it's probably not very clear from the release notes, but 
QuantLib-SWIG 0.8.0 will work with QuantLib 0.8.1.

Later,
	Luigi

Re: [Quantlib-dev] Help! newbie question : autogen

[Richard G. <rgo...@ya...>]

Richard Gomes wrote:

> Hi guys,
> 
> I'd like to call QuantLib from Java. I'm following the steps from
> http://www.skwash.com and I've got errors, probably because the .i files
> do not match the 0.8.1 release.
> 
> As Debian does not have a kind of quantlib-swig package, I decided to
> compile from sources, as described:
> 
> I compiled/installed the latest boost release.
> I just checked out QuantLib from SVN.
> I installed autoconf, automake and m4 packages.
> 
> and then...
> 
> $ sh autogen.sh
> /usr/share/aclocal/libmcrypt.m4:17: warning: underquoted definition of
> AM_PATH_LIBMCRYPT
> /usr/share/aclocal/libmcrypt.m4:17:   run info '(automake)Extending
> aclocal'
> /usr/share/aclocal/libmcrypt.m4:17:   or see
> http://sources.redhat.com/automake/automake.html#Extending-aclocal
> configure.ac:63: error: possibly undefined macro: AC_PROG_LIBTOOL
>       If this token and others are legitimate, please use
>       m4_pattern_allow. See the Autoconf documentation.
> autoreconf: /usr/bin/autoconf failed with exit status: 1
> 
> Oh well... I'm a complete newbie on autoconf/automake. I followed ...
> 
> http://sources.redhat.com/automake/automake.html#Extending-aclocal
> 
> ... but it seems to be a steep learning curve.
> 
> Does anyone have an idea how to solve the error I've got?
> 
> Thank you all in advance for any help.
> 

This is a followup to the group.

Version 0.8.0 has the c++ sources and swig sources in download area at
sourceforge.net.

I'm trying it now.

-- 
Richard Gomes
mailto:rg...@ya...
M:+44(77)9955-6813
H:+44(870)068-8205

[Quantlib-dev] Help! newbie question : autogen

[Richard G. <rgo...@ya...>]

Hi guys,

I'd like to call QuantLib from Java. I'm following the steps from
http://www.skwash.com and I've got errors, probably because the .i files do
not match the 0.8.1 release.

As Debian does not have a kind of quantlib-swig package, I decided to
compile from sources, as described:

I compiled/installed the latest boost release.
I just checked out QuantLib from SVN.
I installed autoconf, automake and m4 packages.

and then...

$ sh autogen.sh
/usr/share/aclocal/libmcrypt.m4:17: warning: underquoted definition of
AM_PATH_LIBMCRYPT
/usr/share/aclocal/libmcrypt.m4:17:   run info '(automake)Extending aclocal'
/usr/share/aclocal/libmcrypt.m4:17:   or see
http://sources.redhat.com/automake/automake.html#Extending-aclocal
configure.ac:63: error: possibly undefined macro: AC_PROG_LIBTOOL
      If this token and others are legitimate, please use m4_pattern_allow.
      See the Autoconf documentation.
autoreconf: /usr/bin/autoconf failed with exit status: 1

Oh well... I'm a complete newbie on autoconf/automake. I followed ... 

http://sources.redhat.com/automake/automake.html#Extending-aclocal

... but it seems to be a steep learning curve.

Does anyone have an idea how to solve the error I've got?

Thank you all in advance for any help.

-- 
Richard Gomes
mailto:rg...@ya...
M:+44(77)9955-6813
H:+44(870)068-8205

[Quantlib-dev] multidimensional integration

[Warren C. <war...@al...>]

Does Quantlib have any multidimensional integration routines?

If not, do the 1-D integration routines lend themselves to multidimensional integration?

Thanks

[Quantlib-dev] Projecting Forward Volatility Surfaces

[David B. <doc...@gm...>]

Greetings All-

   Given a portfolio of american options on one stock, let's suppose that we
wanted to figure out what this portfolio was worth in an arbitrary number of
days, suppose 30.  One could take the stock, let's suppose that we could
project the pdf of the price forward 30 days perfectly with no loss of fat
tails.  (This is another problem in itself.)  No we have this accurate pdf,
let's say composed of 50000 MCMC simulations.  Every one of these values
could be plugged into the options we hold and now and semi-accurate values
could be obtained using whatever we priced the option with in the first
place, baroni-adesi-whaley / whatever.  The only matter now is that we don't
have the slightest idea of what the volatility would be in this arbitrary
number of days.
   How do we begin to solve this problem?  My approach given my own thoughts
and a number of papers, including Peter Carr and Jim Gatheral, is to find
the market implied volatility surface at time zero and figure out using
historical data what the surface should approximately look like in the
future.  This is way easier said than done.  Dupire has a slide-show that
discusses the different types of local volatility and stochastic volatility
surfaces, but his book on volatity only has a short chapter, number 8,
called Dynamics of the Skew / Smile or something to that effect.  He doesn't
really seem to answer the question as to what must be done to solve this
problem accurately.
   So I've raised this question here, because I think this is the only place
where there is a possibility of this problem being solved given the
expertise of the parties involved in similar projects.  Peter Carr addressed
this a little bit with his FMLS model which uses stable distributions to fit
the market smile, but the model has yet to be built for stochastic
volatility.  I contacted Johnathon Nolan on the matter and he told me a good
source to try would be Hugh McCullough's papers.  I tried his papers and it
appears that he hasn't fully solved the matter either.  Peter Carr himself
said that some type of PDE could be used to solve for the speed of the
volatility surface, but from what I gather, this was his form of
conjecturing.  Perhaps not?
   Any advice on how to proceed with this problem?  I am highly interested
in programming this into Quantlib.


Best Wishes,
David Brown

[Quantlib-dev] Binomial Tree Fever

[Mark j. <mar...@gm...>]

Re the trees being buggy, I think the unit test would have picked up
an error that big.

Re trees, in general, I've been doing a lot of playing with them
lately. My main conclusion is that it's not what tree you use for
American puts, but what acceleration techniques you use: i.e.
Richardson extrapolation, smoothing, truncation, that matters.

(expect a paper soon!)

Mark

Re: [Quantlib-dev] Binomial Tree Fever

[John M. <jwm...@ya...>]

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
 <qua...@li...>. 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
 <qua...@li...>. 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);
	}

}