hw5_chengliangg.cpp

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

			for(int j=1;j<=6;j++)
				Gaussian[j] = Path2.getW(i,j);
			double s[7];s[0]=100;
			double TT;
			double SumGeometric=1;
			double SumArithmetic=0;


			for(int k=1;k<=6;k++)
			{
				TT=T*k/6.0;
				rootT = sqrt(TT);
				s[k] = s[0]* exp((r-div-vol*vol/2.0)*TT+ vol*Gaussian[k]);
				SumGeometric*=s[k];
				SumArithmetic+=s[k];
			}

			double GeometricMean=pow(SumGeometric,1.0/6);					// calculate mean
			double ArithmeticMean=SumArithmetic/6.0;

			double tPayoff = stockK - K;
			double GeometricAsianPayoff=GeometricMean-K;
			double ArithmeticAsianPayoff=ArithmeticMean-K;

			if (tPayoff > 0){
				// We will exercise the call option
				callsum += tPayoff;
			} 
			if (GeometricAsianPayoff> 0){
				GeometricCallSum += GeometricAsianPayoff;
			} 
			if (ArithmeticAsianPayoff > 0){
				ArithmeticCallSum += ArithmeticAsianPayoff;
			} 
		}
		VanillaCallPrice = callsum / NumberOfPaths * exp(-r*T)	;
		GeometricAsianCallPrice=GeometricCallSum/NumberOfPaths * exp(-r*T)	;
		ArithmeticAsianCallPrice=ArithmeticCallSum/NumberOfPaths * exp(-r*T)	;
	}

private:
	double T;   // time to maturity 
	double K;   // strike price
	double S;   // spot price
	double vol; // volatility
	double r;   // interest rate
	double div; // dividend
	unsigned long NumberOfPaths; 
};


///////////////////////////////////////////////////////////
//
//        Sobol.h
//
////////////////////////////////////////////////////////////


//--------------------------------------   Sobol.h  ------------------------------------------------------------
class SobolSequence :public RandomBase
{
public:
	SobolSequence(unsigned long Dimensionality);  // default constructor

	virtual RandomBase* clone() const;			 // virtual duplication function
	virtual void GetUniforms(MJArray& variates);	     // generate Sobol sequence
	virtual void Skip(unsigned long numberOfPaths);
	virtual void SetSeed(unsigned long Seed){};
	virtual void Reset();
	virtual	void ResetDimensionality(unsigned long NewDimensionality);

private:
	VecDoub x;

};

//	---------------------------------------------  	Sobol.cpp	--------------------------------------------- //



SobolSequence::SobolSequence(unsigned long Dimensionality):RandomBase(Dimensionality)
{
	sobseq(-1,x);    // initial Sobol sequence
	x.resize(Dimensionality);
}


RandomBase* SobolSequence::clone() const
{
	return new SobolSequence(*this);
}


void SobolSequence::GetUniforms(MJArray& variates)
{

	sobseq(GetDimensionality(),x);			// get d-dimensional Sobol sequence

	for (unsigned long j=0; j< GetDimensionality(); j++)
	{
		variates[j] = x[j];
	}
}


void SobolSequence::Skip(unsigned long numberOfPaths)
{
	MJArray tmp(GetDimensionality());
	for (unsigned long j=0; j < numberOfPaths; j++)      // skip numberofpaths' Sobol sequence
	{
		GetUniforms(tmp);
	}
}


void SobolSequence::Reset()     
{  
	sobseq(-1,x);            //	reset Sobol sequence
}


void SobolSequence::ResetDimensionality(unsigned long NewDimensionality)
{
	RandomBase::ResetDimensionality(NewDimensionality);		//	reset RandomBase class dimensionality
	sobseq(-1,x);											//  reset Sobol sequence
	x.resize(NewDimensionality);							//  reset vector x dimension
}


///////////////////////////////////////////////////////////
//
// Close Form vanilla
//
////////////////////////////////////////////////////////////

