The mathematical reason behind "time-like partial derivatives" in Local Vol Surface, Dupire formula

lballabio / QuantLib

The QuantLib C++ library

http://quantlib.org

Other

5.24k stars 1.78k forks source link

The mathematical reason behind "time-like partial derivatives" in Local Vol Surface, Dupire formula #1124

Closed qiqvo closed 2 years ago

qiqvo commented 3 years ago

Hi all,

I was going through the testing of localvolsurface and arrived at the "decreasing variance at strike" error. I have two comments here:

During the time derivative calculation we do:
```
        Real strikept = strike*dr*dqpt/(drpt*dq);
        Real strikemt = strike*dr*dqmt/(drmt*dq);

        wpt = blackTS_->blackVariance(t+dt, strikept, true);
        wmt = blackTS_->blackVariance(t-dt, strikemt, true);

       dwdt = (wpt-wmt)/(2.0*dt);
```
I don't think this is correct: what is the mathematical reason for shifting the strikes while calculating the time derivative? And even if it has some reason to it, then there should be an additional multiplier --- i.e. a strike derivative wrt time (due to the chain rule). This is a discussion which produced the thing; and commits one and two.
Unfortunately this does not fix my problem, since it might have some financial reason to assume that the blacks variance (= implied vols squared * time) along the "adjusted" strikes should be increasing. Would you know if I could find more info on this: for example, how does this produce an arbitrage; what are the ways to handle such cases?

Thanks very much, and don't hesitate to address me to go somewhere else to look for answers. Best, Oleksandr

boring-cyborg[bot] commented 3 years ago

Thanks for posting! It might take a while before we look at your issue, so don't worry if there seems to be no feedback. We'll get to it.

klausspanderen commented 3 years ago

The implementation follows Gatheral's book "The Volatility Surface", formula 1.10. If you prefer the explicit use of the chain rule then the code can be rewritten into



        Real strikept = strike*dr*dqpt/(drpt*dq);
        Real strikemt = strike*dr*dqmt/(drmt*dq);
        Real dstrikedt = (strikept - strikemt)/(2*dt);

        wpt = blackTS_->blackVariance(t+dt, strike, true);
        wmt = blackTS_->blackVariance(t-dt, strike, true);

        Real wpu = blackTS_->blackVariance(t, strike*(1+1e-6), true);
        Real wmu = blackTS_->blackVariance(t, strike*(1-1e-6), true);

        dwdt = (wpt-wmt)/(2.0*dt) + (wpu-wmu)/(strike*2*1e-6)*dstrikedt;


2. The Dupire formula has a proven track record to exhibit "numerical issues" but also to expose arbitrage conditions in the input surtace (here in particular calendar arbitrage). Often the issues show up in the far ITM or OTM wings. In this case you might even ignore the exception and use a default value, if the Dupire formula fails. If this is not acceptable you might consider to use more evolved numerical methods like Andreas Huge interpolation etc. IMO there is no easy solution to this problem.

qiqvo commented 3 years ago

Hi Klaus,

1.1. It's only partially true about Gatheral's book: I don't see him discussing shifting strikes (as you do in the code) any near 1.10. He defines: w = w(spot, strike, time) and does not mention that strike is time-dependent, is it implied by the context? Maybe it's somewhere in the book? 1.2. Yes, this was exactly my point about the chain rule. I'm not sure if this term is negligible.

Thank you, this is very interesting.

Best, Oleksandr

klausspanderen commented 3 years ago

Hi,

in formula 1.10 it is $\omega=\omega(S_0, y_t, t), \text{with}\ \ y_t=\ln \frac{K}{F_t}$ . If you prefer to take the derivatives with respect to K and t rather than with respect to y_t and t then please have a look into formula 46 in Fabrice's paper or 2.23 in Roel's thesis. In here the numerator is the coming from the chain rule when calculating

$\frac{d\omega}{dt}=\frac{\partial \omega}{\partial t} + \frac{\partial \omega}{\partial K}\frac{\partial K}{\partial y_t}\frac{\partial y_t}{t} = \frac{\partial \omega}{\partial t} + (r_t-q_t)K\frac{\partial \omega}{\partial K}$

and this corresponds to the implementation in QuantLib. If you want to test your implementation I've added another test to Luigi's master, HestonModelTest::testLocalVolFromHestonModel. In this test a local volatility surface is derived from a Heston model and afterwards it is ensured that the local volaitlity option prices match the semi-analytic Heston prices. The test is very sensitive to the questions discussed above.

regards Klaus

qiqvo commented 3 years ago

Hi Klaus,

Thank you for setting up a test! I will try to run it over the weekend, but I already see the reason behind the original question after writing some formulas down. I have another question though with line 101: d2wdy2 = (wp-2.0*w+wm)/(dy*dy); I would say that this actually gives: $\frac{w(K e^{dy}) - 2 w(K e^{dy}) + w(K e^{-dy})}{dy^2} = \frac{w(K e^{dy}) - 2 w(K e^{dy}) + w(K e^{-dy})}{K^2 \Big(e^{dy} - 1\Big)^2}\frac{K^2 \Big(e^{dy} - 1\Big)^2}{dy^2} = \frac{\partial^2 w}{\partial K^2} K^2$ whereas the correct derivative is $\frac{\partial^2 w}{\partial y^2} = \frac{\partial}{\partial K} \Big(\frac{\partial w}{\partial K} \cdot \frac{\partial K}{\partial y} \Big)\cdot \frac{\partial K}{\partial y} = \frac{\partial^2 w}{\partial K^2} K^2 + \frac{\partial w}{\partial K} K.$

Also, I've checked the first derivative wrt y and it seems correct. So the correct line in my opinion should be: d2wdy2 = (wp-2.0*w+wm)/(dy*dy) + dwdy;

I'd rather prefer to see the derivatives wrt to plain strike and time, since then it would be straightforward to substitute them with analytically calculated ones. For example using sabr model. Does it sound like something people do in practice or in research?

Regards, Oleksandr

klausspanderen commented 3 years ago

Hi Oleksandr,

the current line 101 is actually correct with

$y + \Delta y = \ln\frac{K+\Delta K}{F} \Rightarrow K + \Delta K = Fe^{y+\Delta y} = Ke^{\Delta y},\ K -\Delta K = Ke^{-\Delta y} \\[0.5cm] \Rightarrow \frac{\partial^2 \omega(K=Fe^y)}{\partial y^2}\approx \frac{\omega(K e^{\Delta y})-2\omega(K)+\omega(K e^{-\Delta y})}{\Delta y^2}$

I agree, an implementation w.r.t. plain strike and time would be easier to read (but might be a bit slower at runtime). Maybe you want to rewrite the method with derivatives in plain strike and time and raise a pull request.

best regards Klaus

github-actions[bot] commented 3 years ago

This issue was automatically marked as stale because it has been open 60 days with no activity. Remove stale label or comment, or this will be closed in two weeks.

github-actions[bot] commented 2 years ago

This issue was automatically closed because it has been stalled for two weeks with no further activity.