Skip to content

Instantly share code, notes, and snippets.

@tlmaloney

tlmaloney/func_names.cpp

Created December 14, 2011 04:34
Show Gist options
  • Save tlmaloney/1475273 to your computer and use it in GitHub Desktop.
Save tlmaloney/1475273 to your computer and use it in GitHub Desktop.
Download ZIP
Implied volatility for American Options
Raw
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;
}
Raw
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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment