* Test the pricing of a 1D option with the local volatility surface produced by the mixed log-normal model solving both the forward and backwards PDE, and compare
* with the result from the implied volatility surface. Neither solving the forward nor the backwards PDE give great accuracy
*/
@Test
public void shortDatedOptionTest() {
final ConvectionDiffusionPDESolver solver = new ThetaMethodFiniteDifference(0.5, false);
final PDE1DCoefficientsProvider pdeProvider = new PDE1DCoefficientsProvider();
final InitialConditionsProvider initialConProvider = new InitialConditionsProvider();
//set up the mixed log-normal model
final double[] weights = new double[] {0.1, 0.8, 0.1 };
final double[] sigmas = new double[] {0.15, 0.5, 0.9 };
final double[] mus = new double[] {-0.1, 0, 0.1 };
final MultiHorizonMixedLogNormalModelData mlnData = new MultiHorizonMixedLogNormalModelData(weights, sigmas, mus);
final double spot = 100.;
final double r = 0.1;
final ForwardCurve fwdCurve = new ForwardCurve(spot, r);
final double t = 1.0 / 365.;
final double rootT = Math.sqrt(t);
final double fwd = fwdCurve.getForward(t);
BlackVolatilitySurfaceStrike ivs = MixedLogNormalVolatilitySurface.getImpliedVolatilitySurface(fwdCurve, mlnData);
LocalVolatilitySurfaceStrike lvs = MixedLogNormalVolatilitySurface.getLocalVolatilitySurface(fwdCurve, mlnData);
LocalVolatilitySurfaceMoneyness lvsm = LocalVolatilitySurfaceConverter.toMoneynessSurface(lvs, fwdCurve);
// PDEUtilityTools.printSurface("imp vol", ivs.getSurface(), 0, t, 0.9, 1.1);
// DupireLocalVolatilityCalculator dupire = new DupireLocalVolatilityCalculator();
// LocalVolatilitySurfaceMoneyness dlv = dupire.getLocalVolatility(BlackVolatilitySurfaceConverter.toMoneynessSurface(ivs, fwdCurve));
// PDEUtilityTools.printSurface("local vol (dupire)", dlv.getSurface(), 0, t, 0.99, 1.01);
// PDEUtilityTools.printSurface("local vol", lvsm.getSurface(), 0, t, 0.99, 1.01);
//set up for solving the forward PDE
//TODO shunt this setup into its own class
ConvectionDiffusionPDE1DStandardCoefficients pde = pdeProvider.getForwardLocalVol(lvsm);
Function1D<Double, Double> initialCond = initialConProvider.getForwardCallPut(true);
double xL = 0.8;
double xH = 1.2;
BoundaryCondition lower = new NeumannBoundaryCondition(-1.0, xL, true);
BoundaryCondition upper = new NeumannBoundaryCondition(0.0, xH, false);
final MeshingFunction spaceMeshF = new HyperbolicMeshing(xL, xH, 1.0, 200, 0.001);
final MeshingFunction timeMeshF = new ExponentialMeshing(0, t, 50, 4.0);
final MeshingFunction timeMeshB = new DoubleExponentialMeshing(0, t, t / 2, 50, 2.0, -4.0);
final PDEGrid1D grid = new PDEGrid1D(timeMeshF, spaceMeshF);
PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> dbF = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(pde, initialCond, lower, upper, grid);
PDETerminalResults1D res = (PDETerminalResults1D) solver.solve(dbF);
final double minK = Math.exp(-6 * rootT);
final double maxK = Math.exp(6 * rootT);
Map<Double, Double> vols = PDEUtilityTools.priceToImpliedVol(fwdCurve, t, res, minK, maxK, true);
DoubleQuadraticInterpolator1D interpolator = Interpolator1DFactory.DOUBLE_QUADRATIC_INSTANCE;
Interpolator1DDataBundle idb = interpolator.getDataBundle(vols);
//set up for solving backwards PDE
ConvectionDiffusionPDE1DStandardCoefficients pdeB = pdeProvider.getBackwardsLocalVol(t, lvsm);
double sL = xL * spot;
double sH = xH * spot;
final MeshingFunction spaceMeshB = new HyperbolicMeshing(sL, sH, spot, 200, 0.001);
final PDEGrid1D gridB = new PDEGrid1D(timeMeshB, spaceMeshB);
int index = SurfaceArrayUtils.getLowerBoundIndex(gridB.getSpaceNodes(), spot);
double s1 = gridB.getSpaceNode(index);
double s2 = gridB.getSpaceNode(index + 1);
final double w = (s2 - spot) / (s2 - s1);
//solve a separate backwards PDE for each strike
for (int i = 0; i < 10; i++) {
double z = -5 + i;
double k = spot * Math.exp(0.4 * rootT * z);
double x = k / fwd;
double vol = ivs.getVolatility(t, k);
double volFPDE = interpolator.interpolate(idb, x);
boolean isCall = (k >= fwd);
BoundaryCondition lowerB = new NeumannBoundaryCondition(isCall ? 0 : -1, sL, true);
BoundaryCondition upperB = new NeumannBoundaryCondition(isCall ? 1 : 0, sH, false);
Function1D<Double, Double> bkdIC = initialConProvider.getEuropeanPayoff(k, isCall);
PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> dbB = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(pdeB, bkdIC, lowerB, upperB, gridB);
PDEResults1D resB = solver.solve(dbB);
double price1 = resB.getFunctionValue(index);
double price2 = resB.getFunctionValue(index + 1);
double vol1 = BlackFormulaRepository.impliedVolatility(price1, s1, k, t, isCall);
double vol2 = BlackFormulaRepository.impliedVolatility(price2, s2, k, t, isCall);
double volBPDE = w * vol1 + (1 - w) * vol2;