PdeOptionValuation.java

package net.finmath.equities.pricer;

import java.time.LocalDate;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

import org.apache.commons.math3.linear.DecompositionSolver;
import org.apache.commons.math3.linear.LUDecomposition;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;

import net.finmath.equities.marketdata.YieldCurve;
import net.finmath.equities.models.EquityForwardStructure;
import net.finmath.equities.models.FlatVolatilitySurface;
import net.finmath.equities.models.VolatilitySurface;
import net.finmath.equities.pricer.EquityValuationRequest.CalculationRequestType;
import net.finmath.equities.products.EuropeanOption;
import net.finmath.equities.products.Option;
import net.finmath.rootfinder.BisectionSearch;
import net.finmath.rootfinder.SecantMethod;
import net.finmath.time.daycount.DayCountConvention;

/**
 * This class implements a finite difference pricer under a Black-Scholes process or a
 * local volatility process in the presence of Buehler dividends.
 * It supports European and American options. Greeks are calculated inside the grid to the extent possible.
 * The implementation uses a Crank-Nicolson scheme.
 * Payoffs are smoothed using the modified timestepping from Rannacher's 1984 paper.
 * The American exercise feature is priced using the penalty approach from Forsyth's 2001 paper.
 *
 * TODO The linear algebra framework used (apache.commons.math3) is not optimized for the discretized PDE.
 * More performant linear algebra algorithms should be used that take account of
 * the tridiagonal matrix structure of the problem, e.g. direct diagonal operations instead of full-blown
 * matrix multiplication, and the Thomas algorithm instead of LU decomposition for solving equations.
 *
 * @author Andreas Grotz
 */
public class PdeOptionValuation implements OptionValuation {

	private final int timeStepsPerYear;
	private final double spaceMinForwardMultiple;
	private final double spaceMaxForwardMultiple;
	private final int spaceNbOfSteps;
	private final double spaceStepSize;
	private final List<Double> spots;
	private final int spotIndex;
	private final DayCountConvention dayCounter;
	private final boolean isLvPricer;
	private final boolean includeDividendDatesInGrid;

	public PdeOptionValuation(
			double spaceMinForwardMultiple,
			double spaceMaxForwardMultiple,
			int spaceNbPoints,
			final int timeStepsPerYear,
			DayCountConvention dcc,
			final boolean isLvPricer,
			final boolean includeDividendDatesInGrid) {

		assert spaceMinForwardMultiple < 1.0 : "min multiple of forward must be below 1.0";
		assert spaceMaxForwardMultiple > 1.0 : "max multiple of forward must be below 1.0";

		this.timeStepsPerYear = timeStepsPerYear;
		this.dayCounter = dcc;
		this.isLvPricer = isLvPricer;
		this.includeDividendDatesInGrid = includeDividendDatesInGrid;

		// Set up the space grid for the pure volatility process
		double tmpSpaceStepSize = (spaceMaxForwardMultiple - spaceMinForwardMultiple) / spaceNbPoints;
		int tmpSpaceNbPoints = spaceNbPoints;
		List<Double> tmpSpots = new ArrayList<Double>();
		for(int i = 0; i < tmpSpaceNbPoints; i++) {
			tmpSpots.add(spaceMinForwardMultiple + tmpSpaceStepSize * i);
		}

		// The space grid needs to include the forward level 1.0 for the pure volatility process.
		// Hence if necessary, we increase the step size slightly to include it.
		final int lowerBound = Math.abs(Collections.binarySearch(tmpSpots, 1.0)) - 2;
		if(!(tmpSpots.get(lowerBound) == 1.0)) {
			tmpSpaceStepSize += (1.0 - tmpSpots.get(lowerBound)) / lowerBound;
			tmpSpots = new ArrayList<Double>();
			tmpSpaceNbPoints = 0;
			double tmpSpot = 0.0;
			while(tmpSpot < spaceMaxForwardMultiple) {
				tmpSpot = spaceMinForwardMultiple + tmpSpaceStepSize * tmpSpaceNbPoints;
				tmpSpots.add(tmpSpot);
				tmpSpaceNbPoints++;
			}
		}

		this.spaceMinForwardMultiple = spaceMinForwardMultiple;
		this.spaceMaxForwardMultiple = tmpSpots.get(tmpSpots.size() - 1);
		this.spaceNbOfSteps = tmpSpaceNbPoints;
		spots = tmpSpots;
		spaceStepSize = tmpSpaceStepSize;
		spotIndex = lowerBound;
	}

