Is it possible to use autodiff on Monte Carlo Simulation?

I am running a Monte Carlo simulation for a specific function and would like to compute the gradient of this function with respect to several independent variables. I initially assumed that using the autodiff library would be suitable for this purpose. However, when I attempt to differentiate with respect to any particular variable, the gradient consistently evaluates to zero, which is incorrect. I am uncertain about the cause of this issue. Could you provide any insights or suggestions on why this might be occurring?

It would be helpful if you provided a short code sample that can be directly executed and demonstrates your exact problem.

1 Like

Hi @hax10 , please find below my code:

// Include the necessary library for random number generation and autodiff
use rand::SeedableRng;
use rand::rngs::StdRng;
use rand_distr::{StandardNormal, Distribution};
use autodiff::*;  // Import for automatic differentiation

/// Function to simulate the stock price at maturity using Geometric Brownian Motion
fn simulate_gbm(
    s0: FT<f64>,          // Initial stock price (as differentiable type)
    mu: FT<f64>,          // Risk-free interest rate
    sigma: FT<f64>,       // Volatility of the stock
    t: FT<f64>,           // Time to maturity in years
    timestep: f64,        // Timestep for the simulation
    num_simulations: usize, // Number of simulations to run
    seed: u64,            // Seed for random number generator for reproducibility
) -> Vec<FT<f64>> {
    let mut rng = StdRng::seed_from_u64(seed);
    let num_steps = timestep as usize; // Calculate the number of steps
    let dt = t / num_steps as f64;     // Time step size

    let mut final_prices = Vec::with_capacity(num_simulations);

    for _ in 0..num_simulations {
        let mut st = s0; // Start with the initial stock price (as FT<f64>)

        for _ in 0..num_steps {
            let z: f64 = StandardNormal.sample(&mut rng);
            let half_sigma_squared = FT::<f64>::cst(0.5) * sigma * sigma;
            let drift = (mu - half_sigma_squared) * dt;
            let diffusion = sigma * FT::<f64>::sqrt(dt) * FT::<f64>::cst(z);
            st = st * (drift + diffusion).exp(); // Apply drift and diffusion
        }

        final_prices.push(st); // Store the final simulated price for each path
    }

    final_prices
}

/// Function to calculate the price of a European call option using Monte Carlo simulation
fn price_european_call(
    s0: FT<f64>,          // Initial stock price
    k: f64,               // Strike price
    mu: FT<f64>,          // Drift
    r: FT<f64>,           // Risk-free interest rate
    sigma: FT<f64>,       // Volatility of the stock
    t: FT<f64>,           // Time to maturity in years
    timestep: f64,        // Timestep for the simulation
    num_simulations: usize, // Number of simulations to run
    seed: u64,            // Seed for random number generator for reproducibility
) -> FT<f64> {
    let final_prices = simulate_gbm(s0, mu, sigma, t, timestep, num_simulations, seed);
    let mut payoff_sum = FT::<f64>::cst(0.0);

    for final_price in final_prices {
        let payoff = (final_price - k).max(FT::<f64>::cst(0.0)); // Payoff for a call option
        payoff_sum += payoff;
    }

    let average_payoff = payoff_sum / FT::<f64>::cst(num_simulations as f64);
    (-r * t).exp() * average_payoff // Discount the average payoff back to present value
}


/// Function to calculate the gamma of the option through automatic differentiation
fn calculate_gamma(x: &[FT<f64>]) -> FT<f64> {
    let base_s0 = x[0];
    let k = x[1].to_f64().unwrap();
    let mu = x[2];
    let r = x[3];
    let sigma = x[4];
    let t = x[5];
    let timestep = x[6].to_f64().unwrap();
    let num_simulations = x[7].to_f64().unwrap() as usize;
    let seed = x[8].to_f64().unwrap() as u64;

    // Define the price function as a closure
    let price_fn = |vars: &[FT<f64>]| -> FT<f64> {
        let s0 = vars[0];
        let mu = vars[2];
        let r = vars[3];
        let sigma = vars[4];
        let t = vars[5];

        price_european_call(s0, k, mu, r, sigma, t, timestep, num_simulations, seed)
    };

    // Convert x to f64 for differentiation
    let x_f64: Vec<f64> = x.iter().map(|v| v.to_f64().unwrap()).collect();

    // Compute delta as the gradient of the price function
    let delta_fn = |vars: &[FT<f64>]| -> FT<f64> {
        let grad_result = grad(price_fn, &x_f64);
        grad_result[0].into() // Return derivative with respect to s0
    };

    // Calculate gamma as the gradient of the delta function
    let gamma_result = grad(delta_fn, &x_f64);
    gamma_result[0].into() // Return the second derivative (gamma)
}

