A quadratic approximation to American prices

~bernt/gcc_prog/algoritms_v1/algoritms/node24.html

We now discuss an approximation to the option price of an American option on a commodity having a continuous payout, described in (BAW). Their solution technique finds an approximate solution to the differential equation

(7.1)

The case we look at here is an option written on a commodity paying as a continuous payout. Here is the (unknown) formula that determines the price of the contingent claim. For an European option this can be explicitly solved, with the adjusted Black Scholes formula as the solution, but for American options, which may be exercised early, there is no known analytical solution.

To do their approximation, BAW decomposes the American price into the European price and the early exercise premium

Here is the early exercise premium. This will also have to satisfy the partial differential equation. To come up with an approximation BAW transformed equation () into one where the terms involving are neglible, removed these, and ended up with a standard linear homeogenous second order equation, which has a known solution.

The functional form of the approximation is

where

and solves

In implementing this, the only problem is finding the critical value .

This is the classical problem of finding a root of the equation

This is solved using Newton's algorithm for finding the root.

We start by finding a first ``seed'' value . The next estimate of is found

At each step we need to evaluate and its derivative .

Call option

where is the Black Scholes value for commodities.

Similarly, for a put option the approximation is

and we solve for again by Newton's procedure, where now

---

// file approx_am_call.cc
// author: Bernt Arne Oedegaard
// implements the quadratic approximation of Barone-Adesi and Whaley
// described in Journal of Finance, june 87.

#include 
#include          
#include "normdist.h"               // normal distribution
#include "fin_algoritms.h"         // define other option pricing formulas

const double ACCURACY=1.0e-6;

double option_price_american_call_approximated_baw( double S,
						    double X, 
						    double r,
						    double b,
						    double sigma,
						    double time) 
{
    double sigma_sqr = sigma*sigma;
    double time_sqrt = sqrt(time);
    double nn = 2.0*b/sigma_sqr; 
    double m = 2.0*r/sigma_sqr;  
    double K = 1.0-exp(-r*time); 
    double q2 = (-(nn-1)+sqrt(pow((nn-1),2.0)+(4*m/K)))*0.5;

    // seed value from paper
    double q2_inf = 0.5 * ( (-nn-1.0) + sqrt(pow((nn-1),2.0)+4.0*m));
    double S_star_inf = X / (1.0 - 1.0/q2_inf);
    double h2 = -(b*time+2.0*sigma*time_sqrt)*(X/(S_star_inf-X));
    double S_seed = X + (S_star_inf-X)*(1.0-exp(h2));

    int no_iterations=0; // iterate on S to find S_star, using Newton steps
    double Si=S_seed;         
    double g=1;
    double gprime=1.0;
    while ((fabs(g) > ACCURACY) 
	   && (fabs(gprime)>ACCURACY) // to avoid exploding Newton's  
	   && ( no_iterations++<500)
	   && (Si>0.0)) {
	double c  = option_price_european_call_payout(Si,X,r,b,sigma,time);
	double d1 = (log(Si/X)+(b+0.5*sigma_sqr)*time)/(sigma*time_sqrt);
	g=(1.0-1.0/q2)*Si-X-c+(1.0/q2)*Si*exp((b-r)*time)*N(d1);
	gprime=( 1.0-1.0/q2)*(1.0-exp((b-r)*time)*N(d1))
	    +(1.0/q2)*exp((b-r)*time)*n(d1)*(1.0/(sigma*time_sqrt));
	Si=Si-(g/gprime); 
    };
    double S_star = 0;
    if (fabs(g)>ACCURACY) { S_star = S_seed; } // did not converge
    else { S_star = Si; };
    double C=0;
    double c  = option_price_european_call_payout(S,X,r,b,sigma,time);
    if (S>=S_star) {
	C=S-X;
    } 
    else {
	double d1 = (log(S_star/X)+(b+0.5*sigma_sqr)*time)/(sigma*time_sqrt);
	double A2 =  (1.0-exp((b-r)*time)*N(d1))* (S_star/q2); 
	C=c+A2*pow((S/S_star),q2);
    };
    return max(C,c); // know value will never be less than BS value
}