# Integration, the Simpson´s rule and QuantLib

QuantLib provides with numerous mathematical and function-related tools, and it closely works with the Boost library.

Here, we deal with integrals. More specifically, we implement the well-known Simpson’s integral. We keep the things concrete by calculating the value of a call option from its integral representation.

The integral of a function is the area between the curve and the x-axis. Limits can be defined so that we have a definite integral. The integral will thus be computed by dividing the x-axis and the corresponding area between the curve and the x-axis between the upper and the lower bound.

We can write :

where a is the lower limit, and b the upper limit.

The integral representation of a call is

where f(x) is the lognormal density with mean

And standard deviation

The Simpson’s rule gets numerical approximations of definite integrals

Our code uses the boost::bind class, which is a predefined metafunction. Please, refer to the Boost documentation.

The a-lower limit is the strike, and the b-upper limit is an arbitrary higher value. In this case 10 x a.

```#include<ql\quantlib.hpp>
#include <boost/math/distributions.hpp>
#include<math.h>

using namespace QuantLib;

Real callOptionFunction(
Real underlying,
Real strike,
Rate riskFreeInterestRate,
Volatility volatility,
Time tau,
Real x){
Real mean = log(underlying) + (riskFreeInterestRate
- 0.5 * volatility * volatility) * tau ;

Real standardDeviation = volatility * sqrt(tau);

boost::math::lognormal_distribution<>d(mean, standardDeviation);
return (x - strike) * pdf(d,x) * exp (-riskFreeInterestRate * tau);
}

int main(int, char*[]){

//option parameters
Real underlying = 30;
Real strike = 36;
Rate riskFreeInterestRate = 0.06;
Time tau = 0.5; // 0,5 year
Volatility volatility = 0.20;

//Integration parameters : absolute accuracy
//and maximum number of evaluations
Real absoluteAccuracy = 1e-4;
Size maxEvaluations = 1e3;

boost::function<Real(Real)> ptrF;
ptrF = boost::bind(
&callOptionFunction,
underlying,
strike,
riskFreeInterestRate,
volatility,
tau,
_1);

SimpsonIntegral numInt(
absoluteAccuracy,
maxEvaluations);

std::cout << "Call option value : " << numInt(
ptrF,
strike,/*lower limit a*/
10*strike)/*upper limit b*/ << std::endl;

return 0;
}```
Anuncios