/// Main function to demonstrate Greek calculation
fn main() {
    // Parameters for the simulation
    let s0 = 100.0; // Initial stock price
    let k = 100.0;  // Strike price
    let mu = 0.01;  // Drift (expected return)
    let r = 0.01;   // Risk-free interest rate (1%)
    let sigma = 0.2; // Volatility (20%)
    let t = 1.0;    // Time to maturity (1 year)
    let timestep = 365.0; // Number of time steps per year
    let num_simulations = 1000; // Number of Monte Carlo simulations
    let seed = 42;  // Seed for random number generator

  

    // Calculate Gamma
    let gamma = calculate_gamma(&v);
    println!("Gamma: {}", gamma.to_f64().unwrap());
}

I get the result

 Finished `dev` profile [unoptimized + debuginfo] target(s) in 0.01s
     Running `target/debug/bs-option-pricing`
Gamma: 0

Do you know why this might be happening? Any suggestions on avoiding the error, as I am pretty sure it has gamma derivative!

Hi, I'm not able to run your code as-is, because you have not provided the value for v in

    let gamma = calculate_gamma(&v);

I invented a dummy value:

    let gamma = calculate_gamma(&vec![
        autodiff::F::new(1.0, 1.0),
        autodiff::F::new(2.0, 2.0),
        autodiff::F::new(3.0, 3.0),
        autodiff::F::new(4.0, 4.0),
        autodiff::F::new(5.0, 5.0),
        autodiff::F::new(6.0, 6.0),
        autodiff::F::new(7.0, 7.0),
        autodiff::F::new(8.0, 8.0),
        autodiff::F::new(9.0, 9.0),
    ]);

which does reproduce the problem of a zero gradient.

I don't know if this is a problem with the autodiff library though. If you look at the return value from the price_european_call function, it seems to always be zero as well, and this is before you call grad().

I suggest you check your own simulation code to make sure it's doing what you think it's doing. You can add some unit tests to your code to help you.

Hi @hax10 ,

Yeah sorry, didn't check the v value, I edited the code, please note the option value is non-zero, not sure what you mean, I am also able to find first order derivatives, it is the second order derivative that gives me troubles and I am wondering, if for discretized calculation, for instance Monte Carlo Simulation, does autodiff fail to find second order derivative? If so, why is it? I am curious.

// Include the necessary library for random number generation and autodiff
use rand::SeedableRng;
use rand::rngs::StdRng;
use rand_distr::{StandardNormal, Distribution};
use autodiff::*;  // Import for automatic differentiation

/// Function to simulate the stock price at maturity using Geometric Brownian Motion
fn simulate_gbm(
    s0: FT<f64>,          // Initial stock price (as differentiable type)
    mu: FT<f64>,          // Risk-free interest rate
    sigma: FT<f64>,       // Volatility of the stock
    t: FT<f64>,           // Time to maturity in years
    timestep: f64,        // Timestep for the simulation
    num_simulations: usize, // Number of simulations to run
    seed: u64,            // Seed for random number generator for reproducibility
) -> Vec<FT<f64>> {
    let mut rng = StdRng::seed_from_u64(seed);
    let num_steps = timestep as usize; // Calculate the number of steps
    let dt = t / num_steps as f64;     // Time step size

    let mut final_prices = Vec::with_capacity(num_simulations);

    for _ in 0..num_simulations {
        let mut st = s0; // Start with the initial stock price (as FT<f64>)

        for _ in 0..num_steps {
            let z: f64 = StandardNormal.sample(&mut rng);
            let half_sigma_squared = FT::<f64>::cst(0.5) * sigma * sigma;
            let drift = (mu - half_sigma_squared) * dt;
            let diffusion = sigma * FT::<f64>::sqrt(dt) * FT::<f64>::cst(z);
            st = st * (drift + diffusion).exp(); // Apply drift and diffusion
        }

        final_prices.push(st); // Store the final simulated price for each path
    }

    final_prices
}

