final LocalVolatilitySurfaceStrike lvDown = new LocalVolatilitySurfaceStrike(SurfaceShiftFunctionFactory.getShiftedSurface(localVol.getSurface(), -volShift, true));
//first order shifts
final PDEFullResults1D pdeRes = runForwardPDESolver(forwardCurve, localVol, _isCall, _theta, maxT, maxProxyDelta,
_timeSteps, _spaceSteps, _timeGridBunching, _spaceGridBunching, 1.0);
final PDEResults1D pdeResUp = runForwardPDESolver(forwardCurve, lvUp, _isCall, _theta, maxT, maxProxyDelta,
_timeSteps, _spaceSteps, _timeGridBunching, _spaceGridBunching, 1.0);
final PDEResults1D pdeResDown = runForwardPDESolver(forwardCurve, lvDown, _isCall, _theta, maxT, maxProxyDelta,
_timeSteps, _spaceSteps, _timeGridBunching, _spaceGridBunching, 1.0);
//second order shifts
final PDEResults1D pdeResUpUp = runForwardPDESolver(forwardCurve.withFractionalShift(fwdShift), lvUp, _isCall, _theta, maxT, maxProxyDelta,
_timeSteps, _spaceSteps, _timeGridBunching, _spaceGridBunching, 1.0);
final PDEFullResults1D pdeResUpDown = runForwardPDESolver(forwardCurve.withFractionalShift(fwdShift), lvDown, _isCall, _theta, maxT, maxProxyDelta,
_timeSteps, _spaceSteps, _timeGridBunching, _spaceGridBunching, 1.0);
final PDEFullResults1D pdeResDownUp = runForwardPDESolver(forwardCurve.withFractionalShift(-fwdShift), lvUp, _isCall, _theta, maxT, maxProxyDelta,
_timeSteps, _spaceSteps, _timeGridBunching, _spaceGridBunching, 1.0);
final PDEFullResults1D pdeResDownDown = runForwardPDESolver(forwardCurve.withFractionalShift(-fwdShift), lvDown, _isCall, _theta, maxT, maxProxyDelta,
_timeSteps, _spaceSteps, _timeGridBunching, _spaceGridBunching, 1.0);
ps.println("Strike\tBS Vega\tVega\tBS Vanna\tVanna\tBS Vomma\tVomma");
final int n = pdeRes.getNumberSpaceNodes();
for (int i = 0; i < n; i++) {
final double x = pdeRes.getSpaceValue(i);
final double k = x * forward;
final double mPrice = pdeRes.getFunctionValue(i);
try {
final double bsVol = BlackFormulaRepository.impliedVolatility(mPrice, 1.0, x, maxT, _isCall);
final double bsVega = BlackFormulaRepository.vega(forward, k, maxT, bsVol);
final double bsVanna = BlackFormulaRepository.vanna(forward, k, maxT, bsVol);
final double bsVomma = BlackFormulaRepository.vomma(forward, k, maxT, bsVol);
final double modelVega = forward * (pdeResUp.getFunctionValue(i) - pdeResDown.getFunctionValue(i)) / 2 / volShift;
//xVanna is the vanna if the moneyness parameterised local vol surface was invariant to changes in the forward curve
final double xVanna = (pdeResUp.getFunctionValue(i) - pdeResDown.getFunctionValue(i)
- x * (pdeResUp.getFirstSpatialDerivative(i) - pdeResDown.getFirstSpatialDerivative(i))) / 2 / volShift;
//this is the vanna coming purely from deformation of the local volatility surface
final double surfaceVanna = (pdeResUpUp.getFunctionValue(i) + pdeResDownDown.getFunctionValue(i) -
pdeResUpDown.getFunctionValue(i) - pdeResDownUp.getFunctionValue(i)) / 4 / fwdShift / volShift;
final double modelVanna = xVanna + surfaceVanna;
final double modelVomma = forward * (pdeResUp.getFunctionValue(i) + pdeResDown.getFunctionValue(i)
- 2 * pdeRes.getFunctionValue(i)) / volShift / volShift;
ps.println(k + "\t" + bsVega + "\t" + modelVega + "\t" + bsVanna + "\t" + modelVanna + "\t" + bsVomma + "\t" + modelVomma);