/**
* 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.fitting;
import java.util.BitSet;
import org.apache.commons.lang.Validate;
import com.opengamma.analytics.financial.model.volatility.smile.function.SmileModelData;
import com.opengamma.analytics.financial.model.volatility.smile.function.VolatilityFunctionProvider;
import com.opengamma.analytics.math.function.Function1D;
import com.opengamma.analytics.math.linearalgebra.DecompositionFactory;
import com.opengamma.analytics.math.matrix.DoubleMatrix1D;
import com.opengamma.analytics.math.matrix.DoubleMatrix2D;
import com.opengamma.analytics.math.matrix.MatrixAlgebra;
import com.opengamma.analytics.math.matrix.OGMatrixAlgebra;
import com.opengamma.analytics.math.minimization.NonLinearParameterTransforms;
import com.opengamma.analytics.math.minimization.NonLinearTransformFunction;
import com.opengamma.analytics.math.statistics.leastsquare.LeastSquareResults;
import com.opengamma.analytics.math.statistics.leastsquare.LeastSquareResultsWithTransform;
import com.opengamma.analytics.math.statistics.leastsquare.NonLinearLeastSquare;
/**
*
* @param <T> The data for the smile model used
*/
public abstract class SmileModelFitter<T extends SmileModelData> {
private static final MatrixAlgebra MA = new OGMatrixAlgebra();
private static final NonLinearLeastSquare SOLVER = new NonLinearLeastSquare(DecompositionFactory.SV_COLT, MA, 1e-12);
private static final Function1D<DoubleMatrix1D, Boolean> UNCONSTRAINED = new Function1D<DoubleMatrix1D, Boolean>() {
@Override
public Boolean evaluate(final DoubleMatrix1D x) {
return true;
}
};
private final VolatilityFunctionProvider<T> _model;
private final Function1D<T, double[]> _volFunc;
private final Function1D<T, double[][]> _volAdjointFunc;
private final DoubleMatrix1D _marketValues;
private final DoubleMatrix1D _errors;
/**
* Attempts to calibrate a model to the implied volatilities of European vanilla options, by minimising the sum of squares between the
* market and model implied volatilities. All the options must be for the same expiry and (implicitly) on the same underlying.
* @param forward The forward value of the underlying
* @param strikes ordered values of strikes
* @param timeToExpiry The time-to-expiry
* @param impliedVols The market implied volatilities
* @param error The 'measurement' error to apply to the market volatility of a particular option TODO: Review should this be part of EuropeanOptionMarketData?
* @param model VolatilityFunctionProvider
*/
public SmileModelFitter(final double forward, final double[] strikes, final double timeToExpiry, final double[] impliedVols,
final double[] error, final VolatilityFunctionProvider<T> model) {
Validate.notNull(strikes, "null strikes");
Validate.notNull(impliedVols, "null implied vols");
Validate.notNull(error, "null errors");
Validate.notNull(model, "null model");
final int n = strikes.length;
Validate.isTrue(n == impliedVols.length, "vols not the same length as strikes");
Validate.isTrue(n == error.length, "errors not the same length as strikes");
_marketValues = new DoubleMatrix1D(impliedVols);
_errors = new DoubleMatrix1D(error);
_volFunc = model.getVolatilityFunction(forward, strikes, timeToExpiry);
_volAdjointFunc = model.getModelAdjointFunction(forward, strikes, timeToExpiry);
_model = model;
}
/**
* Solve using the default NonLinearParameterTransforms for the concrete implementation
* @param start The first guess at the parameter values
* @return The LeastSquareResults
*/
public LeastSquareResultsWithTransform solve(final DoubleMatrix1D start) {
return solve(start, new BitSet());
}
/**
* Solve using the default NonLinearParameterTransforms for the concrete implementation, but with some parameters fixed to their initial
* values (indicated by fixed)
* @param start The first guess at the parameter values
* @param fixed Indicates which parameters are fixed
* @return The LeastSquareResults
*/
public LeastSquareResultsWithTransform solve(final DoubleMatrix1D start, final BitSet fixed) {
final NonLinearParameterTransforms transform = getTransform(start, fixed);
return solve(start, transform);
}
/**
* Solve using a user supplied NonLinearParameterTransforms
* @param start The first guess at the parameter values
* @param transform Transform from model parameters to fitting parameters, and vice versa
* @return The LeastSquareResults
*/
public LeastSquareResultsWithTransform solve(final DoubleMatrix1D start, final NonLinearParameterTransforms transform) {
final NonLinearTransformFunction transFunc = new NonLinearTransformFunction(getModelValueFunction(), getModelJacobianFunction(), transform);
final LeastSquareResults solRes = SOLVER.solve(_marketValues, _errors, transFunc.getFittingFunction(), transFunc.getFittingJacobian(),
transform.transform(start), getConstraintFunction(transform), getMaximumStep());
return new LeastSquareResultsWithTransform(solRes, transform);
}
protected Function1D<DoubleMatrix1D, DoubleMatrix1D> getModelValueFunction() {
return new Function1D<DoubleMatrix1D, DoubleMatrix1D>() {
@SuppressWarnings("synthetic-access")
@Override
public DoubleMatrix1D evaluate(final DoubleMatrix1D x) {
final T data = toSmileModelData(x);
final double[] res = _volFunc.evaluate(data);
return new DoubleMatrix1D(res);
}
};
}
protected Function1D<DoubleMatrix1D, DoubleMatrix2D> getModelJacobianFunction() {
return new Function1D<DoubleMatrix1D, DoubleMatrix2D>() {
@SuppressWarnings("synthetic-access")
@Override
public DoubleMatrix2D evaluate(final DoubleMatrix1D x) {
final T data = toSmileModelData(x);
//this thing will be (#strikes/vols) x (# model Params)
final double[][] volAdjoint = _volAdjointFunc.evaluate(data);
return new DoubleMatrix2D(volAdjoint);
}
};
}
protected abstract DoubleMatrix1D getMaximumStep();
protected abstract NonLinearParameterTransforms getTransform(final DoubleMatrix1D start);
protected abstract NonLinearParameterTransforms getTransform(final DoubleMatrix1D start, final BitSet fixed);
protected abstract T toSmileModelData(final DoubleMatrix1D modelParameters);
protected Function1D<DoubleMatrix1D, Boolean> getConstraintFunction(@SuppressWarnings("unused") final NonLinearParameterTransforms t) {
return UNCONSTRAINED;
}
public VolatilityFunctionProvider<T> getModel() {
return _model;
}
}