net.finmath.timeseries.models.parametric.DisplacedLognormal Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of finmath-lib Show documentation
Show all versions of finmath-lib Show documentation
finmath lib is a Mathematical Finance Library in Java.
It provides algorithms and methodologies related to mathematical finance.
The newest version!
/*
* (c) Copyright Christian P. Fries, Germany. Contact: [email protected].
*
* Created on 15.07.2012
*/
package net.finmath.timeseries.models.parametric;
import java.io.Serializable;
import java.util.HashMap;
import java.util.Map;
import org.apache.commons.math3.analysis.MultivariateFunction;
import org.apache.commons.math3.optimization.GoalType;
import org.apache.commons.math3.optimization.PointValuePair;
import net.finmath.timeseries.HistoricalSimulationModel;
/**
* Displaced log-normal process with constanst volatility.
*
* This class estimate the process
* \[
* \mathrm{d} \log(X + a) = \frac{\sigma}{b + a} \mathrm{d}W(t)
* \]
* where \( a > -min(X(t_{i}) \) and thus \( X+a > 0 \) and \( b = 1 - -min(X(t_{i}) \) \) and
* \( \sigma \) is a constant.
*
* The choice of b ensures that b+a ≥ 1.
* For a=0 we have a log-normal process with volatility σ/(b + a).
* For a=infinity we have a normal process with volatility σ.
*
* @author Christian Fries
* @version 1.0
*/
public class DisplacedLognormal implements HistoricalSimulationModel {
private final double[] values;
private final double lowerBoundDisplacement;
private final double upperBoundDisplacement = 10000000;
private final int windowIndexStart;
private final int windowIndexEnd;
private final int maxIterations = 1000000;
public DisplacedLognormal(final double[] values) {
this.values = values;
windowIndexStart = 0;
windowIndexEnd = values.length-1;
double valuesMin = Double.MAX_VALUE;
for (int i = windowIndexStart; i <= windowIndexEnd; i++) {
valuesMin = Math.min(values[i], valuesMin);
}
lowerBoundDisplacement = -valuesMin+1;
}
public DisplacedLognormal(final double[] values, final double lowerBoundDisplacement) {
this.values = values;
windowIndexStart = 0;
windowIndexEnd = values.length-1;
double valuesMin = Double.MAX_VALUE;
for (int i = windowIndexStart; i <= windowIndexEnd; i++) {
valuesMin = Math.min(values[i], valuesMin);
}
this.lowerBoundDisplacement = Math.max(-valuesMin+1,lowerBoundDisplacement);
}
public DisplacedLognormal(final double[] values, final int windowIndexStart, final int windowIndexEnd) {
this.values = values;
this.windowIndexStart = windowIndexStart;
this.windowIndexEnd = windowIndexEnd;
double valuesMin = Double.MAX_VALUE;
for (int i = windowIndexStart; i <= windowIndexEnd; i++) {
valuesMin = Math.min(values[i], valuesMin);
}
lowerBoundDisplacement = -valuesMin+1;
}
public DisplacedLognormal(final double[] values, final double lowerBoundDisplacement, final int windowIndexStart, final int windowIndexEnd) {
this.values = values;
this.windowIndexStart = windowIndexStart;
this.windowIndexEnd = windowIndexEnd;
double valuesMin = Double.MAX_VALUE;
for (int i = windowIndexStart; i <= windowIndexEnd; i++) {
valuesMin = Math.min(values[i], valuesMin);
}
this.lowerBoundDisplacement = Math.max(-valuesMin+1,lowerBoundDisplacement);
}
@Override
public HistoricalSimulationModel getCloneWithWindow(final int windowIndexStart, final int windowIndexEnd) {
return new DisplacedLognormal(values, windowIndexStart, windowIndexEnd);
}
public HistoricalSimulationModel getCloneWithWindow(final double lowerBoundDisplacement, final int windowIndexStart, final int windowIndexEnd) {
return new DisplacedLognormal(values, lowerBoundDisplacement, windowIndexStart, windowIndexEnd);
}
public double getLogLikelihoodForParameters(final double omega, final double alpha, final double beta, final double displacement)
{
double logLikelihood = 0.0;
final double volScaling = (1+Math.abs(displacement));
double volSquaredEstimate = 0.0;
for (int i = windowIndexStart+1; i <= windowIndexEnd-1; i++) {
final double eval = volScaling * (Math.log((values[i]+displacement)/(values[i-1]+displacement)));
volSquaredEstimate += eval*eval;
}
volSquaredEstimate /= windowIndexEnd-windowIndexStart;
double eval = volScaling * (Math.log((values[windowIndexStart+1]+displacement)/(values[windowIndexStart+1-1]+displacement)));
for (int i = windowIndexStart+1; i <= windowIndexEnd-1; i++) {
final double evalNext = volScaling * (Math.log((values[i+1]+displacement)/(values[i]+displacement)));
final double volSquared = volSquaredEstimate / volScaling * volScaling; // h = (sigma*)^2, volSquared = (sigma^a)^2
logLikelihood += - Math.log(volSquaredEstimate) - 2 * Math.log((values[i+1]+displacement)/volScaling) - evalNext*evalNext / volSquaredEstimate;
eval = evalNext;
}
logLikelihood += - Math.log(2 * Math.PI) * (windowIndexEnd-windowIndexStart);
logLikelihood *= 0.5;
return logLikelihood;
}
public double getLastResidualForParameters(final double omega, final double alpha, final double beta, final double displacement) {
final double volScaling = (1+Math.abs(displacement));
double h = omega / (1.0 - alpha - beta);
for (int i = windowIndexStart+1; i <= windowIndexEnd; i++) {
final double eval = volScaling * (Math.log((values[i]+displacement)/(values[i-1]+displacement)));
// double eval = volScaling * (values[i]-values[i-1])/(values[i-1]+displacement);
h = omega + alpha * eval * eval + beta * h;
}
return h;
}
public double[] getQuantilPredictionsForParameters(final double omega, final double alpha, final double beta, final double displacement, final double[] quantiles) {
final double[] szenarios = new double[windowIndexEnd-windowIndexStart+1-1];
final double volScaling = (1+Math.abs(displacement));
double volSquaredEstimate = 0.0;
for (int i = windowIndexStart+1; i <= windowIndexEnd-1; i++) {
final double eval = volScaling * (Math.log((values[i]+displacement)/(values[i-1]+displacement)));
volSquaredEstimate += eval*eval;
}
volSquaredEstimate /= windowIndexEnd-windowIndexStart;
double vol = Math.sqrt(volSquaredEstimate) / volScaling;
for (int i = windowIndexStart+1; i <= windowIndexEnd; i++) {
final double y = Math.log((values[i]+displacement)/(values[i-1]+displacement));
// double y = (values[i]-values[i-1])/(values[i-1]+displacement);
szenarios[i-windowIndexStart-1] = y / vol;
final double eval = volScaling * y;
vol = Math.sqrt(volSquaredEstimate) / volScaling;
}
java.util.Arrays.sort(szenarios);
final double[] quantileValues = new double[quantiles.length];
for(int i=0; i getBestParameters() {
return getBestParameters(null);
}
/* (non-Javadoc)
* @see net.finmath.timeseries.HistoricalSimulationModel#getBestParameters(java.util.Map)
*/
@Override
public Map getBestParameters(final Map guess) {
// Create the objective function for the solver
class GARCHMaxLikelihoodFunction implements MultivariateFunction, Serializable {
private static final long serialVersionUID = 7072187082052755854L;
@Override
public double value(final double[] variables) {
/*
* Transform variables: The solver variables are in (-\infty, \infty).
* We transform the variable to the admissible domain for GARCH, that is
* omega > 0, 0 < alpha < 1, 0 < beta < (1-alpha), displacement > lowerBoundDisplacement ??????
* ???? usually for GARCH the restrictions are written like omega > 0, alpha > 0, beta > 0, and alpha + beta < 1
*/
final double omega = Math.exp(variables[0]);
final double mucorr = Math.exp(-Math.exp(-variables[1]));
final double muema = Math.exp(-Math.exp(-variables[2]));
final double beta = mucorr * muema;
final double alpha = mucorr - beta;
// double alpha = 1.0/(1.0+Math.exp(-variables[1]));
// double beta = (1.0-alpha)*1.0/(1.0+Math.exp(-variables[2]));
final double displacementNormed = 1.0/(1.0+Math.exp(-variables[3]));
final double displacement = (upperBoundDisplacement-lowerBoundDisplacement)*displacementNormed+lowerBoundDisplacement;
double logLikelihood = getLogLikelihoodForParameters(omega,alpha,beta,displacement);
// Penalty to prevent solver from hitting the bounds
logLikelihood -= Math.max(1E-30-omega,0)/1E-30;
logLikelihood -= Math.max(1E-30-alpha,0)/1E-30;
logLikelihood -= Math.max((alpha-1)+1E-30,0)/1E-30;
logLikelihood -= Math.max(1E-30-beta,0)/1E-30;
logLikelihood -= Math.max((beta-1)+1E-30,0)/1E-30;
logLikelihood -= Math.max(1E-30-displacementNormed,0)/1E-30;
logLikelihood -= Math.max((displacementNormed-1)+1E-30,0)/1E-30;
return logLikelihood;
}
}
final GARCHMaxLikelihoodFunction objectiveFunction = new GARCHMaxLikelihoodFunction();
// Create a guess for the solver
double guessOmega = 1.0;
double guessAlpha = 0.2;
double guessBeta = 0.2;
double guessDisplacement = (lowerBoundDisplacement + upperBoundDisplacement) / 2.0;
if(guess != null) {
// A guess was provided, use that one
guessOmega = (Double)guess.get("Omega");
guessAlpha = (Double)guess.get("Alpha");
guessBeta = (Double)guess.get("Beta");
guessDisplacement = (Double)guess.get("Displacement");
}
// Constrain guess to admissible range
guessOmega = restrictToOpenSet(guessOmega, 0.0, Double.MAX_VALUE);
guessAlpha = restrictToOpenSet(guessAlpha, 0.0, 1.0);
guessBeta = restrictToOpenSet(guessBeta, 0.0, 1.0-guessAlpha);
guessDisplacement = restrictToOpenSet(guessDisplacement, lowerBoundDisplacement, upperBoundDisplacement);
final double guessMucorr = guessAlpha + guessBeta;
final double guessMuema = guessBeta / (guessAlpha+guessBeta);
// Transform guess to solver coordinates
final double[] guessParameters = new double[4];
guessParameters[0] = Math.log(guessOmega);
guessParameters[1] = -Math.log(-Math.log(guessMucorr));
guessParameters[2] = -Math.log(-Math.log(guessMuema));
guessParameters[3] = -Math.log(1.0/((guessDisplacement-lowerBoundDisplacement)/(upperBoundDisplacement-lowerBoundDisplacement))-1.0);
// Seek optimal parameter configuration
// org.apache.commons.math3.optimization.direct.BOBYQAOptimizer optimizer2 = new org.apache.commons.math3.optimization.direct.BOBYQAOptimizer(6);
final org.apache.commons.math3.optimization.direct.CMAESOptimizer optimizer2 = new org.apache.commons.math3.optimization.direct.CMAESOptimizer();
double[] bestParameters = null;
try {
final PointValuePair result = optimizer2.optimize(
maxIterations,
objectiveFunction,
GoalType.MAXIMIZE,
guessParameters /* start point */
);
bestParameters = result.getPoint();
} catch(final org.apache.commons.math3.exception.MathIllegalStateException e) {
// Retry with new guess. This guess corresponds to omaga=1, alpha=0.5; beta=0.25; displacement=1+lowerBoundDisplacement;
final double[] guessParameters2 = {0.0, 0.0, 0.0, 10.0};
/* PointValuePair result = optimizer2.optimize(
maxIterations,
objectiveFunction,
GoalType.MAXIMIZE,
guessParameters2
);*/
System.out.println("Solver failed");
bestParameters = guessParameters2;//result.getPoint();
}
// Transform parameters to GARCH parameters
final double omega = Math.exp(bestParameters[0]);
final double mucorr = Math.exp(-Math.exp(-bestParameters[1]));
final double muema = Math.exp(-Math.exp(-bestParameters[2]));
final double beta = mucorr * muema;
final double alpha = mucorr - beta;
final double displacementNormed = 1.0/(1.0+Math.exp(-bestParameters[3]));
final double displacement = (upperBoundDisplacement-lowerBoundDisplacement)*displacementNormed+lowerBoundDisplacement;
final double[] quantiles = {0.01, 0.05, 0.5};
final double[] quantileValues = this.getQuantilPredictionsForParameters(omega, alpha, beta, displacement, quantiles);
final Map results = new HashMap<>();
results.put("Omega", omega);
results.put("Alpha", alpha);
results.put("Beta", beta);
results.put("Displacement", displacement);
results.put("Likelihood", this.getLogLikelihoodForParameters(omega, alpha, beta, displacement));
results.put("Vol", Math.sqrt(this.getLastResidualForParameters(omega, alpha, beta, displacement)));
results.put("Quantile=1%", quantileValues[0]);
results.put("Quantile=5%", quantileValues[1]);
results.put("Quantile=50%", quantileValues[2]);
return results;
}
private static double restrictToOpenSet(double value, final double lowerBond, final double upperBound) {
value = Math.max(value, lowerBond * (1.0+Math.signum(lowerBond)*1E-15) + 1E-15);
value = Math.min(value, upperBound * (1.0-Math.signum(upperBound)*1E-15) - 1E-15);
return value;
}
}