Last time in the Black Scholes series I wrote about how to write the Black Scholes algorithm in Java in a maintainable and readable way by decomposing its constituent parts for easy reference. This time I look at modifying the previous implementation to incorporate The Greeks. Though – this time – there’s isn’t any great degree of decomposition to do as each Greek has its own formula and each one is implemented separately anyway.

What are The Greeks and how do they relate to the Black Scholes algorithm? This is best answered by Risk Glossary. Below I merely present a decomposed implementation of Black Scholes and The Greeks. In terms of performance considerations I’ve ensured that no equivalent operations are called more than once. I’ve also eliminated object allocation completely – the only allocation made is of a result array which contains six doubles. However, I have not applied advanced optimisations such as the use of multiplication in place of the pow() function.

I want to add a disclaimer that there isn’t anything very new in this post that isn’t already out there. I’m sure this has been done lots of times out there but I’m posting this for two reasons: 1) I wanted to post my take on the implementation in a decomposed fashion 2) I didn’t really find anything out there implementing the Greeks in Java. I hope it helps others looking for something similar.

### The Greeks in Java

The inputs to the The Greeks formulas are the same as those available to the The Black Scholes algorithm so both can implemented in common scope. Most Greeks have different formulas for call and put options with the exception of gamma and vega which have been implemented common to call and put options. Theta was the most complex formula compared to the rest of them and so I’ve broken that one down into left and right halves. The calculate() method returns a six element double array which contains the following values in order: (1) price (2) delta (3) gamma (4) vega (5) theta (6) rho. If you want to know what the inputs are see my previous article and for the formulas implemented below look here.

[java]

package name.dhruba.black.scholes;

import static java.lang.Math.PI;

import static java.lang.Math.abs;

import static java.lang.Math.exp;

import static java.lang.Math.log;

import static java.lang.Math.pow;

import static java.lang.Math.sqrt;

public enum BlackScholesGreeks2 {

_;

private static final double P = 0.2316419;

private static final double B1 = 0.319381530;

private static final double B2 = -0.356563782;

private static final double B3 = 1.781477937;

private static final double B4 = -1.821255978;

private static final double B5 = 1.330274429;

public static double[] calculate(boolean c,

double s, double k, double r, double t, double v) {

double[] p = new double[6];

double d1 = d1(s, k, r, t, v);

double d2 = d2(s, k, r, t, v);

double sd1 = standardNormalDistribution(d1);

double cd1 = cumulativeDistribution(d1, sd1);

double thetaLeft = -(s * sd1 * v) / (2 * sqrt(t));

if (c) {

double cd2 = cumulativeDistribution(d2);

// price

p[0] = s * cd1 – k * exp(-r * t) * cd2;

// delta

p[1] = cd1;

// theta

double thetaRight = r * k * exp(-r * t) * cd2;

p[4] = thetaLeft – thetaRight;

// rho

p[5] = k * t * exp(-r * t) * cd2;

} else {

double pcd1 = cumulativeDistribution(-d1);

double pcd2 = cumulativeDistribution(-d2);

// price

p[0] = k * exp(-r * t) * pcd2 – s * pcd1;

// delta

p[1] = cd1 – 1;

// theta

double thetaRight = r * k * exp(-r * t) * pcd2;

p[4] = thetaLeft + thetaRight;

// rho

p[5] = -k * t * exp(-r * t) * pcd2;

}

// gamma

p[2] = sd1 / (s * v * sqrt(t));

// vega

p[3] = s * sd1 * sqrt(t);

return p;

}

private static double d1(double s, double k, double r, double t, double v) {

double top = log(s / k) + (r + pow(v, 2) / 2) * t;

double bottom = v * sqrt(t);

return top / bottom;

}

private static double d2(double s, double k, double r, double t, double v) {

return d1(s, k, r, t, v) – v * sqrt(t);

}

public static double cumulativeDistribution(double x) {

return cumulativeDistribution(x, standardNormalDistribution(x));

}

public static double cumulativeDistribution(double x, double sdx) {

double t = 1 / (1 + P * abs(x));

double t1 = B1 * pow(t, 1);

double t2 = B2 * pow(t, 2);

double t3 = B3 * pow(t, 3);

double t4 = B4 * pow(t, 4);

double t5 = B5 * pow(t, 5);

double b = t1 + t2 + t3 + t4 + t5;

double cd = 1 – sdx * b;

return x < 0 ? 1 – cd : cd;

}

public static double standardNormalDistribution(double x) {

double top = exp(-0.5 * pow(x, 2));

double bottom = sqrt(2 * PI);

return top / bottom;

}

}

