// Newton Raphson algorithm for finding critical price Si
double /*@Real*/ Q, LHS, RHS, bi;
double /*@Real*/ forwardSi = Si * dividendDiscount / riskFreeDiscount;
double /*@Real*/ d1 = (Math.log(forwardSi/payoff.strike()) + 0.5*variance) / Math.sqrt(variance);
final CumulativeNormalDistribution cumNormalDist = new CumulativeNormalDistribution();
final double /*@Real*/ K = (riskFreeDiscount!=1.0 ? -2.0*Math.log(riskFreeDiscount)/ (variance*(1.0-riskFreeDiscount)) : 0.0);
final double /*@Real*/ temp = BlackFormula.blackFormula(payoff.optionType(), payoff.strike(), forwardSi, Math.sqrt(variance))*riskFreeDiscount;
switch (payoff.optionType()) {
case Call:
Q = (-(n-1.0) + Math.sqrt(((n-1.0)*(n-1.0)) + 4 * K)) / 2;
LHS = Si - payoff.strike();
RHS = temp + (1 - dividendDiscount * cumNormalDist.op(d1)) * Si / Q;
bi = dividendDiscount * cumNormalDist.op(d1) * (1 - 1/Q)
+ (1 - dividendDiscount * cumNormalDist.derivative(d1) / Math.sqrt(variance)) / Q;
while (Math.abs(LHS - RHS)/payoff.strike() > tolerance) {
Si = (payoff.strike() + RHS - bi * Si) / (1 - bi);
forwardSi = Si * dividendDiscount / riskFreeDiscount;
d1 = (Math.log(forwardSi/payoff.strike())+0.5*variance)/Math.sqrt(variance);
LHS = Si - payoff.strike();
final double /*@Real*/ temp2 = BlackFormula.blackFormula(payoff.optionType(), payoff.strike(), forwardSi, Math.sqrt(variance))*riskFreeDiscount;
RHS = temp2 + (1 - dividendDiscount * cumNormalDist.op(d1)) * Si / Q;
bi = dividendDiscount * cumNormalDist.op(d1) * (1 - 1 / Q)
+ (1 - dividendDiscount * cumNormalDist.derivative(d1) / Math.sqrt(variance)) / Q;
}
break;
case Put:
Q = (-(n-1.0) - Math.sqrt(((n-1.0)*(n-1.0)) + 4 * K)) / 2;
LHS = payoff.strike() - Si;
RHS = temp - (1 - dividendDiscount * cumNormalDist.op(-d1)) * Si / Q;
bi = -dividendDiscount * cumNormalDist.op(-d1) * (1 - 1/Q)
- (1 + dividendDiscount * cumNormalDist.derivative(-d1) / Math.sqrt(variance)) / Q;
while (Math.abs(LHS - RHS)/payoff.strike() > tolerance) {
Si = (payoff.strike() - RHS + bi * Si) / (1 + bi);
forwardSi = Si * dividendDiscount / riskFreeDiscount;
d1 = (Math.log(forwardSi/payoff.strike())+0.5*variance)/Math.sqrt(variance);
LHS = payoff.strike() - Si;
final double /*@Real*/ temp2 = BlackFormula.blackFormula(payoff.optionType(), payoff.strike(), forwardSi, Math.sqrt(variance))*riskFreeDiscount;
RHS = temp2 - (1 - dividendDiscount * cumNormalDist.op(-d1)) * Si / Q;
bi = -dividendDiscount * cumNormalDist.op(-d1) * (1 - 1 / Q)
- (1 + dividendDiscount * cumNormalDist.derivative(-d1) / Math.sqrt(variance)) / Q;
}
break;
default:
throw new LibraryException(UNKNOWN_OPTION_TYPE); // QA:[RG]::verified
}