Figlewski and Gao (1999) set out the Adaptive Mesh Model (AMM) which dramatically reduces nonlinearity error by grafting one or more smaller sections of fine high-resolution lattice onto the subtrate trinomial tree with coarser time and price steps. This approach is highly innovative and piggybacks the Trinomial Lattice previously and introduces some improvement in performance. The AMM approach is sufficiently flexible to accommodate a wide variety of contingent claims. Accuracy increases by several orders of magnitude with no commensurate increase in execution time. Closed-form valuation solutions while useful tend to be arbitrarily limited by accuracy constraints and only exist for a small subset of derivative securities. Practitioners very often default to Binomial and Trinomial lattice frameworks. Alternatively, Monte Carlo is employed. These models are widely used because they are intuitive and accuracy can be controlled by altering computational intensity. Broadie and Detemple (1996) for example employed a 15,000-step binomial model to obtain “True values”.Amin and Khanna (1994) prove that the discrete-time models converge to the corresponding continuous-time models. This is important, as it provides a roadmap for establishing benchmarks. A number of variants of the original CRR trees have been proposed. Boyle (1986) proposed the trinomial tree in which each node generates three branches. The lattice framework is therefore acknowledged as a reliable workhorse and can serve to benchmark other techniques. Unfortunately, convergence can be non-monotonic, and the adapted mesh approach can serve mitigate these problems.
// Adaptative Mesh Method
// Trinomial Tree with one grafted level
// Requires Boost C++ libraries Version 1.54.0
// By Fabrice Douglas Rouah, FRouah.com and Volopta.com
#include <stdio.h>
#include <cstdlib>
#include <math.h>
#include <iostream>
#include <boost/numeric/ublas/matrix.hpp>
using namespace boost::numeric::ublas;
using namespace std;
double AMM(double Spot, double Strike, double T, double r, double v, int n, char OpType, char ExType) {
// Tree parameters
double alpha = r - 0.5*v*v;
double K = T / n;
double h = v*sqrt(3.0*K);
double pu = (v*v*K / h / h + alpha*alpha*K*K / h / h + alpha*K / h) / 2.0;
double pd = (v*v*K / h / h + alpha*alpha*K*K / h / h - alpha*K / h) / 2.0;
double pm = 1.0 - pu - pd;
double logK = log(Strike);
double exp_rT = exp(-r*T / n);
// Compute the Coarse Mesh
matrix<double> A(2 * n + 1, n + 1);
A(n, 0) = log(Spot);
for (int j = 1; j <= n; j++)
for (int i = n - j; i <= n + j; i++)
A(i, j) = A(n, 0) + (n - i)*h;
// Identify Cells for the Fine Mesh
int Locate;
matrix<double> Ind(8, 2);
for (int i = 0; i <= 2 * n - 1; i++)
if ((A(i, n - 1) >= logK) & (logK >= A(i + 1, n - 1))) Locate = i;
for (int i = 0; i <= 7; i++)
Ind(i, 1) = A(Locate - 3 + i, n);
for (int i = 2; i <= 5; i++)
Ind(i, 0) = Ind(i, 1);
// Compute the Fine Mesh 15 x 5, First Column
matrix<double> F(15, 5);
for (int i = 2; i <= 5; i++)
F(2 * i, 0) = Ind(i, 0);
// Compute the Fine Mesh 5 x 15, Last Column
for (int i = 0; i <= 7; i++)
F(2 * i, 4) = Ind(i, 1);
// Impute values for the Fine Mesh, Last Column
for (int i = 0; i <= 6; i++)
F(2 * i + 1, 4) = F(2 * i, 4) - h / 2.0;
// Impute values for the Fine Mesh, Remaining Columns
for (int j = 3; j >= 1; j--) {
for (int i = 4 - j; i <= 10 + j; i++)
if (F(i, j + 1) != 0) F(i, j) = F(i, j + 1);
}
// Convert to Prices
matrix<double> expF(15, 5);
matrix<double> S(2 * n + 1, n + 1);
for (int j = 0; j <= 4; j++) {
for (int i = 0; i <= 14; i++)
if (F(i, j) != 0) expF(i, j) = exp(F(i, j));
}
for (int j = 0; j <= n; j++) {
for (int i = 0; i <= 2 * n; i++)
if (A(i, j) != 0) S(i, j) = exp(A(i, j));
}
// Option Prices for Fine Mesh
matrix<double> FineOp(15, 5);
for (int i = 0; i <= 14; i++) {
if (OpType == 'C') FineOp(i, 4) = max(expF(i, 4) - Strike, 0.0);
else FineOp(i, 4) = max(Strike - expF(i, 4), 0.0);
}
for (int j = 3; j >= 0; j--) {
for (int i = 4 - j; i <= 10 + j; i++) {
switch (ExType) {
case 'E':
FineOp(i, j) = exp_rT*(pu*FineOp(i - 1, j + 1) + pm*FineOp(i, j + 1) + pd*FineOp(i + 1, j + 1));
break;
case 'A':
if (OpType == 'C')
FineOp(i, j) = max(expF(i, j) - Strike, exp_rT*(pu*FineOp(i - 1, j + 1) + pm*FineOp(i, j + 1) + pd*FineOp(i + 1, j + 1)));
else
FineOp(i, j) = max(Strike - expF(i, j), exp_rT*(pu*FineOp(i - 1, j + 1) + pm*FineOp(i, j + 1) + pd*FineOp(i + 1, j + 1)));
break;
}
}
}
// Option Prices for Coarse Mesh Last two Time Steps
matrix<double> Op(2 * n + 1, n + 1);
for (int i = 0; i <= 2 * n; i++) {
switch (OpType) {
case 'C': Op(i, n) = max(S(i, n) - Strike, 0.0); break;
case 'P': Op(i, n) = max(Strike - S(i, n), 0.0); break;
}
}
for (int i = 1; i <= 2 * n - 1; i++) {
switch (ExType) {
case 'E':
Op(i, n - 1) = exp_rT*(pu*Op(i - 1, n) + pm*Op(i, n) + pd*Op(i + 1, n));
break;
case 'A':
if (OpType == 'C')
Op(i, n - 1) = max(S(i, n - 1) - Strike, exp_rT*(pu*Op(i - 1, n) + pm*Op(i, n) + pd*Op(i + 1, n)));
else
Op(i, n - 1) = max(Strike - S(i, n - 1), exp_rT*(pu*Op(i - 1, n) + pm*Op(i, n) + pd*Op(i + 1, n)));
break;
}
}
// Impute Option Prices from Fine Mesh At Second to Last Step
Op(Locate - 1, n - 1) = FineOp(4, 0);
Op(Locate, n - 1) = FineOp(6, 0);
Op(Locate + 1, n - 1) = FineOp(8, 0);
Op(Locate + 2, n - 1) = FineOp(10, 0);
// Obtain remaining prices on Coarse mesh
for (int j = n - 2; j >= 0; j--) {
for (int i = n - j; i <= n + j; i++) {
switch (ExType) {
case 'E':
Op(i, j) = exp_rT*(pu*Op(i - 1, j + 1) + pm*Op(i, j + 1) + pd*Op(i + 1, j + 1));
break;
case 'A':
if (OpType == 'C')
Op(i, j) = max(S(i, j) - Strike, exp_rT*(pu*Op(i - 1, j + 1) + pm*Op(i, j + 1) + pd*Op(i + 1, j + 1)));
else
Op(i, j) = max(Strike - S(i, j), exp_rT*(pu*Op(i - 1, j + 1) + pm*Op(i, j + 1) + pd*Op(i + 1, j + 1)));
break;
}
}
}
// Return the option price
return Op(n, 0);
}
int main() {
double Spot = 100.0; // Spot Price
double Strike = 100.0; // Strike Price
double T = 1; // Years to Maturity
double r = 0.05; // Risk free Rate
double v = 0.2; // Volatility
int n = 1000; // Number of time steps
char OpType = 'C'; // 'P'ut or 'C'all
char ExType = 'E'; // 'E'uropean or 'A'merican
cout << "The AMM price is ";
cout << AMM(Spot, Strike, T, r, v, n, OpType, ExType) << endl;
system("PAUSE");
return 0.0;
}
// AMM Volopta without Turbo
// Adaptative Mesh Method
// Trinomial Tree with one grafted level
//Originally By Fabrice Douglas Rouah, FRouah.com and Volopta.com
//BB I made some small changes
#include <stdio.h>
#include <cstdlib>
#include <math.h>
#include <iostream>
#include <vector>
#include <cmath>
#include <string>
#include <algorithm>
#include<iomanip>
#include <time.h>
//double max(double a, double b)
//{
// return (a < b) ? b : a;
//}
//#include <boost/numeric/ublas/matrix.hpp>
//using namespace boost::numeric::ublas;
//vector<vector<double> > S(n+1, vector<double>(n+1,0));
using namespace std;
double AMM(double Spot, double Strike, double T, double r, double q, double v, int n, char OpType, char ExType) {
// Tree parameters
// was originally double alpha = r - 0.5*v*v;
// double alpha = r - (q + 0.5*v*v);
double alpha = r - (q + 0.5*v*v);
double K = T / n;
double h = v*sqrt(3.0*K);
double pu = (v*v*K / h / h + alpha*alpha*K*K / h / h + alpha*K / h) / 2.0;
double pd = (v*v*K / h / h + alpha*alpha*K*K / h / h - alpha*K / h) / 2.0;
double pm = 1.0 - pu - pd;
double logK = log(Strike);
double exp_rT = exp(-r*T / n);
// Compute the Coarse Mesh
//matrix<double> A(2 * n + 1, n + 1);
vector<vector<double>> A(2*n+1, vector<double>(n + 1, 0));
A[n][0] = log(Spot);
for (int j = 1; j <= n; j++)
for (int i = n - j; i <= n + j; i++)
A[i][j] = A[n][0] + (n - i)*h;
// Identify Cells for the Fine Mesh
int Locate;
//matrix<double> Ind(8, 2);
vector<vector<double>> Ind(8,vector<double>(2,0));
for (int i = 0; i <= 2 * n - 1; i++)
if ((A[i][n - 1] >= logK) & (logK >= A[i + 1][ n - 1])) Locate = i;
for (int i = 0; i <= 7; i++)
Ind[i][1] = A[Locate-3+i][n];
for (int i = 2; i <= 5; i++)
Ind[i][0] = Ind[i][1];
// Compute the Fine Mesh 15 x 5, First Column
//matrix<double> F(15, 5);
vector<vector<double>> F(15,vector<double>(5,0));
for (int i = 2; i <= 5; i++)
F[2 * i][0] = Ind[i][0];
// Compute the Fine Mesh 5 x 15, Last Column
for (int i = 0; i <= 7; i++)
F[2 * i][4] = Ind[i][1];
// Impute values for the Fine Mesh, Last Column
for (int i = 0; i <= 6; i++)
F[2 * i + 1][4] = F[2 * i][4] - h / 2.0;
// Impute values for the Fine Mesh, Remaining Columns
for (int j = 3; j >= 1; j--) {
for (int i = 4 - j; i <= 10 + j; i++)
if (F[i][ j + 1] != 0) F[i][j] = F[i][j + 1];
}
// Convert to Prices
//matrix<double> expF(15, 5);
vector<vector<double> > expF((15), vector<double>(5, 0));
//matrix<double> S(2 * n + 1, n + 1);
vector<vector<double> > S((2*n+1), vector<double>(n + 1, 0));
for (int j = 0; j <= 4; j++) {
for (int i = 0; i <= 14; i++)
if (F[i][j] != 0) expF[i][j] = exp(F[i][j]);
}
for (int j = 0; j <= n; j++) {
for (int i = 0; i <= 2 * n; i++)
if (A[i][j] != 0) S[i][j] = exp(A[i][j]);
}
// Option Prices for Fine Mesh
//matrix<double> FineOp(15, 5);
vector<vector<double> > FineOp((15), vector<double>(5, 0));
for (int i = 0; i <= 14; i++) {
if (OpType == 'C') FineOp[i][4] = max(expF[i][4] - Strike, 0.0);
else FineOp[i][4] = max(Strike - expF[i][4], 0.0);
}
for (int j = 3; j >= 0; j--) {
for (int i = 4 - j; i <= 10 + j; i++) {
switch (ExType) {
case 'E':
FineOp[i][j] = exp_rT*(pu*FineOp[i - 1][j + 1] + pm*FineOp[i][j + 1] + pd*FineOp[i + 1][j + 1]);
break;
case 'A':
if (OpType == 'C')
FineOp[i][j] = max(expF[i][j] - Strike, exp_rT*(pu*FineOp[i - 1] [j + 1] + pm*FineOp[i][j + 1] + pd*FineOp[i + 1][j + 1]));
else
FineOp[i][j] = max(Strike - expF[i][j], exp_rT*(pu*FineOp[i - 1][j + 1] + pm*FineOp[i][j + 1] + pd*FineOp[i + 1][j + 1]));
break;
}
}
}
// Option Prices for Coarse Mesh Last two Time Steps
//matrix<double> Op(2 * n + 1, n + 1);
vector<vector<double> > Op((2*n+1), vector<double>(n+1, 0));
for (int i = 0; i <= 2 * n; i++) {
switch (OpType) {
case 'C': Op[i][n] = max(S[i][n] - Strike, 0.0); break;
case 'P': Op[i][n] = max(Strike - S[i][n], 0.0); break;
}
}
for (int i = 1; i <= 2 * n - 1; i++) {
switch (ExType) {
case 'E':
Op[i][n-1] = exp_rT*(pu*Op[i - 1][n] + pm*Op[i][n] + pd*Op[i + 1][n]);
break;
case 'A':
if (OpType == 'C')
Op[i][n - 1] = max(S[i][n - 1] - Strike, exp_rT*(pu*Op[i - 1][n] + pm*Op[i][n] + pd*Op[i + 1][n]));
else
Op[i][n - 1] = max(Strike - S[i][n - 1], exp_rT*(pu*Op[i - 1][n] + pm*Op[i][n] + pd*Op[i + 1][n]));
break;
}
}
// Impute Option Prices from Fine Mesh At Second to Last Step
Op[Locate - 1][n - 1] = FineOp[4][0];
Op[Locate][n - 1] = FineOp[6][0];
Op[Locate + 1][n - 1] = FineOp[8][0];
Op[Locate + 2][n - 1] = FineOp[10][0];
// Obtain remaining prices on Coarse mesh
for (int j = n - 2; j >= 0; j--) {
for (int i = n - j; i <= n + j; i++) {
switch (ExType) {
case 'E':
Op[i][j] = exp_rT*(pu*Op[i - 1][j + 1] + pm*Op[i][j + 1] + pd*Op[i + 1][j + 1]);
break;
case 'A':
if (OpType == 'C')
Op[i][j] = max(S[i][j] - Strike, exp_rT*(pu*Op[i - 1][j + 1] + pm*Op[i][j + 1] + pd*Op[i + 1][j + 1]));
else
Op[i][j] = max(Strike - (S[i][j]), exp_rT*(pu*Op[i - 1][j + 1] + pm*Op[i][j + 1] + pd*Op[i + 1][j + 1]));
break;
}
}
}
// Return the option price
return Op[n][0];
}
int main() {
clock_t start_time, end_time;
start_time = clock();
double Spot = 100.0; // Spot Price
double Strike = 100.0; // Strike Price
double T = 3; // Years to Maturity
double r = 0.03; // Risk free Rate
double q = 0.07;
double v = 0.2; // Volatility
int n = 5000; // Number of time steps
char OpType = 'C'; // 'P'ut or 'C'all
char ExType = 'A'; // 'E'uropean or 'A'merican
cout << setprecision(10);
cout << "The AMM price is ";
cout << AMM(Spot, Strike, T, r, q, v, n, OpType, ExType) << endl;
end_time = clock();
cout << " " << start_time << " " << end_time << " " << (end_time - start_time) << endl;
cout << " " << start_time << " " << end_time << " " << (end_time - start_time) / (double)CLOCKS_PER_SEC << " seconds" << endl;
//system("PAUSE");
//return 0.0;
}