Finite Difference Methods

[animesh]

Hi,
    I am trying to code a simple Finite Difference method (Crank
Nicholson Scheme). I was taking help from reference manual but it seems
the manual is quite old. There have been lots of changes in the version
I have which is 1.0.1. So after lots of hits and trials I came up with
this errorenous piece of code. Can anyone help me with this. If a sample
FD method is available it will be helpful.

Thanks,

Animesh


#include <ql/quantlib.hpp>
#include <limits.h>
#include <ql/mathconstants.hpp>

#ifdef BOOST_MSVC
/* Uncomment the following lines to unmask floating-point
 exceptions. Warning: unpredictable results can arise...
 
 See http://www.wilmott.com/messageview.cfm?catid=10&threadid=9481
 Is there anyone with a definitive word about this?
 */
// #include <float.h>
// namespace { unsigned int u = _controlfp(_EM_INEXACT, _MCW_EM); }
#endif

#include <boost/timer.hpp>
#include <iostream>
#include <iomanip>
#include <float.h>
#include <math.h>

#define MAX_PATHS 1024
using namespace QuantLib;

#if defined(QL_ENABLE_SESSIONS)
namespace QuantLib {
   
    Integer sessionId() { return 0; }
   
}
#endif
 std::vector <double> paths(MAX_PATHS);
 std::vector <double> variates(MAX_PATHS);
 
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] = std::max(a[i], exercisingValue_[i]);
    }
    private:
    Array exercisingValue_;
};

class BlackScholesOperator : public TridiagonalOperator {
public:

    BlackScholesOperator( double sigma, double nu, Rate r, unsigned int
points, double h)
: TridiagonalOperator( - (sigma*sigma/2.0) * DPlusDMinus(points,h) - nu
* DZero(points,h) + r * TridiagonalOperator::identity(points)
                      ) {}
};

std::vector<Time> stoppingTimes = std::vector<Time>();


int main()
{

    Option::Type type = Option::Call;
    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 = log(underlying);
    double Delta = 4.0*volatility*pow(residualTime,0.5);
    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] = exp(grid[i]);
    }
   
    Array exercisingValue(gridPoints);
    for (unsigned int i=0; i<gridPoints; i++)
        exercisingValue[i] = ((prices[i]-strike)>0 ?
(prices[i]-strike):0.0);
   
    double nu = riskFreeRate - dividendYield - volatility*volatility/2.0;
    BlackScholesOperator L1 = BlackScholesOperator(volatility,
nu,riskFreeRate, gridPoints, h);
   
    BoundaryCondition<TridiagonalOperator>::Side s2 =
BoundaryCondition<TridiagonalOperator>::None;
    Real s1;
   
//    NeumannBC(Real value, Side side);

    NeumannBC bCondition = NeumannBC(s1,s2);
   
   
//    
L1.setLowerBC(BoundaryCondition(BoundaryCondition::Neumann,exercisingValue[1]-exercisingValue[0]));
   
//    
L1.setUpperBC(BoundaryCondition(BoundaryCondition::Neumann,exercisingValue[gridPoints_-1]-exercisingValue[gridPoints_-2]));

   
    OperatorTraits<TridiagonalOperator>::operator_type t11;
//    OperatorTraits<TridiagonalOperator>::bc_type bc1;
   
   
//    typedef Operator operator_type;
//    typedef typename Operator::array_type array_type;
//    typedef BoundaryCondition<operator_type> bc_type;
//    typedef std::vector<boost::shared_ptr<bc_type> > bc_set;
//    typedef StepCondition<array_type> condition_type;
   
   
   
    CrankNicolson<BlackScholesOperator> cn(t11,bCondition);
    unsigned int timeSteps = 365;
    FiniteDifferenceModel<CrankNicolson<TridiagonalOperator> > model =
    FiniteDifferenceModel<CrankNicolson<TridiagonalOperator>
 >(cn,stoppingTimes);
   
   
    Array f = exercisingValue; // initial condition
   
    model.rollback(f, residualTime, 0.0, timeSteps);
    double europeanValue = valueAtCenter(f);
    f = exercisingValue;    // reset
    Handle<StepCondition<Array> > condition(new
ExerciseCondition(exercisingValue));
   
    model.rollback(f, residualTime, 0.0, timeSteps, condition);
    double americanValue = valueAtCenter(f);
   
   
   
   
   
}


------------------------------------------------------------------------------
This SF.net email is sponsored by Sprint
What will you do first with EVO, the first 4G phone?
Visit sprint.com/first -- http://p.sf.net/sfu/sprint-com-first
_______________________________________________
QuantLib-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-users