Bachelier Options Pricer Issue

[broadhed]

Describe the bug
Options priced using the modified Bachelier algorithm are inaccurate against reference implementations.

To Reproduce
Look at the unit tests for ModifiedBachelier which has an answer of 2.513031723793472

Expected behavior
Using CME QuikStrike and PyFENG get a different answer:

>>> import pyfeng as pf
>>> (pf.Norm(sigma=0.2, intr=0.05, divr=0.00, is_fwd=True)).price(strike=100.0, spot=100.0, texp=1, cp=1)
np.float64(0.07589712715905148)  # not 2.513 as in the tests

Additional context
I believe we have solved it, and it explains why there is such a large variance between these values with a large strike price such as in your example. Will follow up with a PR.

[avhz]

The formulas used currently are:

Bachelier Model

$$ C = (S - K) \Phi(d) + \sigma \sqrt{T} \phi(d) $$

$$ P = (K - S) \Phi(-d) + \sigma \sqrt{T} \phi(d) $$

$$ d = \frac{S - K}{ \sigma \sqrt{T} } $$

This does not account for time value of money, so instead:

Modified Bachelier Model

$$ \begin{aligned} C &= S\Phi(d) - K e^{-rT} \Phi(d) + \sigma \sqrt{T} \phi(d) \\ &= (S - K e^{-rT}) \Phi(d) + \sigma \sqrt{T} \phi(d) \\ \end{aligned} $$

$$ \begin{aligned} P &= K e^{-rT} \Phi(-d) - S\Phi(-d) + \sigma \sqrt{T} \phi(d) \\ &= (K e^{-rT} - S)\Phi(-d) + \sigma \sqrt{T} \phi(d) \\ \end{aligned} $$

$$ d = \frac{S - K}{ \sigma \sqrt{T} } $$

In terms of the T-Forward price

I am guessing that PR #268 implements the Bachelier model in terms of the forward price.

The paper here https://arxiv.org/pdf/2104.08686 has the explanation of the forward representation and the relation to the spot representation on page 3.

If you could expand on the formula used in the PR that would be great (where it's from, if forward, etc), that would be helpful.

Some output from Quick Strike would also be great.

[broadhed]

The Bachelier I'm using is from page 14 of the paper you linked:

$(F_{0}-K)N(d_{N}+\sigma_{N}\sqrt{T}n(d_{N}))$.

I believe this brings it into alignment with the other library I mentioned which is an academic implementation that cites the same paper. I'm proposing this PyFENG implementation is correct because it agrees with QuikStrike.

I unfortunately don't think I can share any quantities of QuikStrike data but you should be able to sign up yourself. Once in, I click on the calculator in the upper right, QuikCalculator, switch the model type to Bachelier, and put in the theoretical pricing featured in unit tests for ModifiedBachelier which has an answer of 2.513031723793472:

Instead this yields a Price 0.08, the same as the PyFENG answer, rounded to the significant digits used in the same theoretical option:

>>> import pyfeng as pf
>>> (pf.Norm(sigma=0.2, intr=0.05, divr=0.00, is_fwd=True)).price(strike=100.0, spot=100.0, texp=1, cp=1)
np.float64(0.07589712715905148)  # not 2.513 as in the tests

Perhaps these are different forms and my proposed modification needs to be a separate case? I can't recreate the 2.51 in the current tests using any formulation I've found, so that's why I thought the implementation in ModifiedBachelier was wrong. Another way you could look at it would be investigating the PyFENG notebooks. They have one for the paper you cited here: https://github.com/PyFE/PyfengForPapers/blob/main/ipynb/ChoiEtAl2022-Fut-BachelierModel.ipynb .

[avhz]

The discount factor should be removed in the PR to match the paper.

match self.option_type {
    TypeFlag::Call => ((S - K) * (-r * T).exp()) * n.cdf(d1) + v * T.sqrt() * n.pdf(d1),
    TypeFlag::Put => ((K - S) * (-r * T).exp()) * n.cdf(-d1) + v * T.sqrt() * n.pdf(-d1),
}

should be:

match self.option_type {
    TypeFlag::Call => (S - K) * n.cdf(d1) + v * T.sqrt() * n.pdf(d1),
    TypeFlag::Put => (K - S) * n.cdf(-d1) + v * T.sqrt() * n.pdf(-d1),
}

assuming S represents the forward.

Edit: Note that the representation on p.14 is un-discounted.