Probability ITM formula for options

Question

Given a stock of price price and annual volatility annual_volatility, and given an option with strike price strike and expiry in calendar_days_remaining calendar days, I want to know the probability that it will expire in-the-money.

In other words, I need what you see as "probability ITM" in TOS or InteractiveBrokers.

So far I've found no answer for this on the internet, other than this calculator that thankfully works client side, so we can actually see the code (variable names are edited by me to try to make sense of this):

  price = form.price.value;
  strike = form.strike.value;
  calendar_years_remaining = calendar_days_remaining/365;
  annual_volatility = percent_annual_volatility/100;

  vt = annual_volatility*Math.sqrt(calendar_years_remaining);
  lnpq = Math.log(strike/price);
  d1 = lnpq / vt;

  y = Math.floor(1/(1+.2316419*Math.abs(d1))*100000)/100000;
  z = Math.floor(.3989423*Math.exp(-((d1*d1)/2))*100000)/100000;
  y5 = 1.330274*Math.pow(y,5);
  y4 = 1.821256*Math.pow(y,4);
  y3 = 1.781478*Math.pow(y,3);
  y2 = .356538*Math.pow(y,2);
  y1 = .3193815*y;
  x = 1-z*(y5-y4+y3-y2+y1);
  x = Math.floor(x*100000)/100000;

  if (d1<0) {x=1-x};

  pabove = Math.floor(x*1000)/10; 
  pbelow = Math.floor((1-x)*1000)/10;

Here we have, in pabove, the probability it'll expire above the strike price (and the opposite in pbelow).

My question is: why are there hardcoded numbers (1.330274, .3989423, etc)? What do they represent?

What is the actual formula to compute probability ITM and does it have hardcoded constants like the above script?

Answer 1

The formula with the hard coded constants is a polynomial approximation for the cumulative normal distribution function. This can be found in Abromowitz and Stegun, Handbook of Mathematical Functions, 1972. It is also used in Hull's Options, Futures, and Other Derivatives.

If you have the 4th edition of Hull, the forumulas for this approximation is on p.252.