	@Override
	public EquityValuationResult calculate(
			EquityValuationRequest request,
			EquityForwardStructure forwardStructure,
			YieldCurve discountCurve,
			VolatilitySurface volaSurface) {

		final Map<CalculationRequestType, Double> results = new HashMap<CalculationRequestType, Double>();
		if(request.getCalcsRequested().isEmpty()) {
			return new EquityValuationResult(request, results);
		}

		double price = 0.0;
		if(request.getCalcsRequested().contains(CalculationRequestType.EqDelta)
				|| request.getCalcsRequested().contains(CalculationRequestType.EqGamma)) {

			final double[] spotSensis = getPdeSensis(request.getOption(), forwardStructure, discountCurve, volaSurface);
			price = spotSensis[0];

			if(request.getCalcsRequested().contains(CalculationRequestType.EqDelta)) {
				results.put(CalculationRequestType.EqDelta, spotSensis[1]);
			}
			if(request.getCalcsRequested().contains(CalculationRequestType.EqGamma)) {
				results.put(CalculationRequestType.EqGamma, spotSensis[2]);
			}
		}
		else {
			price = getPrice(request.getOption(), forwardStructure, discountCurve, volaSurface);
		}

		if(request.getCalcsRequested().contains(CalculationRequestType.Price)) {
			results.put(CalculationRequestType.Price, price);
		}

		if(request.getCalcsRequested().contains(CalculationRequestType.EqVega)) {
			final double volShift = 0.0001; // TODO Make part of class members
			final double priceShifted = getPrice(
					request.getOption(),
					forwardStructure,
					discountCurve,
					volaSurface.getShiftedSurface(volShift));
			results.put(CalculationRequestType.EqVega, (priceShifted - price) / volShift);
		}

		return new EquityValuationResult(request, results);
	}

	public double getPrice(
			Option option,
			EquityForwardStructure forwardStructure,
			YieldCurve discountCurve,
			VolatilitySurface volSurface) {
		return evolvePde(option, forwardStructure, discountCurve, volSurface, false)[0];
	}

	public double[] getPdeSensis(
			Option option,
			EquityForwardStructure forwardStructure,
			YieldCurve discountCurve,
			VolatilitySurface volSurface) {
		return evolvePde(option, forwardStructure, discountCurve, volSurface, true);
	}

	public double getVega(
			Option option,
			EquityForwardStructure forwardStructure,
			YieldCurve discountCurve,
			VolatilitySurface volSurface,
			double basePrice,
			double volShift) {

		final double shiftedPrice = getPrice(
				option,
				forwardStructure,
				discountCurve,
				volSurface.getShiftedSurface(volShift));
		return (shiftedPrice - basePrice) / volShift;
	}

	public double getTheta(
			Option option,
			EquityForwardStructure forwardStructure,
			YieldCurve discountCurve,
			VolatilitySurface volSurface,
			double basePrice) {

		final LocalDate valDate = forwardStructure.getValuationDate();
		final LocalDate thetaDate = valDate.plusDays(1);
		final double thetaSpot = forwardStructure.getForward(thetaDate);
		final EquityForwardStructure shiftedFwdStructure = forwardStructure.cloneWithNewSpot(thetaSpot).cloneWithNewDate(thetaDate);
		final double shiftedPrice = getPrice(option, shiftedFwdStructure, discountCurve, volSurface);
		return (shiftedPrice - basePrice) / dayCounter.getDaycountFraction(valDate, thetaDate);
	}

