/**
* Copyright (C) 2012 - present by OpenGamma Inc. and the OpenGamma group of companies
*
* Please see distribution for license.
*/
package com.opengamma.analytics.financial.equity.variance.pricing;
import com.opengamma.analytics.financial.model.finitedifference.BoundaryCondition;
import com.opengamma.analytics.financial.model.finitedifference.ConvectionDiffusionPDE1DCoefficients;
import com.opengamma.analytics.financial.model.finitedifference.ConvectionDiffusionPDE1DStandardCoefficients;
import com.opengamma.analytics.financial.model.finitedifference.ConvectionDiffusionPDESolver;
import com.opengamma.analytics.financial.model.finitedifference.DirichletBoundaryCondition;
import com.opengamma.analytics.financial.model.finitedifference.ExponentialMeshing;
import com.opengamma.analytics.financial.model.finitedifference.HyperbolicMeshing;
import com.opengamma.analytics.financial.model.finitedifference.MeshingFunction;
import com.opengamma.analytics.financial.model.finitedifference.NeumannBoundaryCondition;
import com.opengamma.analytics.financial.model.finitedifference.PDE1DDataBundle;
import com.opengamma.analytics.financial.model.finitedifference.PDEGrid1D;
import com.opengamma.analytics.financial.model.finitedifference.PDETerminalResults1D;
import com.opengamma.analytics.financial.model.finitedifference.ThetaMethodFiniteDifference;
import com.opengamma.analytics.financial.model.finitedifference.applications.InitialConditionsProvider;
import com.opengamma.analytics.financial.model.finitedifference.applications.PDE1DCoefficientsProvider;
import com.opengamma.analytics.financial.model.interestrate.curve.ForwardCurve;
import com.opengamma.analytics.financial.model.interestrate.curve.YieldAndDiscountCurve;
import com.opengamma.analytics.financial.model.volatility.local.LocalVolatilitySurfaceMoneyness;
import com.opengamma.analytics.financial.model.volatility.local.PureLocalVolatilitySurface;
import com.opengamma.analytics.math.FunctionUtils;
import com.opengamma.analytics.math.function.Function1D;
import com.opengamma.util.ArgumentChecker;
/**
* Class to calculate the expected variance (<b>not</b> annualised) of an equity variance swap when a discount curve, affine dividends and a <b>pure</b> local volatility surface is
* specified. See White (2012), Equity Variance Swap with Dividends, for details of the model.
*/
public class EquityVarianceSwapForwardPurePDE {
// TODO The parameters here (in particular those of the meshing functions) have been chosen to minimise the errors in VarianceSwapWithDividendsTest, and will not be the
// best choice of parameters in other situations
/** Default weighting for averaging an explicit and implicit scheme */
private static final double DEFAULT_THETA = 0.5;
/** Sets up the PDE */
private static final PDE1DCoefficientsProvider PDE_PROVIDER = new PDE1DCoefficientsProvider();
/** Sets up the initial conditions */
private static final InitialConditionsProvider INITIAL_COND_PROVIDER = new InitialConditionsProvider();
/** The number of time steps */
private final int _nTimeSteps;
/** The number of space steps */
private final int _nSpaceSteps;
/** The minimum number of time steps between dividends */
private final int _minStepsBetweenDividends;
/** The PDE solver */
private final ConvectionDiffusionPDESolver _solver;
/**
* Sets up the calculator with the default values; 100 time steps, 100 space steps and a theta of 0.5 (i.e. Crank-Nicolson).
*/
public EquityVarianceSwapForwardPurePDE() {
_nTimeSteps = 100;
_nSpaceSteps = 100;
_minStepsBetweenDividends = 5;
_solver = new ThetaMethodFiniteDifference(DEFAULT_THETA, false);
}
/**
* @param nTimeSteps The number of time steps, greater than zero
* @param nSpaceSteps The number of space steps, greater than zero
* @param minStepsBetweenDividends The minimum number of steps to use between dividends, greater than zero
* @param theta The weight used in the finite difference scheme in the range 0 to 1 inclusive, where 0 is fully explicit,
* 1 is fully implicit, and 0.5 gives the Crank-Nicolson scheme.
*/
public EquityVarianceSwapForwardPurePDE(final int nTimeSteps, final int nSpaceSteps, final int minStepsBetweenDividends,
final double theta) {
ArgumentChecker.notNegativeOrZero(nTimeSteps, "number of time steps");
ArgumentChecker.notNegativeOrZero(nSpaceSteps, "number of space steps");
ArgumentChecker.notNegativeOrZero(minStepsBetweenDividends, "minimum number of steps between dividends");
ArgumentChecker.isInRangeInclusive(0, 1, theta);
_nTimeSteps = nTimeSteps;
_nSpaceSteps = nSpaceSteps;
_minStepsBetweenDividends = minStepsBetweenDividends;
_solver = new ThetaMethodFiniteDifference(theta, false);
}
/**
* Computes the expected variance by solving the forward PDE for the price of a call on the pure stock price, then using these computed prices at the expiry and all
* dividend dates, computes the expected variance with and without adjustments for the dividend payments.
* @param spot The current level of the stock or index, greater than zero
* @param discountCurve The risk free interest rate curve, not null
* @param dividends The dividends structure, not null
* @param expiry The expiry of the variance swap, greater than zero
* @param pureLocalVolSurface A <b>pure</b> local volatility surface, not null
* @return The expected variance with and without adjustments for the dividend payments (the former is usually the case for single stock and the latter for indices)
*/
public double[] expectedVariance(final double spot, final YieldAndDiscountCurve discountCurve, final AffineDividends dividends, final double expiry,
final PureLocalVolatilitySurface pureLocalVolSurface) {
ArgumentChecker.isTrue(spot > 0, "spot > 0");
ArgumentChecker.notNull(discountCurve, "discount curve");
ArgumentChecker.notNull(dividends, "dividends");
ArgumentChecker.isTrue(expiry > 0, "expiry > 0");
ArgumentChecker.notNull(pureLocalVolSurface, "pure local volatility surface");
//"convert" to a LocalVolatilitySurfaceMoneyness _ DO NOT interpret this as an actual LocalVolatilitySurfaceMoneyness
final LocalVolatilitySurfaceMoneyness localVolSurface = new LocalVolatilitySurfaceMoneyness(pureLocalVolSurface.getSurface(), new ForwardCurve(1.0));
final ConvectionDiffusionPDE1DStandardCoefficients pde = PDE_PROVIDER.getForwardLocalVol(localVolSurface);
final Function1D<Double, Double> initialCond = INITIAL_COND_PROVIDER.getForwardCallPut(true);
final EquityDividendsCurvesBundle divCurves = new EquityDividendsCurvesBundle(spot, discountCurve, dividends);
final double terminalFwd = divCurves.getF(expiry);
final double logNoDivFwd = Math.log(spot) + discountCurve.getInterestRate(expiry) * expiry;
final int n = dividends.getNumberOfDividends();
final double stepLength = expiry / _nTimeSteps;
int count = 0;
while (n > 0 && count < n && dividends.getTau(count) < expiry) {
count++;
}
final boolean divAtExp = n > 0 && count < n && Double.doubleToLongBits(dividends.getTau(count)) == Double.doubleToLongBits(expiry);
final int nDivsBeforeExpiry = count;
final int[] steps = new int[nDivsBeforeExpiry + 1];
int tSteps = 0;
for (int i = 0; i < nDivsBeforeExpiry; i++) {
steps[i] = (int) Math.max(_minStepsBetweenDividends, (dividends.getTau(i) - (i == 0 ? 0.0 : dividends.getTau(i - 1))) / stepLength);
tSteps += steps[i];
}
steps[nDivsBeforeExpiry] = Math.max(_minStepsBetweenDividends, _nTimeSteps - tSteps);
//common boundary conditions
//TODO the grid involves some magic numbers that should be possible to alter externally by expert users
final double xL = 0.0;
final double xH = 6;
final BoundaryCondition lower = new DirichletBoundaryCondition(1.0, xL);
final BoundaryCondition upper = new NeumannBoundaryCondition(0.0, xH, false);
final MeshingFunction spaceMesh = new HyperbolicMeshing(xL, xH, 1.0, _nSpaceSteps + 1, 0.001); //0.05
final MeshingFunction[] timeMeshes = new MeshingFunction[nDivsBeforeExpiry + 1];
final PDEGrid1D[] grids = new PDEGrid1D[nDivsBeforeExpiry + 1];
if (nDivsBeforeExpiry == 0) {
timeMeshes[0] = new ExponentialMeshing(0, expiry, _nTimeSteps, 5.0);
} else {
timeMeshes[0] = new ExponentialMeshing(0, dividends.getTau(0), steps[0], 5.0);
for (int i = 1; i < nDivsBeforeExpiry; i++) {
timeMeshes[i] = new ExponentialMeshing(dividends.getTau(i - 1), dividends.getTau(i), steps[i], 0.0);
}
timeMeshes[nDivsBeforeExpiry] = new ExponentialMeshing(dividends.getTau(nDivsBeforeExpiry - 1), expiry, steps[nDivsBeforeExpiry], 0.0);
}
final PDETerminalResults1D[] res = new PDETerminalResults1D[nDivsBeforeExpiry + 1];
grids[0] = new PDEGrid1D(timeMeshes[0], spaceMesh);
PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> db = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(pde, initialCond, lower, upper, grids[0]);
res[0] = (PDETerminalResults1D) _solver.solve(db);
for (int i = 1; i <= nDivsBeforeExpiry; i++) {
grids[i] = new PDEGrid1D(timeMeshes[i], spaceMesh);
db = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(pde, res[i - 1].getTerminalResults(), lower, upper, grids[i]);
res[i] = (PDETerminalResults1D) _solver.solve(db);
}
double corrDivAdj = 0;
double uncorrDivAdj = 0;
for (int i = 0; i < nDivsBeforeExpiry; i++) {
final double f = divCurves.getF(dividends.getTau(i));
corrDivAdj += integrate(getCorrectedDividendAdjustmentWeight(i, divCurves, dividends), res[i]) + getCorrectedDividendAdjustment(f, i, dividends);
uncorrDivAdj += integrate(getUncorrectedDividendAdjustmentWeight(i, divCurves, dividends), res[i]) + getUncorrectedDividendAdjustment(f, i, dividends);
}
if (divAtExp) {
corrDivAdj += integrate(getCorrectedDividendAdjustmentWeight(nDivsBeforeExpiry, divCurves, dividends), res[nDivsBeforeExpiry]) +
getCorrectedDividendAdjustment(terminalFwd, nDivsBeforeExpiry, dividends);
uncorrDivAdj += integrate(getUncorrectedDividendAdjustmentWeight(nDivsBeforeExpiry, divCurves, dividends), res[nDivsBeforeExpiry]) +
getUncorrectedDividendAdjustment(terminalFwd, nDivsBeforeExpiry, dividends);
}
final double logContract = integrate(getLogContractWeight(divCurves, expiry), res[nDivsBeforeExpiry]) + Math.log(terminalFwd);
final double rvNoDivs = -2 * (logContract - logNoDivFwd);
final double rvCorrDivs = rvNoDivs + 2 * corrDivAdj;
final double rvUncorrDivs = rvNoDivs + 2 * uncorrDivAdj;
return new double[] {rvCorrDivs, rvUncorrDivs };
}
private double integrate(final Function1D<Double, Double> weightFunc, final PDETerminalResults1D pdeRes) {
final int n = pdeRes.getNumberSpaceNodes();
final double[] x = pdeRes.getGrid().getSpaceNodes();
double sum = 0.0;
for (int i = 1; i < n - 1; i++) {
final double pureCallPrice = pdeRes.getFunctionValue(i);
final double otmPrice = pureCallPrice - (x[i] < 1.0 ? 1 - x[i] : 0);
final double w = weightFunc.evaluate(x[i]);
sum += w * otmPrice * (x[i + 1] - x[i - 1]);
}
sum /= 2.0;
//don't add on the end points;
return sum;
}
private Function1D<Double, Double> getLogContractWeight(final EquityDividendsCurvesBundle divCurves, final double expiry) {
final double f = divCurves.getF(expiry);
final double d = divCurves.getD(expiry);
return new Function1D<Double, Double>() {
@Override
public Double evaluate(final Double x) {
return -FunctionUtils.square((f - d) / ((f - d) * x + d));
}
};
}
private Function1D<Double, Double> getCorrectedDividendAdjustmentWeight(final int dividendIndex, final EquityDividendsCurvesBundle divCurves, final AffineDividends dividends) {
final double tau = dividends.getTau(dividendIndex);
final double f = divCurves.getF(tau);
final double d = divCurves.getD(tau);
final double alpha = dividends.getAlpha(dividendIndex);
final double fMd2 = FunctionUtils.square(f - d);
return new Function1D<Double, Double>() {
@Override
public Double evaluate(final Double x) {
final double s = (f - d) * x + d;
final double sPalpha = s + alpha;
final double ddH = -alpha * (s + sPalpha) / s / s / sPalpha / sPalpha;
final double weight = fMd2 * ddH;
return weight;
}
};
}
private Function1D<Double, Double> getUncorrectedDividendAdjustmentWeight(final int dividendIndex, final EquityDividendsCurvesBundle divCurves, final AffineDividends dividends) {
final double tau = dividends.getTau(dividendIndex);
final double f = divCurves.getF(tau);
final double d = divCurves.getD(tau);
final double alpha = dividends.getAlpha(dividendIndex);
final double beta = dividends.getBeta(dividendIndex);
final double fMd2 = FunctionUtils.square(f - d);
return new Function1D<Double, Double>() {
@Override
public Double evaluate(final Double x) {
final double s = (f - d) * x + d;
final double sPalpha = s + alpha;
final double h = Math.log(s * (1 - beta) / sPalpha);
final double dH = alpha / s / sPalpha;
final double ddH = -alpha * (s + sPalpha) / s / s / sPalpha / sPalpha;
final double weight = fMd2 * ((1 + h) * ddH + dH * dH);
return weight;
}
};
}
private double getCorrectedDividendAdjustment(final double s, final int dividendIndex, final AffineDividends dividends) {
final double alpha = dividends.getAlpha(dividendIndex);
final double beta = dividends.getBeta(dividendIndex);
final double h = Math.log(s * (1 - beta) / (s + alpha));
return h;
}
private double getUncorrectedDividendAdjustment(final double s, final int dividendIndex, final AffineDividends dividends) {
final double alpha = dividends.getAlpha(dividendIndex);
final double beta = dividends.getBeta(dividendIndex);
final double h = Math.log(s * (1 - beta) / (s + alpha));
return h + 0.5 * h * h;
}
}