Delta Gamma hedging

#include <iostream>
#include <string>
#include "normdist.h"
#include <cmath>
#include <fstream>
#include<ctime>
#include<cstdlib>

// We hedge an European call option with a strike price K of 100 on stock S with initial price S0 = 100
// and drift =0.15. The implied volatility sigma1 = 0.2 is at the purchase.
// For the purpose of designing a delta hedge, we 1 year option on the stock with a strike price of 90
// and implied volatility sigma2 = 0.3. The interest rate is 0.05.


using namespace std;
class CallOption {

public:

    double S, X, sigma, r, time;
    double CallOptionPrice(double S, double X,double sigma,double r,double time);
    double GetDelta(double S, double X,double sigma,double r,double time);
};

double CallOption::CallOptionPrice(double S,double X,double sigma,double r,double time) {
  double d1=(log(S/X)+r*time)/(sigma*sqrt(time))+0.5*sigma*sqrt(time);
  double d2=d1-(sigma*sqrt(time));
  double c=S* N(d1) - X*exp(-r*time)*N(d2);

  return c;
}
double CallOption::GetDelta(double S,double X,double sigma,double r,double time){
double d1=(log(S/X)+r*time)/(sigma*sqrt(time))+0.5*sigma*sqrt(time);
return N(d1);
}

double GetOneGaussianByBoxMuller ( ) {
double result , x , y , SizeSquared ;
do {
x = 2.0*rand()/static_cast<double>(RAND_MAX)-1;
y = 2.0*rand()/static_cast<double>(RAND_MAX)-1;
SizeSquared = x*x+y*y;
}
while ( SizeSquared >= 1.0);
result = x*sqrt(-2*log(SizeSquared)/SizeSquared);
return result;
}

int main(){
std::cout.setf(std::ios_base::fixed);
std::cout.setf(std::ios_base::showpoint);
std::cout.precision(4);
std::ofstream fileOutput("Delta_63.csv", std::ios::app);
srand(time(NULL));
    CallOption myOption;
for (int j=1;j<1000;j++)
{
    double S=100;
    double X=100;
    double sigma=0.2;
    double r=0.05;
    double time=1.0;
    double trend=0.15;
    double pnl=0;
    double loan=myOption.CallOptionPrice(S,X,sigma,r,time);
    double deposit=myOption.GetDelta(S,X,sigma,r,time)*S;
    double drift=(trend/63.0-0.5*pow(sigma/sqrt(63),2));

for( double i = 1; i < 63; i++ )
   {
    double Delta0=myOption.GetDelta(S,X,sigma,r,time);
    loan=loan*exp(r/63.0);
    time=time-(1.0/63.0);
    S=S*exp((drift+GetOneGaussianByBoxMuller()*(sigma/sqrt(63))));
    double Delta1=myOption.GetDelta(S,X,sigma,r,time);
    double sell=Delta0*S+(Delta1-Delta0)*S;
    deposit=deposit*exp(r/63.0);
    deposit=deposit+(Delta1-Delta0)*S;
    pnl=myOption.CallOptionPrice(S,X,sigma,r,time)-loan-sell+deposit;
   }
cout<<pnl<<endl;
fileOutput<<pnl<<endl;
}
fileOutput.close();
return 0;
}