Virtual methods in MQL5 – an application on options
[Versiunea romaneasca] [MQLmagazine.com in romana] [English edition]
This is a continuation of our article Intro to Options. Part II : Option valuation and greeks . However, it can be viewed and understood seperately, because it focuses on the programming side.
The option models are exceedingly complicated from the math point of view. The purpose of this article is to present an implementation of the Black-Scholes valuation model, also by introducing virtual methods.
One of the functions that appears over and over again in financial math, also in the Black-Scholes model, is the cummulative distribution function.
The function that describes the bell-shaped form of the distribution is called probability density function. The area that is below it, meaningly the definite integral , is the cummulative distribution function. And for the case with 0 mean and 1 standard deviation, (assuming also no skewness , mezokurtic) this distribution is called normal.
The definite integral, also known as Riemann sum, calculates the area below a function’s graphics by adding infinitezimal rectangles that span below the function and x axis. Obviously, the Riemann calculus procedure is always the same, while the function is something user defined. Both Pascal and C++ allowed function types. A C++ declaration would have looked like this:
typedef double (*SingleVarFunction) (double) ;
double Riemann(SingleVarFunction pFunc,double a, double b, int divisions)
In this C++ case, the Riemann could have received as parameter any function that received one double as parameter and returned double as result. However, this mechanism does not exist in MQL5.
This is where virtual methods come into play. Objects have properties. Properties that are designed to be public can be accessed from outside the objects, and then methods called again. For instance the user can change the radius of a circle and call the redraw method. But in the case of a Riemann sum what has to be passed as a property value is itself a method. User-defined methods that are called by inherited methods are called virtual methods.
This is the MQL5 implementation for Riemann sum:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 | class RiemannSumCalculus { public: virtual double IntegratableFunction(double x) { return (0.0); } double Riemann(double a,double b,uint divisions) { double sum,swap,stepsize,x; long istep; if (a>b) { swap=a; a=b; b=swap; } sum=0.0; if (divisions==0) divisions=1; stepsize=(b-a)/divisions; for (istep=1;istep<=divisions;istep++) { x=a+(istep-0.5)*stepsize; sum=sum+IntegratableFunction(x)*stepsize; } return(sum); } }; |
As you see, the IntegratableFunction method is designed as virtual. The function itself does nothing, returns always a zero. Because the function is not made to be used directly. It is just a place holder for user-defined IntegratableFunction(s) that are to be used by the Riemann method.
Next step is to have the Gauss cummulative distribution.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | class GaussCummulativeFunction : public RiemannSumCalculus { private: double sigma; public: virtual double IntegratableFunction(double x) { return( exp((-x*x)/(2*sigma*sigma)) ); } double CummulativeDistribution(double stddev,double l1,double l2,uint divisions) { sigma=stddev; return( (1/(stddev*sqrt(2*M_PI))) * Riemann(l1,l2,divisions) ); } }; |
As you can see the class is inherited from RiemannSumCalculus. The IntegratableFunction is redefined, and the CummulativeDistribution function calls the inherited Riemann, which, in turn, will call the local IntegratableFunction, instead of the inherited null-value place holder function. And since this is a class, we can create a more easy way to access this, by calling this method from a function:
1 2 3 4 5 | double cum_norm(double x) { GaussCummulativeFunction normaldistribution; return(normaldistribution.CummulativeDistribution(1,-7.0,x,1000000)); } |
Perhaps 1,000,000 divisions seem a bit too much, but it’s necesary to keep results accurate.
What follows is the MQL5 version of the Black-Scholes procedure written by George Levy in his book, “Computational finance using C and C++”, page 76 (with some small modifications):
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 | struct Greeks { double delta; double gamma; double vega; double theta; double rho; }; void black_scholes(double& value, Greeks &greeks, double s0, double x, double sigma, double t, double r, double q, bool iscall, long &iflag) { /* Input parameters: ================= s0 - the current price of the underlying asset x - the strike price sigma - the volatility t - the time to maturity r - the interest rate q - the continuous dividend yield iscall - true if option is a call Output parameters: ================== value - the value of the option greeks[] - the hedge statistics output; iflag - an error indicator */ greeks.delta=0.0; greeks.gamma=0.0; greeks.vega=0.0; greeks.theta=0.0; greeks.rho=0.0; double one=1.0,two=2.0,zero=0.0; double eps,d1,d2,temp,temp1,temp2,np; /* Check if any of the the input arguments are too small */ eps=1.0e-16; if ( (x < eps) || (sigma < eps) || (t < eps) ) { iflag = 2; return; } temp = log(s0/x); d1 = temp+(r-q+(sigma*sigma/two))*t; d1 = d1/(sigma*sqrt(t)); d2 = d1-sigma*sqrt(t); /* evaluate the option price */ if (iscall==true) value = (s0*exp(-q*t)*cum_norm(d1)- x*exp(-r*t)*cum_norm(d2)); else value = (-s0*exp(-q*t)*cum_norm(-d1) + x*exp(-r*t)*cum_norm(-d2)); /* then calculate the Greeks */ temp1 = -d1*d1/two; d2 = d1-sigma*sqrt(t); np = (one/sqrt(two*M_PI)) * exp(temp1); if (iscall==true) { /* a call option */ greeks.delta = (cum_norm(d1))*exp(-q*t); greeks.theta = -s0*exp(-q*t)*np*sigma/(two*sqrt(t)) + q*s0*cum_norm(d1)*exp(-q*t)- r*x*exp(-r*t)*cum_norm(d2); greeks.rho = x*t*exp(-r*t)*cum_norm(d2); } else { /* a put option */ greeks.delta = (cum_norm(d1) - one)*exp(-q*t); greeks.theta = -s0*exp(-q*t)*np*sigma/(two*sqrt(t)) - q*s0*cum_norm(-d1)*exp(-q*t) + r*x*exp(-r*t)*cum_norm(-d2); greeks.rho = -x*t*exp(-r*t)*cum_norm(-d2); } greeks.gamma = np*exp(-q*t)/(s0*sigma*sqrt(t)); greeks.vega = s0*sqrt(t)*np*exp(-q*t); return; } |
And the following , translated from the same book, to get the implied volatility:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 | void implied_volatility(double value, double s0, double x, double& sigma[],double t, double r, double q, bool iscall, long& iflag) { /* Input parameters: ================= value - the current value of the option s0 - the current price of the underlying asset x - the strike price sigma[] - the input bounds on the volatility: sigma[0], the lower bound and, sigma[1],the upper bound t - the time to maturity r - the interest rate q - the continuous dividend yield iscall - true if option is a call; Output parameters: ================== sigma[] - the element sigma[0] contains the estimated implied volatility iflag - an error indicator */ double zero=0.0; double fx, sig1, sig2; double val,tolx; double temp,eps,epsqrt,temp1,v1; long max_iters, i, ind, ir; double c[20]; double sig,vega; Greeks greeks; long done; tolx = eps; eps=1.0e-16; epsqrt = sqrt(eps); if (iscall == true) /* a call option */ temp1 = MathMax(s0*exp(-q*t)-x*exp(-r*t),zero); else /* a put option */ temp1 = MathMax(x*exp(-r*t)-s0*exp(-q*t),zero); v1 = fabs(value-temp1); if (v1 <= epsqrt) { /* the volatility is too small */ iflag = 3; return; } iflag = 0; i = 0; max_iters = 50; done = 0; sig = sigma[0]; /* initial estimate */ val = value; while ((i < max_iters) && (!done)) { /* Newton iteration */ black_scholes(val,greeks,s0,x,sig,t,r,q,iscall,iflag); /* compute the Black-Scholes option value, val */ vega = greeks.vega; /* and vega. */ sig1 = sig - ((val - value)/vega); /* compute the new estimate of sigma using Newton’s method */ if (tolx > fabs((sig1 - sig)/sig1)) { /* check whether the specified accuracy has been reached */ done = 1; } sig = sig1; /* up date sigma */ ++i; } sigma[0] = sig1; /* return the estimate for sigma */ return; } |
Now , keeping these parts as the beginning of our script, we will complete our script to test the Black-Scholes model.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 | //+------------------------------------------------------------------+ //| Script program start function | //+------------------------------------------------------------------+ void OnStart() { double StockPrice=100.00; double SelectedStrike=100.00; double r=0.1; double sigma=0.3; double q=0.06; ObjectsDeleteAll(ChartID(),-1,-1); int s0; Greeks greeks; double value; double t; int T; long iflag; for (T=10;T>=1;T--) { t=0.1*T; black_scholes(value, greeks, StockPrice, SelectedStrike, sigma, t, r, q, true, iflag); Print(DoubleToString(t,1)," :: ",DoubleToString(value,3)," ",DoubleToString(greeks.delta,3)," ",DoubleToString(greeks.gamma,3)," ",DoubleToString(greeks.vega,3)," ",DoubleToString(greeks.theta,3)," ",DoubleToString(greeks.rho,3)); } Print("Time Value Delta Gamma Vega Theta Rho"); Print("Calls"); Print(""); for (T=10;T>=1;T--) { t=0.1*T; black_scholes(value, greeks, StockPrice, SelectedStrike, sigma, t, r, q, false, iflag); Print(DoubleToString(t,1)," :: ",DoubleToString(value,3)," ",DoubleToString(greeks.delta,3)," ",DoubleToString(greeks.gamma,3)," ",DoubleToString(greeks.vega,3)," ",DoubleToString(greeks.theta,3)," ",DoubleToString(greeks.rho,3)); } Print("Time Value Delta Gamma Vega Theta Rho"); Print("Puts"); Print(""); Print("Continuous dividend yield: ",q*100,"%"); Print("Implied volatility: ",sigma*100,"%"); Print("Riskless rate: ",r*100,"%"); Print("Selected strike: ",SelectedStrike); Print("Stock price: ",StockPrice); Print("ASSUMPTIONS: "); } //+------------------------------------------------------------------+ |
And the output should look like this:
Now delete the previous added lines that completed the script, and add the following ones, to test implied volatility function:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 | //+------------------------------------------------------------------+ //| Script program start function | //+------------------------------------------------------------------+ void OnStart() { double S = 10.0; double X = 10.5; double r = 0.1; double sigmat = 0.5; double q = 0.04; Greeks greeks; double sigma[2]; double value, T; long i, ifail, put; double origsigma=sigmat; long flag; ifail = 0; for (i = 5;i >= 1; i--) { T = (double)i*0.5; black_scholes(value,greeks,S,X,sigmat,T,r,q,true,flag); sigma[0] = 0.05; sigma[1] = 1.0; implied_volatility(value,S,X,sigma,T,r,q,true,flag); printf("%8.4f %15.4f %15.4f (%8.4e) ",T,value,sigma[0],fabs(sigmat-sigma[0])); sigmat = sigmat - 0.1; } printf (" Time(years) Option value True sigma Error"); Print(""); Print("Underlying= ",DoubleToString(S,2)); Print("Strike= ",DoubleToString(X,2)); Print("Continuous dividend yield= ",DoubleToString(q*100,2),"%"); Print("Implied volatility= ",DoubleToString(origsigma*100,2),"%"); Print("Riskless rate= ",DoubleToString(r*100,2),"%"); Print("x= ",DoubleToString(X,2)); Print("S= ",DoubleToString(S,2)); Print("ASSUMPTIONS"); return; } |
The script plugs some values into the Black-Scholes model, gets option values, then sets a guess range and calls the implied volatility function, that gives pretty accurate results:
