Example 3: Monte Carlo Tutorial

This example demonstrates how AADC can be applied to Monte Carlo simulations, showing how to efficiently compute both option values and their derivatives using stochastic methods. We’ll estimate the same Black-Scholes option value from Example 2, but using Monte Carlo simulation instead of the analytical formula.

Monte Carlo Method for Option Pricing

To estimate the integral \(f(S(0), K, r, \sigma, T) = \mathbb{E}[e^{-rT}(S(T)-K)^{+}]\) from the previous example, we use the following Monte Carlo procedure:

Generate independent normal variables \(Z_i, i = 1, \ldots, n\)
Simulate stock prices at maturity: \(S_i(T) = S(0)\exp\left(\left(r-\frac{1}{2}\sigma^2\right)T + \sigma \sqrt{T} Z_i\right)\)
Calculate discounted payoffs: \(C_i = e^{-rT}(S_i(T)-K)^{+}\)
Estimate the integral: \(\frac{C_1 + \cdots + C_n}{n}\)

Variable Type Strategy

For Monte Carlo simulations, we use different variable marking strategies based on the parameter roles:

aadc::AADCArgument random_sample_arg(random_sample.markAsInputNoDiff());
aadc::AADCArgument S0_arg(S0.markAsDiff());
aadc::AADCArgument K_arg(K.markAsDiff());
aadc::AADCArgument r_arg(r.markAsDiff());
aadc::AADCArgument vol_arg(vol.markAsDiff());
aadc::AADCArgument T_arg(T.markAsDiff());

Variable Marking Methods

markAsInput(): Use when the variable is a function argument AND you need derivatives with respect to it
markAsDiff(): Use when the variable value is fixed across all simulations, but you need derivatives (model parameters)
markAsInputNoDiff(): Use when the variable changes between simulations but derivatives are not needed (random samples)

Why This Strategy?

Random samples: Marked as InputNoDiff because they vary between simulations but we don’t need derivatives w.r.t. random noise
Model parameters: Marked as Diff because they’re constant during simulation but we need sensitivities (Greeks)

Complete Implementation

void exampleMonteCarlo()
{
    typedef __m256d mmType;

    // Random number generation setup
    std::mt19937_64 gen;  
    std::normal_distribution<> norm_distrib(0.0, 1.0);
    
    // Option parameters
    idouble S0 = 100, K = 110, r = 0.1, vol = 0.25, T = 1, random_sample;
    
    // Simulation parameters
    int num_sim = 1000000;
    int count = aadc::mmSize<mmType>();
    int num_sim_avx = (int)floor(num_sim / count);
       
    aadc::AADCFunctions<mmType> aad_funcs;

    aad_funcs.startRecording();
    
    // Mark variables according to their roles
    aadc::AADCArgument random_sample_arg(random_sample.markAsInputNoDiff());
    aadc::AADCArgument S0_arg(S0.markAsDiff());
    aadc::AADCArgument K_arg(K.markAsDiff());
    aadc::AADCArgument r_arg(r.markAsDiff());
    aadc::AADCArgument vol_arg(vol.markAsDiff());
    aadc::AADCArgument T_arg(T.markAsDiff());

    // Record the Monte Carlo payoff calculation
    idouble S_T = S0 * std::exp((r - vol*vol / 2) * T + 
                                vol * std::sqrt(T) * random_sample);
    idouble price = std::exp(-r*T) * std::max(S_T - K, 0.0);

    aadc::AADCResult price_res(price.markAsOutput());
    aad_funcs.stopRecording();

    std::shared_ptr<aadc::AADCWorkSpace<mmType>> ws(aad_funcs.createWorkSpace());
    
    // Initialize integral accumulators using vector intrinsics
    mmType integral_S0(mmSetConst<mmType>(0));
    mmType integral_K(mmSetConst<mmType>(0)); 
    mmType integral_r(mmSetConst<mmType>(0));
    mmType integral(mmSetConst<mmType>(0));
    
    // Main simulation loop
    for (int i = 0; i < num_sim_avx; ++i) {
        
        // Generate random samples for vectorized execution
        double* _x = ws->valp(random_sample_arg);
        for (int ci = 0; ci < count; ++ci) {
            _x[ci] = norm_distrib(gen);
        }
            
        // Forward pass: compute option values
        aad_funcs.forward(*ws);
        
        // Accumulate function values
        integral = aadc::mmAdd(ws->val(price_res), integral);
        
        // Prepare for derivative computation
        ws->resetDiff();  // TODO: Check. This should not be necessary
        ws->setDiff(price_res, 1.0);

        // Reverse pass: compute derivatives
        aad_funcs.reverse(*ws);        
        
        // Accumulate derivatives
        integral_S0 = mmAdd(ws->diff(S0_arg), integral_S0);
        integral_K = mmAdd(ws->diff(K_arg), integral_K);
        integral_r = mmAdd(ws->diff(r_arg), integral_r);
    }

    // Compute final estimates
    double final_price = aadc::mmSum(integral) / num_sim;
    double delta = aadc::mmSum(integral_S0) / num_sim;
    double strike_sensitivity = aadc::mmSum(integral_K) / num_sim;
    double rho = aadc::mmSum(integral_r) / num_sim;
    
    // Print results
    std::cout << "f(100,110,0.1,0.25,1)=" << final_price
              << ". df/dS0=" << delta 
              << "; df/dK=" << strike_sensitivity 
              << "; df/dr=" << rho << std::endl;
}