double CloseFormPricing(double Expiry,
						double Strike,
						double Spot,
						double Vol,                             //  volatility
						double r,                              //   interest rate
						double d,                                //  dividend rate
						string type)
{
	double variance = Vol*Vol*Expiry;
	double rootVariance = sqrt(variance);
	double itoCorrection = 0.5*variance;
	double d1 = (1/rootVariance)*(log(Spot/Strike)+(r-d)*Expiry + itoCorrection);        //  the d+  in Black Scholes
	double d2 = (1/rootVariance)*(log(Spot/Strike)+(r-d)*Expiry - itoCorrection);        //   the d- in Black Scholes

	double price;
	if (type=="call")          //  If the type is call option
		price = Spot*exp(-d*Expiry)*CumulativeNormal(d1) - exp(-r*Expiry)*Strike*CumulativeNormal(d2);
	else                        //  If the type is put option (default)
		price = exp(-r*Expiry)*Strike*CumulativeNormal(-d2) - Spot*exp(-d*Expiry)*CumulativeNormal(-d1);

	return price;             //  return the result
}


////////////////////////////////////////////////
//
//     Main Program  
//
////////////////////////////////////////////

int main()
{

	double Expiry = 1;               //   Initialize all the variable
	double Strike = 110; 
	double Spot = 100; 
	double Vol = 0.3; 
	double r = 0.05; 
	double d = 0.02;
	double period=6;
	unsigned long NumberOfPaths = 100000;
	unsigned NumberOfDates=6;
	double CallPrice;
	double GeoAsianCall;
	double ArithAsianCall;
	double NumberOfStrata=100;


	////////////////////////////////////////////////
	//
	//        Vanilla
	//
	///////////////////////////////////////////////

	//-------------------------closed-form vanilla--------------------------



	PayOffCall thePayOff(Strike);

	VanillaOption theVanilla(thePayOff, Expiry);

	ParametersConstant VolParam(Vol);
	ParametersConstant rParam(r);
	ParametersConstant dParam(d);

	StatisticsMean gatherer;                   // initialize the statistic object
	ConvergenceTable gathererTwo(gatherer);

	RandomParkMiller generator1(1);      // Park Miller generator object

	//  AntiThetic GenTwo(generator);


	double CF = CloseFormPricing(Expiry, Strike, Spot,Vol, r, d, "call");           // call close form function
	cout<<"Closed-form vanilla call price ="<<CF<<endl;



	//-------------------------parkmiller vanilla--------------------------


	RandomParkMiller generator3(1);             //   Since Joshi's original code use the inverse distribution normal generator, therefore I use the original generator of Park Miller by Joshi
	AntiThetic GenTwo(generator3);    // Enable the antithetic function

	StatisticsMean gatherer3;                  // initialize the statistic object
	ConvergenceTable gathererTwo3(gatherer3);

	SimpleMonteCarlo(theVanilla,
		Spot, 
		VolParam,
		r,
		d,
		NumberOfPaths,
		gathererTwo3,
		GenTwo);

	vector<vector<double> >results =gathererTwo3.GetResultsSoFar();
	double size = results.size() - 1; 
	cout<<"MC vanilla call, Park-Miller uniforms, antithetics = "<<results[size][0]<<endl;


	//-------------------------Stratified vanilla--------------------------


	MCbyStratify stratify(Expiry,Strike,Spot,Vol,r,d,NumberOfPaths);
	stratify.StratifyPricing(NumberOfStrata,NumberOfDates,CallPrice,GeoAsianCall,ArithAsianCall);
	cout<<"MC vanilla call Park-Miller uniforms, stratified =";
	cout<<CallPrice<<endl;




	//-------------------------sobol vanilla--------------------------


	SobolSequence generator6(1);          //  Sobol Method to generate Unif[0,1] RV 

	StatisticsMean gatherer6;                  // initialize the statistic object
	ConvergenceTable gathererTwo6(gatherer6);

	SimpleMonteCarlo(theVanilla,
		Spot, 
		VolParam,
		r,
		d,
		NumberOfPaths,
		gathererTwo6,
		generator6);

	results =gathererTwo6.GetResultsSoFar();          //  get the result

	cout<<"QMC call option price with Sobol sequence = "<<results[size][0]<<endl<<endl;	   // output the result


	////////////////////////////////////////////////
	//
	//       Geometric Asian
	//
	///////////////////////////////////////////////

	//-------------------------closed form geometric asia--------------------------	

	GeometricAsianCF closedform(Spot,Strike,r,d,Vol,Expiry,period);
	cout<<"Closed-form geometric Asian call price =";
	cout<<closedform.option_price_asian_geometric_average_price_call()<<endl;	


	//-------------------------park miller geometric asia--------------------------	


	MJArray times(NumberOfDates);

	for (unsigned long i=0; i < NumberOfDates; i++)
		times[i] = (i+1.0)*Expiry/NumberOfDates;


	PathDependentAsianG theOption2(times, Expiry, thePayOff);

	StatisticsMean gathererB;
	ConvergenceTable gathererTwoB(gathererB);


	RandomParkMiller generatorB(NumberOfDates);

	AntiThetic GenTwoB(generatorB);

	ExoticBSEngine theEngineB(theOption2, rParam, dParam, VolParam, GenTwoB, Spot);

	theEngineB.DoSimulation(gathererTwoB, NumberOfPaths);


	cout <<"MC Geometric Asian call, Park-Miller uniforms, antithetics =";
	results =gathererTwoB.GetResultsSoFar();

	cout << results[results.size()-1][0] ;

	cout << "\n";
	//-------------------------park miller geometric asia--------------------------	

	cout<<"MC geometric Asian call, Park-Miller uniforms, stratified =";
	cout<<GeoAsianCall<<endl;

	//-------------------------sobol geometric asia--------------------------	


	SobolSequence generatorC(6);          //  Sobol Method to generate Unif[0,1] RV 

	StatisticsMean gathererC;                  // initialize the statistic object
	ConvergenceTable gathererTwoC(gathererC);

	ExoticBSEngine theEngineC(theOption2, rParam, dParam, VolParam, generatorC, Spot);
	theEngineC.DoSimulation(gathererTwoC, NumberOfPaths);

	results =gathererTwoC.GetResultsSoFar();          //  get the result

	cout<<"QMC geometric Asian call, Sobol sequence = "<<results[size][0]<<endl<<endl;	   // output the result



	////////////////////////////////////////////////
	//
	//       Arithmetic Asian
	//
	///////////////////////////////////////////////


	//-------------------------parkmiller arithmetic asia--------------------------	

	PathDependentAsian theOption(times, Expiry, thePayOff);

	StatisticsMean gathererA;
	ConvergenceTable gathererTwoA(gathererA);


	RandomParkMiller generatorA(NumberOfDates);

	AntiThetic GenTwoA(generatorA);

	ExoticBSEngine theEngine(theOption, rParam, dParam, VolParam, GenTwoA, Spot);

	theEngine.DoSimulation(gathererTwoA, NumberOfPaths);


	cout <<"MC arithmetic Asian call, Park-Miller uniforms, antithetics= ";
	results =gathererTwoA.GetResultsSoFar();

	cout << results[results.size()-1][0] ;

	cout << "\n";


	//-------------------------stratified arithmetic asia--------------------------	

	cout<<"MC arithmetic Asian call, Park-Miller uniforms, stratified =";
	cout<<ArithAsianCall<<endl;


	//-------------------------sobol arithmetic asia--------------------------	

	SobolSequence generatorD(6);          //  Sobol Method to generate Unif[0,1] RV 

	StatisticsMean gathererD;                  // initialize the statistic object
	ConvergenceTable gathererTwoD(gathererD);

	ExoticBSEngine theEngineD(theOption, rParam, dParam, VolParam, generatorD, Spot);
	theEngineD.DoSimulation(gathererTwoD, NumberOfPaths);

	results =gathererTwoD.GetResultsSoFar();          //  get the result

	cout<<"QMC arithmetic Asian call, Sobol sequence = "<<results[size][0]<<endl;	   // output the result

	return 0;
}