FdDividendAmericanOption

[Graham Miller-7]

So the code at the bottom of this email (which I think maybe someone on
this list pointed me to) produces the following output:

assetSteps=20 price=2.16967 delta=0.51742 vega=7.72539
assetSteps=40 price=2.16968 delta=0.517419 vega=7.72561
assetSteps=60 price=4.75005 delta=-4.39606 vega=131.91

Basically the price, delta and vega for the same option with three
different values for the assetSteps parameter.  The first two sets of
values are reasonable.  The third set (the one with price=4.75005) is
definitely way off.  The only thing I can see (given that I'm new to
QuantLib) is that this is the only one where assetSteps > timeSteps.  
Does anyone know if this should be a problem?

graham

// and the code...

#include <ql/quantlib.hpp>

using namespace std;
using namespace QuantLib;

int main() {
  const unsigned timeSteps = 40;
  vector<unsigned> assetStepsA;
  assetStepsA.push_back(20);
  assetStepsA.push_back(40);
  assetStepsA.push_back(60);

  vector<double> dividendA, dividendTimeA;
  dividendA.push_back(0.025);
  dividendTimeA.push_back(1e-4);

  double under=54.625;
  double strike=55;
  double dividendYield=0;
  double interestRate=0.052706;
  double expirationTime=0.126027;
  double sigma=0.282922;


  for (unsigned i=0; i<assetStepsA.size(); i++) {
    const unsigned assetSteps = assetStepsA[i];

    FdDividendAmericanOption opt(Option::Call, under,
                                 strike,
                                 dividendYield, interestRate,
                                 expirationTime, sigma,
                                 dividendA, dividendTimeA,
                                 timeSteps, assetSteps);


    cout << "assetSteps=" << assetSteps
         << " price=" << opt.value()
         << " delta=" << opt.delta()
         << " vega=" << opt.vega()
         << "\n";
  }

  return 0;
}

[Joseph Wang]

Just looking over the code, what seems to be happening is that with a
large grid
and a short time step, the code puts some grid points at fairly extreme
values,
which result in large dividends.

If I run the code with zero dividend, I still get odd values.  I don't
get any odd
values if I just run FdAmericanOption.

[Joseph Wang]

Still looking but I think I know what the basic trouble is.  Basically there
is a numerical instability if you make the time steps too large.  I
think that the
stepper is missing the Courant condition.

One way of thinking about this is to imagine the simulation as sound waves.
If the time step is longer than the time it takes for a wave to "move"
across
the grid, then the grid elements start getting disconnected from each other.
The consequence is that if you make the grid steps small, you have to also
make the time steps small, so that the grid steps don't get disconnected
from each other.

There are a number of ways to fix the problem.  The first is to just refuse
to do the simulation when the Courant condition fails.  The second is
involves changing the way that the time steps get calculated.

What confuses me is that the Courant condition normally doesn't affect
financial equations since financial equations usually involve "smoothing"
rather than "waves."  But I suspect that the interest rate which involves
"drifting" rather than "smoothing" is causing the finite differencer to
break.

[Joseph Wang]

More data:

* Bermudan options also seems to fail when you run it with small grids.
* The failure appears independent of dividends (i.e. it still fails when
you get rid of dividends).
* The failure appears to be independent of finite differencing models
* The numerical instability doesn't seem to happen when you have values
that are very much
   in the money or out of the money
* As I mentioned before, the code fails when interest rates are being
set to zero which suggests
   that the problem involves calculating the second derivative

The one thing that seems to pop out is the BSM operator which takes as
input the log of the gridSpacing_.  Knowing nothing about the code, it would
seem to me that the BSMOperator is expecting a linear grid, and putting in
a logarithmic grid would appear to do bad things.  It looks like I'm
going to
have to crack open Hull and White and check the derivation from first
principles.

Can anyone who knows more about the code, let me know if I'm on the right
track?  This does seem to be a pretty major bug, unless I'm missing
something.

[Luigi Ballabio]