	private double[] evolvePde(
			Option option,
			EquityForwardStructure forwardStructure,
			YieldCurve discountCurve,
			VolatilitySurface volSurface,
			boolean calculateSensis) {

		// Get data
		final LocalDate valDate = forwardStructure.getValuationDate();
		final LocalDate expiryDate = option.getExpiryDate();
		final double expiryTime = dayCounter.getDaycountFraction(valDate, expiryDate);

		assert !forwardStructure.getValuationDate().isAfter(expiryDate)
		: "Valuation date must not be after option expiry";

		final double impliedVol = volSurface.getVolatility(option.getStrike(), expiryDate, forwardStructure);
		double forward = forwardStructure.getForward(expiryDate);
		double fdf = forwardStructure.getFutureDividendFactor(expiryDate);

		// Build matrices
		final RealMatrix idMatrix = MatrixUtils.createRealIdentityMatrix(spaceNbOfSteps);
		final RealMatrix tridiagMatrix = MatrixUtils.createRealMatrix(spaceNbOfSteps, spaceNbOfSteps);
		final double spaceStepSq = spaceStepSize * spaceStepSize;

		for(int i = 0; i < spaceNbOfSteps; i++) {
			for(int j = 0; j < spaceNbOfSteps; j++) {
				if(i == j) {
					tridiagMatrix.setEntry(i, j, Math.pow(spots.get(i), 2) / spaceStepSq);
				}
				else if(i == j - 1 || i == j + 1) {
					tridiagMatrix.setEntry(i, j, -0.5 * Math.pow(spots.get(i), 2) / spaceStepSq);
				}
				else {
					tridiagMatrix.setEntry(i, j, 0);
				}
			}
		}

		// Set initial values
		RealVector prices = MatrixUtils.createRealVector(new double[spaceNbOfSteps]);
		for(int i = 0; i < spaceNbOfSteps; i++) {
			prices.setEntry(i, option.getPayoff((forward - fdf) * spots.get(i) + fdf));
		}

		// Set time intervals to evolve the PDE (i.e. from dividend to dividend)
		final List<LocalDate> diviDates = forwardStructure.getDividendStream().getDividendDates();
		final List<Double> anchorTimes = new ArrayList<Double>();
		anchorTimes.add(0.0);

		if(includeDividendDatesInGrid) {
			for(final LocalDate date : diviDates) {
				if(date.isAfter(valDate) && date.isBefore(expiryDate)) {
					anchorTimes.add(dayCounter.getDaycountFraction(valDate, date));
				}
			}
		}

		anchorTimes.add(expiryTime);
		anchorTimes.sort(Comparator.comparing(pt -> pt));

		double lastAtmPrice = 0.0;
		double dt = 0.0;

		// Evolve PDE
		for(int a = anchorTimes.size() - 1; a > 0; a--) {

			// Set time steps
			final double timeInterval = anchorTimes.get(a) - anchorTimes.get(a - 1);
			int timeNbOfSteps;
			double timeStepSize;

			if(timeStepsPerYear == 0) {
				// Use optimal ratio of time and space step size
				timeNbOfSteps = (int) Math.ceil(2 * impliedVol * Math.pow(timeInterval, 1.5) / spaceStepSize);
				timeStepSize = timeInterval / timeNbOfSteps;
			}
			else {
				// Use time step size provided externally
				timeNbOfSteps = (int) Math.floor(timeInterval * timeStepsPerYear);
				timeStepSize = timeInterval / timeNbOfSteps;
			}

			final List<Double> times = new ArrayList<Double>();
			for(int i = 0; i <= 4; i++) {
				times.add(anchorTimes.get(a) - i * 0.25 * timeStepSize);
			}
			for(int i = timeNbOfSteps - 2; i >= 0; i--) {
				times.add(anchorTimes.get(a - 1) + i * timeStepSize);
			}

			// Evolve PDE in current time interval
			for(int i = 1; i < times.size(); i++) {

				lastAtmPrice = prices.getEntry(spotIndex);
				dt = times.get(i - 1) - times.get(i);

				double theta = 0.5;
				if(i <= 4) {
					theta = 1.0;
				}
				final double theta1 = 1.0 - theta;
				final double volSq = impliedVol * impliedVol;

				RealMatrix implicitMatrix;
				RealMatrix explicitMatrix;

				if(isLvPricer) {
					implicitMatrix = tridiagMatrix.scalarMultiply(theta * dt);
					explicitMatrix = tridiagMatrix.scalarMultiply(-theta1 * dt);

					final double[] localVol = new double[spaceNbOfSteps];
					for(int s = 0; s < spaceNbOfSteps; s++) {
						final double lv = volSurface.getLocalVolatility(
								Math.log(spots.get(s)),
								times.get(i - 1),
								forwardStructure,
								spaceStepSize,
								dt);
						localVol[s] = lv * lv;
					}

					final RealMatrix volaMatrix = MatrixUtils.createRealDiagonalMatrix(localVol);
					implicitMatrix = volaMatrix.multiply(implicitMatrix);
					explicitMatrix = volaMatrix.multiply(explicitMatrix);
				}
				else {
					implicitMatrix = tridiagMatrix.scalarMultiply(theta * dt * volSq);
					explicitMatrix = tridiagMatrix.scalarMultiply(-theta1 * dt * volSq);
				}

				implicitMatrix = idMatrix.add(implicitMatrix);
				explicitMatrix = idMatrix.add(explicitMatrix);

				if(option.isAmericanOption()) {
					// Use the penalty algorithm from Forsyth's 2001 paper to solve the
					// linear complementary problem for the American exercise feature.
					final double penaltyFactor = 1 / Math.min(timeStepSize * timeStepSize, spaceStepSize * spaceStepSize);

					forward = forwardStructure.getForward(times.get(i));
					fdf = forwardStructure.getFutureDividendFactor(times.get(i));
					final double discountFactor = discountCurve.getForwardDiscountFactor(times.get(i), expiryTime);

					final RealVector payoffs = MatrixUtils.createRealVector(new double[spaceNbOfSteps]);
					final RealMatrix penaltyMatrix = MatrixUtils.createRealMatrix(spaceNbOfSteps, spaceNbOfSteps);

					for(int j = 1; j < spaceNbOfSteps - 1; j++) {
						final double payoff = option.getPayoff((forward - fdf) * spots.get(j) + fdf) / discountFactor;
						payoffs.setEntry(j, payoff);
						penaltyMatrix.setEntry(j, j, prices.getEntry(j) < payoff ? penaltyFactor : 0);
					}

					final RealVector b = explicitMatrix.operate(prices);
					RealVector oldPrices = prices.copy();
					final RealMatrix oldPenaltyMatrix = penaltyMatrix.copy();
					final double tol = 1 / penaltyFactor;

					int iterations = 0;
					while(true) {

						assert iterations++ < 100 : "Penalty algorithm for american exercise did not converge in 100 steps";

						final RealVector c = b.add(penaltyMatrix.operate(payoffs));
						final RealMatrix A = implicitMatrix.add(penaltyMatrix);
						final DecompositionSolver solver = new LUDecomposition(A).getSolver();

						prices = solver.solve(c);

						for(int j = 1; j < spaceNbOfSteps - 1; j++) {
							penaltyMatrix.setEntry(j, j, prices.getEntry(j) < payoffs.getEntry(j) ? penaltyFactor : 0);
						}

						if(penaltyMatrix.equals(oldPenaltyMatrix)
								|| (prices.subtract(oldPrices).getLInfNorm())
								/ Math.max(oldPrices.getLInfNorm(), 1.0) < tol) {
							break;
						}

						oldPrices = prices.copy();
					}
				}
				else {
					// Solve the PDE step directly
					prices = explicitMatrix.operate(prices);
					final DecompositionSolver solver = new LUDecomposition(implicitMatrix).getSolver();
					prices = solver.solve(prices);
				}

				// Set boundary conditions
				prices.setEntry(0, option.getPayoff((forward - fdf) * spaceMinForwardMultiple + fdf));
				prices.setEntry(spaceNbOfSteps - 1, option.getPayoff((forward - fdf) * spaceMaxForwardMultiple + fdf));
			}
		}

		final double discountFactor = discountCurve.getDiscountFactor(expiryDate);
		final double price = discountFactor * prices.getEntry(spotIndex);

		if(calculateSensis) {

			final double dFdX = forwardStructure.getDividendAdjustedStrike(
					forwardStructure.getForward(expiryDate),
					expiryDate);
			final double dFdS = forwardStructure.getGrowthDiscountFactor(valDate, expiryDate);

			final double delta = discountFactor
					* 0.5 * (prices.getEntry(spotIndex + 1) - prices.getEntry(spotIndex - 1))
					/ spaceStepSize * dFdS / dFdX;

			final double gamma = discountFactor
					* (prices.getEntry(spotIndex + 1) + prices.getEntry(spotIndex - 1) - 2 * prices.getEntry(spotIndex))
					/ spaceStepSq * dFdS * dFdS / dFdX / dFdX;

			final double discountFactorTheta = discountCurve.getDiscountFactor(expiryTime - dt);
			final double theta = (discountFactorTheta * lastAtmPrice - price) / dt;

			return new double[] { price, delta, gamma, theta };
		}
		else {
			return new double[] { price, Double.NaN, Double.NaN, Double.NaN };
		}
	}

