LocalVolPDEPricer { public double price(final double s0, final double k, final double r, final double b, final double t, final LocalVolatilitySurfaceStrike locVol, final boolean isCall, final boolean isAmerican, final int spaceNodes, final int timeNodes) { final double q = r - b; final double sigma0 = locVol.getVolatility(t, k); final double mult = Math.exp(6.0 * sigma0 * Math.sqrt(t)); final double sMin = Math.min(0.8 * k, s0 / mult); final double sMax = Math.max(1.25 * k, s0 * mult); final int tBurnNodes = (int) (USE_BURNIN ? Math.max(2, timeNodes * BURNIN_FRACTION) : 0); final double tBurn = USE_BURNIN ? BURNIN_FRACTION * t * t / timeNodes : 0.0; MeshingFunction xMesh = new ExponentialMeshing(sMin, sMax, spaceNodes, 0.0, new double[] {s0, k}); MeshingFunction tMeshBurn = USE_BURNIN ? new ExponentialMeshing(0.0, tBurn, tBurnNodes, 0.0) : null; MeshingFunction tMesh = new ExponentialMeshing(tBurn, t, timeNodes - tBurnNodes, 0.0); PDEGrid1D gridBurn = USE_BURNIN ? new PDEGrid1D(tMeshBurn, xMesh) : null; PDEGrid1D grid = new PDEGrid1D(tMesh, xMesh); final int index = Arrays.binarySearch(grid.getSpaceNodes(), s0); ArgumentChecker.isTrue(index >= 0, "cannot find spot on grid"); ConvectionDiffusionPDE1DStandardCoefficients coef = PDE.getBackwardsLocalVol(r, q, t, locVol); Function1D<Double, Double> payoff = ICP.getEuropeanPayoff(k, isCall); BoundaryCondition lower; BoundaryCondition upper; PDEResults1D res; if (isAmerican) { if (isCall) { lower = new NeumannBoundaryCondition(0.0, sMin, true); upper = new NeumannBoundaryCondition(1.0, sMax, false); } else { lower = new NeumannBoundaryCondition(-1.0, sMin, true); upper = new NeumannBoundaryCondition(0.0, sMax, false); } final Function<Double, Double> func = new Function<Double, Double>() { @Override public Double evaluate(Double... tx) { final double x = tx[1]; return isCall ? Math.max(x - k, 0.0) : Math.max(k - x, 0); } }; final FunctionalDoublesSurface free = new FunctionalDoublesSurface(func); if (USE_BURNIN) { PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> dataBurn = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(coef, payoff, lower, upper, free, gridBurn); PDEResults1D resBurn = INITIAL_SOLVER.solve(dataBurn); PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> data = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(coef, resBurn.getTerminalResults(), lower, upper, free, grid); res = SOLVER.solve(data); } else { PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> data = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(coef, payoff, lower, upper, free, grid); res = SOLVER.solve(data); } } else { if (isCall) { lower = new NeumannBoundaryCondition(0.0, sMin, true); Function1D<Double, Double> upFunc = new Function1D<Double, Double>() { @Override public Double evaluate(Double time) { return Math.exp(-q * time); } }; upper = new NeumannBoundaryCondition(upFunc, sMax, false); } else { Function1D<Double, Double> downFunc = new Function1D<Double, Double>() { @Override public Double evaluate(Double time) { return -Math.exp(-q * time); } }; lower = new NeumannBoundaryCondition(downFunc, sMin, true); upper = new NeumannBoundaryCondition(0.0, sMax, false); } if (USE_BURNIN) { PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> dataBurn = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(coef, payoff, lower, upper, gridBurn); PDEResults1D resBurn = INITIAL_SOLVER.solve(dataBurn); PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> data = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(coef, resBurn.getTerminalResults(), lower, upper, grid); res = SOLVER.solve(data); } else { PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> data = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(coef, payoff, lower, upper, grid); res = SOLVER.solve(data); } } return res.getFunctionValue(index); }  double price(final double s0, final double k, final double r, final double b, final double t, final LocalVolatilitySurfaceStrike locVol, final boolean isCall, final boolean isAmerican,
      final int spaceNodes, final int timeNodes);  }
@Test public void testFlat() { final double t = 1.5; final double spot = 100.0; final double r = 0.05; final double q = 0; final double b = r - q; final double k = 120.0; final double sigma = 0.4; final boolean isCall = true; final ForwardCurve fc = new ForwardCurve(spot, r); final LocalVolatilitySurfaceStrike lv = new LocalVolatilitySurfaceStrike(ConstantDoublesSurface.from(sigma)); final int tN = 60; final int sN = 80 * tN; final double pdePrice = PRICER.price(spot, k, r, b, t, lv, true, false, sN, tN); final double price = Math.exp(-r * t) * BlackFormulaRepository.price(fc.getForward(t), k, t, sigma, isCall); assertEquals(price, pdePrice, 1e-5 * price); }
@Test public void test() { final double t = 1.5; final double spot = 100.0; final double r = 0.1; final double q = 0.07; final double b = r - q; final double k = 120.0; final boolean isCall = true; final ForwardCurve fc = new ForwardCurve(spot, b); final YieldAndDiscountCurve discountCurve = new YieldCurve("test", ConstantDoublesCurve.from(r)); final double[] w = new double[] {0.7, 0.25, 0.05 }; final double[] sigma = new double[] {0.3, 0.6, 1.0 }; final double[] mu = new double[] {0.0, 0.3, -0.5 }; final MultiHorizonMixedLogNormalModelData data = new MultiHorizonMixedLogNormalModelData(w, sigma, mu); final LocalVolatilitySurfaceStrike lv = MixedLogNormalVolatilitySurface.getLocalVolatilitySurface(fc, data); final int tN = 400; final int sN = 2 * tN; final double pdePrice = PRICER.price(spot, k, r, b, t, lv, true, false, sN, tN); final double price = MixedLogNormalVolatilitySurface.getPriceSurface(fc, discountCurve, data).getPrice(t, k); final double vol = MixedLogNormalVolatilitySurface.getImpliedVolatilitySurface(fc, data).getVolatility(t, k); final double price2 = Math.exp(-r * t) * BlackFormulaRepository.price(spot * Math.exp(b * t), k, t, vol, isCall); assertEquals(price, pdePrice, 1e-3 * price); }