CoupledFiniteDifference { public PDEResults1D[] solve(final CoupledPDEDataBundle pdeData1, final CoupledPDEDataBundle pdeData2) { Validate.notNull(pdeData1, "pde1 data"); Validate.notNull(pdeData2, "pde2 data"); final PDEGrid1D grid = pdeData1.getGrid(); ArgumentChecker.isTrue(grid == pdeData2.getGrid(), "grids must be same object"); final ConvectionDiffusionPDE1DCoupledCoefficients coeff1 = pdeData1.getCoefficients(); final ConvectionDiffusionPDE1DCoupledCoefficients coeff2 = pdeData2.getCoefficients(); final double[] initalCond1 = pdeData1.getInitialCondition(); final double[] initalCond2 = pdeData2.getInitialCondition(); final BoundaryCondition lower1 = pdeData1.getLowerBoundary(); final BoundaryCondition lower2 = pdeData2.getLowerBoundary(); final BoundaryCondition upper1 = pdeData1.getUpperBoundary(); final BoundaryCondition upper2 = pdeData2.getUpperBoundary(); final double lambda1 = coeff1.getLambda(); final double lambda2 = coeff2.getLambda(); final int tNodes = grid.getNumTimeNodes(); final int xNodes = grid.getNumSpaceNodes(); double[] f = new double[2 * xNodes]; double[][] full1 = null; double[][] full2 = null; if (_showFullResults) { full1 = new double[tNodes][xNodes]; full2 = new double[tNodes][xNodes]; } final double[] q = new double[2 * xNodes]; final double[][] m = new double[2 * xNodes][2 * xNodes]; double[][] a1 = new double[2][xNodes - 2]; final double[][] a2 = new double[2][xNodes - 2]; double[][] b1 = new double[2][xNodes - 2]; final double[][] b2 = new double[2][xNodes - 2]; double[][] c1 = new double[2][xNodes - 2]; final double[][] c2 = new double[2][xNodes - 2]; double dt, t1, t2, x; double[] x1st, x2nd; System.arraycopy(initalCond1, 0, f, 0, xNodes); System.arraycopy(initalCond2, 0, f, xNodes, xNodes); if (_showFullResults) { if (full1 != null && full2 != null) { full1[0] = Arrays.copyOfRange(f, 0, xNodes); full2[0] = Arrays.copyOfRange(f, xNodes, 2 * xNodes); } } for (int i = 0; i < xNodes - 2; i++) { x = grid.getSpaceNode(i + 1); a1[0][i] = coeff1.getA(0, x); b1[0][i] = coeff1.getB(0, x); c1[0][i] = coeff1.getC(0, x); a1[1][i] = coeff2.getA(0, x); b1[1][i] = coeff2.getB(0, x); c1[1][i] = coeff2.getC(0, x); } final boolean first = true; DecompositionResult decompRes = null; for (int n = 1; n < tNodes; n++) { t1 = grid.getTimeNode(n - 1); t2 = grid.getTimeNode(n); dt = grid.getTimeStep(n - 1); for (int i = 1; i < xNodes - 1; i++) { x = grid.getSpaceNode(i); x1st = grid.getFirstDerivativeCoefficients(i); x2nd = grid.getSecondDerivativeCoefficients(i); q[i] = f[i]; q[i] -= (1 - _theta) * dt * (x2nd[0] * a1[0][i - 1] + x1st[0] * b1[0][i - 1]) * f[i - 1]; q[i] -= (1 - _theta) * dt * (x2nd[1] * a1[0][i - 1] + x1st[1] * b1[0][i - 1] + c1[0][i - 1]) * f[i]; q[i] -= (1 - _theta) * dt * (x2nd[2] * a1[0][i - 1] + x1st[2] * b1[0][i - 1]) * f[i + 1]; q[i] -= (1 - _theta) * dt * lambda1 * f[i + xNodes]; q[xNodes + i] = f[xNodes + i]; q[xNodes + i] -= (1 - _theta) * dt * (x2nd[0] * a1[1][i - 1] + x1st[0] * b1[1][i - 1]) * f[xNodes + i - 1]; q[xNodes + i] -= (1 - _theta) * dt * (x2nd[1] * a1[1][i - 1] + x1st[1] * b1[1][i - 1] + c1[1][i - 1]) * f[xNodes + i]; q[xNodes + i] -= (1 - _theta) * dt * (x2nd[2] * a1[1][i - 1] + x1st[2] * b1[1][i - 1]) * f[xNodes + i + 1]; q[xNodes + i] -= (1 - _theta) * dt * lambda2 * f[i]; a2[0][i - 1] = coeff1.getA(t2, x); b2[0][i - 1] = coeff1.getB(t2, x); c2[0][i - 1] = coeff1.getC(t2, x); a2[1][i - 1] = coeff2.getA(t2, x); b2[1][i - 1] = coeff2.getB(t2, x); c2[1][i - 1] = coeff2.getC(t2, x); m[i][i - 1] = _theta * dt * (x2nd[0] * a2[0][i - 1] + x1st[0] * b2[0][i - 1]); m[i][i] = 1 + _theta * dt * (x2nd[1] * a2[0][i - 1] + x1st[1] * b2[0][i - 1] + c2[0][i - 1]); m[i][i + 1] = _theta * dt * (x2nd[2] * a2[0][i - 1] + x1st[2] * b2[0][i - 1]); m[i][i + xNodes] = dt * _theta * lambda1; m[xNodes + i][xNodes + i - 1] = _theta * dt * (x2nd[0] * a2[1][i - 1] + x1st[0] * b2[1][i - 1]); m[xNodes + i][xNodes + i] = 1 + _theta * dt * (x2nd[1] * a2[1][i - 1] + x1st[1] * b2[1][i - 1] + c2[1][i - 1]); m[xNodes + i][xNodes + i + 1] = _theta * dt * (x2nd[2] * a2[1][i - 1] + x1st[2] * b2[1][i - 1]); m[xNodes + i][i] = dt * _theta * lambda2; } double[] temp = lower1.getLeftMatrixCondition(pdeData1.getCoefficients(), grid, t2); for (int k = 0; k < temp.length; k++) { m[0][k] = temp[k]; } temp = upper1.getLeftMatrixCondition(pdeData1.getCoefficients(), grid, t2); for (int k = 0; k < temp.length; k++) { m[xNodes - 1][xNodes - temp.length + k] = temp[k]; } temp = lower2.getLeftMatrixCondition(pdeData2.getCoefficients(), grid, t2); for (int k = 0; k < temp.length; k++) { m[xNodes][xNodes + k] = temp[k]; } temp = upper2.getLeftMatrixCondition(pdeData2.getCoefficients(), grid, t2); for (int k = 0; k < temp.length; k++) { m[2 * xNodes - 1][2 * xNodes - temp.length + k] = temp[k]; } temp = lower1.getRightMatrixCondition(pdeData1.getCoefficients(), grid, t1); double sum = 0; for (int k = 0; k < temp.length; k++) { sum += temp[k] * f[k]; } q[0] = sum + lower1.getConstant(pdeData1.getCoefficients(), t2); temp = upper1.getRightMatrixCondition(pdeData1.getCoefficients(), grid, t1); sum = 0; for (int k = 0; k < temp.length; k++) { sum += temp[k] * f[xNodes - 1 - k]; } q[xNodes - 1] = sum + upper1.getConstant(pdeData1.getCoefficients(), t2); temp = lower2.getRightMatrixCondition(pdeData2.getCoefficients(), grid, t1); sum = 0; for (int k = 0; k < temp.length; k++) { sum += temp[k] * f[k]; } q[xNodes] = sum + lower2.getConstant(pdeData2.getCoefficients(), t2); temp = upper2.getRightMatrixCondition(pdeData2.getCoefficients(), grid, t1); sum = 0; for (int k = 0; k < temp.length; k++) { sum += temp[k] * f[xNodes - 1 - k]; } q[2 * xNodes - 1] = sum + upper2.getConstant(pdeData2.getCoefficients(), t2); if (first) { final DoubleMatrix2D mM = new DoubleMatrix2D(m); decompRes = DCOMP.evaluate(mM); } f = decompRes.solve(q); a1 = a2; b1 = b2; c1 = c2; if (_showFullResults) { if (full1 != null && full2 != null) { full1[n] = Arrays.copyOfRange(f, 0, xNodes); full2[n] = Arrays.copyOfRange(f, xNodes, 2 * xNodes); } } } final PDEResults1D[] res = new PDEResults1D[2]; if (_showFullResults) { res[0] = new PDEFullResults1D(grid, full1); res[1] = new PDEFullResults1D(grid, full2); } else { final double[] res1 = Arrays.copyOfRange(f, 0, xNodes); final double[] res2 = Arrays.copyOfRange(f, xNodes, 2 * xNodes); res[0] = new PDETerminalResults1D(grid, res1); res[1] = new PDETerminalResults1D(grid, res2); } return res; } CoupledFiniteDifference();  CoupledFiniteDifference(final double theta, final boolean showFullResults); double getTheta(); boolean showFullResults(); Decomposition<?> getDecomposition(); PDEResults1D[] solve(final CoupledPDEDataBundle pdeData1,
      final CoupledPDEDataBundle pdeData2);  }