On 12/30/04 09:06:01, Joseph Wang wrote:
>
> The one thing that seems to pop out is the BSM operator which takes as
> input the log of the gridSpacing_.  Knowing nothing about the code, it  
> would seem to me that the BSMOperator is expecting a linear grid, and  
> putting in a logarithmic grid would appear to do bad things.

Joseph,
        the BSM operator discretizes the equation for x = log(S), therefore  
a logarithmic grid in S is used which corresponds to a linear grid in x.

Also, I very much appreciate that you dived in to get data, as  
unfortunately I don't have the time to dive in myself now---thanks a lot.

If I can suggest a line of investigation, I find interesting that, as you  
said,

> If I run the code with zero dividend, I still get odd values.  I don't  
> get any odd values if I just run FdAmericanOption.

You might get some insight by comparing the two calculations. Are the two  
grids the same when the two classes are passed the same arguments? Are the  
calculated values the same until the first dividend payment is reached?  
What happens to the asset values when the code for applying dividends is  
triggered, even with zero dividends?

Later,
        Luigi

[Joseph Wang]

I'm suspecting that rescaling the terms causes the operator to break.  
Attached is a proposed
fix that calculates BSMOperator from a flexible grid, that seems to
cause the numbers to
be more reasonable......

assetSteps=20 price=2.17038 delta=0.518664 vega=7.74307
assetSteps=40 price=2.16897 delta=0.516824 vega=7.71849
assetSteps=60 price=2.21022 delta=0.534935 vega=7.28234
assetSteps=70 price=2.22188 delta=0.590997 vega=8.0055
assetSteps=80 price=2.20972 delta=0.580363 vega=7.96264
assetSteps=100 price=2.19306 delta=0.560017 vega=7.90149




Index: ql/FiniteDifferences/bsmoperator.cpp
===================================================================
RCS file: /cvsroot/quantlib/QuantLib/ql/FiniteDifferences/bsmoperator.cpp,v
retrieving revision 1.16
diff -u -r1.16 bsmoperator.cpp
--- ql/FiniteDifferences/bsmoperator.cpp 12 May 2004 09:46:11 -0000 1.16
+++ ql/FiniteDifferences/bsmoperator.cpp 31 Dec 2004 06:20:46 -0000
@@ -16,6 +16,7 @@
 */
 
 #include <ql/FiniteDifferences/bsmoperator.hpp>
+#include <ql/Math/array.hpp>
 
 namespace QuantLib {
 
@@ -30,4 +31,22 @@
         setMidRows(pd,pm,pu);
     }
 
+    BSMOperator::BSMOperator(const Array &grid, Rate r,
+                             Rate q, Volatility sigma)
+    : TridiagonalOperator(grid.size()) {
+        Real sigma2 = sigma*sigma;
+        Real nu = r-q-sigma2/2;
+ Size grid_size = grid.size();
+ for (Size i=1; i < grid_size - 1; i++) {
+  Real dxm = grid[i] - grid[i-1];
+  Real dxp = grid[i+1] - grid[i];
+  Real dxtotal = dxm+dxp;
+
+  Real pd = -(sigma2/dxm-nu)/dxtotal;
+  Real pu = -(sigma2/dxp+nu)/dxtotal;
+  Real pm = sigma2/(dxm*dxp)+r;
+  setMidRow(i, pd,pm,pu);
+ }
+    }
+
 }
Index: ql/FiniteDifferences/bsmoperator.hpp
===================================================================
RCS file: /cvsroot/quantlib/QuantLib/ql/FiniteDifferences/bsmoperator.hpp,v
retrieving revision 1.16
diff -u -r1.16 bsmoperator.hpp
--- ql/FiniteDifferences/bsmoperator.hpp 12 May 2004 09:46:11 -0000 1.16
+++ ql/FiniteDifferences/bsmoperator.hpp 31 Dec 2004 06:20:47 -0000
@@ -24,6 +24,8 @@
 
 #include <ql/FiniteDifferences/tridiagonaloperator.hpp>
 
+class Array;
+
 namespace QuantLib {
 
     //! Black-Scholes-Merton differential operator
@@ -32,6 +34,7 @@
       public:
         BSMOperator() {}
         BSMOperator(Size size, Real dx, Rate r, Rate q, Volatility sigma);
+      BSMOperator(const Array &grid, Rate r, Rate q, Volatility sigma);
     };
 
 }
