Java tutorial
/** * Copyright (C) 2016 - present by OpenGamma Inc. and the OpenGamma group of companies * * Please see distribution for license. */ package com.opengamma.strata.pricer.impl.volatility.local; import static com.opengamma.strata.market.curve.interpolator.CurveInterpolators.LINEAR; import static com.opengamma.strata.market.curve.interpolator.CurveInterpolators.TIME_SQUARE; import java.util.ArrayList; import java.util.Arrays; import java.util.List; import java.util.function.Function; import com.google.common.collect.ImmutableList; import com.opengamma.strata.collect.array.DoubleArray; import com.opengamma.strata.collect.array.DoubleMatrix; import com.opengamma.strata.collect.tuple.DoublesPair; import com.opengamma.strata.collect.tuple.Pair; import com.opengamma.strata.market.ValueType; import com.opengamma.strata.market.surface.DefaultSurfaceMetadata; import com.opengamma.strata.market.surface.InterpolatedNodalSurface; import com.opengamma.strata.market.surface.Surface; import com.opengamma.strata.market.surface.SurfaceMetadata; import com.opengamma.strata.market.surface.SurfaceName; import com.opengamma.strata.market.surface.interpolator.GridSurfaceInterpolator; import com.opengamma.strata.market.surface.interpolator.SurfaceInterpolator; import com.opengamma.strata.pricer.fxopt.RecombiningTrinomialTreeData; import com.opengamma.strata.pricer.impl.option.BlackFormulaRepository; import com.opengamma.strata.pricer.impl.option.BlackScholesFormulaRepository; /** * Local volatility calculation based on trinomila tree model. * <p> * Emanuel Derman, Iraj Kani and Neil Chriss, "Implied Trinomial Trees of the Volatility Smile" (1996). */ public class ImpliedTrinomialTreeLocalVolatilityCalculator implements LocalVolatilityCalculator { /** * The default interpolator. */ private static final GridSurfaceInterpolator DEFAULT_INTERPOLATOR = GridSurfaceInterpolator.of(TIME_SQUARE, LINEAR); /** * The number of steps in trinomial tree. */ private final int nSteps; /** * The maximum value of time in trinomial tree. * <p> * The time step in the tree is then given by {@code maxTime/nSteps}. */ private final double maxTime; /** * The interpolator for local volatilities. * <p> * The resulting local volatilities are interpolated by this interpolator along time and spot dimensions. */ private final SurfaceInterpolator interpolator; /** * Creates an instance with default setups. * <p> * The number of time steps is 20, and the tree covers up to 3 years. * The time square linear interpolator is used for time direction, * whereas the linear interpolator is used for spot dimension. * The extrapolation is flat for both the dimensions. */ public ImpliedTrinomialTreeLocalVolatilityCalculator() { this(20, 3d, DEFAULT_INTERPOLATOR); } /** * Creates an instance with the number of steps and maximum time fixed. * <p> * The default interpolators are used: the time square linear interpolator for time direction, * the linear interpolator for spot dimension, and flat extrapolator for both the dimensions. * * @param nSteps the number of steps * @param maxTime the maximum time */ public ImpliedTrinomialTreeLocalVolatilityCalculator(int nSteps, double maxTime) { this(nSteps, maxTime, DEFAULT_INTERPOLATOR); } /** * Creates an instance by specifying the number of steps, maximum time, and 2D interpolator. * * @param nSteps number of steps * @param maxTime the maximum time * @param interpolator the interpolator */ public ImpliedTrinomialTreeLocalVolatilityCalculator(int nSteps, double maxTime, SurfaceInterpolator interpolator) { this.nSteps = nSteps; this.maxTime = maxTime; this.interpolator = interpolator; } //------------------------------------------------------------------------- @Override public InterpolatedNodalSurface localVolatilityFromImpliedVolatility(Surface impliedVolatilitySurface, double spot, Function<Double, Double> interestRate, Function<Double, Double> dividendRate) { Function<DoublesPair, Double> surface = new Function<DoublesPair, Double>() { @Override public Double apply(DoublesPair tk) { return impliedVolatilitySurface.zValue(tk); } }; ImmutableList<double[]> localVolData = calibrate(surface, spot, interestRate, dividendRate).getFirst(); SurfaceMetadata metadata = DefaultSurfaceMetadata.builder().xValueType(ValueType.YEAR_FRACTION) .yValueType(ValueType.STRIKE).zValueType(ValueType.LOCAL_VOLATILITY) .surfaceName(SurfaceName.of("localVol_" + impliedVolatilitySurface.getName())).build(); return InterpolatedNodalSurface.ofUnsorted(metadata, DoubleArray.ofUnsafe(localVolData.get(0)), DoubleArray.ofUnsafe(localVolData.get(1)), DoubleArray.ofUnsafe(localVolData.get(2)), interpolator); } /** * Calibrate trinomial tree to implied volatility surface. * * @param impliedVolatilitySurface the implied volatility surface * @param spot the spot * @param interestRate the interest rate * @param dividendRate the dividend rate * @return the trinomial tree */ public RecombiningTrinomialTreeData calibrateImpliedVolatility( Function<DoublesPair, Double> impliedVolatilitySurface, double spot, Function<Double, Double> interestRate, Function<Double, Double> dividendRate) { return calibrate(impliedVolatilitySurface, spot, interestRate, dividendRate).getSecond(); } @Override public InterpolatedNodalSurface localVolatilityFromPrice(Surface callPriceSurface, double spot, Function<Double, Double> interestRate, Function<Double, Double> dividendRate) { double[][] stateValue = new double[nSteps + 1][]; double[] df = new double[nSteps]; List<DoubleMatrix> probability = new ArrayList<DoubleMatrix>(nSteps); int nTotal = (nSteps - 1) * (nSteps - 1) + 1; double[] timeRes = new double[nTotal]; double[] spotRes = new double[nTotal]; double[] volRes = new double[nTotal]; // uniform grid based on TrigeorgisLatticeSpecification, using reference values double refPrice = callPriceSurface.zValue(maxTime, spot) * Math.exp(interestRate.apply(maxTime) * maxTime); double refForward = spot * Math.exp((interestRate.apply(maxTime) - dividendRate.apply(maxTime)) * maxTime); double refVolatility = BlackFormulaRepository.impliedVolatility(refPrice, refForward, spot, maxTime, true); double dt = maxTime / nSteps; double dx = refVolatility * Math.sqrt(3d * dt); double upFactor = Math.exp(dx); double downFactor = Math.exp(-dx); double[] adSec = new double[2 * nSteps + 1]; double[] assetPrice = new double[2 * nSteps + 1]; for (int i = nSteps; i > -1; --i) { if (i == 0) { resolveFirstLayer(interestRate, dividendRate, nTotal, dt, spot, adSec, assetPrice, timeRes, spotRes, volRes, df, stateValue, probability); } else { double time = dt * i; double zeroRate = interestRate.apply(time); double zeroDividendRate = dividendRate.apply(time); int nNodes = 2 * i + 1; double[] assetPriceLocal = new double[nNodes]; double[] callOptionPrice = new double[nNodes]; double[] putOptionPrice = new double[nNodes]; int position = i - 1; double assetTmp = spot * Math.pow(upFactor, i); // call options for upper half nodes for (int j = nNodes - 1; j > position - 1; --j) { assetPriceLocal[j] = assetTmp; callOptionPrice[j] = callPriceSurface.zValue(time, assetPriceLocal[j]); assetTmp *= downFactor; } // put options for lower half nodes assetTmp = spot * Math.pow(downFactor, i); for (int j = 0; j < position + 2; ++j) { assetPriceLocal[j] = assetTmp; putOptionPrice[j] = callPriceSurface.zValue(time, assetPriceLocal[j]) - spot * Math.exp(-zeroDividendRate * time) + Math.exp(-zeroRate * time) * assetPriceLocal[j]; assetTmp *= upFactor; } resolveLayer(interestRate, dividendRate, i, nTotal, position, dt, zeroRate, zeroDividendRate, callOptionPrice, putOptionPrice, adSec, assetPrice, assetPriceLocal, timeRes, spotRes, volRes, df, stateValue, probability); } } SurfaceMetadata metadata = DefaultSurfaceMetadata.builder().xValueType(ValueType.YEAR_FRACTION) .yValueType(ValueType.STRIKE).zValueType(ValueType.LOCAL_VOLATILITY) .surfaceName(SurfaceName.of("localVol_" + callPriceSurface.getName())).build(); return InterpolatedNodalSurface.ofUnsorted(metadata, DoubleArray.ofUnsafe(timeRes), DoubleArray.ofUnsafe(spotRes), DoubleArray.ofUnsafe(volRes), interpolator); } //----------------------------------------------------------------------- private Pair<ImmutableList<double[]>, RecombiningTrinomialTreeData> calibrate( Function<DoublesPair, Double> impliedVolatilitySurface, double spot, Function<Double, Double> interestRate, Function<Double, Double> dividendRate) { double[][] stateValue = new double[nSteps + 1][]; double[] df = new double[nSteps]; double[] timePrim = new double[nSteps + 1]; List<DoubleMatrix> probability = new ArrayList<DoubleMatrix>(nSteps); int nTotal = (nSteps - 1) * (nSteps - 1) + 1; double[] timeRes = new double[nTotal]; double[] spotRes = new double[nTotal]; double[] volRes = new double[nTotal]; // uniform grid based on TrigeorgisLatticeSpecification double volatility = impliedVolatilitySurface.apply(DoublesPair.of(maxTime, spot)); double dt = maxTime / nSteps; double dx = volatility * Math.sqrt(3d * dt); double upFactor = Math.exp(dx); double downFactor = Math.exp(-dx); double[] adSec = new double[2 * nSteps + 1]; double[] assetPrice = new double[2 * nSteps + 1]; for (int i = nSteps; i > -1; --i) { timePrim[i] = dt * i; if (i == 0) { resolveFirstLayer(interestRate, dividendRate, nTotal, dt, spot, adSec, assetPrice, timeRes, spotRes, volRes, df, stateValue, probability); } else { double zeroRate = interestRate.apply(timePrim[i]); double zeroDividendRate = dividendRate.apply(timePrim[i]); double zeroCostRate = zeroRate - zeroDividendRate; int nNodes = 2 * i + 1; double[] assetPriceLocal = new double[nNodes]; double[] callOptionPrice = new double[nNodes]; double[] putOptionPrice = new double[nNodes]; int position = i - 1; double assetTmp = spot * Math.pow(upFactor, i); // call options for upper half nodes for (int j = nNodes - 1; j > position - 1; --j) { assetPriceLocal[j] = assetTmp; double impliedVol = impliedVolatilitySurface .apply(DoublesPair.of(timePrim[i], assetPriceLocal[j])); callOptionPrice[j] = BlackScholesFormulaRepository.price(spot, assetPriceLocal[j], timePrim[i], impliedVol, zeroRate, zeroCostRate, true); assetTmp *= downFactor; } // put options for lower half nodes assetTmp = spot * Math.pow(downFactor, i); for (int j = 0; j < position + 2; ++j) { assetPriceLocal[j] = assetTmp; double impliedVol = impliedVolatilitySurface .apply(DoublesPair.of(timePrim[i], assetPriceLocal[j])); putOptionPrice[j] = BlackScholesFormulaRepository.price(spot, assetPriceLocal[j], timePrim[i], impliedVol, zeroRate, zeroCostRate, false); assetTmp *= upFactor; } resolveLayer(interestRate, dividendRate, i, nTotal, position, dt, zeroRate, zeroDividendRate, callOptionPrice, putOptionPrice, adSec, assetPrice, assetPriceLocal, timeRes, spotRes, volRes, df, stateValue, probability); } } ImmutableList<double[]> localVolData = ImmutableList.of(timeRes, spotRes, volRes); RecombiningTrinomialTreeData treeData = RecombiningTrinomialTreeData.of(DoubleMatrix.ofUnsafe(stateValue), probability, DoubleArray.ofUnsafe(df), DoubleArray.ofUnsafe(timePrim)); return Pair.of(localVolData, treeData); } // resolve the t=0 layer private void resolveFirstLayer(Function<Double, Double> interestRate, Function<Double, Double> dividendRate, int nTotal, double dt, double spot, double[] adSec, double[] assetPrice, double[] timeRes, double[] spotRes, double[] volRes, double[] df, double[][] stateValue, List<DoubleMatrix> probability) { double discountFactor = Math.exp(-interestRate.apply(dt) * dt); double fwdFactor = Math.exp((interestRate.apply(dt) - dividendRate.apply(dt)) * dt); double upProb = adSec[2] / discountFactor; double midProb = getMiddle(upProb, fwdFactor, spot, assetPrice[0], assetPrice[1], assetPrice[2]); double dwProb = 1d - upProb - midProb; double fwd = spot * fwdFactor; timeRes[nTotal - 1] = dt; spotRes[nTotal - 1] = spot; double var = (dwProb * Math.pow(assetPrice[0] - fwd, 2) + midProb * Math.pow(assetPrice[1] - fwd, 2) + upProb * Math.pow(assetPrice[2] - fwd, 2)) / (fwd * fwd * dt); volRes[nTotal - 1] = Math.sqrt(0.5 * (var + volRes[nTotal - 2] * volRes[nTotal - 2])); probability.add(0, DoubleMatrix.ofUnsafe(new double[][] { { dwProb, midProb, upProb } })); df[0] = discountFactor; stateValue[0] = new double[] { spot }; } // resolve the i-th layer private void resolveLayer(Function<Double, Double> interestRate, Function<Double, Double> dividendRate, int i, int nTotal, int position, double dt, double zeroRate, double zeroDividendRate, double[] callOptionPrice, double[] putOptionPrice, double[] adSec, double[] assetPrice, double[] assetPriceLocal, double[] timeRes, double[] spotRes, double[] volRes, double[] df, double[][] stateValue, List<DoubleMatrix> probability) { int positionLocal = position; int nNodes = callOptionPrice.length; double[] adSecLocal = new double[nNodes]; // AD security prices from call options for (int j = nNodes - 1; j > positionLocal; --j) { adSecLocal[j] = callOptionPrice[j - 1]; for (int k = j + 1; k < nNodes; ++k) { adSecLocal[j] -= (assetPriceLocal[k] - assetPriceLocal[j - 1]) * adSecLocal[k]; } adSecLocal[j] /= (assetPriceLocal[j] - assetPriceLocal[j - 1]); } ++positionLocal; // AD security prices from put options for (int j = 0; j < positionLocal; ++j) { adSecLocal[j] = putOptionPrice[j + 1]; for (int k = 0; k < j; ++k) { adSecLocal[j] -= (assetPriceLocal[j + 1] - assetPriceLocal[k]) * adSecLocal[k]; } adSecLocal[j] /= (assetPriceLocal[j + 1] - assetPriceLocal[j]); } if (i != nSteps) { double time = dt * i; double timeNext = dt * (i - 1); double rate = (zeroRate * time - interestRate.apply(timeNext) * timeNext) / dt; double dividend = (zeroDividendRate * time - dividendRate.apply(timeNext) * timeNext) / dt; double cost = rate - dividend; double discountFactor = Math.exp(-rate * dt); double fwdFactor = Math.exp(cost * dt); double[][] prob = new double[nNodes][3]; // highest node prob[nNodes - 1][2] = adSec[nNodes + 1] / adSecLocal[nNodes - 1] / discountFactor; prob[nNodes - 1][1] = getMiddle(prob[nNodes - 1][2], fwdFactor, assetPriceLocal[nNodes - 1], assetPrice[nNodes - 1], assetPrice[nNodes], assetPrice[nNodes + 1]); prob[nNodes - 1][0] = 1d - prob[nNodes - 1][2] - prob[nNodes - 1][1]; correctProbability(prob[nNodes - 1], fwdFactor, assetPriceLocal[nNodes - 1], assetPrice[nNodes - 1], assetPrice[nNodes], assetPrice[nNodes + 1]); // second highest node prob[nNodes - 2][2] = (adSec[nNodes] / discountFactor - prob[nNodes - 1][1] * adSecLocal[nNodes - 1]) / adSecLocal[nNodes - 2]; prob[nNodes - 2][1] = getMiddle(prob[nNodes - 2][2], fwdFactor, assetPriceLocal[nNodes - 2], assetPrice[nNodes - 2], assetPrice[nNodes - 1], assetPrice[nNodes]); prob[nNodes - 2][0] = 1d - prob[nNodes - 2][2] - prob[nNodes - 2][1]; correctProbability(prob[nNodes - 2], fwdFactor, assetPriceLocal[nNodes - 2], assetPrice[nNodes - 2], assetPrice[nNodes - 1], assetPrice[nNodes]); // subsequent nodes for (int j = nNodes - 3; j > -1; --j) { prob[j][2] = (adSec[j + 2] / discountFactor - prob[j + 2][0] * adSecLocal[j + 2] - prob[j + 1][1] * adSecLocal[j + 1]) / adSecLocal[j]; prob[j][1] = getMiddle(prob[j][2], fwdFactor, assetPriceLocal[j], assetPrice[j], assetPrice[j + 1], assetPrice[j + 2]); prob[j][0] = 1d - prob[j][1] - prob[j][2]; correctProbability(prob[j], fwdFactor, assetPriceLocal[j], assetPrice[j], assetPrice[j + 1], assetPrice[j + 2]); } // local variance int offset = nTotal - i * i - 1; double[] varBare = new double[nNodes]; for (int k = 0; k < nNodes; ++k) { double fwd = assetPriceLocal[k] * fwdFactor; varBare[k] = (prob[k][0] * Math.pow(assetPrice[k] - fwd, 2) + prob[k][1] * Math.pow(assetPrice[k + 1] - fwd, 2) + prob[k][2] * Math.pow(assetPrice[k + 2] - fwd, 2)) / (fwd * fwd * dt); if (varBare[k] < 0d) { throw new IllegalArgumentException("Negative variance"); } } // smoothing for (int k = 0; k < nNodes - 2; ++k) { double var = (k == 0 || k == nNodes - 3) ? (varBare[k] + varBare[k + 1] + varBare[k + 2]) / 3d : (varBare[k - 1] + varBare[k] + varBare[k + 1] + varBare[k + 2] + varBare[k + 3]) / 5d; volRes[offset + k] = i == nSteps - 1 ? Math.sqrt(var) : Math.sqrt(0.5 * (var + volRes[offset - (2 * i - k)] * volRes[offset - (2 * i - k)])); timeRes[offset + k] = dt * (i + 1d); spotRes[offset + k] = assetPriceLocal[k + 1]; } probability.add(0, DoubleMatrix.ofUnsafe(prob)); df[i] = discountFactor; } stateValue[i] = Arrays.copyOf(assetPriceLocal, nNodes); System.arraycopy(adSecLocal, 0, adSec, 0, nNodes); System.arraycopy(assetPriceLocal, 0, assetPrice, 0, nNodes); } private void correctProbability(double[] probability, double factor, double assetBase, double assertPriceLow, double assertPriceMid, double assetPriceHigh) { if (!(probability[2] > 0d && probability[1] > 0d && probability[0] > 0d)) { double fwd = assetBase * factor; if (fwd <= assertPriceMid && fwd > assertPriceLow) { probability[0] = 0.5 * (fwd - assertPriceLow) / (assetPriceHigh - assertPriceLow); probability[2] = 0.5 * ((assetPriceHigh - fwd) / (assetPriceHigh - assertPriceLow) + (assertPriceMid - fwd) / (assertPriceMid - assertPriceLow)); } else if (fwd < assetPriceHigh && fwd > assertPriceMid) { probability[0] = 0.5 * ((fwd - assertPriceMid) / (assetPriceHigh - assertPriceLow) + (fwd - assertPriceLow) / (assetPriceHigh - assertPriceLow)); probability[2] = 0.5 * (assetPriceHigh - fwd) / assetPriceHigh; } probability[1] = 1d - probability[0] - probability[2]; } } private double getMiddle(double upProbability, double factor, double assetBase, double assetPrevDw, double assetPrevMd, double assetPrevUp) { return (factor * assetBase - assetPrevDw - upProbability * (assetPrevUp - assetPrevDw)) / (assetPrevMd - assetPrevDw); } }