/**
* Copyright (C) 2011 - present by OpenGamma Inc. and the OpenGamma group of companies
*
* Please see distribution for license.
*/
package com.opengamma.analytics.financial.interestrate.swaption.provider;
import java.util.Arrays;
import org.apache.commons.lang.Validate;
import com.opengamma.analytics.financial.interestrate.annuity.derivative.AnnuityPaymentFixed;
import com.opengamma.analytics.financial.interestrate.swaption.derivative.SwaptionBermudaFixedIbor;
import com.opengamma.analytics.financial.model.interestrate.HullWhiteOneFactorPiecewiseConstantInterestRateModel;
import com.opengamma.analytics.financial.model.interestrate.definition.HullWhiteOneFactorPiecewiseConstantParameters;
import com.opengamma.analytics.financial.provider.calculator.discounting.CashFlowEquivalentCalculator;
import com.opengamma.analytics.financial.provider.description.interestrate.HullWhiteOneFactorProviderInterface;
import com.opengamma.analytics.financial.provider.description.interestrate.MulticurveProviderInterface;
import com.opengamma.analytics.math.statistics.distribution.NormalDistribution;
import com.opengamma.analytics.math.statistics.distribution.ProbabilityDistribution;
import com.opengamma.util.ArgumentChecker;
import com.opengamma.util.money.Currency;
import com.opengamma.util.money.MultipleCurrencyAmount;
/**
* Method to compute the present value of Bermuda swaptions with the Hull-White one factor model by numerical integration.
* Reference: Henrard, M. Bermudan Swaptions in Gaussian HJM One-Factor Model: Analytical and Numerical Approaches. SSRN, October 2008. Available at SSRN: http://ssrn.com/abstract=1287982
*/
public final class SwaptionBermudaFixedIborHullWhiteNumericalIntegrationMethod {
/**
* The method unique instance.
*/
private static final SwaptionBermudaFixedIborHullWhiteNumericalIntegrationMethod INSTANCE = new SwaptionBermudaFixedIborHullWhiteNumericalIntegrationMethod();
/**
* Return the unique instance of the class.
* @return The instance.
*/
public static SwaptionBermudaFixedIborHullWhiteNumericalIntegrationMethod getInstance() {
return INSTANCE;
}
/**
* Private constructor.
*/
private SwaptionBermudaFixedIborHullWhiteNumericalIntegrationMethod() {
}
// TODO: The number of integration points should be a setting.
/**
* The number of points used in the numerical integration process.
*/
private static final int NB_POINT = 50;
/**
* The cash flow equivalent calculator used in computations.
*/
private static final CashFlowEquivalentCalculator CFEC = CashFlowEquivalentCalculator.getInstance();
/**
* The model used in computations.
*/
private static final HullWhiteOneFactorPiecewiseConstantInterestRateModel MODEL = new HullWhiteOneFactorPiecewiseConstantInterestRateModel();
/**
* The normal distribution implementation.
*/
private static final ProbabilityDistribution<Double> NORMAL = new NormalDistribution(0, 1);
/**
* Computes the present value of the Physical delivery swaption.
* @param swaption The swaption.
* @param hullWhite The Hull-White parameters and the curves.
* @return The present value.
*/
public MultipleCurrencyAmount presentValue(final SwaptionBermudaFixedIbor swaption, final HullWhiteOneFactorProviderInterface hullWhite) {
ArgumentChecker.notNull(swaption, "Swaption");
ArgumentChecker.notNull(hullWhite, "Hull-White provider");
MulticurveProviderInterface multicurves = hullWhite.getMulticurveProvider();
HullWhiteOneFactorPiecewiseConstantParameters parameters = hullWhite.getHullWhiteParameters();
Currency ccy = swaption.getCurrency();
final int nbExpiry = swaption.getExpiryTime().length;
Validate.isTrue(nbExpiry > 1, "At least two expiry dates required for this method");
double tmpdb;
final double[] theta = new double[nbExpiry + 1]; // Extended expiry time (with 0).
theta[0] = 0.0;
System.arraycopy(swaption.getExpiryTime(), 0, theta, 1, nbExpiry);
final AnnuityPaymentFixed[] cashflow = new AnnuityPaymentFixed[nbExpiry];
for (int loopexp = 0; loopexp < nbExpiry; loopexp++) {
cashflow[loopexp] = swaption.getUnderlyingSwap()[loopexp].accept(CFEC, multicurves);
}
final int[] n = new int[nbExpiry];
final double[][][] alpha = new double[nbExpiry][][];
final double[][][] alpha2 = new double[nbExpiry][][]; // alpha^2
for (int loopexp = 0; loopexp < nbExpiry; loopexp++) {
n[loopexp] = cashflow[loopexp].getNumberOfPayments();
alpha[loopexp] = new double[loopexp + 1][];
alpha2[loopexp] = new double[loopexp + 1][];
for (int k = 0; k <= loopexp; k++) {
alpha[loopexp][k] = new double[n[loopexp]];
alpha2[loopexp][k] = new double[n[loopexp]];
for (int l = 0; l < alpha[loopexp][k].length; l++) {
alpha[loopexp][k][l] = MODEL.alpha(parameters, theta[k], theta[k + 1], theta[k + 1], cashflow[loopexp].getNthPayment(l).getPaymentTime());
alpha2[loopexp][k][l] = alpha[loopexp][k][l] * alpha[loopexp][k][l];
}
}
}
final int nbPoint2 = 2 * NB_POINT + 1;
final int[] startInt = new int[nbExpiry - 1];
final int[] endInt = new int[nbExpiry - 1];
for (int i = 1; i < nbExpiry - 1; i++) {
startInt[i] = 0;
endInt[i] = nbPoint2 - 1;
}
startInt[0] = NB_POINT;
endInt[0] = NB_POINT;
final double[][] t = new double[nbExpiry][]; // payment time
final double[][] dfS = new double[nbExpiry][]; // discount factor
final double[] beta = new double[nbExpiry];
final double[][] h = new double[nbExpiry][];
final double[][] sa2 = new double[nbExpiry][];
for (int loopexp = 0; loopexp < nbExpiry; loopexp++) {
beta[loopexp] = MODEL.beta(parameters, theta[loopexp], theta[loopexp + 1]);
t[loopexp] = new double[n[loopexp]];
dfS[loopexp] = new double[n[loopexp]];
h[loopexp] = new double[n[loopexp]];
sa2[loopexp] = new double[n[loopexp]];
for (int loopcf = 0; loopcf < n[loopexp]; loopcf++) {
t[loopexp][loopcf] = cashflow[loopexp].getNthPayment(loopcf).getPaymentTime();
dfS[loopexp][loopcf] = multicurves.getDiscountFactor(ccy, t[loopexp][loopcf]);
h[loopexp][loopcf] = (1 - Math.exp(-parameters.getMeanReversion() * t[loopexp][loopcf])) / parameters.getMeanReversion();
tmpdb = 0.0;
for (int k = 0; k <= loopexp; k++) {
tmpdb += alpha2[loopexp][k][loopcf];
}
sa2[loopexp][loopcf] = tmpdb;
}
}
final double[] discountedCashFlowN = new double[n[nbExpiry - 1]];
for (int loopcf = 0; loopcf < n[nbExpiry - 1]; loopcf++) {
discountedCashFlowN[loopcf] = dfS[nbExpiry - 1][loopcf] * cashflow[nbExpiry - 1].getNthPayment(loopcf).getAmount();
}
final double lambda = MODEL.lambda(discountedCashFlowN, sa2[nbExpiry - 1], h[nbExpiry - 1]);
final double[] betaSort = new double[nbExpiry];
System.arraycopy(beta, 0, betaSort, 0, nbExpiry);
Arrays.sort(betaSort);
final double minbeta = betaSort[0];
final double maxbeta = betaSort[nbExpiry - 1];
final double b = Math.min(10 * minbeta, maxbeta);
final double epsilon = -2.0 / NB_POINT * NORMAL.getInverseCDF(1.0 / (200.0 * NB_POINT)) * b; // <-
final double[] bX = new double[nbPoint2];
for (int looppt = 0; looppt < nbPoint2; looppt++) {
bX[looppt] = -NB_POINT * epsilon + looppt * epsilon;
}
final double[] bX2 = new double[4 * NB_POINT + 1];
for (int looppt = 0; looppt < 4 * NB_POINT + 1; looppt++) {
bX2[looppt] = -2 * NB_POINT * epsilon + looppt * epsilon;
}
final double[] htheta = new double[nbExpiry];
for (int loopexp = 0; loopexp < nbExpiry; loopexp++) {
htheta[loopexp] = (1 - Math.exp(-parameters.getMeanReversion() * theta[loopexp + 1])) / parameters.getMeanReversion();
}
final double[][] vZ = new double[nbExpiry - 1][nbPoint2];
for (int i = nbExpiry - 2; i >= 0; i--) {
for (int looppt = 0; looppt < nbPoint2; looppt++) {
vZ[i][looppt] = Math.exp(bX[looppt] * htheta[i]);
}
}
final double[][] vW = new double[nbExpiry][]; // Swaption
final double[][] vT = new double[nbExpiry - 1][]; // Swap
final double omega = -Math.signum(cashflow[nbExpiry - 1].getNthPayment(0).getAmount());
final double[] kappaL = new double[nbPoint2];
for (int looppt = 0; looppt < nbPoint2; looppt++) {
kappaL[looppt] = (lambda - bX[looppt]) / beta[nbExpiry - 1];
}
final double[] sa2N1 = new double[n[nbExpiry - 1]];
for (int i = 0; i < n[nbExpiry - 1]; i++) {
tmpdb = 0;
for (int k = 0; k <= nbExpiry - 2; k++) {
tmpdb += alpha2[nbExpiry - 1][k][i];
}
sa2N1[i] = tmpdb;
}
vW[nbExpiry - 1] = new double[4 * NB_POINT + 1];
for (int j = 0; j < n[nbExpiry - 1]; j++) {
for (int looppt = 0; looppt < nbPoint2; looppt++) {
vW[nbExpiry - 1][NB_POINT + looppt] += discountedCashFlowN[j] * Math.exp(-sa2N1[j] / 2.0 - h[nbExpiry - 1][j] * bX[looppt])
* NORMAL.getCDF(omega * (kappaL[looppt] + alpha[nbExpiry - 1][nbExpiry - 1][j]));
}
}
for (int looppt = 0; looppt < NB_POINT; looppt++) {
vW[nbExpiry - 1][looppt] = vW[nbExpiry - 1][NB_POINT];
}
for (int looppt = 0; looppt < NB_POINT; looppt++) {
vW[nbExpiry - 1][3 * NB_POINT + 1 + looppt] = vW[nbExpiry - 1][3 * NB_POINT];
}
final double c1sqrt2pi = 1.0 / Math.sqrt(2 * Math.PI);
final double[][] pvcfT = new double[nbExpiry - 1][];
double[] vL; // Left side of intersection
double[] vR; // Right side of intersection
double[][] labc;
double[][] rabc;
final double[][] labcM = new double[3][4 * NB_POINT + 1];
final double[][] rabcM = new double[3][4 * NB_POINT + 1];
final double[] dabc = new double[3];
final int[] indSwap = new int[nbExpiry - 1]; // index of the intersection
double xroot;
final double[][] xN = new double[nbExpiry - 1][nbPoint2];
double ci;
double coi;
int is;
final double[] ncdf0 = new double[nbPoint2];
final double[] ncdf1 = new double[nbPoint2];
final double[] ncdf2 = new double[nbPoint2];
final double[] ncdf0X = new double[nbPoint2 + 1];
final double[] ncdf1X = new double[nbPoint2 + 1];
final double[] ncdf2X = new double[nbPoint2 + 1];
double ncdf0x;
double ncdf1x;
double ncdf2x;
double ncdfinit;
// Main loop for the different expiry dates (except the last one)
for (int i = nbExpiry - 2; i >= 0; i--) {
vW[i] = new double[4 * NB_POINT + 1];
vT[i] = new double[4 * NB_POINT + 1];
// T: swap
pvcfT[i] = new double[n[i]];
for (int j = 0; j < n[i]; j++) {
pvcfT[i][j] = cashflow[i].getNthPayment(j).getAmount() * dfS[i][j];
for (int looppt = 0; looppt < 4 * NB_POINT + 1; looppt++) {
vT[i][looppt] += pvcfT[i][j] * Math.exp(-sa2[i][j] / 2.0 - h[i][j] * bX2[looppt]);
}
}
// Preparation
for (int looppt = 0; looppt < nbPoint2; looppt++) {
xN[i][looppt] = bX[looppt] / beta[i];
}
ci = htheta[i] * beta[i];
coi = Math.exp(ci * ci / 2);
// Left/Right
if (omega < 0) {
vL = vW[i + 1];
vR = vT[i];
} else {
vR = vW[i + 1];
vL = vT[i];
}
indSwap[i] = 0;
while (vL[indSwap[i] + 1] >= vR[indSwap[i] + 1]) {
indSwap[i]++;
}
// Parabola fit
labc = parafit(epsilon / beta[i], vL);
rabc = parafit(epsilon / beta[i], vR);
for (int k = 0; k < 3; k++) {
dabc[k] = labc[k][indSwap[i]] - rabc[k][indSwap[i]];
System.arraycopy(labc[k], 0, labcM[k], 0, indSwap[i] + 1);
System.arraycopy(labc[k], indSwap[i], labcM[k], indSwap[i] + 1, labc[k].length - indSwap[i]);
System.arraycopy(rabc[k], 0, rabcM[k], 0, indSwap[i] + 1);
System.arraycopy(rabc[k], indSwap[i], rabcM[k], indSwap[i] + 1, rabc[k].length - indSwap[i]);
}
for (int looppt = 0; looppt < 4 * NB_POINT + 1; looppt++) {
labcM[1][looppt] = labcM[1][looppt] + labcM[0][looppt] * 2 * ci;
labcM[2][looppt] = labcM[2][looppt] + labcM[1][looppt] * ci - labcM[0][looppt] * ci * ci;
rabcM[1][looppt] = rabcM[1][looppt] + rabcM[0][looppt] * 2 * ci;
rabcM[2][looppt] = rabcM[2][looppt] + rabcM[1][looppt] * ci - rabcM[0][looppt] * ci * ci;
}
xroot = (-dabc[1] - Math.sqrt(dabc[1] * dabc[1] - 4 * dabc[0] * dabc[2])) / (2 * dabc[0]);
ncdfinit = NORMAL.getCDF(xN[i][0]);
for (int looppt = 0; looppt < nbPoint2; looppt++) {
ncdf0[looppt] = NORMAL.getCDF(xN[i][looppt] - ci);
ncdf1[looppt] = -c1sqrt2pi * Math.exp(-(xN[i][looppt] - ci) * (xN[i][looppt] - ci) / 2.0);
ncdf2[looppt] = ncdf1[looppt] * (xN[i][looppt] - ci) + ncdf0[looppt];
}
for (int j = startInt[i]; j <= endInt[i]; j++) {
is = indSwap[i] - j + 1;
// % all L
if (j + 2 * NB_POINT <= indSwap[i]) {
final double[][] xabc = new double[3][2 * NB_POINT];
for (int k = 0; k < 3; k++) {
System.arraycopy(labcM[k], j, xabc[k], 0, 2 * NB_POINT);
}
for (int looppt = 0; looppt < 2 * NB_POINT; looppt++) {
xabc[1][looppt] = xabc[1][looppt] + xabc[0][looppt] * 2 * xN[i][j];
xabc[2][looppt] = xabc[2][looppt] + xabc[1][looppt] * xN[i][j] - xabc[0][looppt] * xN[i][j] * xN[i][j];
}
vW[i][j + NB_POINT] = 0;
vW[i][j + NB_POINT] = vW[i][j + NB_POINT] + coi * ni2ncdf(ncdf2, ncdf1, ncdf0, xabc);
} else if (j < indSwap[i]) {
final double[][] xabc = new double[3][2 * NB_POINT + 1];
tmpdb = xroot - xN[i][j] - ci;
ncdf0x = NORMAL.getCDF(tmpdb);
ncdf1x = -Math.exp(-(tmpdb * tmpdb) / 2) * c1sqrt2pi;
ncdf2x = ncdf1x * tmpdb + ncdf0x;
for (int k = 0; k < 3; k++) {
// System.arraycopy(rabcM[k], j, xabc[k], 0, 2 * _nbPoint + 1); // Swap
System.arraycopy(labcM[k], j, xabc[k], 0, 2 * NB_POINT + 1);
System.arraycopy(rabcM[k], indSwap[i] + 1, xabc[k], indSwap[i] + 1 - j, j + 2 * NB_POINT - indSwap[i]);
}
for (int looppt = 0; looppt < 2 * NB_POINT; looppt++) {
xabc[1][looppt] = xabc[1][looppt] + xabc[0][looppt] * 2 * xN[i][j];
xabc[2][looppt] = xabc[2][looppt] + xabc[1][looppt] * xN[i][j] - xabc[0][looppt] * xN[i][j] * xN[i][j];
}
System.arraycopy(ncdf0, 0, ncdf0X, 0, is);
ncdf0X[is] = ncdf0x;
System.arraycopy(ncdf0, is, ncdf0X, is + 1, ncdf0.length - is);
System.arraycopy(ncdf1, 0, ncdf1X, 0, is);
ncdf1X[is] = ncdf1x;
System.arraycopy(ncdf1, is, ncdf1X, is + 1, ncdf1.length - is);
System.arraycopy(ncdf2, 0, ncdf2X, 0, is);
ncdf2X[is] = ncdf2x;
System.arraycopy(ncdf2, is, ncdf2X, is + 1, ncdf2.length - is);
vW[i][j + NB_POINT] = vL[j] * vZ[i][0] * ncdfinit;
vW[i][j + NB_POINT] += coi * ni2ncdf(ncdf2X, ncdf1X, ncdf0X, xabc);
vW[i][j + NB_POINT] += vR[j + 2 * NB_POINT] * vZ[i][vZ[i].length - 1] * ncdfinit;
} else {
final double[][] xabc = new double[3][2 * NB_POINT];
for (int k = 0; k < 3; k++) {
System.arraycopy(rabcM[k], j + 1, xabc[k], 0, 2 * NB_POINT);
// System.arraycopy(labcM[k], j + 1, xabc[k], 0, 2 * _nbPoint); // Swaption
}
for (int looppt = 0; looppt < 2 * NB_POINT; looppt++) {
xabc[1][looppt] = xabc[1][looppt] + xabc[0][looppt] * 2 * xN[i][j];
xabc[2][looppt] = xabc[2][looppt] + xabc[1][looppt] * xN[i][j] - xabc[0][looppt] * xN[i][j] * xN[i][j];
}
vW[i][j + NB_POINT] = vR[j] * vZ[i][0] * ncdfinit;
vW[i][j + NB_POINT] += coi * ni2ncdf(ncdf2, ncdf1, ncdf0, xabc);
vW[i][j + NB_POINT] += vR[j + 2 * NB_POINT] * vZ[i][vZ[i].length - 1] * ncdfinit;
}
}
for (int j = 0; j < NB_POINT; j++) { // Flat extrapolation
vW[i][j] = vW[i][NB_POINT];
vW[i][3 * NB_POINT + 1 + j] = vW[i][3 * NB_POINT];
}
} // End main loop
return MultipleCurrencyAmount.of(swaption.getUnderlyingSwap()[0].getFixedLeg().getCurrency(), vW[0][2 * NB_POINT] * (swaption.isLong() ? 1.0 : -1.0));
}
/**
* Fit the parabolas.
* @param dx Distance between the x values.
* @param y The y values.
* @return The parabolas coefficients.
*/
private static double[][] parafit(final double dx, final double[] y) {
final int nbPts = y.length;
final int nbStep = (nbPts - 1) / 2;
final double dx2 = dx * dx;
final double[] x1 = new double[nbStep];
final double[] x = new double[nbStep];
for (int i = 0; i < nbStep; i++) {
x1[i] = -nbStep + 2.0 * i;
x[i] = x1[i] * dx;
}
final double[][] abc = new double[3][2 * nbStep];
double tmp;
for (int i = 0; i < nbStep; i++) {
tmp = (y[2 + 2 * i] - 2 * y[1 + 2 * i] + y[2 * i]) / (2 * dx2);
abc[0][2 * i] = tmp;
abc[0][2 * i + 1] = tmp;
tmp = (y[1 + 2 * i] - y[2 * i]) / dx - (2 * x[i] + dx) * abc[0][2 * i];
abc[1][2 * i] = tmp;
abc[1][2 * i + 1] = tmp;
tmp = y[2 * i] - x[i] * x[i] * abc[0][2 * i] - x[i] * abc[1][2 * i];
abc[2][2 * i] = tmp;
abc[2][2 * i + 1] = tmp;
}
return abc;
}
/**
* Numerical integration of the parabolas against the normal distribution.
* @param n2 Second order integrals.
* @param n1 First order integrals.
* @param n0 Order 0 integrals.
* @param abc The parabolas coefficients.
* @return The integral.
*/
private static double ni2ncdf(final double[] n2, final double[] n1, final double[] n0, final double[][] abc) {
double s = 0;
for (int i = 0; i < abc[0].length; i++) {
s += abc[0][i] * (n2[i + 1] - n2[i]);
s += abc[1][i] * (n1[i + 1] - n1[i]);
s += abc[2][i] * (n0[i + 1] - n0[i]);
}
return s;
}
}