Problems with the Finite Difference Framework

Posted by Mario Rometsch on
URL: http://quantlib.414.s1.nabble.com/Problems-with-the-Finite-Difference-Framework-tp4882.html

Hello everybody,
I'm trying to use the FD Framework in QuantLib, but unfortunately the
documentation is out-of-date and I'm failing to understand the
sourcecode (I'm a C++ newbie....). With the example from the tutorial, I
got to this point:

#include <iostream>
#include <ql/quantlib.hpp>

using namespace QuantLib;
using namespace std;

#define QL_MAX(a,b) max(a,b)
#define QL_LOG(a) log(a)
#define QL_EXP(a) exp(a)
#define QL_SQRT(a) sqrt(a)

class BlackScholesOperator : public TridiagonalOperator {
    public:
        BlackScholesOperator(
            double sigma, double nu, // parameters of the
            Rate r, // Black-Scholes equation;
            unsigned int points, // number of discretized points;
            double h) // grid spacing.
        : TridiagonalOperator(
        // build the operator by adding basic ones
        - (sigma*sigma/2.0) * DPlusDMinus(points,h)
        - nu * DZero(points,h)
        + r * TridiagonalOperator::identity(points)
        ) {}
};

class ExerciseCondition : public StepCondition<Array> {
    public:
        ExerciseCondition(const Array& exercisingValue)
            : exercisingValue_(exercisingValue) {}
        void applyTo(Array& a, Time) const {
            for (unsigned int i = 0; i < a.size(); i++) {
                a[i] = QL_MAX(a[i], exercisingValue_[i]);
            }
        }
    private:
        Array exercisingValue_;
};

int main() {
    cout << "Finite Difference Framework" << endl << endl;
    double underlying = 100.0, strike = 95.0;
    Time residualTime = 1.0;
    Rate dividendYield = 0.03, riskFreeRate = 0.05;
    double volatility = 0.10;
   
    unsigned int gridPoints = 101;
    Array grid(gridPoints), prices(gridPoints);
    double x0 = QL_LOG(underlying);
    double Delta = 4.0*volatility*QL_SQRT(residualTime);
    double xMin = x0 - Delta, xMax = x0 + Delta;
    double h = (xMax-xMin)/(gridPoints-1);
   
    for (unsigned int i=0; i<gridPoints; i++) {
        grid[i] = xMin + i*h;
        prices[i] = QL_EXP(grid[i]);
    }
   
    Array exercisingValue(gridPoints);
    for (unsigned int i=0; i<gridPoints; i++) {
        exercisingValue[i] = QL_MAX(prices[i]-strike,0.0);
    }
   
    double nu = riskFreeRate - dividendYield - volatility*volatility/2.0;
    TridiagonalOperator L = BlackScholesOperator(volatility, nu,
        riskFreeRate, gridPoints, h);
    NeumannBC
LowerBC(exercisingValue[1]-exercisingValue[0],NeumannBC::Lower),
             
UpperBC(exercisingValue[gridPoints-1]-exercisingValue[gridPoints-2],NeumannBC::Upper);
   
    LowerBC.applyBeforeApplying(L);
    UpperBC.applyBeforeApplying(L);

    FiniteDifferenceModel<CrankNicolson<TridiagonalOperator> >::bc_set
bcs();
   
    unsigned int timeSteps = 365;
    // build the model - Crank-Nicolson scheme chosen
   
    FiniteDifferenceModel<CrankNicolson<TridiagonalOperator> > model(L,bcs);

and the last line fails. Could someone provide me a working example or
something else?

Thanks Mario