Adjoint Greeks VI – Source Code Transformation Pt. II

I will be talking about Automatic differentiation beyond typedef and operator overloading at this year’s QuantLib User Conference. Registration is free, but seats are limited. The blog post today is sort of a teaser, so do not miss to be in Düsseldorf on November 30th / December 1st.

To recall what has been done to make QuantLib AD-ready: We have two approaches. One is by CompatibL which in essence changes the typedef for Real in the library from double to some active type of an operator overloading tool, like CppAD::AD. Of course this doesn’t work just by relinking the typedef, but in effect it describes well what is going on.

The other approach is to transform the library’s code and introduce template parameters whenever we might potentially be interested in calculating derivatives. You can read more about this by selecting the tag Automatic Differentiation in the sidebar that can be expanded by clicking on the three horizontal lines in the upper right corner when you scroll to the top of this page.

Backward compatibility in the template approach is achieved by typedef’ing the original class names to template instantiations with T = Real, so existing code compiles without any changes. Well let’s say 99% of it, some few modifications have to be applied, but really not many (from my experience with the original examples and the test-suite).

The transformation is not applied to the full library (yet), but the vanilla IRD stuff works, so coupons, yield term structures, volatility termstructures, the GSR model. You can get it (together with some working examples) from my adjoint branch at

Both approaches have clear advantages and disadvantages. One potential problem with the typedef approach is, that this might kill your performance in parts of the code where you do not even want to use AD. For example,

    Matrix_t<T> A(1024, 1024);
    Matrix_t<T> B(1024, 1024);
    Matrix_t<T> C = A * B;

will take 764 ms on my machine with T=double, but 8960 ms with T=CppAD::AD. This seems to be due to compiler optimizations (like using SIMD instructions) that are no longer applicable when switching from native double to the custom AD type. This is a penalty of 11.7x and you do not get anything for that. Note that I do not record any tape or something, I just use the double functionality wrapped in the AD type.

Of course there might be other operator overloading tools that are more efficient. To get a lower bound for the penalty factor I replaced T by a minimal wrapper, just storing a double and an additional native pointer and forwarding the operators +,*,+= in question here to the respective double operations. With this wrapper the penalty is 2.0x. This is much better. The perfect operator overloading tool will probably lie somewhere between 2.0 and 11.7. But even 2.0 is not good, if you can avoid it, is it?

In any case the template approach allows for selective activation of variables where AD is needed and use native instances of the library classes elsewhere. You can also mix both. For example you could do a model calibration with the double-instance of a class, then extract the calibrated parameters and plug them into the AD – version of the same model, which then maybe generates Monte Carlo paths for some CVA calculations in which you want to calculate an adjoint gradient.

Still there are use cases left where even the template approach will not help and this is when you want to apply AD to code sections like the one above, where a lot of dense double calculations take place, enabling the compiler to apply optimizations that are not possible with custom types.

We saw an example in The time equivalent of around 80 NPV calculations (holy cow!) is consumed for one adjoint gradient calculation of a Bermudan swaption in the Hull White model. From theory we would expect a complexity of 4x for the adjoint sweep, in practice we are probably happy with any factor up to 10x. But with 80 NPV calculations you can do a lot of things and even if your gradient vector has 100 – 200 entries, the application of AD would merely result in a speed up of around 2x compared to a naive bump and revalue approach. This would not really justify the whole AD-thing.

One possible way out of this is to write the adjoint code by hand.

Another approach is to use source code transformation. I started with a toy example in one of my last posts.

In the meantime I implemented a Bermudan swaption engine that produces an AD gradient w.r.t. the H and \zeta vectors describing an LGM model and the relevant discount factors of an interest rate curve. For this I extracted the computational core of the Gaussian1dSwaptionEngine (which relies on numerical integration) and rewrote it in Fortran. Then I applied OpenAD/F to the code to generate adjoint code and compiled everything into a library. You can have a look at the project here:

On the QuantLib side I wrote a pricing engine acting as a wrapper for the Fortran core: It can be used just as any other pricing engine, but provides the adjoint gradient vector as an additional result. The gradient is computed w.r.t. the model parameters, but it is easy to tranform this into the usual market variables’ sensitivities. This is done exactly the same way as in the post I mentioned earlier.

