hw5_chengliangg.cpp

蒙特卡罗期权定价 布朗桥方法 分层样本方法

//HW 5 Math 623 
// Chengliang Ge


/*
AntiThetic.cpp
Arrays.cpp,  
ConvergenceTable.cpp, 
ExoticBSEngine.cpp
ExoticEngine.cpp
MCStatistics.cpp
Normals.cpp
Parameters.cpp,
ParkMiller.cpp
PathDependent.cpp
PathDependentAsian.cpp
PayOff3.cpp, 
PayOffBridge.cpp,
Random2.cpp,
SimpleMC8.cpp 
Vanilla3.cpp, 
*/
#include<SimpleMC8.h>
#include<ParkMiller.h>
#include<iostream>
#include<Vanilla3.h>
#include<MCStatistics.h>
#include<ConvergenceTable.h>
#include<AntiThetic.h>
#include <Normals.h>
#include <Random2.h>
#include <cmath>
#include <string>
#include <cstdlib>
#include <vector>
#include <nr3.h>
#include <sobseq.h>
#include<PathDependentAsian.h>
#include<ExoticBSEngine.h>

using namespace std;



///////////////////////////////////////////////////////////
//
//  GeometricAsianCF.H  close-form
//
////////////////////////////////////////////////////////////


class GeometricAsianCF
{
public:
	GeometricAsianCF(double Spot,        // constructor 
		double Strike,
		double r,
		double d,
		double Vol,
		double Expiry,
		double period)
	{
		S=Spot;
		K=Strike;
		rate=r;
		q=d;
		sigma=Vol;
		time=Expiry;
		N=period;
	}

	double option_price_asian_geometric_average_price_call()    //close form fuction
	{
		double deltaT=1/N;
		double sigma_sqr = pow(sigma,2);
		double adj_div_yield=rate-q-sigma_sqr/2.0;
		double V0=exp(-rate*time)*S*exp(0.5*(N+1)*adj_div_yield*deltaT + 1.0/12/N*(N+1)*(2*N+1)*sigma_sqr*deltaT);
		double adj_sigma=sigma*pow(N,-1.5)*pow((1.0/6)*N*(N+1)*(2*N+1),0.5);
		double adj_sigma_sqr = pow(adj_sigma,2);
		double time_sqrt = sqrt(time);
		double d1 = (log(V0/K) + (rate + 0.5*adj_sigma_sqr)*time)/(adj_sigma*time_sqrt);
		double d2 = d1-(adj_sigma*time_sqrt);
		double call_price = V0* CumulativeNormal(d1) - K * exp(-rate*time) * CumulativeNormal(d2);
		return call_price;
	};


private:

	double S;
	double K;
	double rate;
	double q;
	double sigma;
	double time;
	double N;
};




///////////////////////////////////////////////////////////
//
//  PathDependentAsianG_H
//
////////////////////////////////////////////////////////////

class PathDependentAsianG : public PathDependent   //inherite
{
public:

	PathDependentAsianG(const MJArray& LookAtTimes_, 
		double DeliveryTime_,
		const PayOffBridge& ThePayOff_);

	virtual unsigned long MaxNumberOfCashFlows() const;
	virtual MJArray PossibleCashFlowTimes() const;
	virtual unsigned long CashFlows(const MJArray& SpotValues, 
		std::vector<CashFlow>& GeneratedFlows) const;
	virtual ~PathDependentAsianG(){}
	virtual PathDependent* clone() const;

private:

	double DeliveryTime;
	PayOffBridge ThePayOff;
	unsigned long NumberOfTimes;
};


//                    PathDependentAsianG.cpp

PathDependentAsianG::PathDependentAsianG(const MJArray& LookAtTimes_, 
										 double DeliveryTime_,
										 const PayOffBridge& ThePayOff_)
										 :
PathDependent(LookAtTimes_),
DeliveryTime(DeliveryTime_),
ThePayOff(ThePayOff_),
NumberOfTimes(LookAtTimes_.size())
{
}

unsigned long PathDependentAsianG::MaxNumberOfCashFlows() const
{
	return 1UL;
}

MJArray PathDependentAsianG::PossibleCashFlowTimes() const
{
	MJArray tmp(1UL);
	tmp[0] = DeliveryTime;
	return tmp;
}

unsigned long PathDependentAsianG::CashFlows(const MJArray& SpotValues, 
											 std::vector<CashFlow>& GeneratedFlows) const
{
	long double sum =1;                                 // this function is reviesed
	for (int i=0; i<SpotValues.size();i++)              // change the original arithmetic mean to geometric mean
		sum*=SpotValues[i];

	long double p=1.0/NumberOfTimes;
	long double mean = pow(sum,p);

	GeneratedFlows[0].TimeIndex = 0UL;
	GeneratedFlows[0].Amount = ThePayOff(mean);

	return 1UL; 
}

PathDependent* PathDependentAsianG::clone() const
{
	return new PathDependentAsianG(*this);
}

///////////////////////////////////////////////////////////
//
//  BwronianBridgePath.h 
//
////////////////////////////////////////////////////////////

class BrownianBridgePath
{
public:
	BrownianBridgePath(int k, int M,int NumberOfPaths)    //constructor 
	{
		K=k;		//number of strata
		m=M;		// Number of dates
		PathNumber=NumberOfPaths; // Number of path
		for(int i=0;i<=100000;i++)	//initiate the 2 dimension vector W
		{
			vector<double>v;
			W.push_back(v);
		}
		for(int i=0;i<=100000;i++)
			for(int j=0;j<=m;j++)
				W[i].push_back(0);
	}

