final double /* @DiscountFactor */dividendDiscount = process.dividendYield().currentLink().discount(ex.lastDate());
final double /* @DiscountFactor */riskFreeDiscount = process.riskFreeRate().currentLink().discount(ex.lastDate());
final double /* @Real */spot = process.stateVariable().currentLink().value();
QL.require(spot > 0.0, "negative or null underlying given"); // QA:[RG]::verified // TODO: message
final double /* @Real */forwardPrice = spot * dividendDiscount / riskFreeDiscount;
final BlackCalculator black = new BlackCalculator(payoff, forwardPrice, Math.sqrt(variance), riskFreeDiscount);
if (dividendDiscount>=1.0 && payoff.optionType()==Option.Type.Call) {
// early exercise never optimal
r.value = black.value();
greeks.delta = black.delta(spot);
moreGreeks.deltaForward = black.deltaForward();
moreGreeks.elasticity = black.elasticity(spot);
greeks.gamma = black.gamma(spot);
final DayCounter rfdc = process.riskFreeRate().currentLink().dayCounter();
final DayCounter divdc = process.dividendYield().currentLink().dayCounter();
final DayCounter voldc = process.blackVolatility().currentLink().dayCounter();
double /*@Time*/ t = rfdc.yearFraction(process.riskFreeRate().currentLink().referenceDate(), a.exercise.lastDate());
greeks.rho = black.rho(t);
t = divdc.yearFraction(process.dividendYield().currentLink().referenceDate(), a.exercise.lastDate());
greeks.dividendRho = black.dividendRho(t);
t = voldc.yearFraction(process.blackVolatility().currentLink().referenceDate(), a.exercise.lastDate());
greeks.vega = black.vega(t);
greeks.theta = black.theta(spot, t);
moreGreeks.thetaPerDay = black.thetaPerDay(spot, t);
moreGreeks.strikeSensitivity = black.strikeSensitivity();
moreGreeks.itmCashProbability = black.itmCashProbability();
} else {
// early exercise can be optimal
final CumulativeNormalDistribution cumNormalDist = new CumulativeNormalDistribution();
final NormalDistribution normalDist = new NormalDistribution();
final double /*@Real*/ tolerance = 1e-6;
final double /*@Real*/ Sk = BaroneAdesiWhaleyApproximationEngine.criticalPrice(payoff, riskFreeDiscount, dividendDiscount, variance, tolerance);
final double /*@Real*/ forwardSk = Sk * dividendDiscount / riskFreeDiscount;
final double /*@Real*/ alpha = -2.0*Math.log(riskFreeDiscount)/(variance);
final double /*@Real*/ beta = 2.0*Math.log(dividendDiscount/riskFreeDiscount)/
(variance);
final double /*@Real*/ h = 1 - riskFreeDiscount;
double /*@Real*/ phi;
switch (payoff.optionType()) {
case Call:
phi = 1;
break;
case Put:
phi = -1;
break;
default:
throw new LibraryException(UNKNOWN_OPTION_TYPE); // QA:[RG]::verified
}
// TODO: study how zero interest rate could be handled
QL.ensure(h != 0.0 , DIVIDING_BY_ZERO_INTEREST_RATE); // QA:[RG]::verified
final double /* @Real */temp_root = Math.sqrt((beta - 1) * (beta - 1) + (4 * alpha) / h);
final double /* @Real */lambda = (-(beta - 1) + phi * temp_root) / 2;
final double /* @Real */lambda_prime = -phi * alpha / (h * h * temp_root);
final double /* @Real */black_Sk = BlackFormula.blackFormula(payoff.optionType(), payoff.strike(), forwardSk, Math.sqrt(variance)) * riskFreeDiscount;
final double /* @Real */hA = phi * (Sk - payoff.strike()) - black_Sk;
final double /* @Real */d1_Sk = (Math.log(forwardSk / payoff.strike()) + 0.5 * variance) / Math.sqrt(variance);
final double /* @Real */d2_Sk = d1_Sk - Math.sqrt(variance);
final double /* @Real */part1 = forwardSk * normalDist.op(d1_Sk) / (alpha * Math.sqrt(variance));
final double /* @Real */part2 = -phi * forwardSk * cumNormalDist.op(phi * d1_Sk) * Math.log(dividendDiscount) / Math.log(riskFreeDiscount);
final double /* @Real */part3 = +phi * payoff.strike() * cumNormalDist.op(phi * d2_Sk);
final double /* @Real */V_E_h = part1 + part2 + part3;
final double /* @Real */b = (1 - h) * alpha * lambda_prime / (2 * (2 * lambda + beta - 1));
final double /* @Real */c = -((1 - h) * alpha / (2 * lambda + beta - 1)) * (V_E_h / (hA) + 1 / h + lambda_prime / (2 * lambda + beta - 1));
final double /* @Real */temp_spot_ratio = Math.log(spot / Sk);
final double /* @Real */chi = temp_spot_ratio * (b * temp_spot_ratio + c);
if (phi * (Sk - spot) > 0) {
r.value = black.value() + hA * Math.pow((spot / Sk), lambda) / (1 - chi);
} else {
r.value = phi * (spot - payoff.strike());
}
final double /* @Real */temp_chi_prime = (2 * b / spot) * Math.log(spot / Sk);