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.
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.