Question on PathGenerator

[Penschke, Walter]

Hi,

I am currently investigating quantlib's MonteCarlo package.

I tried to simulate a geometric brownian motion process, which I did define
in my own class GeometricBrownianMotionProcess. Basically this class is
derived from StochasticProcess and defines the methods drift and diffusion.
Here my source code:

 <<GeometricBrownianMotionProcess.cpp>>  
<<GeometricBrownianMotionProcess.hpp>>

The stochastic differential equation of this process is:

        dS = mue * S * dt + sigma * S * dz


The problem now is that class PathGenerator from Quantlib (defined in
pathgenerator.hpp) does not produce proper paths for this process. Having
had a look at the next() method of PathGenerator (brownianBridge_ == false)
I realised that the drift and the diffusion arrays in the returned Path are
computed correctly in principal. However the asset_ variable, which tracks
the evolvement of the process, in my case S, within this method is adjusted
as follows per step:

        asset_ *= QL_EXP(next_.value[i])

This is correct if the Path variable (next_) is modelling a process defined
on dS/S, i.e.: on the change of the process, rather than the process itself.
For my GeometricBrownianProcess, which defines the process on dS, the line
above in my (simple) opinion should be (May be that's not quite correct, but
I hope you get the idea.):

        asset_ *= QL_EXP(next_.value[i] / asset_)

As I am modelling S in my path, next_.value[i] does not contain a percentage
change, but the new simulated value itself.

This would mean that I can't reuse the existing PathGenerator class, but I
need to implement my own one which basically only differs in the line above.
I would not fear this extra effort. However, I think this code duplication
is not nice and should be avoided if possible.

Are there any comments on this? Please let me know if you need more
information on this issue.

Any help would be highly appreciated.


wpe

[Penschke, Walter]

I'd like to add the formulas for the drift and diffusion implemented in
GeometricBrownianMotionProcess, because the source files cannot be retrieved
from the mailing list archive:

drift(t,x) = mue * x
diffusion(t,x) = sigma * x

Calls to the methods expectation() and variance()
are delegated to the sub-class StochasticProcess.

wpe

-----Original Message-----
From: [hidden email]
[mailto:[hidden email]]On Behalf Of
Penschke, Walter
Sent: Wednesday, September 01, 2004 5:57 PM
To: [hidden email]
Subject: [Quantlib-users] Question on PathGenerator


Hi,

I am currently investigating quantlib's MonteCarlo package.

I tried to simulate a geometric brownian motion process, which I did define
in my own class GeometricBrownianMotionProcess. Basically this class is
derived from StochasticProcess and defines the methods drift and diffusion.
Here my source code:

 <<GeometricBrownianMotionProcess.cpp>>  
<<GeometricBrownianMotionProcess.hpp>>

The stochastic differential equation of this process is:

        dS = mue * S * dt + sigma * S * dz


The problem now is that class PathGenerator from Quantlib (defined in
pathgenerator.hpp) does not produce proper paths for this process. Having
had a look at the next() method of PathGenerator (brownianBridge_ == false)
I realised that the drift and the diffusion arrays in the returned Path are
computed correctly in principal. However the asset_ variable, which tracks
the evolvement of the process, in my case S, within this method is adjusted
as follows per step:

        asset_ *= QL_EXP(next_.value[i])

This is correct if the Path variable (next_) is modelling a process defined
on dS/S, i.e.: on the change of the process, rather than the process itself.
For my GeometricBrownianProcess, which defines the process on dS, the line
above in my (simple) opinion should be (May be that's not quite correct, but
I hope you get the idea.):

        asset_ *= QL_EXP(next_.value[i] / asset_)

As I am modelling S in my path, next_.value[i] does not contain a percentage
change, but the new simulated value itself.

This would mean that I can't reuse the existing PathGenerator class, but I
need to implement my own one which basically only differs in the line above.
I would not fear this extra effort. However, I think this code duplication
is not nice and should be avoided if possible.

Are there any comments on this? Please let me know if you need more
information on this issue.

Any help would be highly appreciated.


wpe

[Luigi Ballabio-2]

Hello Walter,

On 2004.09.01 17:56, "Penschke, Walter" wrote:

> I tried to simulate a geometric brownian motion process, which I did
> define in my own class GeometricBrownianMotionProcess. Basically this  
> class is derived from StochasticProcess and defines the methods drift  
> and diffusion.
>
> The stochastic differential equation of this process is:
>
> dS = mue * S * dt + sigma * S * dz
>
> The problem now is that class PathGenerator from Quantlib (defined in
> pathgenerator.hpp) does not produce proper paths for this process.
> [...] the asset_ variable, which tracks the evolvement of the  
> process, in my case S, within this method is adjusted as follows per  
> step:
>
> asset_ *= QL_EXP(next_.value[i])
>
> This is correct if the Path variable (next_) is modelling a process
> defined on dS/S, i.e.: on the change of the process, rather than the  
> process itself.

