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

Introduction

Here is the source code for com.opengamma.analytics.financial.model.option.pricing.fourier.IntegralLimitCalculator.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 org.apache.commons.lang.Validate;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

import com.opengamma.analytics.math.ComplexMathUtils;
import com.opengamma.analytics.math.MathException;
import com.opengamma.analytics.math.function.Function1D;
import com.opengamma.analytics.math.number.ComplexNumber;
import com.opengamma.analytics.math.rootfinding.BracketRoot;
import com.opengamma.analytics.math.rootfinding.BrentSingleRootFinder;
import com.opengamma.analytics.math.rootfinding.RealSingleRootFinder;

/**
 * A calculator to determine the upper limit of the Fourier integral for a
 * characteristic function $\phi$.
 * <p>
 * The upper limit is found by determining the root of the function:
 * $$
 * \begin{align*}
 * f(x) = \ln\left(\left|\phi(x - i(1 + \alpha))\right|\right)
 * \end{align*}
 * $$
 * where $\alpha$ is the contour (which is parallel to the real axis and
 * shifted down by $1 + \alpha$) over which to integrate.
 * 
 */
public class IntegralLimitCalculator {
    private static Logger s_log = LoggerFactory.getLogger(IntegralLimitCalculator.class);
    private static BracketRoot s_bracketRoot = new BracketRoot();
    private static final RealSingleRootFinder s_root = new BrentSingleRootFinder(1e-1);

    /**
     * 
     * @param psi The characteristic function, not null
     * @param alpha The value of $\alpha$, not 0 or -1
     * @param tol The tolerance for the root
     * @return The root
     */
    public double solve(final Function1D<ComplexNumber, ComplexNumber> psi, final double alpha, final double tol) {
        Validate.notNull(psi, "psi null");
        Validate.isTrue(alpha != 0.0 && alpha != -1.0, "alpha cannot be -1 or 0");
        Validate.isTrue(tol > 0.0, "need tol > 0");

        final double k = Math.log(tol)
                + Math.log(ComplexMathUtils.mod(psi.evaluate(new ComplexNumber(0.0, -(1 + alpha)))));
        final Function1D<Double, Double> f = new Function1D<Double, Double>() {
            @Override
            public Double evaluate(final Double x) {
                final ComplexNumber z = new ComplexNumber(x, -(1 + alpha));
                return Math.log(ComplexMathUtils.mod(psi.evaluate(z))) - k;
            }
        };
        double[] range = null;
        try {
            range = s_bracketRoot.getBracketedPoints(f, 0.0, 200.0);
        } catch (MathException e) {
            s_log.warn("Could not find integral limit. Using default of 500");
            return 500.0;
        }
        return s_root.getRoot(f, range[0], range[1]);
    }

}