Newtons Method to find implied volatility

[riotburn]

I am trying to use Newtons method to find implied volatility of a Black Scholes Call. So given all the other variables, I need to find what volatility makes the black scholes call equal to the current call price. I'm running what I have below but keep getting a random string returned. I tested the black function and that is working right. Any thoughts?

http://en.wikipedia.org/wiki/Newton%27s_method
http://en.wikipedia.org/wiki/Black_scholes

// black_scholesC.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

double cum_dist_normal(double t);
double black(double S, double K, double T, double vol, double r, double q);
double vega(double S, double K, double T, double vol, double r, double q);

int _tmain(int argc, _TCHAR* argv[])
{
	double C = 2.50;
	double S = 30;
	double K = 30;
	double T = 180;
	double q = 0.01;
	double r = 0.03;
	double tol = 10E-6;
	double vol = 0.25;
	double xnew = vol;
	double xold = vol - 1.0;
	
	while (abs(xnew-xold) > tol)
	{
		xold = xnew;
		xnew = xnew - (black(S, K, T, xnew, r, q)/vega(S, K, T, xnew, r, q));
	}
	cout << xnew << endl;

	return 0;
}

double cum_dist_normal(double t)
{
	const double PI = 3.141592;
	double z = abs(t);
	double y = (1.0/(1.0 + 0.2316419*z));
	double a1 = 0.31938150;
	double a2 = -0.356563782;
	double a3 = 1.781477937;
	double a4 = -1.821255978;
	double a5 = 1.330274429;
	double m = 1.0 - ((exp(-(pow(t, 2.0)/2.0))*(a1*y + a2*pow(y, 2) + a3*pow(y, 3) + a4*pow(y, 4) + a5*pow(y, 5)))/sqrt(2*PI));
	double nm;
	if (t>0)
		nm = m;
	else
		nm = 1.0 - m;
	return nm;
}

double black(double S, double K, double T, double vol, double r, double q)
{
	double d1 = (log(S/K) + (r - q + (pow(vol, 2)/2)*T))/(vol*sqrt(T));
	double d2 = d1 - vol*sqrt(T);
	double cblack = S*exp(-(q*T))*cum_dist_normal(d1) - (K*exp(-(r*T))*cum_dist_normal(d2));
	return cblack;
}

double vega(double S, double K, double T, double vol, double r, double q)
{
	const double PI = 3.141592;
	double d1 = (log(S/K) + (r - q + (pow(vol, 2)/2)*T))/(vol*sqrt(T));
	double vegaC = (1.0/sqrt(2*PI))*(S*exp(-q*T))*(exp(-(pow(d1, 2))))*sqrt(T);
	return vegaC;
}

[HealBrains]

So you're saying your maths is all working a-ok, but when it comes to printing the string it's printing garbage?

What is an example of the output you're getting from 'xnew'?

[riotburn]

1.#QNAN

Well the math function is right but I am not so sure about the while function in the main. I was following a pseudocode for that.

[HealBrains]

Source: http://steve.hollasch.net/cgindex/coding/ieeefloat.html

The value NaN (Not a Number) is used to represent a value that does not represent a real number...

...A QNaN is a NaN with the most significant fraction bit set. QNaN's propagate freely through most arithmetic operations. These values pop out of an operation when the result is not mathematically defined. ...

...Semantically, QNaN's denote indeterminate operations, while SNaN's denote invalid operations.

----

Source: http://en.wikipedia.org/wiki/NaN

Ways to create NaNs:

The divisions 0/0, ∞/∞, ∞/-∞, -∞/∞, and -∞/-∞
The multiplications 0×∞ and 0×-∞
The power 1∞
The additions ∞ + (-∞), (-∞) + ∞ and equivalent subtractions.
Real operations with complex results:
The square root of a negative number
The logarithm of a negative number
The tangent of an odd multiple of 90 degrees (or π/2 radians)
The inverse sine or cosine of a number which is less than -1 or greater than +1.

I think the bold one is the one that is causing your problem. But I could be horribly wrong. (line 31). What numbers are being returned by your functions? Have you tried playing with the values you're passing into the function to get an expected result, then work back towards where you were to create the error and isolate it?

Chances are I'm just regurgitating stuff you've already seen, but perhaps someone here will see something that will help solve your problem.

Edit: Gosh, I feel like such a hack - just copy/pasting stuff from the net.

[riotburn]

I changed this initial vol = .50, and now I get a value. Not the right number but a value. I think I have to play around with it make sure I have the syntax. I am looking to for the program to find what vol plugged into the black function will make it equal C, which is 2.50.

[riotburn]

Found it. Left out one thing below. Thanks!

while (abs(xnew-xold) > tol)
{
xold = xnew;
xnew = xnew - ((black(S, K, T, xnew, r, q) - C)/vega(S, K, T, xnew, r, q));
}