Yes. I've been aware of this for some time, but I still haven't found  
the time to try and remove the assumption. Right now, PathGenerator can  
only be used with a Black-Scholes process, which is hardly a happy  
state of affairs--as there are in the library a couple of other  
processes which suffer the same problem you found.

> For my GeometricBrownianProcess, which defines the process on dS, the
> line above in my (simple) opinion should be (May be that's not quite
> correct, but I hope you get the idea.):
>
> asset_ *= QL_EXP(next_.value[i] / asset_)

or more simply,

        asset_ += next_.value[i];

> This would mean that I can't reuse the existing PathGenerator class,
> but I need to implement my own one which basically only differs in  
> the line above. I think this code duplication is not nice and should  
> be avoided if possible.

I agree. A new method should be added to StochasticProcess, so that the  
above can be written as:

        asset_ = process_->apply(asset_, next_.value[i])

where the process would take care of updating the asset value in the  
proper way.


Then again, an ambiguity would remain---as of today, the generated path  
stores asset variations, and not asset values.  This implies that given  
a path, it is not possible to use it properly without knowing whether  
the variations are absolute or relative, i.e., without knowing which  
process was used to generate it.  The Monte Carlo pricers in the  
library work under the implicit assumption that the paths they're given  
were generated by using a Black-Scholes process--they would fail if a  
process were used which modeled dS instead of dS/S.

It would make me much happier if the actual asset values along the path  
were stored instead.

Thoughts?

Later,
        Luigi

[Neil P Firth]

Hi,

I have also been having a look at this lately. As you say, it will
become clearer what to do when we can simulate other processes. I have
implemented gamma and poisson random number generators (from MC Methods in
Financial Engineering - Glasserman), with the aim of simulating Jump
Diffusion and Variance Gamma processes. Hopefully it should be possible to
use the European path pricer, etc, to value Euro options on those
processes. Then we can have a PathGenerator base class and
GBMPathGenerator, Merton76PathGenerator, and VarianceGammaPathGenerator
sub classes.

Storing the drift and diffusion values at each step in the path is going
to have to change for jump processes. Storing the drift and diffusion
values separately makes it easy to do antithetic variates, so that will
need some thought.

Thoughts? Any volunteers to write general test cases for random variates
generated from different distributions (normal, poisson, gamma)?

Neil

---------------------------------------------------
  Neil Firth
  Brasenose College Oxford OX1 4AJ United Kingdom
  Office: 01865 280616
  [hidden email]
  http://www.maths.ox.ac.uk/~firth
---------------------------------------------------

[Penschke, Walter]

Hi Luigi,

I would say we should refactor this a little bit. I think your idea with
creating a new evolve() method, which is part of a StochasticProcess is
valid.

As you lined out, we could then call this process-method from within the
next()
method of the PathGenerator. To me such an evolve() method would be a
natural
thing for a process. The parameter would be path[i].

Who decides wether to do this refactoring?
How could I possibly contribute to this extension?


wpe

-----Original Message-----
From: Luigi Ballabio [mailto:[hidden email]]
Sent: Thursday, September 02, 2004 11:22 AM
To: Penschke, Walter
Cc: [hidden email]
Subject: Re: [Quantlib-users] Question on PathGenerator



Hello Walter,

On 2004.09.01 17:56, "Penschke, Walter" wrote:

> I tried to simulate a geometric brownian motion process, which I did
> define in my own class GeometricBrownianMotionProcess. Basically this  
> class is derived from StochasticProcess and defines the methods drift  
> and diffusion.
>
> The stochastic differential equation of this process is:
>
> dS = mue * S * dt + sigma * S * dz
>
> The problem now is that class PathGenerator from Quantlib (defined in
> pathgenerator.hpp) does not produce proper paths for this process.
> [...] the asset_ variable, which tracks the evolvement of the  
> process, in my case S, within this method is adjusted as follows per  
> step:
>
> asset_ *= QL_EXP(next_.value[i])
>
> This is correct if the Path variable (next_) is modelling a process
> defined on dS/S, i.e.: on the change of the process, rather than the  
> process itself.

Yes. I've been aware of this for some time, but I still haven't found  
the time to try and remove the assumption. Right now, PathGenerator can  
only be used with a Black-Scholes process, which is hardly a happy  
state of affairs--as there are in the library a couple of other  
processes which suffer the same problem you found.