	void GetVanillaPath()   //stratifid sampling  m=1 scenerio for vanilla
	{
		RandomParkMiller generator1(1);  
		generator1.ResetDimensionality(1);
		MJArray VariateArray1(1);

		int PathPerStrata=PathNumber/K;

		for(int k=0;k<PathPerStrata;k++)
		{
			for(int i=1;i<=K;i++)
			{
				generator1.GetUniforms(VariateArray1);
				double U=VariateArray1[0];
				double V=(i-1+U)/K;
				W[i+K*k][1]=InverseCumulativeNormal(V);
			}
		}
	}


	void GetAsianPath()   //stratifid sampling   m>1 scenerio for asian and other options
	{

		RandomParkMiller generator1(1);  
		generator1.ResetDimensionality(1);
		MJArray VariateArray1(1);

		RandomParkMiller generator2(1);  
		generator2.ResetDimensionality(1);
		MJArray VariateArray2(1);

		int PathPerStrata=PathNumber/K;

		for(int i=1;i<=K;i++)
		{
			for(int e=1;e<=PathPerStrata;e++)
			{
				generator1.GetUniforms(VariateArray1);
				double U=VariateArray1[0];
				double V=(i-1+U)/K;
				W[(i-1)*PathPerStrata+e][m]=InverseCumulativeNormal(V);
			}

			for(int j=1;j<=m-1;j++)
			{

				double tj=j/(double)m;
				double tjminus1=j/(double)m-1.0/m;   
				double a=(1-tj)/(1-tjminus1);
				double b=(tj-tjminus1)/(1-tjminus1);
				double c=sqrt((1-tj)*(tj-tjminus1)/(1-tjminus1));
				for(int e=1;e<=PathPerStrata;e++)
				{
					generator2.GetGaussians(VariateArray2);
					double Z=VariateArray2[0];
					W[(i-1)*PathPerStrata+e][j]=a*W[(i-1)*PathPerStrata+e][j-1]+b*W[(i-1)*PathPerStrata+e][m]+c*Z;
				}
			}
		}
	}

	double getW(int i, int j)
	{return W[i][j];}

private:

	int K;
	int m;
	int PathNumber;
	vector<vector<double>>W;
};



///////////////////////////////////////////////////
//
//                mcsimple.h/.cpp
//
/////////////////////////////////////////////////

#define mcsimple_H

#include <Vanilla3.h>
#include <Parameters.h>
#include <Random2.h>
#include <MCStatistics.h>
#include <cmath>
#include <Arrays.h>

void SimpleMonteCarlo6(const VanillaOption& TheOption, 
					   double Spot,
					   const Parameters& Vol, 
					   double r,
					   double d,
					   unsigned long NumberOfPaths,
					   StatisticsMC& gatherer,
					   RandomBase& generator);

// the basic math functions should be in namespace std but aren't in VCPP6
#if !defined(_MSC_VER)
using namespace std;
#endif

void SimpleMonteCarlo(const VanillaOption& TheOption, 
					  double Spot, 
					  const Parameters& Vol, 
					  double r, 
					  double d,
					  unsigned long NumberOfPaths,
					  StatisticsMC& gatherer,
					  RandomBase& generator)
{
	generator.ResetDimensionality(1);

	double Expiry = TheOption.GetExpiry();
	double variance = Vol.IntegralSquare(0,Expiry);
	double rootVariance = sqrt(variance);
	double itoCorrection = -0.5*variance;       // here is revised to fit the dividend
	double movedSpot = Spot*(exp((r-d)*Expiry) +itoCorrection);

	double thisSpot;
	double discounting = exp(-r*Expiry);

	MJArray VariateArray(1);

	for (unsigned long i=0; i < NumberOfPaths; i++)
	{
		generator.GetGaussians(VariateArray);
		thisSpot = movedSpot*exp( rootVariance*VariateArray[0]);
		double thisPayOff = TheOption.OptionPayOff(thisSpot);
		gatherer.DumpOneResult(thisPayOff*discounting);
	}
	return;
}


///////////////////////////////////////////////////////////
//
//  MCbyStratify_H     use stratify method to do the MC
//
////////////////////////////////////////////////////////////

class MCbyStratify
{

public:

	MCbyStratify(double Expiry,     //constructor
				 double Strike,
				 double Spot,
				 double Vol,
				 double rate,
				 double d, 
				 double PathNumber)
	{
		T=Expiry;   // time to maturity 
		K=Strike;   // strike price
		S=Spot;   // spot price
		vol=Vol; // volatility
		r=rate;   // interest rate
		div=d; // dividend
		NumberOfPaths=PathNumber;
	} 

	void StratifyPricing(   double NumStra,         // pricing fuction
							double NumDates,
							double& VanillaCallPrice,
							double& GeometricAsianCallPrice,
							double& ArithmeticAsianCallPrice)
	{	
		double callsum =0;    
		double GeometricCallSum=0;
		double ArithmeticCallSum=0;
		double stockK;  // The simulated stock price at maturity
		double rootT = sqrt(T);
		// Calculate the fixed part of the price.
		double base = S*exp(r-div - vol*vol/2)* T; 
		BrownianBridgePath Path1(NumStra,NumDates,NumberOfPaths);
		Path1.GetVanillaPath();                                           //get vanilla path

		BrownianBridgePath Path2(NumStra,NumDates,NumberOfPaths);
		Path2.GetAsianPath();											// get asian path

		for (unsigned long i=1; i <= NumberOfPaths; i++)
		{
			double tGaussian = Path1.getW(i,1);

			stockK =S* exp((r-div-vol*vol/2.0)*T+ vol*tGaussian);

			double Gaussian[7];