com.opengamma.analytics.financial.model.option.pricing.fourier.IntegratedCIRTimeChangeCharacteristicExponent.java Source code

Introduction

Here is the source code for com.opengamma.analytics.financial.model.option.pricing.fourier.IntegratedCIRTimeChangeCharacteristicExponent.java

Source

/**
 * 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.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.subtract;
import static com.opengamma.analytics.math.number.ComplexNumber.I;

import org.apache.commons.lang.NotImplementedException;

import com.opengamma.analytics.math.TrigonometricFunctionUtils;
import com.opengamma.analytics.math.function.Function1D;
import com.opengamma.analytics.math.number.ComplexNumber;

/**
 * The Cox-Ingersoll-Ross process is a mean-reverting positive process, with SDE:
 * $$
 * \begin{align*}
 * dy_t = \kappa(\theta - y_t)dt + \lambda\sqrt{y_t}dW_t
 * \end{align*}
 * $$
 * and characteristic exponent
 * $$
 * \begin{align*}
 * \psi(u, t; \kappa, \theta, \lambda) &= \frac{2\kappa\theta}{\lambda^2}\left[ 
 * \frac{\kappa t}{2} - \ln\left(\cosh\left(\frac{\gamma t}{2}\right) + \frac{\kappa}{\gamma}\sinh\left(\frac{\gamma t}{2}\right)\right)
 * + \frac{2iu}{\kappa + \gamma \coth\left(\frac{\gamma t}{2}\right)}\right]\\
 * \text{where}\\
 * \gamma &= \sqrt{\kappa^2 - 2 \lambda^2 iu}
 * \end{align*}
 * $$
 */
public class IntegratedCIRTimeChangeCharacteristicExponent implements StocasticClockCharcteristicExponent {
    private final double _kappa;
    private final double _theta;
    private final double _lambda;
    private final double _alphaMax;

    /**
     * 
     * @param kappa The mean-reverting speed
     * @param theta The mean 
     * @param lambda The volatility
     */
    public IntegratedCIRTimeChangeCharacteristicExponent(final double kappa, final double theta,
            final double lambda) {
        _kappa = kappa;
        _theta = theta;
        _lambda = lambda;
        _alphaMax = _kappa * _kappa / 2 / _lambda / _lambda;
    }

    @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) {
        if (u.getReal() == 0.0 && u.getImaginary() == 0.0) {
            return new ComplexNumber(0.0);
        }

        final ComplexNumber ui = multiply(I, u);

        //handle small lambda properly 
        if (2 * mod(u) * _lambda * _lambda / _kappa / _kappa < 1e-6) {
            final double d = _theta * t + (1 - _theta) * (1 - Math.exp(-_kappa * t)) / _kappa;
            return multiply(d, ui);
        }

        ComplexNumber temp = subtract(_kappa * _kappa, multiply(2 * _lambda * _lambda, ui));
        final ComplexNumber gamma = sqrt(temp);
        final ComplexNumber gammaHalfT = multiply(gamma, t / 2.0);
        temp = divide(multiply(2, ui), add(_kappa, divide(gamma, TrigonometricFunctionUtils.tanh(gammaHalfT))));
        final ComplexNumber kappaOverGamma = divide(_kappa, gamma);
        final double power = 2 * _kappa * _theta / _lambda / _lambda;
        final ComplexNumber res = add(
                multiply(power, subtract(_kappa * t / 2, getLogCoshSinh(gammaHalfT, kappaOverGamma))), temp);
        return res;
    }

    // ln(cosh(a) + bsinh(a)
    private ComplexNumber getLogCoshSinh(final ComplexNumber a, final ComplexNumber b) {
        final ComplexNumber temp = add(TrigonometricFunctionUtils.cosh(a),
                multiply(b, TrigonometricFunctionUtils.sinh(a)));
        return log(temp);
    }

    /**
     * 
     * @return $\frac{\kappa^2}{2\lambda^2}$
     */
    @Override
    public double getLargestAlpha() {
        return _alphaMax;
    }

    /**
     * 
     * @return $-\infty$
     */
    @Override
    public double getSmallestAlpha() {
        return Double.NEGATIVE_INFINITY;
    }

    /**
     * Gets the mean-reverting speed.
     * @return kappa
     */
    public double getKappa() {
        return _kappa;
    }

    /**
     * Gets the mean.
     * @return theta
     */
    public double getTheta() {
        return _theta;
    }

    /**
     * Gets the volatility.
     * @return lambda
     */
    public double getLambda() {
        return _lambda;
    }

    @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(_lambda);
        result = prime * result + (int) (temp ^ (temp >>> 32));
        temp = Double.doubleToLongBits(_theta);
        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 IntegratedCIRTimeChangeCharacteristicExponent other = (IntegratedCIRTimeChangeCharacteristicExponent) obj;
        if (Double.doubleToLongBits(_kappa) != Double.doubleToLongBits(other._kappa)) {
            return false;
        }
        if (Double.doubleToLongBits(_lambda) != Double.doubleToLongBits(other._lambda)) {
            return false;
        }
        return Double.doubleToLongBits(_theta) == Double.doubleToLongBits(other._theta);
    }

    @Override
    public ComplexNumber[] getCharacteristicExponentAdjoint(ComplexNumber u, double t) {
        throw new NotImplementedException();
    }

    @Override
    public Function1D<ComplexNumber, ComplexNumber[]> getAdjointFunction(double t) {
        throw new NotImplementedException();
    }

}