> For my GeometricBrownianProcess, which defines the process on dS, the
> line above in my (simple) opinion should be (May be that's not quite
> correct, but I hope you get the idea.):
>
> asset_ *= QL_EXP(next_.value[i] / asset_)

or more simply,

        asset_ += next_.value[i];

> This would mean that I can't reuse the existing PathGenerator class,
> but I need to implement my own one which basically only differs in  
> the line above. I think this code duplication is not nice and should  
> be avoided if possible.

I agree. A new method should be added to StochasticProcess, so that the  
above can be written as:

        asset_ = process_->apply(asset_, next_.value[i])

where the process would take care of updating the asset value in the  
proper way.


Then again, an ambiguity would remain---as of today, the generated path  
stores asset variations, and not asset values.  This implies that given  
a path, it is not possible to use it properly without knowing whether  
the variations are absolute or relative, i.e., without knowing which  
process was used to generate it.  The Monte Carlo pricers in the  
library work under the implicit assumption that the paths they're given  
were generated by using a Black-Scholes process--they would fail if a  
process were used which modeled dS instead of dS/S.

It would make me much happier if the actual asset values along the path  
were stored instead.

Thoughts?

Later,
        Luigi

[Ferdinando M. Ametrano-3]

Hi Walter

>Who decides wether to do this refactoring?
>How could I possibly contribute to this extension?

the best way is to contribute a patch against the current CVS version of
QuantLib.

Check out QuantLib anonymously (see http://www.quantlib.org/cvs.shtml),
refactor the code, make a diff and send it to the list, or even better
submit it using out Patch Tracker
(http://sourceforge.net/tracker/?group_id=12740&atid=312740)

thank you

ciao -- Nando

[Luigi Ballabio-2]

Hi Walter,

On 2004.09.02 17:15, "Penschke, Walter" wrote:
> Who decides wether to do this refactoring?

Nando and I.  In this case, I'm for refactoring so he better stop me  
quick if he disagrees :)

> How could I possibly contribute to this extension?

As Nando pointed up, send me a patch---but before starting to code, you  
might want to email me privately so that we settle on the interface  
beforehand.

With regard to the second part of the refactoring, i.e., storing the  
asset values in the path---it might be harder to complete since it  
would affect a number of instruments and pricing engines and it should  
be done backward-compatibly. Of course it's ok if you want to tackle  
the full thing---but otherwise, just modify the Path class, send me the  
patch, and I'll retrofit it to the library.

Later,
        Luigi

[Luigi Ballabio-2]

On 2004.09.02 13:19, Neil P Firth wrote:
> Storing the drift and diffusion values separately makes it easy to do  
> antithetic variates

Hmm, not really.  If drift and diffusion are dependent on the value of  
the asset, it doesn't hold that the antithetic of (drift+diffusion) is  
(drift-diffusion)---or rather, it just holds for the first step. At the  
second step, the same random number will give different diffusion  
values as the asset value will be different for the path and its  
antithetic.

Later,
        Luigi

[Ferdinando M. Ametrano-3]

Hi Luigi

At 07:02 PM 9/2/2004, Luigi Ballabio wrote:
>On 2004.09.02 13:19, Neil P Firth wrote:
>>Storing the drift and diffusion values separately makes it easy to do
>>antithetic variates
>
>Hmm, not really.  If drift and diffusion are dependent on the value of
>the asset, it doesn't hold that the antithetic of (drift+diffusion) is
>(drift-diffusion)

yes, you're right. Anyway it is so common to have drift and diffusion not
dependent on the value of the asset: the point about simulating dS/S = mu
dt + sig dZ instead of dS = mu S dt + sig S dZ is indeed to have drift and
diffusion not dependent on S.

Of course this approach fails with local vol and other models, still it
might be worthwhile to allow for storing drift and diffusion separately.

ciao -- Nando

[Ferdinando M. Ametrano-3]

Hi all

Luigi:
>>Who decides wether to do this refactoring?
>
>Nando and I.  In this case, I'm for refactoring so he better stop me
>quick if he disagrees :)
I agree with the refactoring.

>With regard to the second part of the refactoring, i.e., storing the
>asset values in the path---it might be harder to complete since it
>would affect a number of instruments and pricing engines and it should
>be done backward-compatibly.

it's not just backward-compatibility. Calculating the asset values even
when you don't need them it is a waste of exp() computations which is going
to slow down the simulation.

We need more flexibility about stimulating S or log S, about (not) storing
drift and diffusion separately, etc.

Any hint by the MC gurus out there would be appreciated :)

ciao -- Nando

[Ferdinando M. Ametrano-3]

I wrote:
>We need more flexibility about stimulating S or log S

of course just don't worry about sTimulating S... :)

