# (C++) Monte Carlo pricer for SABR model to test Hagan / Paulot formulas

I'm trying to test the so-called Hagan formula (p.6 of this paper) and the Paulot formula, order 1 only (eq. (43) p.19 of this paper. For this, i'm trying to use both Euler and Milstein scheme described here (p.9, eq. (3.1) and (3.2)) to price a call option, but the results seem not very consistent, so i'm asking myself if my code's right...

That is my C++ function:

double MC_SABR_price(const int& num_sims, const int& num_intervals, const double& F_0, const double& K, const double& alpha, const double& beta, const double& rho, const double& nu, const double& r, const double& T)
{
double dt = T / num_intervals;
double F[num_intervals];
double V[num_intervals];
F[0] = F_0;
V[0] = alpha;

double payoff_sum = 0.0;
for (int i=1; i<num_sims; i++)
{
for (int j=1; j<num_intervals; j++)
{
double Z1 = NormalSimulation();
double Z2 = NormalSimulation();
V[j] = V[j-1] * exp((nu * sqrt(dt) * Z1) - (0.5 * nu * nu * dt));
F[j] = F[j-1] + (V[j-1] * pow(F[j-1], beta) * sqrt(dt) * ((rho * Z1) + (sqrt(1 - (rho * rho)) * Z2)));
F[j] = max(F[j], 0.0);
}
payoff_sum += max(F[num_intervals-1] - K, 0.0);
}
return (payoff_sum / num_sims) * exp(-r*T);
}


With those parameters:

double num_sims = 100000;   // Number of simulated asset paths
double num_intervals = 1000;  // Number of intervals for the asset path to be sampled

double F_0 = 5.0;       // Initial forward price
vector<double> K(10);
for (int i=0; i<K.size(); i++) { K[i] = 1.0 + i; }

double r = 0.0;         // Risk-free rate
double T = 2.5;         // One year until expiry

double alpha = 0.3;     // Initial volatility
double beta = 0.7;      // Elasticity
double rho = -0.5;      // Correlation of asset and volatility
double nu = 0.4;        // "Vol of vol"


These are the results I get:

As you can see nothing coincide...

• A few comments about the above: 1) I would just simulate the SDE for V also using Euler, rather than solving the SDE exactly as you've done (mostly for consistency). 2) You enforce F = max(F, 0), but it seems possible (probably that) F is strictly positive. Perhaps if F becomes negative just resample Z2, as this approximate solution may have very difference positivity properties to your SDE solution. 3) For your price estimates you use Monte Carlo but don't compute the standard deviation. Compute this standard deviation and this will give you your error estimates. – oliversm Dec 12 '18 at 10:29
• Just final note on point (3) that I mentioned, if that for only 100'000 samples a relative error of ~0.01 seems optimistic. You can calculate what sort of relative error you can expect to have for such a small sample size. My expectation is that these relative errors you measure and what you might expect will explain the "not very consistent" results you are obtaining. – oliversm Dec 12 '18 at 10:32