Movimiento Browniano Geométrico en C++

Estoy tratando de desarrollar un pequeño calculador de opciones utilizando la API de win32.
Para hacer eso, utilizo la simulación de Monte Carlo para calcular el precio de una opción de compra, pero hay un error en mi simulación y no veo dónde está.
¿Alguien puede decirme por qué los dos precios son diferentes?

El precio de una opción de compra con la fórmula de Black y Scholes es de 6.84, pero el que da la simulación de Monte Carlo es de 7.54.
(No hay error en el precio calculado con la fórmula de Black y Scholes)

Este es el código:

std::vector<double> vecteur_pas(double T) {
    std::vector<double> pas;
    pas.push_back(0);
    double x = T / nb_pas;
    for (int i = 1; i <= nb_pas; i++) {
        pas.push_back(pas[i-1] + x);
    }
    return pas;
}

std::vector<double> NormalRnd() {
    std::vector<double> brow;
    brow.push_back(0);
    double unif_1, unif_2;
    for (int i = 0; i < nb_pas; i++) {
        unif_1 = (double(rand() + 0.5)) / (double(RAND_MAX) + 1.0);
        unif_2 = (double(rand() + 0.5)) / (double(RAND_MAX) + 1.0);
        brow.push_back(sqrt(-2 * log(unif_1)) * cos(2 * M_PI * unif_2));
    }
    return brow;
}

std::vector<double> MBG(double S, double mu, double vol, double T) {
    std::vector<double> traj;
    traj.push_back(S);
    std::vector<double> b =NormalRnd();
    std::vector<double> pas = vecteur_pas(T);
    double drift = 0.0;
    double diffusion = 0.0;
    for (int i = 1; i <= nb_pas; i++) {
        drift = (mu - 0.5 * pow(vol, 2)) * (pas[i]-pas[i-1]);
        diffusion = vol * b[i] * sqrt(pas[i] - pas[i - 1]);
        traj.push_back(traj[i - 1] * exp(drift + diffusion));
    }
    return traj;
}

La función MBG se llama en un bucle después de esto:

for (int i = 0; i < 100000; i++)
{
    if ((i % 1000) == 0)
    {
        SendDlgItemMessage(Calcul, IDE_PB, PBM_STEPIT, 0, 0);
    }
    vector<double> proc_prix = MBG(actif.S, actif.r, actif.v, actif.T);
    double last_value = proc_prix[proc_prix.size() - 1];
    Prime = Prime + std::max(last_value - actif.K, 0.0);
}
Prime = Prime / 100000;

Este es el resultado.