Java tutorial
/** * Copyright (C) 2013 - present by OpenGamma Inc. and the OpenGamma group of companies * * Please see distribution for license. */ package com.opengamma.analytics.math.interpolation; import java.util.ArrayList; import java.util.Arrays; import org.apache.commons.lang.NotImplementedException; import com.opengamma.analytics.math.interpolation.data.Interpolator1DPiecewisePoynomialWithExtraKnotsDataBundle; import com.opengamma.analytics.math.matrix.DoubleMatrix1D; import com.opengamma.analytics.math.matrix.DoubleMatrix2D; import com.opengamma.util.ArgumentChecker; import com.opengamma.util.ParallelArrayBinarySort; /** * Monotone Convex Interpolation based on * P. S. Hagan, G. West, "interpolation Methods for Curve Construction" * Applied Mathematical Finance, Vol. 13, No. 2, 89129, June 2006 * * Given a data set (time) and (spot rates)*(time), {t_i, r_i*t_i}, derive forward rate curve, f(t)=\frac{\partial r(t) t}{\partial t}, by "interpolateFwds" method * or derive a curve of (spot rates) * (time), r(t) * t, by "interpolate" method by applying the interpolation to forward rates f_i estimated from spot rates r_i * When we apply this spline to interest rates, DO INCLUDE THE TRIVIAL POINT (t_0,r_0*t_0) = (0,0). In this case r_0 = 0 is automatically assumed. * Note that f(t_i) = (original) f_i does NOT necessarily hold due to forward modification for ensuring positivity of the curve */ public class MonotoneConvexSplineInterpolator extends PiecewisePolynomialInterpolator { private double[] _time; private double[] _spotRates; /** * */ public MonotoneConvexSplineInterpolator() { _time = null; _spotRates = null; } /** * Determine r(t)t = \int _{xValues_0}^{x} f(s) ds for t >= min{xValues} * Extrapolation by a linear function in the region t > max{xValues}. To employ this extrapolation, use interpolate methods in this class. * @param xValues Data t_i * @param yValues Data r_i*t_i * @return PiecewisePolynomialResult for r(t)t */ @Override public PiecewisePolynomialResult interpolate(final double[] xValues, final double[] yValues) { ArgumentChecker.notNull(xValues, "xValues"); ArgumentChecker.notNull(yValues, "yValues"); ArgumentChecker.isTrue(xValues.length == yValues.length, " xValues length = yValues length"); ArgumentChecker.isTrue(xValues.length > 1, "Data points should be more than 1"); final int nDataPts = xValues.length; for (int i = 0; i < nDataPts; ++i) { ArgumentChecker.isFalse(Double.isNaN(xValues[i]), "xData containing NaN"); ArgumentChecker.isFalse(Double.isInfinite(xValues[i]), "xData containing Infinity"); ArgumentChecker.isFalse(Double.isNaN(yValues[i]), "yData containing NaN"); ArgumentChecker.isFalse(Double.isInfinite(yValues[i]), "yData containing Infinity"); } for (int i = 0; i < nDataPts; ++i) { for (int j = i + 1; j < nDataPts; ++j) { ArgumentChecker.isFalse(xValues[i] == xValues[j], "xValues should be distinct"); } } for (int i = 0; i < nDataPts; ++i) { if (xValues[i] == 0.) { ArgumentChecker.isTrue(yValues[i] == 0., "r_i * t_i = 0 if t_i =0"); } } double[] spotTmp = new double[nDataPts]; for (int i = 0; i < nDataPts; ++i) { spotTmp[i] = xValues[i] == 0. ? 0. : yValues[i] / xValues[i]; } _time = Arrays.copyOf(xValues, nDataPts); _spotRates = Arrays.copyOf(spotTmp, nDataPts); ParallelArrayBinarySort.parallelBinarySort(_time, _spotRates); final DoubleMatrix2D coefMatrix = solve(_time, _spotRates); final DoubleMatrix2D coefMatrixIntegrate = integration(_time, coefMatrix.getData()); for (int i = 0; i < coefMatrixIntegrate.getNumberOfRows(); ++i) { for (int j = 0; j < coefMatrixIntegrate.getNumberOfColumns(); ++j) { ArgumentChecker.isFalse(Double.isNaN(coefMatrixIntegrate.getData()[i][j]), "Too large input"); ArgumentChecker.isFalse(Double.isInfinite(coefMatrixIntegrate.getData()[i][j]), "Too large input"); } } return new PiecewisePolynomialResult(new DoubleMatrix1D(_time), coefMatrixIntegrate, coefMatrixIntegrate.getNumberOfColumns(), 1); } /** * Since Monotone Convex spline method introduces extra knots in some cases and the number of knots depends on yValues, * this multidimensional method can not be supported * @param xValues * @param yValuesMatrix Multidimensional y values * @return Error is returned */ @Override public PiecewisePolynomialResult interpolate(final double[] xValues, final double[][] yValuesMatrix) { throw new IllegalArgumentException("Method with multidimensional yValues is not supported"); } @Override public double interpolate(final double[] xValues, final double[] yValues, final double x) { final PiecewisePolynomialResult result = interpolate(xValues, yValues); final DoubleMatrix2D coefsMatrixIntegrate = result.getCoefMatrix(); final int nKnots = coefsMatrixIntegrate.getNumberOfRows() + 1; final double[] knots = result.getKnots().getData(); int indicator = 0; if (x <= knots[1]) { indicator = 0; } else { for (int i = 1; i < nKnots - 1; ++i) { if (knots[i] < x) { indicator = i; } } } final double[] coefs = coefsMatrixIntegrate.getRowVector(indicator).getData(); final double res = getValue(coefs, x, knots[indicator]); ArgumentChecker.isFalse(Double.isInfinite(res), "Too large/small data values or xKey"); ArgumentChecker.isFalse(Double.isNaN(res), "Too large/small data values or xKey"); return res; } @Override public DoubleMatrix1D interpolate(final double[] xValues, final double[] yValues, final double[] x) { ArgumentChecker.notNull(x, "x"); final int keyLength = x.length; double[] res = new double[keyLength]; final PiecewisePolynomialResult result = interpolate(xValues, yValues); final DoubleMatrix2D coefsMatrixIntegrate = result.getCoefMatrix(); final int nKnots = coefsMatrixIntegrate.getNumberOfRows() + 1; final double[] knots = result.getKnots().getData(); for (int j = 0; j < keyLength; ++j) { int indicator = 0; if (x[j] <= knots[1]) { indicator = 0; } else { for (int i = 1; i < nKnots - 1; ++i) { if (knots[i] < x[j]) { indicator = i; } } } final double[] coefs = coefsMatrixIntegrate.getRowVector(indicator).getData(); res[j] = getValue(coefs, x[j], knots[indicator]); ArgumentChecker.isFalse(Double.isInfinite(res[j]), "Too large/small data values or xKey"); ArgumentChecker.isFalse(Double.isNaN(res[j]), "Too large/small data values or xKey"); } return new DoubleMatrix1D(res); } @Override public DoubleMatrix2D interpolate(final double[] xValues, final double[] yValues, final double[][] xMatrix) { ArgumentChecker.notNull(xMatrix, "xMatrix"); final int keyLength = xMatrix[0].length; final int keyDim = xMatrix.length; double[][] res = new double[keyDim][keyLength]; final PiecewisePolynomialResult result = interpolate(xValues, yValues); final DoubleMatrix2D coefsMatrixIntegrate = result.getCoefMatrix(); final int nKnots = coefsMatrixIntegrate.getNumberOfRows() + 1; final double[] knots = result.getKnots().getData(); for (int j = 0; j < keyDim; ++j) { for (int k = 0; k < keyLength; ++k) { int indicator = 0; if (xMatrix[j][k] <= knots[1]) { indicator = 0; } else { for (int i = 1; i < nKnots - 1; ++i) { if (knots[i] < xMatrix[j][k]) { indicator = i; } } } final double[] coefs = coefsMatrixIntegrate.getRowVector(indicator).getData(); res[j][k] = getValue(coefs, xMatrix[j][k], knots[indicator]); ArgumentChecker.isFalse(Double.isInfinite(res[j][k]), "Too large input"); ArgumentChecker.isFalse(Double.isNaN(res[j][k]), "Too large input"); } } return new DoubleMatrix2D(res); } /** * Since this interpolation method introduces new breakpoints in certain cases, {@link PiecewisePolynomialResultsWithSensitivity} is not well-defined * Instead the node sensitivity is computed in {@link MonotoneConvexSplineInterpolator1D} via {@link Interpolator1DPiecewisePoynomialWithExtraKnotsDataBundle} * @param xValues * @param yValues * @return NotImplementedException */ @Override public PiecewisePolynomialResultsWithSensitivity interpolateWithSensitivity(final double[] xValues, final double[] yValues) { throw new NotImplementedException(); } /** * Determine f(t) = \frac{\partial r(t) t}{\partial t} * @param xValues Data t_i * @param yValues Data r(t_i) * @return PiecewisePolynomialResult for f(t) */ public PiecewisePolynomialResult interpolateFwds(final double[] xValues, final double[] yValues) { ArgumentChecker.notNull(xValues, "xValues"); ArgumentChecker.notNull(yValues, "yValues"); ArgumentChecker.isTrue(xValues.length == yValues.length, " xValues length = yValues length"); ArgumentChecker.isTrue(xValues.length > 1, "Data points should be more than 1"); final int nDataPts = xValues.length; for (int i = 0; i < nDataPts; ++i) { ArgumentChecker.isFalse(Double.isNaN(xValues[i]), "xData containing NaN"); ArgumentChecker.isFalse(Double.isInfinite(xValues[i]), "xData containing Infinity"); ArgumentChecker.isFalse(Double.isNaN(yValues[i]), "yData containing NaN"); ArgumentChecker.isFalse(Double.isInfinite(yValues[i]), "yData containing Infinity"); } for (int i = 0; i < nDataPts; ++i) { for (int j = i + 1; j < nDataPts; ++j) { ArgumentChecker.isFalse(xValues[i] == xValues[j], "xValues should be distinct"); } } for (int i = 0; i < nDataPts; ++i) { if (xValues[i] == 0.) { ArgumentChecker.isTrue(yValues[i] == 0., "r_i * t_i = 0 if t_i =0"); } } double[] spotTmp = new double[nDataPts]; for (int i = 0; i < nDataPts; ++i) { spotTmp[i] = xValues[i] == 0. ? 0. : yValues[i] / xValues[i]; } _time = Arrays.copyOf(xValues, nDataPts); _spotRates = Arrays.copyOf(spotTmp, nDataPts); ParallelArrayBinarySort.parallelBinarySort(_time, _spotRates); final DoubleMatrix2D coefMatrix = solve(_time, _spotRates); for (int i = 0; i < coefMatrix.getNumberOfRows(); ++i) { for (int j = 0; j < coefMatrix.getNumberOfColumns(); ++j) { ArgumentChecker.isFalse(Double.isNaN(coefMatrix.getData()[i][j]), "Too large input"); ArgumentChecker.isFalse(Double.isInfinite(coefMatrix.getData()[i][j]), "Too large input"); } } return new PiecewisePolynomialResult(new DoubleMatrix1D(_time), coefMatrix, coefMatrix.getNumberOfColumns(), 1); } /** * Compute the value of f(t) = \frac{\partial r(t) t}{\partial t} at t=x * @param xValues Data r_i * @param yValues Data r(t_i) * @param x Key larger than min{r_i} * @return f(x) */ public double interpolateFwds(final double[] xValues, final double[] yValues, final double x) { final PiecewisePolynomialResult result = interpolateFwds(xValues, yValues); final DoubleMatrix2D coefsMatrix = result.getCoefMatrix(); final int nKnots = coefsMatrix.getNumberOfRows() + 1; final double[] knots = result.getKnots().getData(); int indicator = 0; if (x <= knots[1]) { indicator = 0; } else { for (int i = 1; i < nKnots - 1; ++i) { if (knots[i] < x) { indicator = i; } } } final double[] coefs = coefsMatrix.getRowVector(indicator).getData(); final double res = getValue(coefs, x, knots[indicator]); ArgumentChecker.isFalse(Double.isInfinite(res), "Too large/small data values or xKey"); ArgumentChecker.isFalse(Double.isNaN(res), "Too large/small data values or xKey"); return res; } /** * Compute the values of f(t) * @param xValues Data t_i * @param yValues Data r(t_i) * @param x Set of xKey * @return r(x) */ public DoubleMatrix1D interpolateFwds(final double[] xValues, final double[] yValues, final double[] x) { ArgumentChecker.notNull(x, "x"); final PiecewisePolynomialResult result = interpolateFwds(xValues, yValues); final DoubleMatrix2D coefsMatrix = result.getCoefMatrix(); final int nKnots = coefsMatrix.getNumberOfRows() + 1; final double[] knots = result.getKnots().getData(); final int keyLength = x.length; double[] res = new double[keyLength]; for (int j = 0; j < keyLength; ++j) { int indicator = 0; if (x[j] <= knots[1]) { indicator = 0; } else { for (int i = 1; i < nKnots - 1; ++i) { if (knots[i] < x[j]) { indicator = i; } } } final double[] coefs = coefsMatrix.getRowVector(indicator).getData(); res[j] = getValue(coefs, x[j], knots[indicator]); } return new DoubleMatrix1D(res); } /** * Compute the values of f(t) * @param xValues Data t_i * @param yValues Data r(t_i) * @param xMatrix Set of xKey * @return r(x) */ public DoubleMatrix2D interpolateFwds(final double[] xValues, final double[] yValues, final double[][] xMatrix) { ArgumentChecker.notNull(xMatrix, "xMatrix"); final int keyLength = xMatrix[0].length; final int keyDim = xMatrix.length; double[][] res = new double[keyDim][keyLength]; for (int i = 0; i < keyDim; ++i) { res[i] = interpolateFwds(xValues, yValues, xMatrix[i]).getData(); } return new DoubleMatrix2D(res); } /** * Derive r(t) * t from f(t) * @param knots * @param coefMatrix * @return coefmatrix of r(t) * t */ private DoubleMatrix2D integration(final double[] knots, final double[][] coefMatrix) { final int nCoefs = coefMatrix[0].length + 1; final int nKnots = knots.length; double[][] res = new double[nKnots][nCoefs]; double sum = _spotRates[0] * _time[0]; for (int i = 0; i < nKnots - 1; ++i) { res[i][0] = coefMatrix[i][0] / 3.; res[i][1] = coefMatrix[i][1] / 2.; res[i][2] = coefMatrix[i][2]; res[i][3] = sum; sum = getValue(res[i], knots[i + 1], knots[i]); } res[nKnots - 1][0] = 0.; res[nKnots - 1][1] = 0.; res[nKnots - 1][2] = getValue(coefMatrix[nKnots - 2], knots[nKnots - 1], knots[nKnots - 2]); res[nKnots - 1][3] = sum; return new DoubleMatrix2D(res); } /** * @param xValues X values of data * @param yValues Y values of data * @return Coefficient matrix whose i-th row vector is {a3, a2, a1, a0} of f(x) = a3 * (x-x_i)^3 + a2 * (x-x_i)^2 +... for the i-th interval */ private DoubleMatrix2D solve(final double[] time, final double[] spotRates) { final int nDataPts = time.length; final double[] discFwds = discFwdsFinder(time, spotRates); double[] fwds = fwdsFinder(time, discFwds); ArrayList<double[]> coefsList = new ArrayList<>(); ArrayList<Double> knots = new ArrayList<>(); for (int i = 0; i < nDataPts - 1; ++i) { final double gValue0 = fwds[i] - discFwds[i]; final double gValue1 = fwds[i + 1] - discFwds[i]; final double gDiff0 = -4. * gValue0 - 2. * gValue1; final double gDiff1 = 2. * gValue0 + 4. * gValue1; final double interval = time[i + 1] - time[i]; final double shift = discFwds[i]; if (Math.abs(gValue0) <= 1e-13) { final double[] coefs = new double[] { 0., 0., shift, }; knots.add(time[i]); coefsList.add(coefs); } else { if (Math.abs(gValue1) <= 1e-13) { final double[] coefs = new double[] { 0., 0., shift, }; knots.add(time[i]); coefsList.add(coefs); } else { if ((gValue0 > 0. && gValue1 > 0.) | (gValue0 < 0. && gValue1 < 0.)) { final double eta = gValue1 / (gValue1 + gValue0); final double cst0 = (gValue0 + 2. * gValue1) * gValue0 / gValue1; final double cst1 = (gValue1 + 2. * gValue0) * gValue1 / gValue0; final double newKnot = time[i] + interval * eta; final double[] coefs1 = new double[] { cst0 / gValue1 * (gValue0 + gValue1) / interval / interval, -2. * cst0 / interval, gValue0 + shift }; final double[] coefs2 = new double[] { cst1 * (gValue0 + gValue1) / gValue0 / interval / interval, 0., -gValue0 * gValue1 / (gValue1 + gValue0) + shift }; knots.add(time[i]); knots.add(newKnot); coefsList.add(coefs1); coefsList.add(coefs2); } else { if ((gDiff0 >= 0. && gDiff1 >= 0.) | (gDiff0 <= 0. && gDiff1 <= 0.)) { final double[] coefs = new double[] { (3. * gValue0 + 3. * gValue1) / interval / interval, (-4. * gValue0 - 2. * gValue1) / interval, gValue0 + shift }; knots.add(time[i]); coefsList.add(coefs); } else { if ((gValue0 < 0. && gValue1 >= -2. * gValue0) | (gValue0 > 0. && gValue1 <= -2. * gValue0)) { final double eta = (gValue1 + 2. * gValue0) / (gValue1 - gValue0); final double newKnot = time[i] + interval * eta; final double cst = (gValue1 - gValue0) / 3. / gValue0; final double[] coefs1 = new double[] { 0., 0., gValue0 + shift }; final double[] coefs2 = new double[] { cst * cst * (gValue1 - gValue0) / interval / interval, 0., gValue0 + shift }; knots.add(time[i]); knots.add(newKnot); coefsList.add(coefs1); coefsList.add(coefs2); } else { final double eta = 3. * gValue1 / (gValue1 - gValue0); final double newKnot = time[i] + interval * eta; final double cst = (gValue0 - gValue1) / 3. / gValue1; final double[] coefs1 = new double[] { cst * cst * (gValue0 - gValue1) / interval / interval, 2. * cst * (gValue0 - gValue1) / interval, gValue0 + shift }; final double[] coefs2 = { 0, 0, gValue1 + shift }; knots.add(time[i]); knots.add(newKnot); coefsList.add(coefs1); coefsList.add(coefs2); } } } } } } knots.add(time[nDataPts - 1]); final int nKnots = knots.size(); _time = new double[nKnots]; for (int i = 0; i < nKnots; ++i) { _time[i] = knots.get(i); } double[][] res = new double[nKnots - 1][3]; for (int i = 0; i < nKnots - 1; ++i) { res[i] = coefsList.get(i); } return new DoubleMatrix2D(res); } /** * @param time * @param spotRates * @return Discrete forwards */ private double[] discFwdsFinder(final double[] time, final double[] spotRates) { final int nDataPts = time.length; double[] res = new double[nDataPts - 1]; for (int i = 0; i < nDataPts - 1; ++i) { res[i] = (spotRates[i + 1] * time[i + 1] - spotRates[i] * time[i]) / (time[i + 1] - time[i]); } return res; } /** * @param time * @param discFwds Discrete forwards * @return Forwards */ private double[] fwdsFinder(final double[] time, final double[] discFwds) { final int nDataPts = time.length; double[] res = new double[nDataPts]; for (int i = 1; i < nDataPts - 1; ++i) { res[i] = (time[i] - time[i - 1]) * discFwds[i] / (time[i + 1] - time[i - 1]) + (time[i + 1] - time[i]) * discFwds[i - 1] / (time[i + 1] - time[i - 1]); } res[0] = 1.5 * discFwds[0] - 0.5 * res[1]; res[nDataPts - 1] = 1.5 * discFwds[nDataPts - 2] - 0.5 * res[nDataPts - 2]; return fwdsModifier(discFwds, res); } /** * Modify forwards such that positivity holds * @param discFwds Discrete forwards * @param fwds Forwards * @return Modified forwards */ private double[] fwdsModifier(final double[] discFwds, final double[] fwds) { final int length = fwds.length; double[] res = new double[length]; res[0] = boundFunc(0., fwds[0], 2. * discFwds[0]); for (int i = 1; i < length - 1; ++i) { res[i] = boundFunc(0., fwds[i], Math.min(2. * discFwds[i - 1], 2. * discFwds[i])); } res[length - 1] = boundFunc(0., fwds[length - 1], 2. * discFwds[length - 2]); return res; } private double boundFunc(final double a, final double b, final double c) { return Math.min(Math.max(a, b), c); } }