Testing the engine on a 10y Bermudan swaption with yearly calls (using 49 grid points per expiry for the numerical integration) yields the following results (you can find the example code here)

  • single pricing on a Intel(R) Core(TM) i7-2760QM CPU @ 2.40GHz, single threaded, non-transformed code: 4.2 ms
  • pricing + gradient (which is 105 dimensional in this example): 25.6 ms
  • transformation of gradient w.r.t. model parameters: 6.2 ms
  • adjoint calculation multiple: 6.1x (7.6x including gradient transformation)

Looks good, no? In particular, because no extra optimization of the input or transformed code was necessary to achieve the good multiple of 6.1, it just came out of OpenAD/F this way.

As a side note when implementing it as described here you can use the power of AD even from the original / “classic” QuantLib just through specialized pricing engines (and by linking against the Fortran support library of course). Or you can mix operator overloading in the template – QuantLib with source code transformation applied to critical sections of the code which slow down too much with operator overloading. Maybe this hybrid solution is the way to go.

More on this and other topics in Düsseldorf.

Adjoint Greeks VI – Source Code Transformation Pt. II

Dynamics II

This is a sequel to the recent post on smile dynamics, generalizing things a bit. It’s just a stream of thoughts, nothing that is actually backed by own experience or anything. Hopefully it makes sense nevertheless.

The PL of a interest rate derivative portfolio reacts on movements in market data like broker quotes for FRAs, Euribor futures, swaps, implied volatilities for caps / floors and swaptions. These quotes are processed into rate curves and volatility surfaces, which are then fed into pricing models which value the portfolio and thereby produce its PL.

A trader is interested in how his portfolio reacts on market movements, i.e. in its greeks, which are by definition partial derivatives of the portfolio value by the single market data inputs. So a portfolio may for example be particularly sensitive on the 10y and 15y Euribor 6M swap quotes and on the 5y/15y and 10y/10y swaption volatilities. The trader would see a large interest rate delta on these two buckets of the 6M curve and a large vega in these two cells of the vega matrix for swaptions.

The input variables are by far not independent. For example if the 10y swap quote rises by 10 basis points, it is likely that the 15y swap quote rises too, and often even roughly by the same amount. Similarly, movements of the cells of implied volatility surfaces are correlated, in particular if their distance is not too big. Even more, an interest rate move will possibly affect the implied volatility, too, in a systematic way: If rates rise, typically lognormal volatilities will decrease, so to roughly keep the product of the forward level and the lognormal volatility (which is by definition the so called (approximate) basispoint volatility) constant. The dependency is obviously dependent on the nature of the volatility, so it will be quite different for normal implied volatilities.

Furthermore we have technical dependencies. Say you work with sensitivites on zero rates instead of market quotes. Then a movement in the EONIA zero rates will affect the Euribor 6M forward curve, because this curve is bootstrapped using the EONIA curve as a discounting curve. So a sensitivity on the EONIA curve might come from a direct dependency of a deal on the EONIA curve, or also from an indirect dependency via forward curves.

So the abstract picture is that you have tons of correlated data points as an input, a very long gradient vector describing the sensitivity of the PL to the input, and maybe also a big Hessian matrix describing second order effects.

Wouldn’t it be reasonable to apply a principal component analysis first, select the most important risk drivers and only manage the sensitivities on these drivers ? Conceptually this is nothing more that a suitable rotation of coordinates such that the new coordinates are independent. Once this is done you can sort the dimensions by the ratio of the total variance of the input data they carry and focus on the most important ones. Which may be only 5 to 10 maybe instead of hundreds on input data points.

Hedging would be be more efficient, because you could hedge the same amount of PL variance with fewer hedge positions. The main risks drivers could be summarized in terms of only a few numbers for management reporting. Sensitivity limits can be put on these main drivers. PL could be explained accurately by movements and sensitivites of only a few quantities. Stress tests could be formulated directly in the new coordinates, yielding more realistic total scenarios, since natural co-movements of other variables not directly stressed are automatically accounted for (e.g. huge rate shifts would imply appropriate adjustments to volatilities and keep them realisitic).

If you think about it, the kind of interest rate delta you see for swaptions and caps / floors this way would include a historically calibrated smile dynamics. Which you could compare to the dynamics we discussed last time. For example.

Also, typically the new coordinates allow for a natural interpretation like parallel movements of curves, rotations and so on, so the view wouldn’t necessarily be more abstract than the original one.

What is needed is the coordinate transformation layer, ideally directly in the front office system, so that it can be used in real-time. Plus a data history has to be collected and maintained, so that the principal components can be updated on a regular basis.

Dynamics II