I get a confusing mistake in using const reference of AReal

[xuruilong100]

Describe the bug

To get Greeks of exotic options, I'm trying to use xad in Monte Carlo simulation. But I get a confusing mistake in using const reference of AReal. To show it, this is a demo below:

To Reproduce

This is my demo, I implement the example in book The Complete Guide to Option Pricing Formulas (2nd edition)(P51), the correct vega is 18.5027.

Using const auto& volatility = model_->volatility; will lead a wrong answer, but you can get a correct number by using auto volatility = model_->volatility;

It is really confusing, I can't use const reference of an AReal object.

#include <XAD/XAD.hpp>
#include <iostream>
#include <memory>
#include <random>
#include <vector>

using Real = double;
using Tape = xad::Tape<Real>;
using AD = xad::AReal<Real, 1>;

enum class OptionType
{
    Call,
    Put
};

struct BlackScholesModel
{
    AD s0;
    Real dividendYield;
    Real riskFreeRate;
    AD volatility;
};

class Simulator
{
  public:
    virtual ~Simulator() = default;
    virtual Real discountFactor() const = 0;
    virtual std::vector<AD> simulate() = 0;
    const std::vector<AD*>& parameters() const { return parameters_; }

  protected:
    std::vector<AD*> parameters_;
};

class BlackScholesSimulator : public Simulator
{
  public:
    BlackScholesSimulator(std::shared_ptr<BlackScholesModel> model, Real dt, std::size_t steps)
        : model_(std::move(model)), dt_(dt), steps_(steps)
    {
        parameters_.push_back(&model_->s0);
        parameters_.push_back(&model_->volatility);
    }

    Real discountFactor() const override { return std::exp(-model_->riskFreeRate * dt_ * steps_); }
    std::vector<AD> simulate() override;

  private:
    std::shared_ptr<BlackScholesModel> model_;
    Real dt_;
    std::size_t steps_;
    std::mt19937 rng_{42};
    std::normal_distribution<Real> norm_{0.0, 1.0};
};

std::vector<AD> BlackScholesSimulator::simulate()
{
    std::vector<AD> path(steps_ + 1);
    path.front() = xad::log(model_->s0);

    const auto& dividendYield = model_->dividendYield;
    const auto& riskFreeRate = model_->riskFreeRate;
    // const auto& volatility = model_->volatility;  // wrong vega
    auto volatility = model_->volatility;  // correct vega

    for (std::size_t i = 1; i <= steps_; ++i)
    {
        Real dw = norm_(rng_);
        AD drift = (riskFreeRate - dividendYield - 0.5 * volatility * volatility) * dt_;
        AD diffusion = volatility * std::sqrt(dt_) * dw;
        AD delta = drift + diffusion;

        path[i] = path[i - 1] + delta;
    }

    std::transform(path.begin(), path.end(), path.begin(), [](const AD& x) { return xad::exp(x); });

    return path;
}

class PathPricer
{
  public:
    virtual ~PathPricer() = default;
    virtual AD evaluate(const std::vector<AD>& path, Real discountFactor) const = 0;
};

class EuropeanPathPricer : public PathPricer
{
  public:
    EuropeanPathPricer(OptionType optionType, Real strike)
        : optionType_(optionType), strike_(strike)
    {
    }

    AD evaluate(const std::vector<AD>& path, Real discountFactor) const override
    {
        AD payoff;
        if (optionType_ == OptionType::Call)
        {
            payoff = xad::max(path.back() - strike_, 0.0);
        }
        else  // Put
        {
            payoff = xad::max(strike_ - path.back(), 0.0);
        }
        return payoff * discountFactor;
    }

  private:
    OptionType optionType_;
    Real strike_;
};

