List of usage examples for java.lang Math exp
@HotSpotIntrinsicCandidate public static double exp(double a)
From source file:ch.unil.genescore.vegas.FarebrotherExperiment.java
/** * ruben.cpp:// ww w . jav a2s . c om * Algorithm AS 204 Appl. Statist. (1984) Vol. 33, No.3 * ruben evaluates the probability that a positive definite quadratic form in Normal variates is less than a given value */ private void ruben() { ifault_ = -41; // I added this initialization --daniel //void ruben(double *lambda, int *mult, double *delta, int *n, double *c, double *mode, int *maxit, double *eps, double *dnsty, int *ifault, double *res) { // Initialized at 0 in the R code double dnsty = 0.0; int i, k, m, j; //double pnorm(double q, double mean, double sd, int lower_tail, int log_p); double ao, aoinv, z, bbeta, eps2, hold, hold2, sum, sum1, dans, lans, pans, prbty, prbtyAgain, tol; double[] gamma = new double[lambda_.length]; double[] theta = new double[lambda_.length]; double[] a = new double[maxit_]; double[] b = new double[maxit_]; if ((lambda_.length < 1) || (q_ <= 0) || (maxit_ < 1) || (eps_ <= 0.0)) { //res_ = -2.0; ifault_ = 2; return; } tol = -200.0; // Preliminaries sum = lambda_[0]; bbeta = sum; for (i = 1; i <= lambda_.length; i++) { hold = lambda_[i - 1]; if ((hold <= 0.0) || (h_[i - 1] < 1) || (delta_[i - 1] < 0.0)) { res_ = -7.1;//new BigDecimal(-71.0); ifault_ = -i; return; } if (bbeta > hold) bbeta = hold; // calcul du max des lambdas if (sum < hold) sum = hold; // calcul du min des lambdas } if (mode_ > 0.0) { // if ((2.0/(1.0/bbeta+1.0/sum))>1.8*sum) bbeta = sum; // comme dans NAG : methode avec betaA bbeta = mode_ * bbeta; } else { bbeta = 2.0 / (1.0 / bbeta + 1.0 / sum); // methode avec betaB } k = 0; sum = 1.0; sum1 = 0.0; for (i = 1; i <= lambda_.length; i++) { hold = bbeta / lambda_[i - 1]; gamma[i - 1] = 1.0 - hold; sum = sum * Math.pow(hold, h_[i - 1]); //???? pas sur .. sum1 = sum1 + delta_[i - 1]; k = k + h_[i - 1]; theta[i - 1] = 1.0; } ao = Math.exp(0.5 * (Math.log(sum) - sum1)); if (ao <= 0.0) { res_ = 0.0;//new BigDecimal(0.0); dnsty = 0.0; ifault_ = 1; // returns after this (the rest of the function is in the else) } else { // evaluate probability and density of chi-squared on k degrees of freedom. The constant 0.22579135264473 is ln(sqrt(pi/2)) z = q_ / bbeta; if ((k % 2) == 0) { // k est un entier donc on regarde si k est divisible par 2: k == (k/2)*k i = 2; lans = -0.5 * z; dans = Math.exp(lans); pans = 1.0 - dans; } else { i = 1; lans = -0.5 * (z + Math.log(z)) - 0.22579135264473; dans = Math.exp(lans); //pans = pnorm(Math.sqrt(z),0.0,1.0,1,0) - pnorm(-Math.sqrt(z),0.0,1.0,1,0); pans = normal_.cumulativeProbability(Math.sqrt(z)) - normal_.cumulativeProbability(-Math.sqrt(z)); } k = k - 2; for (j = i; j <= k; j = j + 2) { if (lans < tol) { lans = lans + Math.log(z / (double) j); dans = Math.exp(lans); } else { dans = dans * z / (double) j; } pans = pans - dans; } // evaluate successive terms of expansion prbty = pans; prbtyAgain = 1 - ao * pans; BigDecimal prbtyBig = new BigDecimal(prbtyAgain); dnsty = dans; eps2 = eps_ / ao; aoinv = 1.0 / ao; sum = aoinv - 1.0; for (m = 1; m <= maxit_; m++) { sum1 = 0.0; for (i = 1; i <= lambda_.length; i++) { hold = theta[i - 1]; hold2 = hold * gamma[i - 1]; theta[i - 1] = hold2; sum1 = sum1 + hold2 * h_[i - 1] + m * delta_[i - 1] * (hold - hold2); } sum1 = 0.5 * sum1; b[m - 1] = sum1; for (i = m - 1; i >= 1; i--) { sum1 = sum1 + b[i - 1] * a[m - i - 1]; } sum1 = sum1 / (double) m; a[m - 1] = sum1; k = k + 2; if (lans < tol) { lans = lans + Math.log(z / (double) k); dans = Math.exp(lans); } else { dans = dans * z / (double) k; } pans = pans - dans; sum = sum - sum1; dnsty = dnsty + dans * sum1; sum1 = pans * sum1; prbty = prbty + sum1; prbtyBig.subtract(new BigDecimal(sum1)); prbtyAgain = prbtyAgain - ao * sum1; if (prbty < (-aoinv)) { res_ = -3.0;//new BigDecimal(-3.0); //res_ = -3.0; ifault_ = 3; return; } if (Math.abs(pans * sum) < eps2) { if (Math.abs(sum1) < eps2) { //if (Math.abs(sum) < eps2) { ifault_ = 0; // this is overwritten below (ifault_=4) -- I now changed the code below // I COMMENTED THIS SO WE CAN See IF THERE WAS CONVERGENCE OR NOT -daniel //m = maxit_+1; // and WHY would you do that? break; //} } } } if (m > maxit_) // I ADDED THIS IF, OTHERWISE IT MAKES NO SENSE -daniel ifault_ = 4; // Check if I understood correctly assert ifault_ == 0 || ifault_ == 4; dnsty = ao * dnsty / (bbeta + bbeta); prbty = ao * prbty; //prbtyAgain=ao*prbtyAgain + (1-ao); //prbtyBig.multiply(new BigDecimal(ao)); //prbtyBig.add(new BigDecimal(1.0)); //prbtyBig.subtract(new BigDecimal(ao)); //prtbyBig.subtract(ao); // With my edits above, this now makes a bit mores sense if (prbty < 0.0 || prbty > 1.0) { ifault_ = ifault_ + 5; // why the ... would they write it like this? I.e., ifault_ = 9 } else { if (dnsty < 0.0) ifault_ = ifault_ + 6; } res_ = prbtyAgain; } return; }
From source file:common.Utilities.java
public static List<Map<Integer, Double>> getNormalizedMaps(List<UserData> userLines, boolean resource) { List<Map<Integer, Integer>> maps = (resource ? getResMaps(userLines) : getUserMaps(userLines)); List<Map<Integer, Double>> relMaps = new ArrayList<Map<Integer, Double>>(); for (Map<Integer, Integer> m : maps) { double denom = getMapDenom(m); Map<Integer, Double> relM = new LinkedHashMap<Integer, Double>(); for (Map.Entry<Integer, Integer> entry : m.entrySet()) { relM.put(entry.getKey(), Math.exp((double) entry.getValue().intValue()) / denom); }// www . j a va 2 s .co m relMaps.add(relM); } return relMaps; }
From source file:br.ufrgs.enq.jcosmo.test.VLEdiagramsGAMESS.java
public JPanel calcAcAcidOctanol() throws Exception { double T = 293.15; setLayout(new BorderLayout()); COSMOSACDataBase db = COSMOSACDataBase.getInstance(); COSMOSACCompound c1 = db.getComp("acetic-acid"); COSMOSACCompound c2 = db.getComp("1-octanol"); double[] cavityVolume = new double[2]; cavityVolume[0] = c1.Vcosmo;//from w w w . j a va2 s.c om cavityVolume[1] = c2.Vcosmo; double[][] sigma = new double[2][]; sigma[0] = c1.sigma; sigma[1] = c2.sigma; SigmaProfileGenerator c1Sigma = new SigmaProfileGenerator(SigmaProfileGenerator.FileType.GAMESS, c1.name.toLowerCase() + ".log"); SigmaProfileGenerator c2Sigma = new SigmaProfileGenerator(SigmaProfileGenerator.FileType.GAMESS, c2.name.toLowerCase() + ".log"); sigma[0] = c1Sigma.getSigmaProfile(); sigma[1] = c2Sigma.getSigmaProfile(); COSMOSAC cosmosac = new COSMOSAC(); cosmosac.setParameters(cavityVolume, c1.charge, sigma); cosmosac.setIgnoreSG(true); cosmosac.setSigmaHB(COSMOSAC.SIGMAHB * 1.2); cosmosac.setTemperature(T); double[] x1 = new double[n]; double[] x2 = new double[n]; double[] gamma1 = new double[n]; double[] gamma2 = new double[n]; double[] z = new double[2]; double[] lnGamma = new double[2]; z[0] = 0.00; int j = 0; while (z[0] < 1.0001) { z[1] = 1 - z[0]; x1[j] = z[0]; x2[j] = z[1]; cosmosac.setComposition(z); cosmosac.activityCoefficient(lnGamma); gamma1[j] = Math.exp(lnGamma[0]); gamma2[j] = Math.exp(lnGamma[1]); z[0] += 0.05; j++; } double[] Psat = new double[2]; Psat[0] = 1600; Psat[1] = 10; double[] P = calcPx(x1, x2, gamma1, gamma2, Psat); double[] y1 = calcY(x1, gamma1, Psat, P); XYPlot plot1; XYSeriesCollection dataset = new XYSeriesCollection(); XYSeries liq = new XYSeries("liquid"); XYSeries vap = new XYSeries("vapor"); XYSeries raoult = new XYSeries("Raoult's Law"); for (int i = 0; i < n; i++) { liq.add(x1[i], P[i]); vap.add(y1[i], P[i]); } raoult.add(0, Psat[1]); raoult.add(1, Psat[0]); dataset.addSeries(liq); dataset.addSeries(vap); dataset.addSeries(raoult); JFreeChart chart = ChartFactory.createXYLineChart(null, "Mole Fraction: x1, y1", "P/KPa", null, PlotOrientation.VERTICAL, true, true, false); plot1 = (XYPlot) chart.getPlot(); plot1.getDomainAxis().setRange(new Range(0.0, 1.0)); plot1.setDataset(dataset); ChartPanel chartPanel = new ChartPanel(chart); JPanel jp1 = new JPanel(new BorderLayout()); jp1.add(chartPanel); add(jp1, BorderLayout.CENTER); setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE); setSize(400, 500); return jp1; }
From source file:hivemall.utils.math.StatsUtils.java
/** * @param mu1 mean vector of the first normal distribution * @param sigma1 covariance matrix of the first normal distribution * @param mu2 mean vector of the second normal distribution * @param sigma2 covariance matrix of the second normal distribution * @return the Hellinger distance between two multivariate normal distributions * @link https://en.wikipedia.org/wiki/Hellinger_distance#Examples *///from w w w .j a v a 2 s . c om public static double hellingerDistance(@Nonnull final RealVector mu1, @Nonnull final RealMatrix sigma1, @Nonnull final RealVector mu2, @Nonnull final RealMatrix sigma2) { RealVector muSub = mu1.subtract(mu2); RealMatrix sigmaMean = sigma1.add(sigma2).scalarMultiply(0.5d); LUDecomposition LUsigmaMean = new LUDecomposition(sigmaMean); double denominator = Math.sqrt(LUsigmaMean.getDeterminant()); if (denominator == 0.d) { return 1.d; // avoid divide by zero } RealMatrix sigmaMeanInv = LUsigmaMean.getSolver().getInverse(); // has inverse iff det != 0 double sigma1Det = MatrixUtils.det(sigma1); double sigma2Det = MatrixUtils.det(sigma2); double numerator = Math.pow(sigma1Det, 0.25d) * Math.pow(sigma2Det, 0.25d) * Math.exp(-0.125d * sigmaMeanInv.preMultiply(muSub).dotProduct(muSub)); return 1.d - numerator / denominator; }
From source file:ArrayUtil.java
public static double sumexp(double[] d) { double r = 0; for (int i = 0; i < d.length; i++) r += Math.exp(d[i]); return r;/*from www.j a v a 2 s.c om*/ }
From source file:com.opengamma.analytics.financial.model.finitedifference.MarkovChainApprox.java
private double getI2(final double t) { final double t2 = t * t; final double lt = _lambda * t; if (lt < 1e-7) { return _probState1; }/*from ww w. j a va 2 s . co m*/ final double a = _theta; final double a2 = a * a; final double b = _lambda; final double b2 = b * b; final double p = _probState1; final double temp1 = a2 * b2 + (2 * a * b * p - 4 * a2 * b + 2 * a * b) / t + (-4 * a * p + 2 * p + 6 * a2 - 4 * a) / t2; final double temp2 = (2 * a * b * p - 2 * b * p - 2 * a2 * b + 2 * a * b) / t + (4 * a * p - 2 * p - 6 * a2 + 4 * a) / t2; return (temp1 + temp2 * Math.exp(-b * t)) / b2; // double a = _theta; // double b = _theta * (_probState1 - _theta); // double c = (_probState1 - _theta) * (1 - _theta); // double d = _theta * (1 - _theta); // return (b + d) * t / _lambda + a * t * t / 2 + (c - b - d) / _lambda / _lambda + Math.exp(-_lambda * t) * (-c * t * _lambda - c + b + d) / _lambda / _lambda; }
From source file:com.opengamma.analytics.financial.interestrate.inflation.method.InflationMarketModelConvexityAdjustmentForCoupon.java
/** * Computes the convexity adjustment for year on year inflation coupon with an interpolated index. * @param coupon The year on year coupon. * @param inflationConvexity The inflation provider. * @return The convexity adjustment.//from www .j a v a 2 s .co m */ public double getYearOnYearConvexityAdjustment(final CouponInflationYearOnYearInterpolation coupon, final InflationConvexityAdjustmentProviderInterface inflationConvexity) { Validate.notNull(coupon, "Coupon"); Validate.notNull(inflationConvexity, "Inflation"); final double firstFixingTime = coupon.getWeightStart() * coupon.getReferenceStartTime()[0] + (1 - coupon.getWeightStart()) * coupon.getReferenceStartTime()[1]; final double secondFixingTime = coupon.getWeightEnd() * coupon.getReferenceEndTime()[0] + (1 - coupon.getWeightEnd()) * coupon.getReferenceEndTime()[1]; final double firstNaturalPaymentTime = coupon.getNaturalPaymentStartTime(); final double secondNaturalPaymentTime = coupon.getNaturalPaymentEndTime(); final double paymentTime = coupon.getPaymentTime(); final double volatilityStart = inflationConvexity.getInflationConvexityAdjustmentParameters() .getPriceIndexAtmVolatility()[0]; final double volatilityEnd = inflationConvexity.getInflationConvexityAdjustmentParameters() .getPriceIndexAtmVolatility()[1]; final double correlationInflation = inflationConvexity.getInflationConvexityAdjustmentParameters() .getPriceIndexCorrelation().getZValue(firstFixingTime, secondFixingTime); final double correlationInflationRateStart = inflationConvexity.getInflationConvexityAdjustmentParameters() .getPriceIndexRateCorrelation().getYValue(firstFixingTime); final double correlationInflationRateEnd = inflationConvexity.getInflationConvexityAdjustmentParameters() .getPriceIndexRateCorrelation().getYValue(secondFixingTime); final double volBondForwardStart = getVolBondForward(firstNaturalPaymentTime, paymentTime, inflationConvexity); final double volBondForwardEnd = getVolBondForward(secondNaturalPaymentTime, paymentTime, inflationConvexity); final double adjustment = volatilityStart * (volatilityStart - volatilityEnd * correlationInflation - volBondForwardStart * correlationInflationRateStart) * firstNaturalPaymentTime + volatilityEnd * volBondForwardEnd * correlationInflationRateEnd * secondNaturalPaymentTime; return Math.exp(adjustment); }
From source file:etomica.models.co2.PNGCPM.java
/** * This returns the pairwise-additive portion of the GCPM potential for a * pair of atoms (dispersion + fixed-charge electrostatics) *///from w w w. ja va 2 s .com public double getNonPolarizationEnergy(IMoleculeList molecules) { IAtomList atoms1 = molecules.getMolecule(0).getChildList(); IAtomList atoms2 = molecules.getMolecule(1).getChildList(); IVectorMutable C1r = atoms1.getAtom(0).getPosition(); IVectorMutable C2r = atoms2.getAtom(0).getPosition(); work.Ev1Mv2(C1r, C2r); shift.Ea1Tv1(-1, work); boundary.nearestImage(work); shift.PE(work); final boolean zeroShift = shift.squared() < 0.1; double r2 = work.squared(); double sum = 0; if (zeroShift) { for (int i = 0; i < atoms1.getAtomCount(); i++) { for (int j = 0; j < atoms2.getAtomCount(); j++) { GCPMAgent pairAgent = getPairAgent(atoms1.getAtom(i).getType(), atoms2.getAtom(j).getType()); double epsilon = pairAgent.epsilon; double r = Double.NaN; if (epsilon > 0) { double sigma = pairAgent.sigma; double gamma = pairAgent.gamma; r2 = atoms1.getAtom(i).getPosition().Mv1Squared(atoms2.getAtom(j).getPosition()); r = Math.sqrt(r2); double rOverSigma = r / sigma; double sigma2OverR2 = 1 / (rOverSigma * rOverSigma); if (1 / sigma2OverR2 < coreFac) return Double.POSITIVE_INFINITY; double sixOverGamma = 6 / gamma; sum += epsilon / (1 - sixOverGamma) * (sixOverGamma * Math.exp(gamma * (1 - rOverSigma)) - sigma2OverR2 * sigma2OverR2 * sigma2OverR2); } double charge2 = pairAgent.charge2; if (charge2 != 0) { double tau = pairAgent.tau; if (Double.isNaN(r)) { r2 = atoms1.getAtom(i).getPosition().Mv1Squared(atoms2.getAtom(j).getPosition()); r = Math.sqrt(r2); } sum += charge2 / r * (1 - org.apache.commons.math3.special.Erf.erfc(Math.sqrt(r2) / (2 * tau))); } } } } else { for (int i = 0; i < atoms1.getAtomCount(); i++) { for (int j = 0; j < atoms2.getAtomCount(); j++) { GCPMAgent pairAgent = getPairAgent(atoms1.getAtom(i).getType(), atoms2.getAtom(j).getType()); double epsilon = pairAgent.epsilon; double r = Double.NaN; if (epsilon > 0) { double sigma = pairAgent.sigma; IVector r1 = atoms1.getAtom(i).getPosition(); shift.PE(r1); r2 = atoms2.getAtom(j).getPosition().Mv1Squared(shift); shift.ME(r1); r = Math.sqrt(r2); double gamma = pairAgent.gamma; double rOverSigma = r / sigma; double sigma2OverR2 = 1 / (rOverSigma * rOverSigma); if (1 / sigma2OverR2 < coreFac) return Double.POSITIVE_INFINITY; double sixOverGamma = 6 / gamma; sum += epsilon / (1 - sixOverGamma) * (sixOverGamma * Math.exp(gamma * (1 - rOverSigma)) - sigma2OverR2 * sigma2OverR2 * sigma2OverR2);//exp-6 potential(Udisp) } double charge2 = pairAgent.charge2; if (charge2 != 0) { double tau = pairAgent.tau; if (Double.isNaN(r)) { IVector r1 = atoms1.getAtom(i).getPosition(); shift.PE(r1); r2 = atoms2.getAtom(j).getPosition().Mv1Squared(shift); shift.ME(r1); r = Math.sqrt(r2); } sum += charge2 / r * (1 - org.apache.commons.math3.special.Erf.erfc(Math.sqrt(r2) / (2 * tau))); } } } } return sum; }
From source file:com.opengamma.analytics.financial.interestrate.payments.method.CapFloorCMSHullWhiteApproximationMethod.java
public CurrencyAmount presentValue(final CapFloorCMS cms, final HullWhiteOneFactorPiecewiseConstantDataBundle hwData) { Validate.notNull(cms);//from w ww . ja v a 2 s . co m Validate.notNull(hwData); final double expiryTime = cms.getFixingTime(); final SwapFixedCoupon<? extends Payment> swap = cms.getUnderlyingSwap(); final double dfPayment = hwData.getCurve(swap.getFirstLeg().getDiscountCurve()) .getDiscountFactor(cms.getPaymentTime()); final int nbFixed = cms.getUnderlyingSwap().getFixedLeg().getNumberOfPayments(); final double[] alphaFixed = new double[nbFixed]; final double[] dfFixed = new double[nbFixed]; final double[] discountedCashFlowFixed = new double[nbFixed]; for (int loopcf = 0; loopcf < nbFixed; loopcf++) { alphaFixed[loopcf] = MODEL.alpha(hwData.getHullWhiteParameter(), 0.0, expiryTime, expiryTime, swap.getFixedLeg().getNthPayment(loopcf).getPaymentTime()); dfFixed[loopcf] = hwData.getCurve(swap.getFixedLeg().getNthPayment(loopcf).getFundingCurveName()) .getDiscountFactor(swap.getFixedLeg().getNthPayment(loopcf).getPaymentTime()); discountedCashFlowFixed[loopcf] = dfFixed[loopcf] * swap.getFixedLeg().getNthPayment(loopcf).getPaymentYearFraction() * swap.getFixedLeg().getNthPayment(loopcf).getNotional(); } final AnnuityPaymentFixed cfeIbor = swap.getSecondLeg().accept(CFEC, hwData); final double[] alphaIbor = new double[cfeIbor.getNumberOfPayments()]; final double[] dfIbor = new double[cfeIbor.getNumberOfPayments()]; final double[] discountedCashFlowIbor = new double[cfeIbor.getNumberOfPayments()]; for (int loopcf = 0; loopcf < cfeIbor.getNumberOfPayments(); loopcf++) { alphaIbor[loopcf] = MODEL.alpha(hwData.getHullWhiteParameter(), 0.0, expiryTime, expiryTime, cfeIbor.getNthPayment(loopcf).getPaymentTime()); dfIbor[loopcf] = hwData.getCurve(cfeIbor.getDiscountCurve()) .getDiscountFactor(cfeIbor.getNthPayment(loopcf).getPaymentTime()); discountedCashFlowIbor[loopcf] = dfIbor[loopcf] * cfeIbor.getNthPayment(loopcf).getAmount(); } final double alphaPayment = MODEL.alpha(hwData.getHullWhiteParameter(), 0.0, expiryTime, expiryTime, cms.getPaymentTime()); final double x0 = -alphaPayment; final double a0 = MODEL.swapRate(x0, discountedCashFlowFixed, alphaFixed, discountedCashFlowIbor, alphaIbor) - cms.getStrike(); final double a1 = MODEL.swapRateDx1(x0, discountedCashFlowFixed, alphaFixed, discountedCashFlowIbor, alphaIbor); final double a2 = MODEL.swapRateDx2(x0, discountedCashFlowFixed, alphaFixed, discountedCashFlowIbor, alphaIbor); // AnnuityPaymentFixed cfe = CFEC.visit(swap.withCoupon(cms.getStrike()), hwData); // double[] alpha = new double[cfe.getNumberOfPayments()]; // double[] df = new double[cfe.getNumberOfPayments()]; // double[] discountedCashFlow = new double[cfe.getNumberOfPayments()]; // for (int loopcf = 0; loopcf < cfe.getNumberOfPayments(); loopcf++) { // alpha[loopcf] = MODEL.alpha(hwData.getHullWhiteParameter(), 0.0, expiryTime, expiryTime, cfe.getNthPayment(loopcf).getPaymentTime()); // df[loopcf] = hwData.getCurve(cfe.getDiscountCurve()).getDiscountFactor(cfe.getNthPayment(loopcf).getPaymentTime()); // discountedCashFlow[loopcf] = df[loopcf] * cfe.getNthPayment(loopcf).getAmount(); // } // double kappaTest = MODEL.kappa(discountedCashFlow, alpha); final double kappa = -a0 / a1 - alphaPayment; // approximation final double kappatilde = kappa + alphaPayment; final double omega = (cms.isCap() ? 1.0 : -1.0); final double s2pi = 1.0 / Math.sqrt(2.0 * Math.PI); double pv = omega * (a0 + a2 / 2) * NORMAL.getCDF(-omega * kappatilde) + s2pi * Math.exp(-kappatilde * kappatilde / 2) * (a1 + a2 * kappatilde); pv *= dfPayment * cms.getNotional() * cms.getPaymentYearFraction(); return CurrencyAmount.of(cms.getCurrency(), pv); }
From source file:com.opengamma.analytics.financial.interestrate.swaption.method.SwaptionBermudaFixedIborHullWhiteNumericalIntegrationMethod.java
/** * Computes the present value of the Physical delivery swaption. * @param swaption The swaption.// w w w . ja v a 2 s . co m * @param hwData The Hull-White parameters and the curves. * @return The present value. */ public CurrencyAmount presentValue(final SwaptionBermudaFixedIbor swaption, final HullWhiteOneFactorPiecewiseConstantDataBundle hwData) { Validate.notNull(swaption); Validate.notNull(hwData); int nbExpiry = swaption.getExpiryTime().length; Validate.isTrue(nbExpiry > 1, "At least two expiry dates required for this method"); double tmpdb; YieldAndDiscountCurve discountingCurve = hwData .getCurve(swaption.getUnderlyingSwap()[0].getFirstLeg().getDiscountCurve()); double[] theta = new double[nbExpiry + 1]; // Extended expiry time (with 0). theta[0] = 0.0; System.arraycopy(swaption.getExpiryTime(), 0, theta, 1, nbExpiry); AnnuityPaymentFixed[] cashflow = new AnnuityPaymentFixed[nbExpiry]; for (int loopexp = 0; loopexp < nbExpiry; loopexp++) { cashflow[loopexp] = CFEC.visit(swaption.getUnderlyingSwap()[loopexp], hwData); } int[] n = new int[nbExpiry]; double[][][] alpha = new double[nbExpiry][][]; double[][][] alpha2 = new double[nbExpiry][][]; // alpha^2 for (int loopexp = 0; loopexp < nbExpiry; loopexp++) { n[loopexp] = cashflow[loopexp].getNumberOfPayments(); alpha[loopexp] = new double[loopexp + 1][]; alpha2[loopexp] = new double[loopexp + 1][]; for (int k = 0; k <= loopexp; k++) { alpha[loopexp][k] = new double[n[loopexp]]; alpha2[loopexp][k] = new double[n[loopexp]]; for (int l = 0; l < alpha[loopexp][k].length; l++) { alpha[loopexp][k][l] = MODEL.alpha(hwData.getHullWhiteParameter(), theta[k], theta[k + 1], theta[k + 1], cashflow[loopexp].getNthPayment(l).getPaymentTime()); alpha2[loopexp][k][l] = alpha[loopexp][k][l] * alpha[loopexp][k][l]; } } } int nbPoint2 = 2 * NB_POINT + 1; int[] startInt = new int[nbExpiry - 1]; int[] endInt = new int[nbExpiry - 1]; for (int i = 1; i < nbExpiry - 1; i++) { startInt[i] = 0; endInt[i] = nbPoint2 - 1; } startInt[0] = NB_POINT; endInt[0] = NB_POINT; double[][] t = new double[nbExpiry][]; // payment time double[][] dfS = new double[nbExpiry][]; // discount factor double[] beta = new double[nbExpiry]; double[][] h = new double[nbExpiry][]; double[][] sa2 = new double[nbExpiry][]; for (int loopexp = 0; loopexp < nbExpiry; loopexp++) { beta[loopexp] = MODEL.beta(hwData.getHullWhiteParameter(), theta[loopexp], theta[loopexp + 1]); t[loopexp] = new double[n[loopexp]]; dfS[loopexp] = new double[n[loopexp]]; h[loopexp] = new double[n[loopexp]]; sa2[loopexp] = new double[n[loopexp]]; for (int loopcf = 0; loopcf < n[loopexp]; loopcf++) { t[loopexp][loopcf] = cashflow[loopexp].getNthPayment(loopcf).getPaymentTime(); dfS[loopexp][loopcf] = discountingCurve.getDiscountFactor(t[loopexp][loopcf]); h[loopexp][loopcf] = (1 - Math.exp(-hwData.getHullWhiteParameter().getMeanReversion() * t[loopexp][loopcf])) / hwData.getHullWhiteParameter().getMeanReversion(); tmpdb = 0.0; for (int k = 0; k <= loopexp; k++) { tmpdb += alpha2[loopexp][k][loopcf]; } sa2[loopexp][loopcf] = tmpdb; } } double[] discountedCashFlowN = new double[n[nbExpiry - 1]]; for (int loopcf = 0; loopcf < n[nbExpiry - 1]; loopcf++) { discountedCashFlowN[loopcf] = dfS[nbExpiry - 1][loopcf] * cashflow[nbExpiry - 1].getNthPayment(loopcf).getAmount(); } double lambda = MODEL.lambda(discountedCashFlowN, sa2[nbExpiry - 1], h[nbExpiry - 1]); double[] betaSort = new double[nbExpiry]; System.arraycopy(beta, 0, betaSort, 0, nbExpiry); Arrays.sort(betaSort); double minbeta = betaSort[0]; double maxbeta = betaSort[nbExpiry - 1]; double b = Math.min(10 * minbeta, maxbeta); double epsilon = -2.0 / NB_POINT * NORMAL.getInverseCDF(1.0 / (200.0 * NB_POINT)) * b; // <- double[] bX = new double[nbPoint2]; for (int looppt = 0; looppt < nbPoint2; looppt++) { bX[looppt] = -NB_POINT * epsilon + looppt * epsilon; } double[] bX2 = new double[4 * NB_POINT + 1]; for (int looppt = 0; looppt < 4 * NB_POINT + 1; looppt++) { bX2[looppt] = -2 * NB_POINT * epsilon + looppt * epsilon; } double[] htheta = new double[nbExpiry]; for (int loopexp = 0; loopexp < nbExpiry; loopexp++) { htheta[loopexp] = (1 - Math.exp(-hwData.getHullWhiteParameter().getMeanReversion() * theta[loopexp + 1])) / hwData.getHullWhiteParameter().getMeanReversion(); } double[][] vZ = new double[nbExpiry - 1][nbPoint2]; for (int i = nbExpiry - 2; i >= 0; i--) { for (int looppt = 0; looppt < nbPoint2; looppt++) { vZ[i][looppt] = Math.exp(bX[looppt] * htheta[i]); } } double[][] vW = new double[nbExpiry][]; // Swaption double[][] vT = new double[nbExpiry - 1][]; // Swap double omega = -Math.signum(cashflow[nbExpiry - 1].getNthPayment(0).getAmount()); double[] kappaL = new double[nbPoint2]; for (int looppt = 0; looppt < nbPoint2; looppt++) { kappaL[looppt] = (lambda - bX[looppt]) / beta[nbExpiry - 1]; } double[] sa2N1 = new double[n[nbExpiry - 1]]; for (int i = 0; i < n[nbExpiry - 1]; i++) { tmpdb = 0; for (int k = 0; k <= nbExpiry - 2; k++) { tmpdb += alpha2[nbExpiry - 1][k][i]; } sa2N1[i] = tmpdb; } vW[nbExpiry - 1] = new double[4 * NB_POINT + 1]; for (int j = 0; j < n[nbExpiry - 1]; j++) { for (int looppt = 0; looppt < nbPoint2; looppt++) { vW[nbExpiry - 1][NB_POINT + looppt] += discountedCashFlowN[j] * Math.exp(-sa2N1[j] / 2.0 - h[nbExpiry - 1][j] * bX[looppt]) * NORMAL.getCDF(omega * (kappaL[looppt] + alpha[nbExpiry - 1][nbExpiry - 1][j])); } } for (int looppt = 0; looppt < NB_POINT; looppt++) { vW[nbExpiry - 1][looppt] = vW[nbExpiry - 1][NB_POINT]; } for (int looppt = 0; looppt < NB_POINT; looppt++) { vW[nbExpiry - 1][3 * NB_POINT + 1 + looppt] = vW[nbExpiry - 1][3 * NB_POINT]; } double c1sqrt2pi = 1.0 / Math.sqrt(2 * Math.PI); double[][] pvcfT = new double[nbExpiry - 1][]; double[] vL; // Left side of intersection double[] vR; // Right side of intersection double[][] labc; double[][] rabc; double[][] labcM = new double[3][4 * NB_POINT + 1]; double[][] rabcM = new double[3][4 * NB_POINT + 1]; double[] dabc = new double[3]; int[] indSwap = new int[nbExpiry - 1]; // index of the intersection double xroot; double[][] xN = new double[nbExpiry - 1][nbPoint2]; double ci; double coi; int is; double[] ncdf0 = new double[nbPoint2]; double[] ncdf1 = new double[nbPoint2]; double[] ncdf2 = new double[nbPoint2]; double[] ncdf0X = new double[nbPoint2 + 1]; double[] ncdf1X = new double[nbPoint2 + 1]; double[] ncdf2X = new double[nbPoint2 + 1]; double ncdf0x; double ncdf1x; double ncdf2x; double ncdfinit; // Main loop for the different expiry dates (except the last one) for (int i = nbExpiry - 2; i >= 0; i--) { vW[i] = new double[4 * NB_POINT + 1]; vT[i] = new double[4 * NB_POINT + 1]; // T: swap pvcfT[i] = new double[n[i]]; for (int j = 0; j < n[i]; j++) { pvcfT[i][j] = cashflow[i].getNthPayment(j).getAmount() * dfS[i][j]; for (int looppt = 0; looppt < 4 * NB_POINT + 1; looppt++) { vT[i][looppt] += pvcfT[i][j] * Math.exp(-sa2[i][j] / 2.0 - h[i][j] * bX2[looppt]); } } // Preparation for (int looppt = 0; looppt < nbPoint2; looppt++) { xN[i][looppt] = bX[looppt] / beta[i]; } ci = htheta[i] * beta[i]; coi = Math.exp(ci * ci / 2); // Left/Right if (omega < 0) { vL = vW[i + 1]; vR = vT[i]; } else { vR = vW[i + 1]; vL = vT[i]; } indSwap[i] = 0; while (vL[indSwap[i] + 1] >= vR[indSwap[i] + 1]) { indSwap[i]++; } // Parabola fit labc = parafit(epsilon / beta[i], vL); rabc = parafit(epsilon / beta[i], vR); for (int k = 0; k < 3; k++) { dabc[k] = labc[k][indSwap[i]] - rabc[k][indSwap[i]]; System.arraycopy(labc[k], 0, labcM[k], 0, indSwap[i] + 1); System.arraycopy(labc[k], indSwap[i], labcM[k], indSwap[i] + 1, labc[k].length - indSwap[i]); System.arraycopy(rabc[k], 0, rabcM[k], 0, indSwap[i] + 1); System.arraycopy(rabc[k], indSwap[i], rabcM[k], indSwap[i] + 1, rabc[k].length - indSwap[i]); } for (int looppt = 0; looppt < 4 * NB_POINT + 1; looppt++) { labcM[1][looppt] = labcM[1][looppt] + labcM[0][looppt] * 2 * ci; labcM[2][looppt] = labcM[2][looppt] + labcM[1][looppt] * ci - labcM[0][looppt] * ci * ci; rabcM[1][looppt] = rabcM[1][looppt] + rabcM[0][looppt] * 2 * ci; rabcM[2][looppt] = rabcM[2][looppt] + rabcM[1][looppt] * ci - rabcM[0][looppt] * ci * ci; } xroot = (-dabc[1] - Math.sqrt(dabc[1] * dabc[1] - 4 * dabc[0] * dabc[2])) / (2 * dabc[0]); ncdfinit = NORMAL.getCDF(xN[i][0]); for (int looppt = 0; looppt < nbPoint2; looppt++) { ncdf0[looppt] = NORMAL.getCDF(xN[i][looppt] - ci); ncdf1[looppt] = -c1sqrt2pi * Math.exp(-(xN[i][looppt] - ci) * (xN[i][looppt] - ci) / 2.0); ncdf2[looppt] = ncdf1[looppt] * (xN[i][looppt] - ci) + ncdf0[looppt]; } for (int j = startInt[i]; j <= endInt[i]; j++) { is = indSwap[i] - j + 1; // % all L if (j + 2 * NB_POINT <= indSwap[i]) { double[][] xabc = new double[3][2 * NB_POINT]; for (int k = 0; k < 3; k++) { System.arraycopy(labcM[k], j, xabc[k], 0, 2 * NB_POINT); } for (int looppt = 0; looppt < 2 * NB_POINT; looppt++) { xabc[1][looppt] = xabc[1][looppt] + xabc[0][looppt] * 2 * xN[i][j]; xabc[2][looppt] = xabc[2][looppt] + xabc[1][looppt] * xN[i][j] - xabc[0][looppt] * xN[i][j] * xN[i][j]; } vW[i][j + NB_POINT] = 0; vW[i][j + NB_POINT] = vW[i][j + NB_POINT] + coi * ni2ncdf(ncdf2, ncdf1, ncdf0, xabc); } else if (j < indSwap[i]) { double[][] xabc = new double[3][2 * NB_POINT + 1]; tmpdb = xroot - xN[i][j] - ci; ncdf0x = NORMAL.getCDF(tmpdb); ncdf1x = -Math.exp(-(tmpdb * tmpdb) / 2) * c1sqrt2pi; ncdf2x = ncdf1x * tmpdb + ncdf0x; for (int k = 0; k < 3; k++) { // System.arraycopy(rabcM[k], j, xabc[k], 0, 2 * _nbPoint + 1); // Swap System.arraycopy(labcM[k], j, xabc[k], 0, 2 * NB_POINT + 1); System.arraycopy(rabcM[k], indSwap[i] + 1, xabc[k], indSwap[i] + 1 - j, j + 2 * NB_POINT - indSwap[i]); } for (int looppt = 0; looppt < 2 * NB_POINT; looppt++) { xabc[1][looppt] = xabc[1][looppt] + xabc[0][looppt] * 2 * xN[i][j]; xabc[2][looppt] = xabc[2][looppt] + xabc[1][looppt] * xN[i][j] - xabc[0][looppt] * xN[i][j] * xN[i][j]; } System.arraycopy(ncdf0, 0, ncdf0X, 0, is); ncdf0X[is] = ncdf0x; System.arraycopy(ncdf0, is, ncdf0X, is + 1, ncdf0.length - is); System.arraycopy(ncdf1, 0, ncdf1X, 0, is); ncdf1X[is] = ncdf1x; System.arraycopy(ncdf1, is, ncdf1X, is + 1, ncdf1.length - is); System.arraycopy(ncdf2, 0, ncdf2X, 0, is); ncdf2X[is] = ncdf2x; System.arraycopy(ncdf2, is, ncdf2X, is + 1, ncdf2.length - is); vW[i][j + NB_POINT] = vL[j] * vZ[i][0] * ncdfinit; vW[i][j + NB_POINT] += coi * ni2ncdf(ncdf2X, ncdf1X, ncdf0X, xabc); vW[i][j + NB_POINT] += vR[j + 2 * NB_POINT] * vZ[i][vZ[i].length - 1] * ncdfinit; } else { double[][] xabc = new double[3][2 * NB_POINT]; for (int k = 0; k < 3; k++) { System.arraycopy(rabcM[k], j + 1, xabc[k], 0, 2 * NB_POINT); // System.arraycopy(labcM[k], j + 1, xabc[k], 0, 2 * _nbPoint); // Swaption } for (int looppt = 0; looppt < 2 * NB_POINT; looppt++) { xabc[1][looppt] = xabc[1][looppt] + xabc[0][looppt] * 2 * xN[i][j]; xabc[2][looppt] = xabc[2][looppt] + xabc[1][looppt] * xN[i][j] - xabc[0][looppt] * xN[i][j] * xN[i][j]; } vW[i][j + NB_POINT] = vR[j] * vZ[i][0] * ncdfinit; vW[i][j + NB_POINT] += coi * ni2ncdf(ncdf2, ncdf1, ncdf0, xabc); vW[i][j + NB_POINT] += vR[j + 2 * NB_POINT] * vZ[i][vZ[i].length - 1] * ncdfinit; } } for (int j = 0; j < NB_POINT; j++) { // Flat extrapolation vW[i][j] = vW[i][NB_POINT]; vW[i][3 * NB_POINT + 1 + j] = vW[i][3 * NB_POINT]; } } // End main loop return CurrencyAmount.of(swaption.getUnderlyingSwap()[0].getFixedLeg().getCurrency(), vW[0][2 * NB_POINT] * (swaption.isLong() ? 1.0 : -1.0)); }