Index: ql/Pricers/fdbsmoption.cpp
===================================================================
RCS file: /cvsroot/quantlib/QuantLib/ql/Pricers/fdbsmoption.cpp,v
retrieving revision 1.21
diff -u -r1.21 fdbsmoption.cpp
--- ql/Pricers/fdbsmoption.cpp 13 May 2004 11:41:46 -0000 1.21
+++ ql/Pricers/fdbsmoption.cpp 31 Dec 2004 06:20:47 -0000
@@ -93,7 +93,7 @@
     }
 
     void FdBsmOption::initializeOperator() const {
-        finiteDifferenceOperator_ = BSMOperator(gridPoints_, gridLogSpacing_,
+        finiteDifferenceOperator_ = BSMOperator(grid_,
                                                 riskFreeRate_, dividendYield_,
                                                 volatility_);
 

[Joseph Wang]

I think the basic problem is that when you turn x = log(S) into an operator,
dx = dS / S and you start getting extra terms in the difference equation
that
aren't being put into BSMOperator.  One could go and work out what the
operator should be with a logarithmic grid, but I suspect that the second
derivative is going to add a lot of nasty terms.

[Graham Miller-7]

Thanks, Joseph for your patch.

I applied it to my code base. Then I thought I'd look at price
convergence versus the number of assetSteps--just for fun--and I came
across an interesting characteristic.  I modified the code I posted
earlier to use every value of assetSteps between 1 and 300.  I attached
a graph of the results.  It's a little hard to tell from the graph, but
the values with an odd number of assetSteps quickly converge to about
2.17. The values with an even number of asset steps diverge a fair
amount before also converging on about 2.17.

Any thoughts?

graham

Joseph Wang wrote:

> I think the basic problem is that when you turn x = log(S) into an
> operator,
> dx = dS / S and you start getting extra terms in the difference
> equation that
> aren't being put into BSMOperator.  One could go and work out what the
> operator should be with a logarithmic grid, but I suspect that the second
> derivative is going to add a lot of nasty terms.
>
>
>

[Luigi Ballabio]

On 12/31/04 08:16:31, Joseph Wang wrote:
> I think the basic problem is that when you turn x = log(S) into an  
> operator, dx = dS / S and you start getting extra terms in the difference  
> equation  that aren't being put into BSMOperator.  One could go and work  
> out what the operator should be with a logarithmic grid, but I suspect  
> that the second derivative is going to add a lot of nasty terms.

Joseph,
        no, the BSM operator is already the correct discretization for the  
logarithm of S. It discretizes in space the equation

df/dt = - sigma^2/2 * d^2f/dx^2 - nu * df/dx + r*f

that one gets from changing the variable from S to x = log S in the BS  
equation. At the point of discretization, it is only incidental that x  
comes from a logarithm; it is simply the variable being discretized, and  
building a "logarithmic" grid is just a way (albeit a roundabout one) to  
build an equally spaced grid for x. In fact, the name "grid_" in the code  
is a misnomer; it is not the grid { x_1, ..., x_N } used for the operator,
but rather the set { S_1, ..., S_N } of prices corresponding to the grid  
points. The actual grid is { log(S) for S in grid_ }, and if you print it  
out, you'll see that it is correctly spaced.  In short, it is beyond me how  
your patch provides results that look sensible :)

I had naively hoped to give the problem some time during the weekend, but I  
found that having a wife and a couple of kids seriously cuts into one's  
debugging time :)
All I was able to put together was the small patch below, which does not  
completely remove the instability, but at least triggers it only above 250  
asset steps for the example code given.

diff -r1.28 fddividendoption.cpp
74c74
<         setGridLimits(center_ + dividends_[step], dates_[step]);
---
>         setGridLimits(center_ + dividends_[step], residualTime_);

Basically, the code tried to be smart (probably too much for its own good)  
and resize the grid based on the time to the dividend payment. Now it  
resizes it based on the time to maturity, so that it doesn't end up with a  
grid too different from the one it started with.