@Test public void testNoCoupling() { final CoupledFiniteDifference solver = new CoupledFiniteDifference(0.5, false); final double lambda12 = 0.0; final double lambda21 = 0.0; final double p0 = 0.5; TwoStateMarkovChainDataBundle chainData = new TwoStateMarkovChainDataBundle(VOL1, VOL2, lambda12, lambda21, p0); ConvectionDiffusionPDE1DCoupledCoefficients[] pdeData = PDE_DATA_PROVIDER.getCoupledBackwardsPair(FORWARD, T, chainData); final int timeNodes = 20; final int spaceNodes = 150; final double lowerMoneyness = 0.4; final double upperMoneyness = 2.5; final MeshingFunction timeMesh = new ExponentialMeshing(0, T, timeNodes, 5.0); final MeshingFunction spaceMesh = new HyperbolicMeshing(LOWER.getLevel(), UPPER.getLevel(), OPTION.getStrike(), spaceNodes, 0.1); final double[] timeGrid = new double[timeNodes]; for (int n = 0; n < timeNodes; n++) { timeGrid[n] = timeMesh.evaluate(n); } final double[] spaceGrid = new double[spaceNodes]; for (int i = 0; i < spaceNodes; i++) { spaceGrid[i] = spaceMesh.evaluate(i); } final PDEGrid1D grid = new PDEGrid1D(timeGrid, spaceGrid); CoupledPDEDataBundle d1 = new CoupledPDEDataBundle(pdeData[0], INITIAL_COND, LOWER, UPPER, grid); CoupledPDEDataBundle d2 = new CoupledPDEDataBundle(pdeData[1], INITIAL_COND, LOWER, UPPER, grid); final PDEResults1D[] res = solver.solve(d1, d2); final double df = YIELD_CURVE.getDiscountFactor(T); final int n = res[0].getNumberSpaceNodes(); for (int i = 0; i < n; i++) { final double spot = res[0].getSpaceValue(i); final double price1 = res[0].getFunctionValue(i); final double price2 = res[1].getFunctionValue(i); final double moneyness = spot / OPTION.getStrike(); final BlackFunctionData data = new BlackFunctionData(spot / df, df, 0.0); double impVol1; try { impVol1 = BLACK_IMPLIED_VOL.getImpliedVolatility(data, OPTION, price1); } catch (final Exception e) { impVol1 = 0.0; } double impVol2; try { impVol2 = BLACK_IMPLIED_VOL.getImpliedVolatility(data, OPTION, price2); } catch (final Exception e) { impVol2 = 0.0; } if (moneyness >= lowerMoneyness && moneyness <= upperMoneyness) { assertEquals(VOL1, impVol1, 1e-3); assertEquals(VOL2, impVol2, 5e-3); } } }
@Test public void testDegenerate() { final CoupledFiniteDifference solver = new CoupledFiniteDifference(); final double lambda12 = 0.2; final double lambda21 = 0.5; TwoStateMarkovChainDataBundle chainData = new TwoStateMarkovChainDataBundle(VOL1, VOL1, lambda12, lambda21, 0.5); ConvectionDiffusionPDE1DCoupledCoefficients[] pdeData = PDE_DATA_PROVIDER.getCoupledBackwardsPair(FORWARD, T, chainData); final int timeNodes = 20; final int spaceNodes = 150; final double lowerMoneyness = 0.4; final double upperMoneyness = 2.5; final MeshingFunction timeMesh = new ExponentialMeshing(0, T, timeNodes, 5.0); final MeshingFunction spaceMesh = new HyperbolicMeshing(LOWER.getLevel(), UPPER.getLevel(), OPTION.getStrike(), spaceNodes, 0.1); final double[] timeGrid = new double[timeNodes]; for (int n = 0; n < timeNodes; n++) { timeGrid[n] = timeMesh.evaluate(n); } final double[] spaceGrid = new double[spaceNodes]; for (int i = 0; i < spaceNodes; i++) { spaceGrid[i] = spaceMesh.evaluate(i); } final PDEGrid1D grid = new PDEGrid1D(timeGrid, spaceGrid); final CoupledPDEDataBundle d1 = new CoupledPDEDataBundle(pdeData[0], INITIAL_COND, LOWER, UPPER, grid); final CoupledPDEDataBundle d2 = new CoupledPDEDataBundle(pdeData[1], INITIAL_COND, LOWER, UPPER, grid); final PDEResults1D[] res = solver.solve(d1, d2); final double df = YIELD_CURVE.getDiscountFactor(T); final int n = res[0].getNumberSpaceNodes(); for (int i = 0; i < n; i++) { final double spot = res[0].getSpaceValue(i); final double price1 = res[0].getFunctionValue(i); final double price2 = res[1].getFunctionValue(i); final double moneyness = spot / OPTION.getStrike(); final BlackFunctionData data = new BlackFunctionData(spot / df, df, 0.0); double impVol1; try { impVol1 = BLACK_IMPLIED_VOL.getImpliedVolatility(data, OPTION, price1); } catch (final Exception e) { impVol1 = 0.0; } double impVol2; try { impVol2 = BLACK_IMPLIED_VOL.getImpliedVolatility(data, OPTION, price2); } catch (final Exception e) { impVol2 = 0.0; } if (moneyness >= lowerMoneyness && moneyness <= upperMoneyness) { assertEquals(impVol1, impVol2, 1e-8); assertEquals(VOL1, impVol1, 2e-3); } } }
@Test public void testSmile() { final CoupledFiniteDifference solver = new CoupledFiniteDifference(0.5, false); final double lambda12 = 0.2; final double lambda21 = 2.0; TwoStateMarkovChainDataBundle chainData = new TwoStateMarkovChainDataBundle(VOL1, VOL2, lambda12, lambda21, 1.0); ConvectionDiffusionPDE1DCoupledCoefficients[] pdeData = PDE_DATA_PROVIDER.getCoupledBackwardsPair(FORWARD, T, chainData); final int timeNodes = 20; final int spaceNodes = 150; final double lowerMoneyness = 0.3; final double upperMoneyness = 3.0; final MarkovChain mc = new MarkovChain(VOL1, VOL2, lambda12, lambda21, 1.0); final double[] mcSims = mc.simulate(T, 10000); final MeshingFunction timeMesh = new ExponentialMeshing(0, T, timeNodes, 0); final MeshingFunction spaceMesh = new ExponentialMeshing(LOWER.getLevel(), UPPER.getLevel(), spaceNodes, 0.0); final double[] timeGrid = new double[timeNodes]; for (int n = 0; n < timeNodes; n++) { timeGrid[n] = timeMesh.evaluate(n); } final double[] spaceGrid = new double[spaceNodes]; for (int i = 0; i < spaceNodes; i++) { spaceGrid[i] = spaceMesh.evaluate(i); } final PDEGrid1D grid = new PDEGrid1D(timeGrid, spaceGrid); final CoupledPDEDataBundle d1 = new CoupledPDEDataBundle(pdeData[0], INITIAL_COND, LOWER, UPPER, grid); final CoupledPDEDataBundle d2 = new CoupledPDEDataBundle(pdeData[1], INITIAL_COND, LOWER, UPPER, grid); final PDEResults1D[] res = solver.solve(d1, d2); final double df = YIELD_CURVE.getDiscountFactor(T); final int n = res[0].getNumberSpaceNodes(); for (int i = 0; i < n; i++) { final double spot = res[0].getSpaceValue(i); final double price1 = res[0].getFunctionValue(i); final double price2 = res[1].getFunctionValue(i); final double delta = res[0].getFirstSpatialDerivative(i); final double gamma = res[0].getSecondSpatialDerivative(i); final double moneyness = spot / OPTION.getStrike(); final BlackFunctionData data = new BlackFunctionData(spot / df, df, 0.0); final double mc_price = mc.price(spot / df, df, STRIKE, T, mcSims); double impVol1; try { impVol1 = BLACK_IMPLIED_VOL.getImpliedVolatility(data, OPTION, price1); } catch (final Exception e) { impVol1 = 0.0; } double impVol2; try { impVol2 = BLACK_IMPLIED_VOL.getImpliedVolatility(data, OPTION, price2); } catch (final Exception e) { impVol2 = 0.0; } double impVol_mc; try { impVol_mc = BLACK_IMPLIED_VOL.getImpliedVolatility(data, OPTION, mc_price); } catch (final Exception e) { impVol_mc = 0.0; } if (moneyness >= lowerMoneyness && moneyness <= upperMoneyness) { assertEquals(impVol_mc, impVol1, 1e-2); } } }