	public double getImpliedVolatility(
			Option option,
			EquityForwardStructure forwardStructure,
			YieldCurve discountCurve,
			double price) {

		double initialGuess = 0.25;
		final double forward = forwardStructure.getForward(option.getExpiryDate());

		// Use analytic pricer as initial guess for Europeans and OTM Americans
		// Use two bisection steps for ITM Americans
		if(option.isAmericanOption() && option.getPayoff(forward) > 0.0) {

			final BisectionSearch bisectionSolver = new BisectionSearch(0.00001, 1.0);
			for(int i = 0; i < 3; i++) {
				final double currentVol = bisectionSolver.getNextPoint();
				final double currentPrice = getPrice(
						option,
						forwardStructure,
						discountCurve,
						new FlatVolatilitySurface(currentVol));

				bisectionSolver.setValue(currentPrice - price);
			}
			initialGuess = bisectionSolver.getBestPoint();
		}
		else {
			final AnalyticOptionValuation anaPricer = new AnalyticOptionValuation(dayCounter);
			Option testOption;

			if(option.isAmericanOption()) {
				testOption = new EuropeanOption(option.getExpiryDate(), option.getStrike(), option.isCallOption());
			}
			else {
				testOption = option;
			}

			initialGuess = anaPricer.getImpliedVolatility(testOption, forwardStructure, discountCurve, price);
		}

		// Solve for implied vol
		final SecantMethod solver = new SecantMethod(initialGuess, initialGuess * 1.01);
		while(solver.getAccuracy() / price > 1e-3 && !solver.isDone()) {

			final double currentVol = solver.getNextPoint();
			final double currentPrice = getPrice(
					option,
					forwardStructure,
					discountCurve,
					new FlatVolatilitySurface(currentVol));

			solver.setValue(currentPrice - price);
		}

		return Math.abs(solver.getBestPoint()); // Note that the PDE only uses sigma^2
	}
}