Java tutorial
/** * Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies * * Please see distribution for license. */ package com.opengamma.analytics.financial.model.option.pricing.fourier; import static com.opengamma.analytics.math.ComplexMathUtils.add; import static com.opengamma.analytics.math.ComplexMathUtils.divide; import static com.opengamma.analytics.math.ComplexMathUtils.exp; import static com.opengamma.analytics.math.ComplexMathUtils.log; import static com.opengamma.analytics.math.ComplexMathUtils.mod; import static com.opengamma.analytics.math.ComplexMathUtils.multiply; import static com.opengamma.analytics.math.ComplexMathUtils.sqrt; import static com.opengamma.analytics.math.ComplexMathUtils.square; import static com.opengamma.analytics.math.ComplexMathUtils.subtract; import static com.opengamma.analytics.math.number.ComplexNumber.I; import static com.opengamma.analytics.math.number.ComplexNumber.ZERO; import org.apache.commons.lang.Validate; import com.opengamma.analytics.financial.model.volatility.smile.function.HestonModelData; import com.opengamma.analytics.math.function.Function1D; import com.opengamma.analytics.math.number.ComplexNumber; /** * The Heston stochastic volatility model is defined as: * $$ * \begin{align*} * dS_t &= \mu S_t dt + \sqrt{V_t}S_t dW_S(t)\\ * dV_t &= \kappa(\theta - V_0)dt + \omega\sqrt{V_t} dW_V(t)\\ * \text{with}\\ * dW_S(t) dW_V(t) = \rho dt * \end{align*} * $$ * This class represents the characteristic function of the Heston model: * $$ * \begin{align*} * \phi(u) &= e^{C(t, u) + D(t, u)V_0 + iu\ln F}\\ * \text{where}\\ * C(t, u) &= \frac{\kappa\theta}{\omega^2} \left((\kappa - \rho \omega ui + d(u))t - 2\ln\left(\frac{c(u)e^{d(u)t} - 1}{c(u) - 1}\right) \right)\\ * D(t, u) &= \frac{\kappa - \rho \omega ui + d(u)}{\omega^2}\left(\frac{e^{d(u)t} - 1}{c(u)e^{d(u)t} - 1}\right)\\ * c(u) &= \frac{\kappa - \rho \omega ui + d(u)}{\kappa - \rho \omega ui - d(u)}\\ * \text{and}\\ * d(u) &= \sqrt{(\rho \omega ui - \kappa)^2 + iu\omega^2 + \omega^2 u^2} * \end{align*} * $$ */ public class HestonCharacteristicExponent implements MartingaleCharacteristicExponent { private final double _kappa; private final double _theta; private final double _vol0; private final double _omega; private final double _rho; private final double _alphaMin; private final double _alphaMax; /** * * @param kappa mean-reverting speed * @param theta mean-reverting level * @param vol0 starting value of volatility * @param omega volatility-of-volatility * @param rho correlation between the spot process and the volatility process */ public HestonCharacteristicExponent(final double kappa, final double theta, final double vol0, final double omega, final double rho) { _kappa = kappa; _theta = theta; _vol0 = vol0; _omega = omega; _rho = rho; final double t1 = _omega - 2 * _kappa * _rho; final double rhoStar = 1 - _rho * _rho; final double root = Math.sqrt(t1 * t1 + 4 * _kappa * _kappa * rhoStar); _alphaMin = (t1 - root) / _omega / rhoStar - 1; _alphaMax = (t1 + root) / _omega / rhoStar + 1; } public HestonCharacteristicExponent(final HestonModelData data) { Validate.notNull(data, "null data"); _kappa = data.getKappa(); _theta = data.getTheta(); _vol0 = data.getVol0(); _omega = data.getOmega(); _rho = data.getRho(); final double t1 = _omega - 2 * _kappa * _rho; final double rhoStar = 1 - _rho * _rho; final double root = Math.sqrt(t1 * t1 + 4 * _kappa * _kappa * rhoStar); _alphaMin = (t1 - root) / _omega / rhoStar - 1; _alphaMax = (t1 + root) / _omega / rhoStar + 1; } @Override public Function1D<ComplexNumber, ComplexNumber> getFunction(final double t) { return new Function1D<ComplexNumber, ComplexNumber>() { @Override public ComplexNumber evaluate(final ComplexNumber u) { return getValue(u, t); } }; } @Override public ComplexNumber getValue(ComplexNumber u, double t) { // that u = 0 gives zero is true for any characteristic function, that u = -i gives zero is because this is already mean corrected if (u.getReal() == 0.0 && (u.getImaginary() == 0.0 || u.getImaginary() == -1.0)) { return ZERO; } //non-stochastic vol limit if (_omega == 0.0 || mod(multiply(multiply(_omega / _kappa, u), add(I, u))) < 1e-6) { final ComplexNumber z = multiply(u, add(I, u)); if (_kappa * t < 1e-6) { return multiply(-_vol0 / 2 * t, z); } final double var = _theta * t + (_vol0 - _theta) * (1 - Math.exp(-_kappa * t)) / _kappa; return multiply(-var / 2, z); } final ComplexNumber c = getC(u, t); final ComplexNumber dv0 = multiply(_vol0, getD(u, t)); return add(c, dv0); } @Override public Function1D<ComplexNumber, ComplexNumber[]> getAdjointFunction(final double t) { return new Function1D<ComplexNumber, ComplexNumber[]>() { @Override public ComplexNumber[] evaluate(final ComplexNumber u) { return getCharacteristicExponentAdjoint(u, t); } }; } private ComplexNumber[] forwardSweep(final ComplexNumber u, final double t) { ComplexNumber[] w = new ComplexNumber[19]; w[0] = new ComplexNumber(_kappa * _theta / _omega / _omega); w[1] = multiply(u, new ComplexNumber(0, _rho * _omega)); w[2] = subtract(w[1], _kappa); w[3] = square(w[2]); w[4] = multiply(u, new ComplexNumber(0, _omega * _omega)); w[5] = square(multiply(u, _omega)); w[6] = add(w[3], w[4], w[5]); w[7] = sqrt(w[6]); w[8] = subtract(w[7], w[2]); w[9] = multiply(-1.0, add(w[7], w[2])); w[10] = divide(w[8], w[9]); w[11] = multiply(t, w[9]); w[12] = exp(multiply(w[7], -t)); w[13] = divide(subtract(w[10], w[12]), subtract(w[10], 1)); w[14] = log(w[13]); w[15] = divide(subtract(1, w[12]), subtract(w[10], w[12])); w[16] = multiply(w[0], subtract(w[11], multiply(2, w[14]))); w[17] = divide(multiply(w[8], w[15]), _omega * _omega); w[18] = add(w[16], multiply(_vol0, w[17])); return w; } private ComplexNumber[] backwardsSweep(final ComplexNumber[] w, final double t) { final double oneOverOmega2 = 1 / _omega / _omega; final ComplexNumber[] wBar = new ComplexNumber[19]; wBar[18] = new ComplexNumber(1.0); //formal start wBar[17] = new ComplexNumber(_vol0); //Suppressing the multiple by wBar18 wBar[16] = wBar[18]; wBar[15] = multiply(oneOverOmega2, w[8], wBar[17]); wBar[14] = multiply(-2, w[0], wBar[16]); wBar[13] = divide(wBar[14], w[13]); ComplexNumber temp1 = subtract(1.0, w[10]); ComplexNumber temp2 = subtract(w[10], w[12]); ComplexNumber temp3 = square(temp2); wBar[12] = add(multiply(wBar[15], divide(temp1, temp3)), divide(wBar[13], temp1)); wBar[11] = multiply(w[0], wBar[16]); wBar[10] = add(multiply(divide(subtract(w[12], 1), temp3), wBar[15]), multiply(divide(subtract(w[12], 1), square(temp1)), wBar[13])); wBar[9] = subtract(multiply(t, wBar[11]), multiply(divide(w[8], square(w[9])), wBar[10])); wBar[8] = add(multiply(oneOverOmega2, w[15], wBar[17]), divide(wBar[10], w[9])); //subtract(add(multiply(wBar12, multiply(-t, w12)), wBar8), wBar9); wBar[7] = subtract(wBar[8], add(multiply(t, w[12], wBar[12]), wBar[9])); wBar[6] = divide(wBar[7], multiply(2, w[7])); wBar[5] = wBar[6]; wBar[4] = wBar[6]; wBar[3] = wBar[6]; wBar[2] = subtract(multiply(2, w[2], wBar[3]), add(wBar[8], wBar[9])); wBar[1] = wBar[2]; wBar[0] = multiply(subtract(w[11], multiply(2, w[14])), wBar[16]); return wBar; } public ComplexNumber[] getCharacteristicExponentAdjointDebug(final ComplexNumber u, final double t) { ComplexNumber[] res = new ComplexNumber[6]; ComplexNumber[] w = forwardSweep(u, t); ComplexNumber[] wBar = backwardsSweep(w, t); res[0] = w[18]; res[1] = subtract(multiply(_theta / _omega / _omega, wBar[0]), wBar[2]); res[2] = multiply(_kappa / _omega / _omega, wBar[0]); res[3] = multiply(w[17], wBar[18]); res[4] = multiply(1 / _omega, add(multiply(-2, w[0], wBar[0]), multiply(w[1], wBar[1]), multiply(2, w[4], wBar[4]), multiply(2, w[5], wBar[5]), multiply(-2, w[17], wBar[17]))); res[5] = multiply(1 / _rho, w[1], wBar[1]); return res; } @Override public ComplexNumber[] getCharacteristicExponentAdjoint(final ComplexNumber u, final double t) { ComplexNumber[] res = new ComplexNumber[6]; if (u.getReal() == 0.0 && (u.getImaginary() == 0.0 || u.getImaginary() == -1.0)) { for (int i = 0; i < 6; i++) { res[i] = ZERO; } return res; } //non-stochastic vol limit if (_omega == 0.0 || mod(multiply(_omega / _kappa, u, add(I, u))) < 1e-6) { final ComplexNumber z = multiply(u, add(I, u)); // //TODO calculate the omega -> 0 sensitivity without resorting to this hack // HestonCharacteristicExponent ceTemp = this.withOmega(1.1e-6 * _kappa / mod(z)); // ComplexNumber[] temp = ceTemp.getCharacteristicExponentAdjoint(u, t); final double var; if (_kappa * t < 1e-6) { var = _vol0 * t + (_vol0 - _theta) * _kappa * t * t / 2; res[1] = multiply(-0.5 * (_vol0 - _theta) * t * t / 2, z); res[2] = multiply(_kappa * t * t / 4, z); res[3] = multiply(-0.5 * (t + _kappa * t * t / 2), z); res[4] = ZERO; //TODO this is wrong res[5] = ZERO; } else { final double expKappaT = Math.exp(-_kappa * t); var = _theta * t + (_vol0 - _theta) * (1 - expKappaT) / _kappa; res[1] = multiply(-0.5 * (_vol0 - _theta) * (expKappaT * (1 + t * _kappa) - 1) / _kappa / _kappa, z); res[2] = multiply(-0.5 * (t - (1 - expKappaT) / _kappa), z); res[3] = multiply(-0.5 * (1 - expKappaT) / _kappa, z); res[4] = ZERO; //TODO this is wrong res[5] = ZERO; } res[0] = multiply(-var / 2, z); return res; } final double oneOverOmega2 = 1 / _omega / _omega; final double w0 = _kappa * _theta * oneOverOmega2; final ComplexNumber w1 = multiply(u, new ComplexNumber(0, _rho * _omega)); //w1 final ComplexNumber w2 = subtract(w1, _kappa); //w2 final ComplexNumber w3 = square(w2); //w3 final ComplexNumber w4 = multiply(u, new ComplexNumber(0, _omega * _omega)); final ComplexNumber w5 = square(multiply(u, _omega)); final ComplexNumber w6 = add(w3, w4, w5); final ComplexNumber w7 = sqrt(w6); final ComplexNumber w8 = subtract(w7, w2); final ComplexNumber w9 = multiply(-1.0, add(w7, w2)); final ComplexNumber w10 = divide(w8, w9); final ComplexNumber w11 = multiply(t, w9); final ComplexNumber w12 = exp(multiply(w7, -t)); final ComplexNumber w13 = divide(subtract(w10, w12), subtract(w10, 1)); final ComplexNumber w14 = log(w13); final ComplexNumber w15 = divide(subtract(1, w12), subtract(w10, w12)); final ComplexNumber w16 = multiply(w0, subtract(w11, multiply(2, w14))); final ComplexNumber w17 = multiply(oneOverOmega2, w8, w15); final ComplexNumber w18 = add(w16, multiply(_vol0, w17)); res[0] = w18; //value of CE //backwards sweep final ComplexNumber wBar16 = new ComplexNumber(1.0); final ComplexNumber wBar17 = new ComplexNumber(_vol0); final ComplexNumber wBar15 = multiply(oneOverOmega2, wBar17, w8); final ComplexNumber wBar14 = multiply(-2 * w0, wBar16); final ComplexNumber wBar13 = divide(wBar14, w13); final ComplexNumber wBar12a = multiply(wBar15, divide(subtract(1, w10), square(subtract(w10, w12)))); final ComplexNumber wBar12b = multiply(wBar13, divide(1.0, subtract(1.0, w10))); final ComplexNumber wBar12 = add(wBar12a, wBar12b); final ComplexNumber wBar11 = multiply(w0, wBar16); final ComplexNumber wBar10a = multiply(wBar15, divide(w15, subtract(w12, w10))); final ComplexNumber wBar10b = multiply(wBar13, divide(subtract(w12, 1.0), square(subtract(w10, 1.0)))); final ComplexNumber wBar10 = add(wBar10a, wBar10b); final ComplexNumber wBar9a = multiply(t, wBar11); final ComplexNumber wBar9b = multiply(-1.0, wBar10, divide(w10, w9)); final ComplexNumber wBar9 = add(wBar9a, wBar9b); final ComplexNumber wBar8a = multiply(oneOverOmega2, wBar17, w15); final ComplexNumber wBar8b = divide(wBar10, w9); final ComplexNumber wBar8 = add(wBar8a, wBar8b); final ComplexNumber wBar7 = subtract(add(multiply(wBar12, multiply(-t, w12)), wBar8), wBar9); final ComplexNumber wBar6 = multiply(wBar7, divide(0.5, w7)); final ComplexNumber wBar5 = wBar6; final ComplexNumber wBar4 = wBar6; final ComplexNumber wBar3 = wBar6; final ComplexNumber wBar2 = subtract(multiply(wBar3, multiply(2.0, w2)), add(wBar8, wBar9)); final ComplexNumber wBar1 = wBar2; final ComplexNumber wBar0 = multiply(wBar16, subtract(w11, multiply(2, w14))); res[1] = subtract(multiply(wBar0, _theta / _omega / _omega), wBar2); //kappaBar res[2] = multiply(wBar0, _kappa / _omega / _omega); //thetaBar res[3] = w17; //vol0 final ComplexNumber omegaBar1 = multiply(-2 / _omega, w17, wBar17); final ComplexNumber omegaBar2 = multiply(-2.0 * w0 / _omega, wBar0); final ComplexNumber omegaBar3 = multiply(1 / _omega, w1, wBar1); final ComplexNumber omegaBar4 = multiply(2 / _omega, w4, wBar4); final ComplexNumber omegaBar5 = multiply(2 / _omega, w5, wBar5); res[4] = add(omegaBar1, omegaBar2, omegaBar3, omegaBar4, omegaBar5); res[5] = multiply(_omega, u, I, wBar1); return res; } private ComplexNumber getC(final ComplexNumber u, final double t) { final ComplexNumber c1 = multiply(u, new ComplexNumber(0, _rho * _omega)); final ComplexNumber c = c(u); final ComplexNumber d = d(u); final ComplexNumber c3 = multiply(t, subtract(_kappa, add(d, c1))); final ComplexNumber e = exp(multiply(d, -t)); ComplexNumber c4 = divide(subtract(c, e), subtract(c, 1)); c4 = log(c4); return multiply(_kappa * _theta / _omega / _omega, subtract(c3, multiply(2, c4))); } private ComplexNumber getD(final ComplexNumber u, final double t) { final ComplexNumber c1 = multiply(u, new ComplexNumber(0, _rho * _omega)); final ComplexNumber c = c(u); final ComplexNumber d = d(u); final ComplexNumber c3 = add(_kappa, subtract(d, c1)); final ComplexNumber e = exp(multiply(d, -t)); final ComplexNumber c4 = divide(subtract(1, e), subtract(c, e)); return divide(multiply(c3, c4), _omega * _omega); } private ComplexNumber c(final ComplexNumber u) { final ComplexNumber c1 = multiply(u, new ComplexNumber(0, _rho * _omega)); final ComplexNumber c2 = d(u); final ComplexNumber num = add(_kappa, subtract(c2, c1)); final ComplexNumber denom = subtract(_kappa, add(c2, c1)); return divide(num, denom); } private ComplexNumber d(final ComplexNumber u) { ComplexNumber c1 = subtract(multiply(u, new ComplexNumber(0, _rho * _omega)), _kappa); c1 = multiply(c1, c1); final ComplexNumber c2 = multiply(u, new ComplexNumber(0, _omega * _omega)); ComplexNumber c3 = multiply(u, _omega); c3 = multiply(c3, c3); final ComplexNumber c4 = add(add(c1, c2), c3); return multiply(1, sqrt(c4)); } /** * The maximum allowable value of $alpha$ which is given by * $$ * \begin{align*} * t &= \omega - 2 \kappa \rho \\ * \rho^* &= 1 - \rho^2\\ * r &= \sqrt{t^2 + 4\kappa^2\rho^*}\\ * \alpha_{max} &= \frac{t + r}{\omega(\rho^* - 1)} * \end{align*} * $$ * * @return $alpha_{max}$ */ @Override public double getLargestAlpha() { return _alphaMax; } /** The maximum allowable value of $alpha$ which is given by * $$ * \begin{align*} * t &= \omega - 2 \kappa \rho \\ * \rho^* &= 1 - \rho^2\\ * r &= \sqrt{t^2 + 4\kappa^2\rho^*}\\ * \alpha_{min} &= \frac{t - r}{\omega(\rho^* - 1)} * \end{align*} * $$ * @return $alpha_{min}$ */ @Override public double getSmallestAlpha() { return _alphaMin; } /** * Gets the mean-reverting speed. * @return kappa */ public double getKappa() { return _kappa; } /** * Gets the mean-reverting level. * @return theta */ public double getTheta() { return _theta; } /** * Gets the initial volatility * @return initial volatility */ public double getVol0() { return _vol0; } /** * Gets the volatility-of-volatility. * @return omega */ public double getOmega() { return _omega; } /** * Gets the correlation. * @return rho */ public double getRho() { return _rho; } public HestonCharacteristicExponent withKappa(final double kappa) { return new HestonCharacteristicExponent(kappa, _theta, _vol0, _omega, _rho); } public HestonCharacteristicExponent withTheta(final double theta) { return new HestonCharacteristicExponent(_kappa, theta, _vol0, _omega, _rho); } public HestonCharacteristicExponent withVol0(final double vol0) { return new HestonCharacteristicExponent(_kappa, _theta, vol0, _omega, _rho); } public HestonCharacteristicExponent withOmega(final double omega) { return new HestonCharacteristicExponent(_kappa, _theta, _vol0, omega, _rho); } public HestonCharacteristicExponent withRho(final double rho) { return new HestonCharacteristicExponent(_kappa, _theta, _vol0, _omega, rho); } @Override public int hashCode() { final int prime = 31; int result = 1; long temp; temp = Double.doubleToLongBits(_kappa); result = prime * result + (int) (temp ^ (temp >>> 32)); temp = Double.doubleToLongBits(_omega); result = prime * result + (int) (temp ^ (temp >>> 32)); temp = Double.doubleToLongBits(_rho); result = prime * result + (int) (temp ^ (temp >>> 32)); temp = Double.doubleToLongBits(_theta); result = prime * result + (int) (temp ^ (temp >>> 32)); temp = Double.doubleToLongBits(_vol0); result = prime * result + (int) (temp ^ (temp >>> 32)); return result; } @Override public boolean equals(final Object obj) { if (this == obj) { return true; } if (obj == null) { return false; } if (getClass() != obj.getClass()) { return false; } final HestonCharacteristicExponent other = (HestonCharacteristicExponent) obj; if (Double.doubleToLongBits(_kappa) != Double.doubleToLongBits(other._kappa)) { return false; } if (Double.doubleToLongBits(_omega) != Double.doubleToLongBits(other._omega)) { return false; } if (Double.doubleToLongBits(_rho) != Double.doubleToLongBits(other._rho)) { return false; } if (Double.doubleToLongBits(_theta) != Double.doubleToLongBits(other._theta)) { return false; } return Double.doubleToLongBits(_vol0) == Double.doubleToLongBits(other._vol0); } }