[/java]

### Testing the implementation

I’ve written a simple test based on values assumed to be correct from an online calculator. Rather oddly for some of the values although my answers matched with that online calculator they differ from those on Wolfram Alpha very slightly. I don’t know why Wolfram Alpha is produced different values. If you know let me know.

[java]

package name.dhruba.black.scholes;

import static org.junit.Assert.assertEquals;

import org.junit.Test;

public class TestBlackScholesGreeks {

@Test

public void testGreeks() {

boolean c;

double s, k, r, t, v;

double[] p;

c = true;

s = 56.25;

k = 55;

r = 0.0285;

t = 0.34;

v = 0.28;

p = BlackScholesGreeks2.calculate(c, s, k, r, t, v);

assertEquals(4.561, round(p[0], 3), 0);

assertEquals(0.610, round(p[1], 3), 0);

assertEquals(0.042, round(p[2], 3), 0);

assertEquals(12.587, round(p[3], 3), 0);

assertEquals(-6.030, round(p[4], 3), 0);

assertEquals(10.110, round(p[5], 3), 0);

c = false;

p = BlackScholesGreeks2.calculate(c, s, k, r, t, v);

assertEquals(2.781, round(p[0], 3), 0);

assertEquals(-0.390, round(p[1], 3), 0);

assertEquals(0.042, round(p[2], 3), 0);

assertEquals(12.587, round(p[3], 3), 0);

assertEquals(-4.478, round(p[4], 3), 0);

assertEquals(-8.409, round(p[5], 3), 0);

}

static double round(double d, int places) {

int factor = (int) Math.pow(10, places);

return (double) Math.round(d * factor) / factor;

}

}

[/java]

### Oddities

One thing I noticed is that for the price formula for a put option reordering certain operations, even though the overall operation was equivalent, produced different digits towards the end of the value; in other words, towards the final decimal places. The following two operations, for example, although equivalent, produce different digits for the final few decimal places. If anyone knows what’s going on here do let me know.

p[0] = pcd2 * k * exp(-r * t) - pcd1 * s; p[0] = k * exp(-r * t) * pcd2 - s * pcd1;

Did this help you or do you have any improvements or fixes? Let me know in the comments.

Is your theta correct?

I tried to implement this using Apple Sept 660 Strike (S=660.59, K=660, T=10/255, R=0.01, V=28%).

Time is 10 days, with 255 trading days in the year.

Your implementation gives me the following:

d=0.52 g=0.011 t=-189.357 v=52.121 r=12.889

Checking against ThinkOrSwim, Delta and Gamma look good. Theta is looks like it is off by several magnitudes.

I used the following only calculator to double check

http://www.volatilitytrading.net/black_scholes_calculator.htm

This return a theta of -0.666.

Did I make a mistake or is the theta incorrect?

Thanks for the great detailed blog post

by definition, theta is time decay over a year, what TOS shows is daily theta, divide by 365 and it matches (decay happens on weekends too, so I think T should be 10/365 also)

Tried it with the OCT 13 AAPL 435 Call. After dividing by 365, code above gave -0,122. TOS show-0.09 and OptionsHouse shows -0.14.

For something closer, APR4 13 AAPL 430 call, calc is -.51, OH is -.61 and TOS is -.45

Close enough for my usual contract size (10-20), much closer than OH (who it looks like is dividing by 255 instead of 365).

Order of operations affecting the last digit of the output is normal. Floating point numbers have finite precision, so the result of most calculations must be rounded to the nearest number within the set of all values a double can take. Essentially unless you take very careful steps to cancel it out, there is a pseudorandom noise component of +-(half the difference between adjacent values at the current magnitude) introduced in every calculation. Changing the order of operations changes the intermediate values, which changes the rounding performed, which changes the errors. This is usually confined to the last digit or two in a simple operation but can be magnified into larger errors by long calculation chains.

Thanks for your posting. One question if you have time…

Why are all the implementations (including this one) based on implied volatility, when implied volatility must either be estimated or inferred from this formula. wouldn’t it make more sense to use the option market price as an input and include implied volatility along with the other greeks as an output? I can’t imagine noone has thought of this, but I don’t see anywhere a reason why this isn’t the norm.

Thanks,

Jonathan