final double dx = Math.log(factor);
final double r = 0.05;
final double q = 0.01;
final double sigma = 0.5;
final BSMOperator ref = new BSMOperator(grid.size(), dx, r, q, sigma);
QL.info("BSMOperator reference diagonals: \n");
outputDiagonals(ref);
final DayCounter dc = new Actual360();
final Date today = Date.todaysDate();
Date exercise = today.clone();
exercise = exercise.add(new Period(2,TimeUnit.Years));
final double residualTime = dc.yearFraction(today, exercise);
final SimpleQuote spot = new SimpleQuote(0.0);
final YieldTermStructure qTS = Utilities.flatRate(today, q, dc);
final YieldTermStructure rTS = Utilities.flatRate(today, r, dc);
final BlackVolTermStructure volTS = Utilities.flatVol(today, sigma, dc);
final GeneralizedBlackScholesProcess stochProcess = new GeneralizedBlackScholesProcess(
new Handle<Quote>(spot),
new Handle<YieldTermStructure>(qTS),
new Handle<YieldTermStructure>(rTS),
new Handle<BlackVolTermStructure>(volTS));
final BSMOperator op1 = new BSMOperator(grid, stochProcess, residualTime);
QL.info("BSMOperator diagonals: \n");
outputDiagonals(op1);
final BSMTermOperator op2 = new BSMTermOperator(grid, stochProcess, residualTime);
QL.info("PdeOperator diagonals: \n");
outputDiagonals(op2);
final double tolerance = 1.0e-6;
Array lderror = ref.lowerDiagonal().clone();
lderror.subAssign(op1.lowerDiagonal());
Array derror = ref.diagonal().clone();
derror.subAssign(op1.diagonal());
Array uderror = ref.upperDiagonal().clone();
uderror.subAssign(op1.upperDiagonal());
for (i=2; i<grid.size()-2; i++) {
QL.info("lderror(" + i + ") = "+ Math.abs(lderror.get(i)) + " tolerance " + tolerance + " \n");
QL.info("derror(" + i + ") = "+ Math.abs(derror.get(i)) + " tolerance " + tolerance + " \n");
QL.info("uderror(" + i + ") = "+ Math.abs(uderror.get(i)) + " tolerance " + tolerance + " \n");
if (Math.abs(lderror.get(i)) > tolerance ||
Math.abs(derror.get(i)) > tolerance ||
Math.abs(uderror.get(i)) > tolerance) {
fail("inconsistency between BSM operators:\n"
+ Integer.toString(i) + " row:\n"
+ "expected: " + ref.lowerDiagonal().get(i) + ", " + ref.diagonal().get(i) + ", " + ref.upperDiagonal().get(i) + "\n"
+ "calculated: " + op1.lowerDiagonal().get(i) + ", " + op1.diagonal().get(i) + ", " + op1.upperDiagonal().get(i));
}
}
lderror = ref.lowerDiagonal();
lderror.subAssign(op2.lowerDiagonal());