Java tutorial
/** * Copyright (C) 2011 - present by OpenGamma Inc. and the OpenGamma group of companies * * Please see distribution for license. */ package com.opengamma.analytics.financial.model.volatility.smile.function; import static com.opengamma.analytics.math.interpolation.Interpolator1DFactory.DOUBLE_QUADRATIC; import static com.opengamma.analytics.math.interpolation.Interpolator1DFactory.FLAT_EXTRAPOLATOR; import java.util.Arrays; import org.apache.commons.lang.Validate; import com.opengamma.analytics.financial.model.option.pricing.analytic.formula.EuropeanVanillaOption; import com.opengamma.analytics.financial.model.option.pricing.fourier.FFTModelGreeks; import com.opengamma.analytics.financial.model.option.pricing.fourier.FFTPricer; import com.opengamma.analytics.financial.model.option.pricing.fourier.HestonCharacteristicExponent; import com.opengamma.analytics.financial.model.option.pricing.fourier.MartingaleCharacteristicExponent; import com.opengamma.analytics.financial.model.volatility.BlackFormulaRepository; import com.opengamma.analytics.math.function.Function1D; import com.opengamma.analytics.math.interpolation.CombinedInterpolatorExtrapolatorFactory; import com.opengamma.analytics.math.interpolation.Interpolator1D; import com.opengamma.analytics.math.interpolation.data.Interpolator1DDataBundle; import com.opengamma.lang.annotation.ExternalFunction; /** * */ public class HestonVolatilityFunction extends VolatilityFunctionProvider<HestonModelData> { /** The FFT pricer */ private static final FFTPricer FFT_PRICER = new FFTPricer(); /** The default interpolator */ private static final Interpolator1D DEFAULT_INTERPOLATOR1D = CombinedInterpolatorExtrapolatorFactory .getInterpolator(DOUBLE_QUADRATIC, FLAT_EXTRAPOLATOR, FLAT_EXTRAPOLATOR); /** The default limit of sigma */ private static final double DEFAULT_LIMIT_SIGMA = 0.3; /** The default limit of alpha */ private static final double DEFAULT_ALPHA = -0.5; /** The limit of sigma */ private final double _limitSigma; /** Alpha */ private final double _alpha; /** The limit tolerance */ private final double _limitTolerance; /** The interpolator */ private final Interpolator1D _interpolator; /** * Default constructor setting sigma, alpha, the limit tolerance and the interpolator to the default values */ public HestonVolatilityFunction() { _limitSigma = DEFAULT_LIMIT_SIGMA; _alpha = DEFAULT_ALPHA; _limitTolerance = 1e-8; _interpolator = DEFAULT_INTERPOLATOR1D; } /** * {@inheritDoc} * Only use this for testing. If you have a set of options with the same expiry but different strikes, use #getVolatilitySetFunction */ @Override public Function1D<HestonModelData, Double> getVolatilityFunction(final EuropeanVanillaOption option, final double forward) { final Function1D<HestonModelData, double[]> func = getVolatilityFunction(forward, new double[] { option.getStrike() }, option.getTimeToExpiry()); return new Function1D<HestonModelData, Double>() { @Override public Double evaluate(final HestonModelData x) { return func.evaluate(x)[0]; } }; } @Override public Function1D<HestonModelData, double[]> getVolatilityFunction(final double forward, final double[] strikes, final double timeToExpiry) { final int n = strikes.length; final double lowestStrike = strikes[0]; final double highestStrike = strikes[n - 1]; return new Function1D<HestonModelData, double[]>() { @SuppressWarnings("synthetic-access") @Override public double[] evaluate(final HestonModelData x) { final MartingaleCharacteristicExponent ce = new HestonCharacteristicExponent(x); //TODO calculations relating to the FFT setup are made each call, even though they will be very similar (depends on Characteristic // Exponent). Maybe worth calculating a typical setup, outside of this function final double[][] strikeNPrice = FFT_PRICER.price(forward, 1.0, timeToExpiry, true, ce, lowestStrike, highestStrike, n, _limitSigma, _alpha, _limitTolerance); final int m = strikeNPrice.length; final double[] k = new double[m]; final double[] vol = new double[m]; int count = 0; for (int i = 0; i < m; i++) { final double strike = strikeNPrice[i][0]; final double price = strikeNPrice[i][1]; if (price > 0.0) { double impVol; try { impVol = BlackFormulaRepository.impliedVolatility(price, forward, strike, timeToExpiry, true); k[count] = strike; vol[count] = impVol; count++; } catch (final IllegalArgumentException e) { //impVol = BlackFormulaRepository.impliedVolatility(price, forward, strike, timeToExpiry, true); } } } final double[] res = new double[n]; if (count == 0) { //i.e. every single price is invalid, which could happen with extreme parameters. All we can do without stopping the // fitter, is return zero vols. for (int i = 0; i < n; i++) { res[i] = 0.0; } } else { double[] validStrikes = new double[count]; double[] validVols = new double[count]; if (count == m) { validStrikes = k; validVols = vol; } else { validStrikes = Arrays.copyOfRange(k, 0, count); validVols = Arrays.copyOfRange(vol, 0, count); } final Interpolator1DDataBundle dataBundle = _interpolator .getDataBundleFromSortedArrays(validStrikes, validVols); for (int i = 0; i < n; i++) { res[i] = _interpolator.interpolate(dataBundle, strikes[i]); } } return res; } }; } /** * Calculates the volatility given Heston model parameters, market data and option data * @param forward The forward * @param strike The strike * @param timeToExpiry The time to expiry * @param kappa kappa * @param theta theta * @param vol0 initial volatility * @param omega omega * @param rho rho * @return The volatility */ @ExternalFunction public double getVolatility(final double forward, final double strike, final double timeToExpiry, final double kappa, final double theta, final double vol0, final double omega, final double rho) { final Function1D<HestonModelData, Double> func = getVolatilityFunction( new EuropeanVanillaOption(strike, timeToExpiry, true), forward); final HestonModelData data = new HestonModelData(kappa, theta, vol0, omega, rho); return func.evaluate(data); } /** * Calculates the volatility given Heston model parameters, market data and an array of strikes * @param forward The forward * @param strikes The strikes * @param timeToExpiry The time to expiry * @param kappa kappa * @param theta theta * @param vol0 initial volatility * @param omega omega * @param rho rho * @return The volatility */ @ExternalFunction public double[] getVolatilitySet(final double forward, final double[] strikes, final double timeToExpiry, final double kappa, final double theta, final double vol0, final double omega, final double rho) { final Function1D<HestonModelData, double[]> func = getVolatilityFunction(forward, strikes, timeToExpiry); final HestonModelData data = new HestonModelData(kappa, theta, vol0, omega, rho); return func.evaluate(data); } @Override public Function1D<HestonModelData, double[]> getVolatilityAdjointFunction(final EuropeanVanillaOption option, final double forward) { final Function1D<HestonModelData, double[][]> func = getVolatilityAdjointFunction(forward, new double[] { option.getStrike() }, option.getTimeToExpiry()); return new Function1D<HestonModelData, double[]>() { @Override public double[] evaluate(final HestonModelData x) { final double[][] temp = func.evaluate(x); Validate.isTrue(temp.length == 1); return temp[0]; } }; } @Override public Function1D<HestonModelData, double[][]> getVolatilityAdjointFunction(final double forward, final double[] strikes, final double timeToExpiry) { final FFTModelGreeks greekCal = new FFTModelGreeks(); final int n = strikes.length; final double lowestStrike = strikes[0]; final double highestStrike = strikes[n - 1]; final double[][] nodeSense = new double[n][]; return new Function1D<HestonModelData, double[][]>() { @SuppressWarnings("synthetic-access") @Override public double[][] evaluate(final HestonModelData x) { final MartingaleCharacteristicExponent ce = new HestonCharacteristicExponent(x); final double[][] greeks = greekCal.getGreeks(forward, 1.0, timeToExpiry, true, ce, lowestStrike, highestStrike, n, _limitSigma, _alpha, _limitTolerance); //1st array is strikes and the second is prices (which we don't need) final double[] k = greeks[0]; final double[] prices = greeks[1]; final int m = k.length; final double[] vols = new double[m]; final double[] vega = new double[m]; for (int i = 0; i < m; i++) { vols[i] = BlackFormulaRepository.impliedVolatility(prices[i], forward, k[i], timeToExpiry, true); } for (int i = 0; i < m; i++) { vega[i] = BlackFormulaRepository.vega(forward, k[i], timeToExpiry, vols[i]); } final Interpolator1DDataBundle dataBundle = _interpolator.getDataBundleFromSortedArrays(k, vols); for (int i = 0; i < n; i++) { nodeSense[i] = _interpolator.getNodeSensitivitiesForValue(dataBundle, strikes[i]); } final int p = greeks.length; final double[][] volSense = new double[p - 2][m]; for (int index = 0; index < p - 2; index++) { for (int i = 0; i < m; i++) { volSense[index][i] = greeks[index + 2][i] / vega[i]; //TODO here is where vega = 0 -> infinity } } //fake the price, forward, and strike sense since we don't used them final double[][] res = new double[n][p + 1]; for (int index = 0; index < p - 2; index++) { final double[] temp = volSense[index]; for (int i = 0; i < n; i++) { final double[] tns = nodeSense[i]; double sum = 0.0; for (int j = 0; j < m; j++) { sum += temp[j] * tns[j]; } res[i][index + 3] = sum; } } return res; } }; } @Override public Function1D<HestonModelData, double[]> getModelAdjointFunction(final EuropeanVanillaOption option, final double forward) { final Function1D<HestonModelData, double[][]> func = getModelAdjointFunction(forward, new double[] { option.getStrike() }, option.getTimeToExpiry()); return new Function1D<HestonModelData, double[]>() { @Override public double[] evaluate(final HestonModelData x) { final double[][] temp = func.evaluate(x); Validate.isTrue(temp.length == 1); return temp[0]; } }; } @Override public Function1D<HestonModelData, double[][]> getModelAdjointFunction(final double forward, final double[] strikes, final double timeToExpiry) { final FFTModelGreeks greekCal = new FFTModelGreeks(); final int n = strikes.length; final double lowestStrike = strikes[0]; final double highestStrike = strikes[n - 1]; final double[][] nodeSense = new double[n][]; return new Function1D<HestonModelData, double[][]>() { @SuppressWarnings("synthetic-access") @Override public double[][] evaluate(final HestonModelData x) { final MartingaleCharacteristicExponent ce = new HestonCharacteristicExponent(x); final double[][] greeks = greekCal.getGreeks(forward, 1.0, timeToExpiry, true, ce, lowestStrike, highestStrike, n, _limitSigma, _alpha, _limitTolerance); //1st array is strikes and the second is prices (which we don't need) final double[] k = greeks[0]; final double[] prices = greeks[1]; final int m = k.length; final double[] kTemp = new double[m]; final double[] vols = new double[m]; final double[] vega = new double[m]; final double[][] modelGreeks = new double[5][m]; int count = 0; for (int i = 0; i < m; i++) { double impVol; try { impVol = BlackFormulaRepository.impliedVolatility(prices[i], forward, k[i], timeToExpiry, true); vols[count] = impVol; vega[count] = BlackFormulaRepository.vega(forward, k[i], timeToExpiry, impVol); kTemp[count] = k[i]; for (int j = 0; j < 5; j++) { modelGreeks[j][count] = greeks[j + 2][i]; } count++; } catch (final IllegalArgumentException e) { //do nothing } } double[] validStrikes = new double[count]; double[] validVols = new double[count]; double[] validVegas = new double[count]; double[][] validModelGreeks = new double[5][count]; if (count == m) { validStrikes = kTemp; validVols = vols; validVegas = vega; validModelGreeks = modelGreeks; } else { validStrikes = Arrays.copyOfRange(k, 0, count); validVols = Arrays.copyOfRange(vols, 0, count); validVegas = Arrays.copyOfRange(vega, 0, count); for (int j = 0; j < 5; j++) { validModelGreeks[j] = Arrays.copyOfRange(modelGreeks[j], 0, count); } } final Interpolator1DDataBundle dataBundle = _interpolator .getDataBundleFromSortedArrays(validStrikes, validVols); for (int i = 0; i < n; i++) { nodeSense[i] = _interpolator.getNodeSensitivitiesForValue(dataBundle, strikes[i]); } final int p = modelGreeks.length; final double[][] volSense = new double[p][count]; for (int index = 0; index < p; index++) { for (int i = 0; i < count; i++) { volSense[index][i] = validModelGreeks[index][i] / validVegas[i]; } } final double[][] res = new double[n][p]; for (int index = 0; index < p; index++) { final double[] temp = volSense[index]; for (int i = 0; i < n; i++) { final double[] tns = nodeSense[i]; double sum = 0.0; for (int j = 0; j < count; j++) { sum += temp[j] * tns[j]; } res[i][index] = sum; } } return res; } }; } @Override public int hashCode() { return toString().hashCode(); } @Override public boolean equals(final Object obj) { if (obj == null) { return false; } if (this == obj) { return true; } if (getClass() != obj.getClass()) { return false; } return true; } @Override public String toString() { return "Heston"; } }