ciao -- Nando

[Luigi Ballabio-2]

On 2004.09.03 10:11, Ferdinando Ametrano wrote:
> it's not just backward-compatibility. Calculating the asset values  
> even when you don't need them it is a waste of exp() computations  
> which is going to slow down the simulation.

Right, but we're calculating them already---right now, storing them  
would not slow down the calculation. I agree that it could prevent  
making it faster, which is an issue we might want to explore :)

My proposal:

 - have the StochasticProcess know whether it's working on S or log(S);
  this is done by implementing its interface correctly, where
  'correctly' is defined as:

  - when the process models S, drift(t,s), diffusion(t,s) and such
    mean the obvious thing;

  - when the process models log(S), drift(t,x) means drift(t,log(s));
    same for diffusion, variance and such;

- at this point, the stochastic process interface is transparent with
  respect to its proper underlying and the path generator just works
  by doing:

  x = ...; // (1) initial value
  for (...) {
      drift = dt*proc->drift(t,x);
      diffusion = rnd*sqrt(proc->variance(t,x,dt);
      x += (drift+diffusion);
      ... // (2) store the value?
  }

  where the exp have disappeared (of course, if the Black-Scholes
  process uses a local volatility surface, it will do the
  exponentiation inside drift() and/or variance(), so it doesn't gain
  anything--but other processes might.)  We still have to specify in
  (1) and (2).

- at (1) we have the problem of giving the above a proper initial
  value x0, to be deduced from the initial asset value; this is done
  by calling proc->underlying(s0), which is the identity function when
  S is modeled and the log function when log(S) is;

- at (2) we'd like to store asset values, but we don't want to
  calculate them needlessly. TimeGrid comes to the rescue: it has the
  possibility to specify that some time points are mandatory stops,
  while others are just intermediate steps.  At the mandatory times,
  we'll need the asset value anyway (they're fixing dates, exercise
  dates...) so that we don't save anything by not exponentiating--it
  would only be done later.  The intermediate points are only used
  for building a more accurate path when drift and diffusion are
  asset-dependent, but the asset value at those points is not needed.

  Therefore, we calculate and store in the path the asset values at
  mandatory times. Actually, I would go a bit further: if we pass the
  generator a time grid with 100 points and 10 mandatory times, it
  would return a path with 10 points, each of them being the
  cumulative result of the intermediate ones.  The values at mandatory
  times would be calculated by calling proc->observable(x), which is
  the identity function when S is modeled and the exp function when
  log(S) is (no, I didn't call it proc->assetValue---although I've
  been using 'asset value' until now for clarity, we might be modeling
  the short rate by using the dynamics returned by a calibrated
  Hull-White model.)

Thoughts?

Later,
        Luigi

[ML-21]

> - at (2) we'd like to store asset values, but we don't want to
>   calculate them needlessly. TimeGrid comes to the rescue: it has the
>   possibility to specify that some time points are mandatory stops,
>   while others are just intermediate steps.

As a newbie to QuantLib and hence at the risk of making an unreasonable
suggestion: what about defining a new "Epoch" class, which is basically the
same as the "Time" class, but augmented by a bitfield "Attributes"? The
individual bits in this bitfield would specify things like "mandatory time",
"save asset price", etc. In "TimeGrid" the std::vector<Time> would have to
be replaced by std::vector<Epoch>. Sorry if my idea is not really helpful.

-Mario

[Luigi Ballabio-2]

On 2004.09.03 15:52, ML wrote:
> As a newbie to QuantLib and hence at the risk of making an
> unreasonable suggestion: what about defining a new "Epoch" class,  
> which is basically the same as the "Time" class, but augmented by a  
> bitfield "Attributes"?

Mario,
        it wouldn't be unreasonable---the only thing is, Time is not a  
class. It's a typedef to double, and I'd keep it that way for the time  
being...

Later,
        Luigi

[Penschke, Walter]

Ok, will do.

wpe

-----Original Message-----
From: Ferdinando Ametrano [mailto:[hidden email]]
Sent: Thursday, September 02, 2004 5:32 PM
To: Penschke, Walter
Cc: [hidden email]; '[hidden email]'
Subject: RE: [Quantlib-users] Question on PathGenerator


Hi Walter

>Who decides wether to do this refactoring?
>How could I possibly contribute to this extension?

the best way is to contribute a patch against the current CVS version of
QuantLib.

Check out QuantLib anonymously (see http://www.quantlib.org/cvs.shtml),
refactor the code, make a diff and send it to the list, or even better
submit it using out Patch Tracker
(http://sourceforge.net/tracker/?group_id=12740&atid=312740)

thank you

ciao -- Nando