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.option; import com.google.common.math.DoubleMath; import com.opengamma.strata.basics.value.ValueDerivatives; import com.opengamma.strata.collect.ArgChecker; import com.opengamma.strata.collect.array.DoubleArray; import com.opengamma.strata.math.impl.statistics.distribution.NormalDistribution; import com.opengamma.strata.math.impl.statistics.distribution.ProbabilityDistribution; import com.opengamma.strata.product.option.SimpleConstantContinuousBarrier; /** * The price function to compute the price of barrier option in the Black world. * Reference: E. G. Haug (2007) The complete guide to Option Pricing Formulas. Mc Graw Hill. Section 4.17.1. */ public class BlackBarrierPriceFormulaRepository { /** * The normal distribution implementation used in the pricing. */ private static final ProbabilityDistribution<Double> NORMAL = new NormalDistribution(0, 1); /** * Small parameter. */ private static final double SMALL = 1.0e-6; /** * Computes the price of a barrier option. * * @param spot the spot * @param strike the strike * @param timeToExpiry the time to expiry * @param costOfCarry the cost of carry * @param rate the interest rate * @param lognormalVol the lognormal volatility * @param isCall true if call, false otherwise * @param barrier the barrier * @return the price */ public double price(double spot, double strike, double timeToExpiry, double costOfCarry, double rate, double lognormalVol, boolean isCall, SimpleConstantContinuousBarrier barrier) { ArgChecker.notNull(barrier, "barrier"); boolean isKnockIn = barrier.getKnockType().isKnockIn(); boolean isDown = barrier.getBarrierType().isDown(); double h = barrier.getBarrierLevel(); ArgChecker.isFalse(isDown && spot <= barrier.getBarrierLevel(), "The Data is not consistent with an alive barrier (DOWN and spot<=barrier)."); ArgChecker.isFalse(!isDown && spot >= barrier.getBarrierLevel(), "The Data is not consistent with an alive barrier (UP and spot>=barrier)."); int phi = isCall ? 1 : -1; double eta = isDown ? 1 : -1; double df1 = Math.exp(timeToExpiry * (costOfCarry - rate)); double df2 = Math.exp(-rate * timeToExpiry); double sigmaSq = lognormalVol * lognormalVol; double sigmaT = lognormalVol * Math.sqrt(timeToExpiry); if (DoubleMath.fuzzyEquals(Math.min(timeToExpiry, sigmaSq), 0d, SMALL)) { if (isKnockIn) { return 0d; } double dscFwd = df1 * spot; double dscStr = df2 * strike; return isCall ? (dscFwd >= dscStr ? dscFwd - dscStr : 0d) : (dscStr >= dscFwd ? dscStr - dscFwd : 0d); } double mu = (costOfCarry - 0.5 * sigmaSq) / sigmaSq; double m1 = sigmaT * (1 + mu); double x1 = Math.log(spot / strike) / sigmaT + m1; double x2 = Math.log(spot / h) / sigmaT + m1; double y1 = Math.log(h * h / spot / strike) / sigmaT + m1; double y2 = Math.log(h / spot) / sigmaT + m1; double xA = getA(spot, strike, df1, df2, x1, sigmaT, phi); double xB = getA(spot, strike, df1, df2, x2, sigmaT, phi); double xC = getC(spot, strike, df1, df2, y1, sigmaT, h, mu, phi, eta); double xD = getC(spot, strike, df1, df2, y2, sigmaT, h, mu, phi, eta); if (isKnockIn) { // KnockIn if (isDown) { if (isCall) { return strike > h ? xC : xA - xB + xD; } return strike > h ? xB - xC + xD : xA; } if (isCall) { return strike > h ? xA : xB - xC + xD; } return strike > h ? xA - xB + xD : xC; } // KnockOut if (isDown) { if (isCall) { return strike > h ? xA - xC : xB - xD; } return strike > h ? xA - xB + xC - xD : 0d; } if (isCall) { return strike > h ? 0d : xA - xB + xC - xD; } return strike > h ? xB - xD : xA - xC; } /** * Computes the price and derivatives of a barrier option. * * The derivatives are [0] spot, [1] strike, [2] rate, [3] cost-of-carry, [4] volatility, [5] timeToExpiry, [6] spot twice * * @param spot the spot * @param strike the strike * @param timeToExpiry the time to expiry * @param costOfCarry the cost of carry * @param rate the interest rate * @param lognormalVol the lognormal volatility * @param isCall true if call, false otherwise * @param barrier the barrier * @return the price and derivatives */ public ValueDerivatives priceAdjoint(double spot, double strike, double timeToExpiry, double costOfCarry, double rate, double lognormalVol, boolean isCall, SimpleConstantContinuousBarrier barrier) { ArgChecker.notNull(barrier, "barrier"); double[] derivatives = new double[7]; boolean isKnockIn = barrier.getKnockType().isKnockIn(); boolean isDown = barrier.getBarrierType().isDown(); double h = barrier.getBarrierLevel(); ArgChecker.isFalse(isDown && spot <= barrier.getBarrierLevel(), "The Data is not consistent with an alive barrier (DOWN and spot<=barrier)."); ArgChecker.isFalse(!isDown && spot >= barrier.getBarrierLevel(), "The Data is not consistent with an alive barrier (UP and spot>=barrier)."); int phi = isCall ? 1 : -1; double eta = isDown ? 1 : -1; double df1 = Math.exp(timeToExpiry * (costOfCarry - rate)); double df2 = Math.exp(-rate * timeToExpiry); double lognormalVolSq = lognormalVol * lognormalVol; double lognormalVolT = lognormalVol * Math.sqrt(timeToExpiry); if (DoubleMath.fuzzyEquals(Math.min(timeToExpiry, lognormalVolSq), 0d, SMALL)) { if (isKnockIn) { return ValueDerivatives.of(0d, DoubleArray.filled(7)); } double dscFwd = df1 * spot; double dscStr = df2 * strike; double price = 0d; if (isCall) { if (dscFwd >= dscStr) { price = dscFwd - dscStr; derivatives[0] = df1; derivatives[1] = -df2; derivatives[2] = -timeToExpiry * price; derivatives[3] = timeToExpiry * dscFwd; derivatives[5] = (costOfCarry - rate) * dscFwd + rate * dscStr; } } else { if (dscStr >= dscFwd) { price = dscStr - dscFwd; derivatives[0] = -df1; derivatives[1] = df2; derivatives[2] = -timeToExpiry * price; derivatives[3] = -timeToExpiry * dscFwd; derivatives[5] = -rate * dscStr - (costOfCarry - rate) * dscFwd; } } return ValueDerivatives.of(price, DoubleArray.ofUnsafe(derivatives)); } double mu = (costOfCarry - 0.5 * lognormalVolSq) / lognormalVolSq; double m1 = lognormalVolT * (1 + mu); double x1 = Math.log(spot / strike) / lognormalVolT + m1; double x2 = Math.log(spot / h) / lognormalVolT + m1; double y1 = Math.log(h * h / spot / strike) / lognormalVolT + m1; double y2 = Math.log(h / spot) / lognormalVolT + m1; double[] aDerivFirst = new double[6]; double[][] aDerivSecond = new double[2][2]; double xA = getAAdjoint(spot, strike, df1, df2, x1, lognormalVolT, phi, aDerivFirst, aDerivSecond); double[] bDerivFirst = new double[6]; double[][] bDerivSecond = new double[2][2]; double xB = getAAdjoint(spot, strike, df1, df2, x2, lognormalVolT, phi, bDerivFirst, bDerivSecond); double[] cDerivFirst = new double[7]; double[][] cDerivSecond = new double[2][2]; double xC = getCAdjoint(spot, strike, df1, df2, y1, lognormalVolT, h, mu, phi, eta, cDerivFirst, cDerivSecond); double[] dDerivFirst = new double[7]; double[][] dDerivSecond = new double[2][2]; double xD = getCAdjoint(spot, strike, df1, df2, y2, lognormalVolT, h, mu, phi, eta, dDerivFirst, dDerivSecond); double xDBar = 0d; double xCBar = 0d; double xBBar = 0d; double xABar = 0d; double price; if (isKnockIn) { // IN start if (isDown) { // DOWN start if (isCall) { // Call start if (strike > h) { xCBar = 1d; price = xC; } else { xABar = 1d; xBBar = -1d; xDBar = 1d; price = xA - xB + xD; } } else { // Put start if (strike > h) { xBBar = 1d; xCBar = -1d; xDBar = 1d; price = xB - xC + xD; } else { xABar = 1d; price = xA; } } // DOWN end } else { // UP start if (isCall) { if (strike > h) { xABar = 1d; price = xA; } else { xBBar = 1d; xCBar = -1d; xDBar = 1d; price = xB - xC + xD; } } else { if (strike > h) { xABar = 1d; xBBar = -1d; xDBar = 1d; price = xA - xB + xD; } else { xCBar = 1d; price = xC; } } // UP end } // IN end } else { // OUT start if (isDown) { // DOWN start if (isCall) { // CALL start if (strike > h) { xABar = 1d; xCBar = -1d; price = xA - xC; } else { xBBar = 1d; xDBar = -1d; price = xB - xD; } } else { // PUT start if (strike > h) { xABar = 1d; xBBar = -1d; xCBar = 1d; xDBar = -1d; price = xA - xB + xC - xD; } else { price = 0d; } // PUT end } // DOWN end } else { // UP start if (isCall) { if (strike > h) { price = 0d; } else { xABar = 1d; xBBar = -1d; xCBar = 1d; xDBar = -1d; price = xA - xB + xC - xD; } } else { if (strike > h) { xBBar = 1d; xDBar = -1d; price = xB - xD; } else { xABar = 1d; xCBar = -1d; price = xA - xC; } // PUT end } // UP end } // OUT end } double dxyds = 1d / spot / lognormalVolT; double x1Bar = aDerivFirst[4] * xABar; double x2Bar = bDerivFirst[4] * xBBar; double y1Bar = cDerivFirst[4] * xCBar; double y2Bar = dDerivFirst[4] * xDBar; double m1Bar = x1Bar + x2Bar + y1Bar + y2Bar; double muBar = cDerivFirst[6] * xCBar + dDerivFirst[6] * xDBar + lognormalVolT * m1Bar; double lognormalVolTBar = aDerivFirst[5] * xABar + bDerivFirst[5] * xBBar + cDerivFirst[5] * xCBar + dDerivFirst[5] * xDBar - Math.log(h / spot) / (lognormalVolT * lognormalVolT) * y2Bar - Math.log(h * h / spot / strike) / (lognormalVolT * lognormalVolT) * y1Bar - Math.log(spot / h) / (lognormalVolT * lognormalVolT) * x2Bar - Math.log(spot / strike) / (lognormalVolT * lognormalVolT) * x1Bar + (1d + mu) * m1Bar; double lognormalVolSqBar = -costOfCarry / (lognormalVolSq * lognormalVolSq) * muBar; double df2Bar = aDerivFirst[3] * xABar + bDerivFirst[3] * xBBar + cDerivFirst[3] * xCBar + dDerivFirst[3] * xDBar; double df1Bar = aDerivFirst[2] * xABar + bDerivFirst[2] * xBBar + cDerivFirst[2] * xCBar + dDerivFirst[2] * xDBar; derivatives[0] = aDerivFirst[0] * xABar + bDerivFirst[0] * xBBar + cDerivFirst[0] * xCBar + dDerivFirst[0] * xDBar + x1Bar * dxyds + x2Bar * dxyds - y1Bar * dxyds - y2Bar * dxyds; derivatives[1] = aDerivFirst[1] * xABar + bDerivFirst[1] * xBBar + cDerivFirst[1] * xCBar + dDerivFirst[1] * xDBar - (x1Bar + y1Bar) / strike / lognormalVolT; derivatives[2] = -timeToExpiry * (df1 * df1Bar + df2 * df2Bar); derivatives[3] = timeToExpiry * df1 * df1Bar + muBar / lognormalVolSq; derivatives[4] = 2d * lognormalVol * lognormalVolSqBar + Math.sqrt(timeToExpiry) * lognormalVolTBar; derivatives[5] = (costOfCarry - rate) * df1 * df1Bar - rate * df2 * df2Bar + lognormalVolTBar * lognormalVolT * 0.5 / timeToExpiry; derivatives[6] = aDerivSecond[0][0] * xABar + bDerivSecond[0][0] * xBBar + cDerivSecond[0][0] * xCBar + dDerivSecond[0][0] * xDBar + 2d * xABar * aDerivSecond[0][1] * dxyds + 2d * xBBar * bDerivSecond[0][1] * dxyds - 2d * xCBar * cDerivSecond[0][1] * dxyds - 2d * xDBar * dDerivSecond[0][1] * dxyds + xABar * aDerivSecond[1][1] * dxyds * dxyds + xBBar * bDerivSecond[1][1] * dxyds * dxyds + xCBar * cDerivSecond[1][1] * dxyds * dxyds + xDBar * dDerivSecond[1][1] * dxyds * dxyds - x1Bar * dxyds / spot - x2Bar * dxyds / spot + y1Bar * dxyds / spot + y2Bar * dxyds / spot; return ValueDerivatives.of(price, DoubleArray.ofUnsafe(derivatives)); } //------------------------------------------------------------------------- private double getA(double s, double k, double df1, double df2, double x, double lognormalVolT, double phi) { return phi * (s * df1 * NORMAL.getCDF(phi * x) - k * df2 * NORMAL.getCDF(phi * (x - lognormalVolT))); } private double getC(double s, double k, double df1, double df2, double y, double lognormalVolT, double h, double mu, double phi, double eta) { return phi * (s * df1 * Math.pow(h / s, 2d * (mu + 1d)) * NORMAL.getCDF(eta * y) - k * df2 * Math.pow(h / s, 2d * mu) * NORMAL.getCDF(eta * (y - lognormalVolT))); } //------------------------------------------------------------------------- // The first derivatives are [0] s, [1] k, [2] df1, [3] df2, [4] x, [5] lognormalVolT // The second derivatives are [0][0] s twice, [0][1] s and x, [1][0] s and x , [1][1] x twice private double getAAdjoint(double s, double k, double df1, double df2, double x, double lognormalVolT, double phi, double[] firstderivatives, double[][] secondderivatives) { // Forward sweep double n1 = NORMAL.getCDF(phi * x); double n2 = NORMAL.getCDF(phi * (x - lognormalVolT)); double a = phi * (s * df1 * n1 - k * df2 * n2); // Backward sweep double n2Bar = phi * -k * df2; double n1Bar = phi * s * df1; firstderivatives[0] = phi * df1 * n1; firstderivatives[1] = phi * -df2 * n2; firstderivatives[2] = phi * s * n1; firstderivatives[3] = phi * -k * n2; double n1df = NORMAL.getPDF(x); double n2df = NORMAL.getPDF(x - lognormalVolT); firstderivatives[4] = n1df * phi * n1Bar + n2df * phi * n2Bar; firstderivatives[5] = n2df * -phi * n2Bar; secondderivatives[0][0] = 0d; secondderivatives[1][0] = n1df * df1; secondderivatives[0][1] = secondderivatives[1][0]; secondderivatives[1][1] = -x * n1df * phi * n1Bar - (x - lognormalVolT) * n2df * phi * n2Bar; return a; } // The first derivatives are [0] s, [1] k, [2] df1, [3] df2, [4] y, [5] lognormalVolT, [6] mu. // The second derivatives are [0][0] s twice, [0][1] s and y, [1][0] s and y , [1][1] y twice. private double getCAdjoint(double s, double k, double df1, double df2, double y, double lognormalVolT, double h, double mu, double phi, double eta, double[] firstDerivatives, double[][] secondDerivatives) { // Forward sweep double n1 = NORMAL.getCDF(eta * y); double n2 = NORMAL.getCDF(eta * (y - lognormalVolT)); double hsMu1 = Math.pow(h / s, 2d * (mu + 1d)); double hsMu = Math.pow(h / s, 2d * mu); double c = phi * (s * df1 * hsMu1 * n1 - k * df2 * hsMu * n2); // Backward sweep double n1df = NORMAL.getPDF(y); double n2df = NORMAL.getPDF(y - lognormalVolT); double hsMuBar = phi * -k * df2 * n2; double hsMu1Bar = phi * s * df1 * n1; double n2Bar = phi * -k * df2 * hsMu; double n1Bar = phi * s * df1 * hsMu1; firstDerivatives[0] = phi * df1 * hsMu1 * n1 - 2d * mu * hsMu / s * hsMuBar - 2d * (mu + 1d) * hsMu1 / s * hsMu1Bar; // s firstDerivatives[1] = phi * -df2 * hsMu * n2; // k firstDerivatives[2] = phi * s * hsMu1 * n1; // df1 firstDerivatives[3] = phi * -k * hsMu * n2; // df2 firstDerivatives[4] = n1df * eta * n1Bar + n2df * eta * n2Bar; // y firstDerivatives[5] = -n2df * eta * n2Bar; // lognormalVolT firstDerivatives[6] = 2d * Math.log(h / s) * hsMu * hsMuBar + 2d * Math.log(h / s) * hsMu1 * hsMu1Bar; // mu secondDerivatives[0][0] = -2d * (mu + 1d) * phi * df1 * hsMu1 * n1 / s + hsMu * hsMuBar * 2d * mu * (2d * mu + 1d) / (s * s) - hsMu1 * hsMu1Bar * 2d * (mu + 1d) / (s * s) + hsMu1 * hsMu1Bar * 2d * (mu + 1d) * (2d * mu + 3d) / (s * s); secondDerivatives[0][1] = hsMu1 * n1df * phi * df1 * eta - 2d * mu * hsMu / s * phi * -k * df2 * n2df * eta - hsMu1 * n1df * 2d * (mu + 1d) * phi * df1 * eta; secondDerivatives[1][0] = secondDerivatives[0][1]; secondDerivatives[1][1] = -y * n1df * eta * n1Bar - (y - lognormalVolT) * n2df * eta * n2Bar; return c; } }