Extending black-scholes example about fastad HOT 9 CLOSED

Screw letting it rest :joy_cat: Very nice:

// [[Rcpp::export]]
Rcpp::NumericMatrix black_scholes(double spot = 105,   		// current spot
                                  double strike = 100, 		// strike
                                  double vol = 5,               // annual vol (500% !!)
                                  double r = 1.25 / 100,        // short-term rate
                                  double tau = 30.0 / 365) {	// time remaining
    ad::Var<double> Delta(spot);
    ad::Var<double> Vol(vol);
    ad::Var<double> Rho(r);
    ad::Var<double> Theta(tau);
    std::vector<ad::Var<double>> cache;

    Rcpp::NumericMatrix M(2,5);
    Rcpp::rownames(M) = Rcpp::CharacterVector::create("call", "put");
    Rcpp::colnames(M) = Rcpp::CharacterVector::create("value", "delta", "vega", "rho", "theta");

    auto call_expr = ad::bind(black_scholes_option_price<option_type::call>(Delta, strike, Vol, Rho, Theta, cache));
    double call_price = ad::autodiff(call_expr);
    M.row(0) = Rcpp::NumericVector::create(call_price, Delta.get_adj(), Vol.get_adj(), Rho.get_adj(), Theta.get_adj());

    // reset adjoints before differentiating again
    Delta.reset_adj();
    Vol.reset_adj();
    Rho.reset_adj();
    Theta.reset_adj();
    for (auto& c : cache) c.reset_adj();

    auto put_expr = ad::bind(black_scholes_option_price<option_type::put>(Delta, strike, Vol, Rho, Theta, cache));
    double put_price = ad::autodiff(put_expr);
    M.row(1) = Rcpp::NumericVector::create(put_price, Delta.get_adj(), Vol.get_adj(), Rho.get_adj(), Theta.get_adj());

    return(M);
}

from fastad.

Hm I tried playing with online calculators:

https://www.mystockoptions.com/black-scholes.cfm
https://goodcalculators.com/black-scholes-calculator/

Fastad seems to be very close to these two calculators and RQuantLib is slightly off from all 3 😮

from fastad.

Ah I see that's why - I see that you're calling:

EuropeanOption("call", 105, 100, 1.25/100, 1.25/100, 30/365, 5)

where dividendYield is set to 1.25/100. Then, I believe the example I wrote is then equivalent to dividendYield=0.

But also, sweet that it's all fine anyways!

from fastad.

Just looked at the code! ptr_pack_t takes in two arguments where the first pointer points to the beginning location of a buffer where fastad will store the node values and second pointer for the adjoints. So you only need 2 buffers that have the total size needed for the whole graph. .bind_cache_size() returns the expression graph-specific amount of space needed.
So https://github.com/eddelbuettel/rcppfastad/blob/6d6730fb25a5f8861b362a052388507bc0498545/src/black_scholes.cpp#L73-L77
needs to change to:

    auto size_pack_c = call_expr.bind_cache_size();
    std::vector<double> value_buf(size_pack_c(0));
    std::vector<double> adjoint_buf(size_pack_c(1));
    call_expr.bind_cache({value_buf.data(), adjoint_buf.data()});

Or if you want to avoid thinking about constructing these buffers yourself, you can remove those 4 lines and just do:

auto call_expr = ad::bind(black_scholes_option_price<option_type::call>(Spot, strike, Vol, tau, r, cache));

I've wrapped this logic of constructing buffers and attaching it to the expression into ad::bind.

Also, feel free to reuse the same buffers for expressions. They are always overwritten and not updated, so you don't have to reset anything like how you have to for ad::Var adjoint buffers.

from fastad.

Lovely. I saw some of the higher-level expression logic in the other examples. Makes sense to do it here too. So some refreshing of the Black-Scholes example may be due -- it for example also doesn't set a seed.

I'll poke some more.

from fastad.

Awesomesauce:

edd@rob:~/git/rcppfastad(master)$ r -lRcppFastAD -pe'black_scholes()'
       value     delta    vega      rho   theta
call 56.5136  0.773818 9.05493  2.03321 275.730
put  51.4109 -0.226182 9.05493 -6.17753 274.481
edd@rob:~/git/rcppfastad(master)$ 

Thanks for the quick help. If you want I can send a PR back to update your example.

Also looks good compare to RQuantLib with known analytical values:

> EuropeanOption("call", 105, 100, 1.25/100, 1.25/100, 30/365, 5)
Concise summary of valuation for EuropeanOption 
    value     delta     gamma      vega     theta       rho    divRho 
  56.4301    0.7728    0.0020    9.0505 -274.5820    2.0313   -6.6694 
> EuropeanOption("put", 105, 100, 1.25/100, 1.25/100, 30/365, 5)
Concise summary of valuation for EuropeanOption 
    value     delta     gamma      vega     theta       rho    divRho 
  51.4352   -0.2262    0.0020    9.0505 -274.6444   -6.1794    1.9518 
> 

I'll let what I have now settle for a day or so; I'll then look into you suggested ad::bind() use.

from fastad.

Yay! Thanks for working this out! Please feel free to submit a PR with the updated Black-Scholes example!

Also looks good compare to RQuantLib with known analytical values:

Hm.. the differences look large for some reason. Is the reason simply that the atomic functions like exp and erf are not exact?

from fastad.

I meant to ask you about that too. The gaps get bigger for different parameters i.e. shorter time to decay, more normal vol (500% is a little ... much, literally) as you can see if you hit any of the available Black-Scholes calculators (as for example the QuantLib one we have in RQuantLib). I haven't really poked around.

from fastad.

There is a technicality that sometimes comes into play: whether a dividend rate is included or not. QuantLib is pretty well respected and tested so with right parameterisation it should give the right numbers. But we're likely good here.

Also of note is that practicioners usually only care about three or so digits and we a definitely good here.

from fastad.