Key Implementation Details

Vectorization Strategy

The code adjusts the number of simulations based on the vector width:

int num_sim = 1000000;
int count = aadc::mmSize<mmType>();          // 4 for AVX2, 8 for AVX512
int num_sim_avx = (int)floor(num_sim / count);

This ensures we process exactly count random samples simultaneously in each iteration.

Vector Intrinsics Usage

AADC provides vector operation utilities:

mmSetConst<mmType>(value): Create vector with all elements set to value
mmAdd(a, b): Element-wise vector addition
mmSum(vector): Sum all elements in the vector

These intrinsics are more efficient than manual loops:

// Efficient vector operation
integral_S0 = mmAdd(ws->diff(S0_arg), integral_S0);

// Equivalent but slower loop
for (int ci = 0; ci < count; ++ci) {
    partialSum += ws->diffp(S0_arg)[ci];
}

Memory Management Pattern

The simulation follows this pattern for each iteration:

Generate random samples directly into workspace memory
Execute forward pass to compute option values
Accumulate function values for final average
Reset adjoints to prepare for derivative computation
Set unit adjoint seed for the output
Execute reverse pass to compute all derivatives
Accumulate derivative values for final averages

Expected Output

f(100,110,0.1,0.25,1)=10.1724. df/dS0=0.55731; df/dK=-0.414169; df/dr=45.5586

Comparison with Analytical Results

Comparing with the Black-Scholes analytical result from Example 2:

Method	Option Value	Delta	Strike Sens.	Rho
Analytical	10.1601	0.557155	-0.41414	45.5554
Monte Carlo	10.1724	0.55731	-0.414169	45.5586
Difference	+0.0123	+0.000155	-0.000029	+0.0032

The Monte Carlo estimates are very close to the analytical values, with differences due to simulation noise that decreases with more samples.

Performance Advantages

AADC Benefits for Monte Carlo

Single kernel compilation: Record the payoff calculation once, execute millions of times
Vectorized execution: Process 4 (AVX2) or 8 (AVX512) paths simultaneously
Automatic differentiation: Compute all Greeks simultaneously without additional simulation
JIT optimization: Compiled kernel runs much faster than interpreted loops

Scaling Considerations

Memory efficiency: Fixed workspace size regardless of simulation count
CPU utilization: Vector instructions maximize processor throughput
Derivative efficiency: All Greeks computed in single reverse pass per iteration

Extensions and Variations

Increasing Accuracy

int num_sim = 10000000;  // More samples for higher accuracy

Multiple Assets

// Extend to basket options with correlation
idouble S1_T = S1_0 * std::exp((r - vol1*vol1/2)*T + vol1*std::sqrt(T)*Z1);
idouble S2_T = S2_0 * std::exp((r - vol2*vol2/2)*T + vol2*std::sqrt(T)*Z2);
idouble payoff = std::exp(-r*T) * std::max(S1_T + S2_T - K, 0.0);

Path-Dependent Options

// Asian option example
idouble running_average = 0.0;
for (int step = 0; step < num_steps; ++step) {
    // Generate path...
    running_average += S_t / num_steps;
}
idouble payoff = std::exp(-r*T) * std::max(running_average - K, 0.0);

Next Steps

This example demonstrates AADC’s effectiveness for stochastic simulations. The next tutorial will show how to compute full Jacobian matrices for systems with multiple inputs and outputs.