Calibrating Bass Local Volatility with Multiple Marginals

# Bass Local Volatility (Bass construction with multi-marginals) Calibration Problem


# Bass Local Volatility (Bass construction with multi-marginals) Calibration Problem












I'm trying to implement the Bass Local Volatility model (Bass construction with multi-marginals) but have encountered a problem.

References:

- Antoine Conze and Henry-Labordere, Bass Construction with Multi-Marginals: Lightspeed Computation in a New Local Volatility Model: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3853085

- Antoine Conze and Henry-Labordere, A new fast local volatility model: https://www.risk.net/media/download/1079736/download

I will refer to formulas from the second article.

The code that I found and adapted is from: https://github.com/igudav/Bass-Local-Volatility

Problem: The main problem is convergence of numerical calibration algorithm for exact case - lognormal model.

Numerical calibration algorithm:

- Initial guess $F^{(0)}_{W_{T_i}}$: use solution of the linearised fixed-point equation (equation (12) in the article): $$F^{-1}(u) = \sqrt{\frac{\Delta}{2}} \int_{\frac{1}{2}}^{u} dy \sqrt{\frac{G_2'(y)}{\int_{0}^{y} (G_1(z) - G_2(z)) dz}},$$ where $\Delta := T_{2}-T_{1}$ , $G_{i}(y) := F^{-1}_{\mu_{i}}(y)$, $i=1,2$.

- Iterate until convergence (in $L^\infty$ - norm) the equation: $$F^{(p+1)}_{W_{T_i}} = \mathcal{A}_i F^{(p)}_{W_{T_i}},$$ where $\mathcal{A}F := F_{\mu_{1}}\circ\big(K_{T_{2}-T_{1}}\star\big(F^{-1}_{\mu_{2}}\circ(K_{T_{2}-T_{1}}\star F)\big)\big)$, convolution $\star$ with heat kernel $K_{t}(x):=e^{-\frac{x^{2}}{2t}}/\sqrt{2\pi t}$ and $F_{\mu_i}$ is $T_i$-marginal cdf.

- All one-dimensional functions are stored and evaluated as linear or spline interpolations, and convolution with heat kernel is done using a Gauss-Hermite quadrature.

In lognormal case the exact solution is $F_{W_{T_1}} = N(\cdot / \sqrt{T_1})$. For interpolation scipy.interpolate.PchipInterpolator is used. Results showing the problem for tenors $T_1=2.0$, $T_2=3.0$ and $T_1=2.0$, $T_2=7.0$:

Figure 1.

Figure 2.

Debug Firstly I thought problem is implementation of solution (equation (12)) of the linearised fixed-point equation (equation (4) in the article). I checked implementation by substituting the solution into linearised fixed-point equation and it appears to be correct ($G_{i}(y) := F^{-1}_{\mu_{i}}(y)$):

Figure 3.

I also checked how the exact solution transforms after N successive applications of the operator $A$ , expecting the equality $F=AF$ (Theorem 1.) to be satisfied. I found an error accumulation after successive application of convolution and interpolation of the obtained result, but even if this is what causes the error, I don't understand how it can be fixed

Figure 4.

Below is a minimal code snippet to reproduce the issue:

```
import numpy as np  
  
from scipy.stats import norm  
  
from typing import Union, List, Callable  
from numpy.typing import NDArray  
from numpy.polynomial.hermite import hermgauss  
from scipy.interpolate import PchipInterpolator, CubicSpline  
  
FloatOrVectorType = Union[float, List[float], NDArray[float]]  
FloatVectorType = Union[List[float], NDArray[float]]  
EPS = np.finfo(float).eps  
  
  
def convolutionWithGaussHermiteQuadrature(  
        t: float,  
        func: Callable[[NDArray], NDArray],  
) -> Callable[[NDArray], NDArray]:  
    nodes, weights = hermgauss(61)  
    def f(x: FloatOrVectorType):  
        shape = (-1, *[1] * x.ndim)  
        newVariable = x[None] - np.sqrt(2 * t) * nodes.reshape(shape)    
        return 1 / np.sqrt(np.pi) * np.sum(  
            weights.reshape(shape) * func(newVariable),  
            axis=0  
        )  
  
    return f  
  
  
class LogNormalMarginal:  
  
    def __init__(self, sigma: float, tenor: float):  
        self._sigma = sigma  
        self._tenor = tenor  
  
    @property  
    def tenor(self):  
        return self._tenor  
  
    @property  
    def sigma(self):  
        return self._sigma  
  
    def inverseCdf(self, u: FloatOrVectorType) -> FloatOrVectorType:  
        return self._inverseCdf(np.clip(u, a_min=EPS, a_max=1 - EPS))  
  
    def cdf(self, x: FloatOrVectorType) -> FloatOrVectorType:  
        return np.clip(self._cdf(x), a_min=EPS, a_max=1 - EPS)  
  
    def pdf(self, x: FloatOrVectorType) -> FloatOrVectorType:  
        return 1 / (x * self._sigma * np.sqrt(self.tenor)) \  
            * norm.pdf(  
                (np.log(x) + 0.5 * self._sigma ** 2 * self.tenor) \  
                / self._sigma / np.sqrt(self.tenor)  
            )  
  
    def _cdf(self, x: FloatOrVectorType) -> FloatOrVectorType:  
        return norm.cdf(  
            (np.log(x) + 0.5 * self._sigma ** 2 * self.tenor) \  
            / self._sigma / np.sqrt(self.tenor)  
        )  
  
    def _inverseCdf(self, u: FloatOrVectorType) -> FloatOrVectorType:  
        return np.exp(  
            - 0.5 * self._sigma ** 2 * self.tenor \  
            + self._sigma * np.sqrt(self.tenor) * norm.ppf(u)  
        )  
  
    def integralOfInverseCdf(self, u: FloatOrVectorType) -> FloatOrVectorType:  
        d_1 = (-np.log(self._inverseCdf(u)) + 0.5 * self._sigma ** 2 * self.tenor) \  
            / self._sigma / np.sqrt(self.tenor)  
        return norm.cdf(-d_1)  
  
    def derivativeOfInverseCdf(self, u: FloatOrVectorType) -> FloatOrVectorType:  
        return 1 / self.pdf(self._inverseCdf(u))  
  
  
class FixedPointEquation:  
  
    @staticmethod  
    def getInitialIterationAsSolutionOfLinearisedFixedPointEquation(  
            marginal1: LogNormalMarginal,  
            marginal2: LogNormalMarginal,  
            gridPoints: int = 2001,  
            bounds: float = 5.  
    ) -> PchipInterpolator:  
        # eq. (12)  
        uGrid = norm.cdf(np.linspace(-bounds, bounds, gridPoints))  
  
        integrand = np.sqrt(  
            marginal2.derivativeOfInverseCdf(uGrid) / (  
                    marginal1.integralOfInverseCdf(uGrid)  
                    - marginal2.integralOfInverseCdf(uGrid) + EPS  
            )  
        )  
        integral = CubicSpline(uGrid, integrand).antiderivative()  
        inverseCdfValues = np.sqrt(  
            (marginal2.tenor - marginal1.tenor) / 2  
        ) * (integral(uGrid) - integral(1 / 2))  
  
        return PchipInterpolator(inverseCdfValues, uGrid)  
  
    @staticmethod  
    def getMappingFunction(  
            solution: PchipInterpolator,  
            marginal1: LogNormalMarginal,  
            marginal2: LogNormalMarginal  
    ) -> Callable[[FloatVectorType], FloatVectorType]:  
        # eq. (3) 
        internalConvolution = \  
            convolutionWithGaussHermiteQuadrature(  
                t=marginal2.tenor - marginal1.tenor,  
                func=solution  
        )  
  
        def applySecondMarginalInverseCdf(x):  
            u = internalConvolution(x)  
            return marginal2.inverseCdf(u)  
  
        externalConvolution = \  
            convolutionWithGaussHermiteQuadrature(  
                t=marginal2.tenor - marginal1.tenor,  
                func=applySecondMarginalInverseCdf  
        )  
        return externalConvolution  
  
    @classmethod  
    def applyOperatorA(  
            cls,  
            solution: PchipInterpolator,  
            marginal1: LogNormalMarginal,  
            marginal2: LogNormalMarginal  
    ) -> Callable[[FloatVectorType], FloatVectorType]:  
  
        def applyFirstMarginalCdf(x):  
            mappingValues = cls.getMappingFunction(  
                solution=solution,  
                marginal1=marginal1,  
                marginal2=marginal2  
            )(x)  
            return marginal1.cdf(mappingValues)  
  
        return applyFirstMarginalCdf  
  
    @staticmethod  
    def LInfinityNorm(  
            sequence1: FloatVectorType,  
            sequence2: FloatVectorType  
    ) -> float:  
        return np.max(np.abs(sequence1 - sequence2))  
  
    @classmethod  
    def solveFixedPointEquation(  
            cls,  
            marginal1: LogNormalMarginal,  
            marginal2: LogNormalMarginal,  
            maxIter: int = 61,  
            tol: float = 1e-4  
    ):  
        wGrid = np.linspace(-5, 5, 2001) * np.sqrt(marginal1.tenor)  
        solutionIteration = \  
            cls.getInitialIterationAsSolutionOfLinearisedFixedPointEquation(  
                marginal1=marginal1,  
                marginal2=marginal2,  
        )  
  
        for iterationNumber in range(maxIter):  
            solutionNextIteration = PchipInterpolator(  
                x=wGrid,  
                y=cls.applyOperatorA(  
                    solution=solutionIteration,  
                    marginal1=marginal1,  
                    marginal2=marginal2  
                )(wGrid)  
            )  
  
            lInfty = cls.LInfinityNorm(  
                sequence1=solutionNextIteration(wGrid),  
                sequence2=solutionIteration(wGrid)  
            )  
            solutionIteration = solutionNextIteration  
            if lInfty < tol:  
                break  
            if iterationNumber == (maxIter - 1):  
                raise Exception(  
                    f'Convergence error: current tol: {tol}, reach {lInfty}, iter: {iterationNumber}'  
                )  
  
        print(f'Current tol: {tol}, reach: {lInfty}, iter: {iterationNumber}')  
        return solutionIteration  
  
  
if __name__ == '__main__':  
    # Let S_0 = 1.  
    import matplotlib  
    import matplotlib.pyplot as plt  
  
    expiries = [2., 3.]  
    vols = [0.2] * 2  
  
    marginal1, marginal2 = [  
        LogNormalMarginal(sigma=sigma, tenor=T)  
        for sigma, T in zip(vols, expiries)  
    ]  
  
    exactSolutionOfFixedPointEq = lambda x: norm.cdf(x / np.sqrt(marginal1.tenor))  
  
    numericalSolution = FixedPointEquation().solveFixedPointEquation(  
        marginal1=marginal1,  
        marginal2=marginal2,  
        maxIter=61,  
        tol=2e-4  
    )  
    x = np.linspace(-5 * np.sqrt(marginal1.tenor), 5 * np.sqrt(marginal1.tenor), 1000)  
    print(np.abs(numericalSolution(x) - exactSolutionOfFixedPointEq(x)).max())  
  
    print(numericalSolution(0.))  
  
    _, (generalAxis, diffAxis) = plt.subplots(nrows=1, ncols=2, figsize=(15, 8), dpi=120)  
    generalAxis.plot(x, numericalSolution(x), label='numericalSolution')  
    generalAxis.plot(x, exactSolutionOfFixedPointEq(x), label='exact')  
    generalAxis.set_title(  
        f"Solutions of fixed point equation, marginal tenors = {expiries}"  
    )  
    generalAxis.legend()  
    generalAxis.grid()  
  
    diffAxis.plot(x, numericalSolution(x) - exactSolutionOfFixedPointEq(x), label='diff')  
    diffAxis.legend()  
    diffAxis.set_title('Difference between the numerical solution and the exact')  
    diffAxis.grid()  
    plt.show()
```

I tried change interpolation methods, grids, verify analytical formulas and double check it by numerical integration and taking the derivative, but I'm completely stuck and don't know how to fix the error, any suggestions would be greatly appreciated.