/// Function to calculate the price of a European call option using Monte Carlo simulation
fn price_european_call(
    s0: FT<f64>,          // Initial stock price
    k: f64,               // Strike price
    mu: FT<f64>,          // Drift
    r: FT<f64>,           // Risk-free interest rate
    sigma: FT<f64>,       // Volatility of the stock
    t: FT<f64>,           // Time to maturity in years
    timestep: f64,        // Timestep for the simulation
    num_simulations: usize, // Number of simulations to run
    seed: u64,            // Seed for random number generator for reproducibility
) -> FT<f64> {
    let final_prices = simulate_gbm(s0, mu, sigma, t, timestep, num_simulations, seed);
    let mut payoff_sum = FT::<f64>::cst(0.0);

    for final_price in final_prices {
        let payoff = (final_price - k).max(FT::<f64>::cst(0.0)); // Payoff for a call option
        payoff_sum += payoff;
    }

    let average_payoff = payoff_sum / FT::<f64>::cst(num_simulations as f64);
    (-r * t).exp() * average_payoff // Discount the average payoff back to present value
}


/// Function to calculate the gamma of the option through automatic differentiation
fn calculate_gamma(x: &[FT<f64>]) -> FT<f64> {
    let base_s0 = x[0];
    let k = x[1].to_f64().unwrap();
    let mu = x[2];
    let r = x[3];
    let sigma = x[4];
    let t = x[5];
    let timestep = x[6].to_f64().unwrap();
    let num_simulations = x[7].to_f64().unwrap() as usize;
    let seed = x[8].to_f64().unwrap() as u64;

    // Define the price function as a closure
    let price_fn = |vars: &[FT<f64>]| -> FT<f64> {
        let s0 = vars[0];
        let mu = vars[2];
        let r = vars[3];
        let sigma = vars[4];
        let t = vars[5];

        price_european_call(s0, k, mu, r, sigma, t, timestep, num_simulations, seed)
    };

    // Convert x to f64 for differentiation
    let x_f64: Vec<f64> = x.iter().map(|v| v.to_f64().unwrap()).collect();
    let price: FT<f64> = price_fn(&x);
    println!("Price of the option: {}", price.to_f64().unwrap());
    // Compute delta as the gradient of the price function
    let delta_fn = |vars: &[FT<f64>]| -> FT<f64> {
        let grad_result = grad(price_fn, &x_f64);
        grad_result[0].into() // Return derivative with respect to s0
    };

    // Calculate gamma as the gradient of the delta function
    let gamma_result = grad(delta_fn, &x_f64);
    gamma_result[0].into() // Return the second derivative (gamma)
}

/// Main function to demonstrate Greek calculation
fn main() {
    // Parameters for the simulation
    let s0 = 100.0; // Initial stock price
    let k = 100.0;  // Strike price
    let mu = 0.01;  // Drift (expected return)
    let r = 0.01;   // Risk-free interest rate (1%)
    let sigma = 0.2; // Volatility (20%)
    let t = 1.0;    // Time to maturity (1 year)
    let timestep = 365.0; // Number of time steps per year
    let num_simulations = 1000; // Number of Monte Carlo simulations
    let seed = 42;  // Seed for random number generator

    // Setup for differentiation
    let v_gamma = vec![
        F::var(s0), // Differentiate with respect to s0
        F::cst(k),
        F::var(mu),
        F::var(r),      // Mark as a variable for differentiation
        F::var(sigma),  // Mark as a variable for differentiation
        F::var(t),      // Mark as a variable for differentiation
        F::cst(timestep),
        F::cst(num_simulations as f64),
        F::cst(seed as f64),
    ];

    // Calculate Gamma
    let gamma = calculate_gamma(&v_gamma);
    println!("Gamma: {}", gamma.to_f64().unwrap());
}

The results are:

 Finished `dev` profile [unoptimized + debuginfo] target(s) in 0.01s
     Running `target/debug/bs-option`
Price of the option: 7.975722396851464
Gamma: 0

Maybe you should use autograd - Rust

AFAIK, There is no interface allow autodiff calculates high order grad

Ah really? Curious now! I will give it a go.