std::vector<Real> MCSimulation(const std::shared_ptr<Simulator>& simulator,
                               const std::shared_ptr<PathPricer>& pathPricer, std::size_t samples)
{
    const auto& parameters = simulator->parameters();
    std::vector<Real> results(parameters.size() + 1, 0.0);

    Tape tape;
    Real df = simulator->discountFactor();

    for (std::size_t i = 0; i < samples; ++i)
    {
        tape.clearAll();
        for (auto p : parameters)
        {
            tape.registerInput(*p);
        }
        tape.newRecording();

        auto path = simulator->simulate();
        AD price = pathPricer->evaluate(path, df);

        tape.registerOutput(price);
        xad::derivative(price) = 1.0;
        tape.computeAdjoints();

        results[0] += xad::value(price);
        for (std::size_t j = 0; j < parameters.size(); ++j)
        {
            results[j + 1] += xad::derivative(*parameters[j]);
        }
    }

    for (auto& r : results)
    {
        r /= static_cast<Real>(samples);
    }

    return results;
}

int main()
{
    Real maturity = 0.75;
    Real spot = 55;
    Real dividendYield = 3.55 / 100.0;
    Real riskFreeRate = 10.5 / 100.0;
    Real volatility = 30 / 100.0;
    Real strike = 60;

    auto model = std::make_shared<BlackScholesModel>();
    model->s0 = spot;
    model->dividendYield = dividendYield;
    model->riskFreeRate = riskFreeRate;
    model->volatility = volatility;

    std::size_t steps = 10;
    Real dt = maturity / static_cast<Real>(steps);
    auto simulator = std::make_shared<BlackScholesSimulator>(model, dt, steps);

    auto opt = OptionType::Call;
    auto pathPricer = std::make_shared<EuropeanPathPricer>(opt, strike);

    std::size_t samples = 10000 * 50;

    auto results = MCSimulation(simulator, pathPricer, samples);

    std::cout << "Option Price: " << results[0] << "\n";
    std::cout << "Delta: " << results[1] << "\n";
    std::cout << "Vega: " << results[2] << "\n";
}

Environment:

• OS: Windows 11
• Compiler: VC++ 2022
• IDE: Visual Studio 2022
• XAD version: 1.8.0

[boring-cyborg[bot]]

Thanks for opening your first issue here. We'll look into it and get back in due time.

[auto-differentiation-dev]

Hi @xuruilong100 , thanks for your report. We've analysed your code and believe this is what is happening:

In the execution of the first Monte-Carlo scenario, you are registering the volatility input within the model shared pointer with the tape. This gives it a slot on the tape and stores it in the AD variable (now both the tape and the variable "know each other"). Before the next scenario, you call tape.clearAll(), which makes the tape completely forget about all variables ever registered. However, when you call registerInput for the same volatility variable (when taken by const reference), the variable already has a slot stored in it and therefore does not register again - but the tape "forgot" about it.

If you take a copy of the input variable rather than using a const reference, this copy depends on a registered variable (a variable with a slot in it), and therefore it registers itself on the tape. Then everything is fine. This is more of a coincidence.

You have 2 clean options to fix your code:

1. Make a copy of your inputs for every scenario. You could for example have BlackScholesModel<T>, initialise with double, and copy into a BlackScholesModel<AD> before every scenario for pathwise derivatives. It's effectively the same as your working version where you take a copy first, but less confusing and safer (and also includes s0).
2. Call the registerInput() only once for the AD inputs you have, before the Monte-Carlo loop, followed by newRecording(), and save the position in the tape at that point. Then use resetTo(position), and clearDerivatives() in the Monte-Carlo paths before every loop. The resetTo makes sure that the variables and operations registered up to that position stay on the tape, while everything after is cleared, and clearDerivatives is making sure that the derivatives are all reset.

[xuruilong100]

Thank you @auto-differentiation-dev , I want to know witch option has better performance?

[auto-differentiation-dev]

We doubt there will be a measurable difference between the approaches. The only difference is the copy of the model parameters in option 1, but compared to the computation that is very small. The second option avoids that copy though, so it has a few less operations.