com.opengamma.strata.math.impl.statistics.distribution.BivariateNormalDistribution Maven / Gradle / Ivy
Show all versions of strata-math Show documentation
/*
* Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
*
* Please see distribution for license.
*/
package com.opengamma.strata.math.impl.statistics.distribution;
import com.opengamma.strata.collect.ArgChecker;
/**
* The bivariate normal distribution is a continuous probability distribution
* of two variables, $x$ and $y$, with cdf
* $$
* \begin{align*}
* M(x, y, \rho) = \frac{1}{2\pi\sqrt{1 - \rho^2}}\int_{-\infty}^x\int_{-\infty}^{y} e^{\frac{-(X^2 - 2\rho XY + Y^2)}{2(1 - \rho^2)}} dX dY
* \end{align*}
* $$
* where $\rho$ is the correlation between $x$ and $y$.
*
* The implementation of the cdf is taken from "Better Approximations to Cumulative Normal Functions", West
* (link).
*/
public class BivariateNormalDistribution implements ProbabilityDistribution {
private static final ProbabilityDistribution NORMAL = new NormalDistribution(0, 1);
private static final double TWO_PI = 2 * Math.PI;
private static final double[] X = new double[] {0.04691008, 0.23076534, 0.5, 0.76923466, 0.95308992};
private static final double[] Y = new double[] {0.018854042, 0.038088059, 0.0452707394, 0.038088059, 0.018854042};
/**
* Calculates CDF.
*
* @param x The parameters for the function, $(x, y, \rho$, with $-1 \geq \rho \geq 1$, not null
* @return The cdf
*/
@Override
public double getCDF(double[] x) {
ArgChecker.notNull(x, "x");
ArgChecker.isTrue(x.length == 3, "Need a, b and rho values");
ArgChecker.isTrue(x[2] >= -1 && x[2] <= 1, "Correlation must be >= -1 and <= 1");
double a = x[0];
double b = x[1];
double rho = x[2];
if (a == Double.POSITIVE_INFINITY || b == Double.POSITIVE_INFINITY) {
return 1;
}
if (a == Double.NEGATIVE_INFINITY || b == Double.NEGATIVE_INFINITY) {
return 0;
}
double sumSq = (a * a + b * b) / 2.;
double rho1, rho2, rho3, ab, absDiff, h5, c, d, mult = 0, rho3Sq, eab, e, result;
if (Math.abs(rho) >= 0.7) {
rho1 = 1 - rho * rho;
rho2 = Math.sqrt(rho1);
if (rho < 0) {
b *= -1;
}
ab = a * b;
eab = Math.exp(-ab / 2.);
if (Math.abs(rho) < 1) {
absDiff = Math.abs(a - b);
h5 = absDiff * absDiff / 2.;
absDiff = absDiff / rho2;
c = 0.5 - ab / 8.;
d = 3. - 2. * c * h5;
mult = 0.13298076 * absDiff * d * (1 - NORMAL.getCDF(absDiff)) - Math.exp(-h5 / rho1) * (d + c * rho1) * 0.053051647;
for (int i = 0; i < 5; i++) {
rho3 = rho2 * X[i];
rho3Sq = rho3 * rho3;
rho1 = Math.sqrt(1 - rho3Sq);
if (eab == 0) {
e = 0;
} else {
e = Math.exp(-ab / (1 + rho1)) / rho1 / eab;
}
mult = mult - Y[i] * Math.exp(-h5 / rho3Sq) * (e - 1 - c * rho3Sq);
}
}
double corr = Double.isNaN(mult) ? 0. : mult * rho2 * eab;
result = corr + NORMAL.getCDF(Math.min(a, b));
if (rho < 0) {
result = NORMAL.getCDF(a) - result;
}
return result;
}
ab = a * b;
if (rho != 0) {
for (int i = 0; i < 5; i++) {
rho3 = rho * X[i];
rho1 = 1 - rho3 * rho3;
mult = mult + Y[i] * Math.exp((rho3 * ab - sumSq) / rho1) / Math.sqrt(rho1);
}
}
double corr = Double.isNaN(mult) ? 0. : rho * mult;
return NORMAL.getCDF(a) * NORMAL.getCDF(b) + corr;
}
/**
* {@inheritDoc}
* @return Not supported
* @throws UnsupportedOperationException always
*/
@Override
public double getInverseCDF(double[] p) {
throw new UnsupportedOperationException();
}
/**
* Calculates PDF.
*
* @param x The parameters for the function, $(x, y, \rho$, with $-1 \geq \rho \geq 1$, not null
* @return The pdf
*/
@Override
public double getPDF(double[] x) {
ArgChecker.notNull(x, "x");
ArgChecker.isTrue(x.length == 3, "Need a, b and rho values");
ArgChecker.isTrue(x[2] >= -1 && x[2] <= 1, "Correlation must be >= -1 and <= 1");
double denom = 1 - x[2] * x[2];
return Math.exp(-(x[0] * x[0] - 2 * x[2] * x[0] * x[1] + x[1] * x[1]) / (2 * denom)) / (TWO_PI * Math.sqrt(denom));
}
/**
* {@inheritDoc}
* @return Not supported
* @throws UnsupportedOperationException always
*/
@Override
public double nextRandom() {
throw new UnsupportedOperationException();
}
}