Implied volatility for American Options

func_names.cpp

/*
 * func_names.cpp
 *
 *  Created on: Dec 1, 2011
 *      Author: tlmaloney
 */

#include "func_names.h"

// boundary condition along tau_final
double f(double x, double strike, double a, bool call)
{
	if (call)
		return double(strike) * exp(a * x) * max(exp(x) - 1, 0.0);
	else
		return double(strike) * exp(a * x) * max(1 - exp(x), 0.0);
}

// boundary condition along x_left
double g_left(double t, double strike, double dividend, double rate, double sigma, double a,
		double b, double x_left, double x_right, bool call)
{
	if (call)
		return 0;
	else
		return strike * exp(a * x_left + b * t)
				* (exp(-2 * rate * t / pow(sigma, 2))
						- exp(x_left - 2.0 * dividend * t / pow(sigma, 2)));
}

// boundary condition along x_right
double g_right(double t, double strike, double dividend, double rate, double sigma, double a,
		double b, double x_left, double x_right, bool call)
{
	if (call)
		return strike * exp(a * x_right + b * t)
				* (exp(x_right - 2 * dividend * t / pow(sigma, 2))
						- exp(-2.0 * rate * t / pow(sigma, 2)));
	else
		return 0;
}

// boundary condition along x_left
double g_left_amer(double t, double strike, double dividend, double rate, double sigma, double a,
		double b, double x_left, double x_right, bool call)
{
	if (call)
		return 0;
	else
		return strike
				* max(
						exp(x_left) - 1,
						exp(a * x_left + b * t)
								* (exp(-2 * rate * t / pow(sigma, 2))
										- exp(x_left - 2.0 * dividend * t / pow(sigma, 2))));
}

// boundary condition along x_right
double g_right_amer(double t, double strike, double dividend, double rate, double sigma, double a,
		double b, double x_left, double x_right, bool call)
{
	if (call)
		return strike
				* max(
						1 - exp(x_right),
						exp(a * x_right + b * t)
								* (exp(x_right - 2 * dividend * t / pow(sigma, 2))
										- exp(-2.0 * rate * t / pow(sigma, 2))));
	else
		return 0;
}

// maximum pointwise error
double pointwise_error(Matrix &u, double a, double x_compute, double b, double tau_final, int M,
		int N_left, double t, double S, double K, double T, double vol, double intRate,
		double divRate, bool call)
{
	double error;
	double V_exact, V_approx;

	if (call)
		V_exact = black_scholes(t, S, K, T, vol, intRate, divRate, true);
	else
		V_exact = black_scholes(t, S, K, T, vol, intRate, divRate, false);

	V_approx = exp(-a * x_compute - b * tau_final) * u(M, N_left);

	error = fabs(V_approx - V_exact);

	return error;
}

// root-mean-squared error
double rms_error(Matrix &u, int N, double x_left, double dx, double tau_final, double a, double b,
		int M, double t, double S, double K, double T, double vol, double intRate, double divRate,
		bool call)
{
	double sum = 0;
	double V_approx, V_exact;
	double x_k;

	for (int k = 0; k < N + 1; k++)
	{
		x_k = x_left + k * dx;
		V_approx = exp(-a * x_k - b * tau_final) * u(M, k);

		if (call)
			V_exact = black_scholes(t, K * exp(x_k), K, T, vol, intRate, divRate, true);
		else
			V_exact = black_scholes(t, K * exp(x_k), K, T, vol, intRate, divRate, false);

		if (V_exact > 0.0001)
			sum += pow(V_approx - V_exact, 2) / pow(V_exact, 2);
	}

	sum /= N + 1;

	sum = sqrt(sum);

	return sum;
}

main.cpp

#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>
#include <cstdlib>
#include <string>
#include <cstring>
#include <sstream>
#include <fstream>
#include "Matrix.h"
#include "pseudocodes.h"
#include "matrix_codes.h"
#include "func_names.h"
using namespace std;

vector<double> V_diff(double sigma, bool isCall);