But the sad truth is, the whole thing should be rewritten. I had a look at  
the code in the past few days (I had'n written it myself) and it's a tangle  
of methods, each one trying to outsmart the user and correct the number of  
grid points, the grid limits and whatever else based on what was thought to  
be a "safe" choice of parameters. Not that it can't be useful, but the user  
should at least have the choice to say "heck no, just do as I said."

Needless to say, at this time I wouldn't use the FD dividend options where  
actual money is involved---at least, not without cross-checking the results  
with some other pricer or some other source.

Later,
        Luigi

[Joseph Wang]

Yes, you are correct.  The patch that I gave seems stable, but it solves
the wrong problem :-(
It took me a bit of algebra to convince myself of this.

There are two possibilities that I see:

1) that there is a really simple error (signs maybe) - I'm starting to
doubt this more and more, because I don't think that a simple error
would have escaped all of the eyes that are staring at this.  The other
issue is that isn't not obvious how a sign error would cause the
equation to be stable at low grid sizes.

2) a much nastier problem is the possibility that there is some subtle
PDE interaction going on.  Diffusion equations generally don't have
these problems, but the discretization adds in an advection term, and
this may be causing things to break.  If this is the problem then one
would expect the differencing to become unstable at delta t > nu * delta
x.  Let me try messing with the difference equations to see what I can do.

Luigi Ballabio wrote:

> On 12/31/04 08:16:31, Joseph Wang wrote:
>
>> I think the basic problem is that when you turn x = log(S) into an  
>> operator, dx = dS / S and you start getting extra terms in the
>> difference  equation  that aren't being put into BSMOperator.  One
>> could go and work  out what the operator should be with a logarithmic
>> grid, but I suspect  that the second derivative is going to add a lot
>> of nasty terms.
>
>
> Joseph,
>     no, the BSM operator is already the correct discretization for
> the  logarithm of S. It discretizes in space the equation
>
> df/dt = - sigma^2/2 * d^2f/dx^2 - nu * df/dx + r*f
>
> that one gets from changing the variable from S to x = log S in the
> BS  equation. At the point of discretization, it is only incidental
> that x  comes from a logarithm; it is simply the variable being
> discretized, and  building a "logarithmic" grid is just a way (albeit
> a roundabout one) to  build an equally spaced grid for x. In fact, the
> name "grid_" in the code  is a misnomer; it is not the grid { x_1,
> ..., x_N } used for the operator,
> but rather the set { S_1, ..., S_N } of prices corresponding to the
> grid  points. The actual grid is { log(S) for S in grid_ }, and if you
> print it  out, you'll see that it is correctly spaced.  In short, it
> is beyond me how  your patch provides results that look sensible :)
>
> I had naively hoped to give the problem some time during the weekend,
> but I  found that having a wife and a couple of kids seriously cuts
> into one's  debugging time :)
> All I was able to put together was the small patch below, which does
> not  completely remove the instability, but at least triggers it only
> above 250  asset steps for the example code given.
>
> diff -r1.28 fddividendoption.cpp
> 74c74
> <         setGridLimits(center_ + dividends_[step], dates_[step]);
> ---
>
>>         setGridLimits(center_ + dividends_[step], residualTime_);
>
>
> Basically, the code tried to be smart (probably too much for its own
> good)  and resize the grid based on the time to the dividend payment.
> Now it  resizes it based on the time to maturity, so that it doesn't
> end up with a  grid too different from the one it started with.
>
> But the sad truth is, the whole thing should be rewritten. I had a
> look at  the code in the past few days (I had'n written it myself) and
> it's a tangle  of methods, each one trying to outsmart the user and
> correct the number of  grid points, the grid limits and whatever else
> based on what was thought to  be a "safe" choice of parameters. Not
> that it can't be useful, but the user  should at least have the choice
> to say "heck no, just do as I said."
>
> Needless to say, at this time I wouldn't use the FD dividend options
> where  actual money is involved---at least, not without cross-checking
> the results  with some other pricer or some other source.
>
> Later,
>     Luigi
>