int main()
{
	// Set output formatting
	cout.setf(ios::showpoint);
	cout.setf(ios::fixed, ios::floatfield);
	cout.precision(15);

	bool isCall = false;
	double sigma_0 = 0.1;
	double sigma_1 = 0.4;

	double xOldest;
	double xNew = sigma_1;
	double xOld = sigma_0;
	double TOL_APPROX = 1e-15;
	double TOL_CONSEC = 1e-15;
	int iter = 1;

	cout << "sigma" << "\t" << "V_approx" << "\t"
			<< "Error" << "\t" << "x_left"
			<< "\t" << "x_right" << "\t" << "tau_final"
			<< "\t" << "N_left" << "\t" << "N_right"
			<< "\t" << "N" << endl;
	cout << xOld << "\t" << V_diff(xOld, isCall).at(1) << "\t"
			<< fabs(V_diff(xOld, isCall).at(0)) << "\t" << V_diff(xOld, isCall).at(2)
			<< "\t" << V_diff(xOld, isCall).at(3) << "\t" << V_diff(xOld, isCall).at(4)
			<< "\t" << (int) V_diff(xOld, isCall).at(5) << "\t" << (int) V_diff(xOld, isCall).at(6)
			<< "\t" << (int) V_diff(xOld, isCall).at(7) << endl;
	cout << xNew << "\t" << V_diff(xNew, isCall).at(1) << "\t"
			<< fabs(V_diff(xNew, isCall).at(0)) << "\t" << V_diff(xNew, isCall).at(2)
			<< "\t" << V_diff(xNew, isCall).at(3) << "\t" << V_diff(xNew, isCall).at(4)
			<< "\t" << (int) V_diff(xNew, isCall).at(5) << "\t" << (int) V_diff(xNew, isCall).at(6)
			<< "\t" << (int) V_diff(xNew, isCall).at(7) << endl;

	while (fabs(V_diff(xNew, isCall).at(0)) > TOL_APPROX || fabs(xNew - xOld) > TOL_CONSEC)
	{

		xOldest = xOld;
		xOld = xNew;
		xNew = xOld
				- V_diff(xOld, isCall).at(0) * (xOld - xOldest)
						/ (V_diff(xOld, isCall).at(0) - V_diff(xOldest, isCall).at(0));

		++iter;
		if (iter < 11)
		{
			cout << xNew << "\t" << V_diff(xNew, isCall).at(1) << "\t"
					<< fabs(V_diff(xNew, isCall).at(0)) << "\t" << V_diff(xNew, isCall).at(2)
					<< "\t" << V_diff(xNew, isCall).at(3) << "\t" << V_diff(xNew, isCall).at(4)
					<< "\t" << (int) V_diff(xNew, isCall).at(5) << "\t" << (int) V_diff(xNew, isCall).at(6)
					<< "\t" << (int) V_diff(xNew, isCall).at(7) << endl;
		}
	}

	return 0;
}

vector<double> V_diff(double sigma, bool isCall)
{
	double stock = 38.;
	double strike = 40.;
	double rate = 0.04;
	double maturity = 1.0 / 12.0;
	double dividend = 0.01;

	double a = (rate - dividend) / (sigma * sigma) - 0.5;
	double b = (a + 1.) * (a + 1.) + 2. * dividend / (sigma * sigma);

	/*************************************************************************/
	/* Forward Euler with M = 16 and alpha = 0.4 */
	/*************************************************************************/

	/*************************************************************************/
	/* Computational domain on the interval [0, tau_final] */
	/*************************************************************************/

	// American option -- x_compute on the grid, temporary left and right boundaries
	// like HW 9
	int M = 16;
	double alpha = 0.4;

	// x_compute
	double x_compute = log(stock / strike);

	// tau_final
	double tau_final = maturity * sigma * sigma / 2.0;

	// dtau and dx
	double dtau = tau_final / static_cast<double>(M);
	double dx = sqrt(dtau / alpha);

	// calculate x_left_tilde
	double x_left = x_compute + (rate - dividend - sigma * sigma / 2.) * maturity
			- 3. * sigma * sqrt(maturity);
	// calculate x_right-tilde
	double x_right = x_compute + (rate - dividend - sigma * sigma / 2.) * maturity
			+ 3. * sigma * sqrt(maturity);

	int N_left = ceil((x_compute - x_left) / dx);
	int N_right = ceil((x_right - x_compute) / dx);
	int N = N_left + N_right;

	x_left = x_compute - N_left * dx;
	x_right = x_compute + N_right * dx;

	/*************************************************************************/
	/* PDE to solve on the interval [0, tau_final] */
	/*************************************************************************/

	// u matrix with solution
	Matrix u(M + 1, N + 1);

	// fill in boundaries
	// t = 0
	for (int i = 0; i <= N; ++i)
	{
		u(0, i) = f(x_left + i * dx, strike, a, isCall);
	}
	// boundaries x = x_left and x = x_right
	for (int i = 0; i <= M; ++i)
	{
		u(i, 0) = g_left_amer(i * dtau, strike, dividend, rate, sigma, a, b, x_left, x_right,
				isCall);
		u(i, N) = g_right_amer(i * dtau, strike, dividend, rate, sigma, a, b, x_left, x_right,
				isCall);
	}

	/*************************************************************************/
	// Forward Euler
	for (int m = 0; m < M; m++)
		for (int n = 1; n < N; n++)
		{
			//american specifics
			double x = x_left + n * dx;
			double t = m * dtau;
			double early_ex_premium =
					isCall ?
							strike * exp(a * x + b * t) * max(exp(x) - 1., 0.) :
							strike * exp(a * x + b * t) * max(1. - exp(x), 0.);

			double u_eu = alpha * u(m, n - 1) + (1. - 2. * alpha) * u(m, n) + alpha * u(m, n + 1);
			u(m + 1, n) = max(u_eu, early_ex_premium);
		}
	/*************************************************************************/

	/*************************************************************************/
	/* Record values */
	/*************************************************************************/

	// calculate V_approx
	double V_approx = exp(-a * x_compute - b * tau_final) * u(M, N_left);
	double V_amer = 2.45;

	// ouput
	vector<double> output(10, 0.0);
	output.at(0) = V_approx - V_amer;
	output.at(1) = V_approx;
	output.at(2) = x_left;
	output.at(3) = x_right;
	output.at(4) = tau_final;
	output.at(5) = N_left;
	output.at(6) = N_right;
	output.at(7) = N;
	return output;
}