ABSTRACT

Consistently fitting vanilla option surface is an important issue in derivative modelling. In this paper, we consider three different models: local and stochastic volatility, local correlation, hybrid local volatility with stochastic rates, and address their exact, nonparametric calibration. This calibration process requires solving a nonlinear partial integro-differential equation. A modified alternating direction implicit algorithm is used, and its theoretical and numerical analysis is performed.

Keywords:

Local Stochastic Volatility, Calibration, Derivative Pricing, Partial Integro-Differential Equations

1. Introduction

One of the most important challenge for real-life applications of a model to derivatives trading is the issue of calibration. Similar to common situations in many areas of physics and engineering, once a model has been suggested, its parameters have to be estimated using external data. In the case of derivative modelling, those data are the liquid (tradable) options, generally called vanilla products.

It is well-known since the pioneering work of [1] and its celebrated extension by [2] that the knowledge of market data such as the prices of vanilla options across all strikes and maturities is equivalent to the knowledge of the risk-neutral marginals of the underlying stock distribution. Here, we are interested in appli- cations of this result to three different cases.

The first one is quite classic, and is the starting point of our work: local and stochastic volatility models. Such models are very useful in practice, since they offer both the flexibility and realistic dynamics of stochastic volatility, and the exact calibration properties of local volatility. The problem of calibrating local and stochastic volatility models has been dealt with for a while now, for instance by [3] . However, practitioners seem to agree that the stability of its resolution becomes uncertain when the volatility’s volatility is too large. We shall see that indeed some kind of instability appears, and offer explanations to the phenome- non. The second case we focus on is the correlation between assets. Empirical measures give a certain set of results. However, when modelling a basket on multiple underlyings, a problem occurs. If one uses local volatility models for each underlying and correlates their brownian motions using the empirical correlation, the basket obtained will not reproduce the vanillas quoted on the market. This raises significant issues when hedging products on multiple under- lyings. One of the solution for this problem is the known ‘local correlation’ approach: the correlation matrix for the n underlyings is deformed using a parameter, function of the time and the basket level. Here, we use that approach to obtain a calibration equation for the basket, relatively similar to the one appearing in the local and stochastic volatility model, and then numerically solve said equation in a two-underlyings framework.

The last topic we shall be interested in are interest rates, we study a hybrid model: local volatility with stochastic rates. Using a partial differential equation approach similar to the local correlation and the local and stochastic volatility, we write a calibration equation for the vanillas of this hybrid model, solve it and verify the accuracy of the fit.

The general form of our calibration equations is nonlinear partial and integro- differential. For their resolution, we chose to adapt the alternating direction implicit scheme (very efficient to solve classic linear second order parabolic equations, [4] ). Being in a nonlinear non-local framework, many questions arise. Is it relevant to use ADI algorithms to solve the equations stemming from our calibration problems? How should we deal with the nonlocal term? Is the finite difference scheme we chose consistent, and what is the order of the truncation error? Can we detect an instability in certain cases? Is it possible to explain it?

The aim of this work is to address and at least partially answer these questions. The paper is organized as follows. In Section 2, we quickly present the case of the local and stochastic volatility model. Section 3 is devoted to the local correlation, its calibration equation and the fit we obtain in the case of a basket on two underlyings. In Section 4, we do the same thing in the stochastic rates frame. Finally, in Sections 5 and 6, we adress the questions raised previously concerning the ADI algorithm used for the resolution, Section 7 is a brief conclusion.

2. Local and Stochastic Volatility Models

2.1. Partial Integro-Differential Equation for the Calibration of LSV Models

The diffusion model is assumed to be the following

$\frac{\text{d}{S}_t}{S_t}=r\left(t\right)\text{d}t+a\left(t,S_t\right)b\left(y_t\right)\text{d}{W}_t^1$ (1)

$\text{d}{y}_t=\mu \left(t,y_t\right)\text{d}t+\xi \left(t,y_t\right)\text{d}{W}_t^2$ (2)

$\left(S_t\mathrm{,}t\ge 0\right)$ is the stock price process and $\left(y_t\mathrm{,}t\ge 0\right)$ the stochastic component of the volatility. The function b simply transforms that factor into a proper volatility. a is the local volatility part of the model, exactly as in Dupire’s formula, its value shall be specified depending on the aimed vanillas. $\xi$ is the volatility of the volatility factor (commonly called “vovol”) and $\mu$ is a drift term. $W^1$ and $W^2$ are one-dimensional standard brownian motions with correlation $\rho$ .

Let us now consider a surface of vanilla prices $C\left(T\mathrm{,}K\right)$ and the corres- ponding Local Volatility ${\sigma }_D$ . Under suitable regularity and ellipticity assum- ptions, the following proposition can be proved

Proposition 1. The diffusion model defined by (1-2) has a density $p\left(t\mathrm{,}S\mathrm{,}y\right)$ with respect to Lebesgue’s measure. Moreover, if the model fits the surface of vanillas $C\left(T\mathrm{,}K\right)$ then necessarily

$a^2\left(t\mathrm{,}S\right)={\sigma }_D^2\left(t\mathrm{,}S\right)\frac{{\displaystyle {\int }_{ℝ}}p\left(t\mathrm{,}S\mathrm{,}y\right)\text{d}y}{{\displaystyle {\int }_{ℝ}}{b}^2\left(y\right)p\left(t\mathrm{,}S\mathrm{,}y\right)\text{d}y}$ (3)

Proof. The exact assumptions and the existence proof can be found in [5] . Here, the main concern is the calibration result. Let us assume that the model fits exactly the surface C. Letting $\left(S_0\mathrm{,}{y}_0\right)$ denote the initial state of the system, the joint density $p\left(t\mathrm{,}S\mathrm{,}y\right)$ of the couple $\left(S_t\mathrm{,}{y}_t\right)$ verifies Kolmogorov for- ward equation

$\frac{\partial p}{\partial t}-\frac{{\partial }^2}{\partial S^2}\left(\frac{1}{2}{a}^2{b}^2{S}^{2}p\right)-\frac{{\partial }^2}{\partial S\partial y}\left(\rho ab\xi Sp\right)-\frac{{\partial }^2}{\partial y^2}\left(\frac{1}{2}{\xi }^{2}p\right)+\frac{\partial }{\partial S}\left(rSp\right)+\frac{\partial }{\partial y}\left(\mu p\right)=0$

$p(0,S,y)={\delta }_{S_0,y_0}$

Applying Fubini, let $q\left(t,S\right)={\displaystyle {\int }_{ℝ}}p\left(t,S,y\right)\text{d}y$ be the first marginal density of our couple. It is possible to integrate the previous equation and obtain

$\frac{\partial q}{\partial t}-\frac{{\partial }^2}{\partial S^2}\left(\frac{1}{2}{a}^2{S}^2\left({\displaystyle {\int }_{ℝ}}{b}^2\left(y\right)p\left(t,S,y\right)\text{d}y\right)\right)+\frac{\partial }{\partial S}\left(rSq\right)=0$

$q\left(0,S\right)={\delta }_{S_0}$

In the case of a local volatility model ( $b=1$ and $a={\sigma }_D$ ), the density $q_D$ of the Spot process solves the equation

$\frac{\partial q_D}{\partial t}-\frac{{\partial }^2}{\partial S^2}\left(\frac{1}{2}{\sigma }_D^2{S}^2{q}_D\right)+\frac{\partial }{\partial S}\left(rS{q}_D\right)=0$

$q_D\left(0,S\right)={\delta }_{S_0}$

The vanillas of the LSV model being perfectly fitted, we have $q=q_D$ . Iden- tifying the terms in the two last formulas gives

$a^2\left(t,S\right)={\sigma }_D^2\left(t,S\right)\frac{q}{{\displaystyle {\int }_{ℝ}}{b}^{2}p\text{d}y}={\sigma }_D^2\left(t,S\right)\frac{{\displaystyle {\int }_{ℝ}}p\left(t,S,y\right)\text{d}y}{{\displaystyle {\int }_{ℝ}}{b}^2\left(y\right)p\left(t,S,y\right)\text{d}y}$

Using this proposition, and reintroducing the value of a in Kolmogorov forward equation, the joint density $p\left(t\mathrm{,}S\mathrm{,}y\right)$ is then solution of the nonlinear partial integro-differential equation

$\begin{array}{l}\frac{\partial p}{\partial t}-\frac{{\partial }^2}{\partial S^2}\left(\frac{1}{2}{\sigma }_D^2{b}^2{S}^2\frac{{\displaystyle {\int }_{ℝ}}p\text{d}y}{{\displaystyle {\int }_{ℝ}}{b}^{2}p\text{d}y}p\right)-\frac{{\partial }^2}{\partial S\partial y}\left(\rho {\sigma }_{D}b\xi S{\left(\frac{{\displaystyle {\int }_{ℝ}}p\text{d}y}{{\displaystyle {\int }_{ℝ}}{b}^{2}p\text{d}y}\right)}^{\frac{1}{2}}p\right)\\ -\frac{{\partial }^2}{\partial y^2}\left(\frac{1}{2}{\xi }^{2}p\right)+\frac{\partial }{\partial S}\left(rSp\right)+\frac{\partial }{\partial y}\left(\mu p\right)=0\end{array}$ (4)

$p\left(0,S,y\right)={\delta }_{S_0,y_0}$ (5)

There is thus equivalence between the existence of a model (1-2) that calibra- tes the vanillas C and the existence of a solution p to the pide (4).

Remark 2.1. The quotient $\frac{{\displaystyle {\int }_{ℝ}}{b}^{2}p\text{d}y}{{\displaystyle {\int }_{ℝ}}p\text{d}y}$ is nothing but the conditional expecta-

tion of the volatility squared, knowing the spot process. This result is not original in itself (by applying the theorem from [6] for instance), the partial differential equation method however is unusual, and will be used on the other models as well.

The theoretical study of Equations (4) and (5) can be found in [7] . Existence of solutions is proved under strong assumptions (especially on b, which must be sufficiently close to a constant). The general resolution remains an open pro- blem.

2.2. Numerical Results

It seems to be well-known among practitioners, that instabilities occur in their calibration when the volatility’s volatility (in the notations, function $\xi$ ) is too large. This seems to confirm the theoretical limitations met trying to prove the global existence of a solution: when the function b oscillates too much (a change of scale in the factor $y_t$ clearly shows the equivalence between a b that moves a lot and a large $\xi$ ), the resolution of the equation is not guaranteed anymore. To assess these statements, we considered our problem from a practical viewpoint.

In this section, the calibration that stems from solving the partial differential Equation (4) is studied, for two stochastic volatility models: a lognormal one and a “Cox-Ingersoll-Ross” process. The details of the algorithm used for the resolu- tion and a study of the instabilities will be treated later, in Sections 5 and 6.

2.2.1. Lognormal Volatility

Starting with a simple mean reverting model for the volatility factor, the func- tion b is chosen as an exponential

$\frac{\text{d}{S}_t}{S_t}=r\left(t\right)\text{d}t+a\left(t,S_t\right)\exp \left(y_t\right)\text{d}{W}_t^1$ (6)

$\text{d}{y}_t=\kappa \left(\delta -y_t\right)\text{d}t+\gamma \text{d}{W}_t^2$ (7)

with

$a^2\left(t,S\right)={\sigma }_D^2\left(t,S\right)\frac{{\displaystyle {\int }_{ℝ}}p\left(t,S,y\right)\text{d}y}{{\displaystyle {\int }_{ℝ}}\exp \left(2y\right)p\left(t,S,y\right)\text{d}y}$

Equation (4) is solved using the functions we just chose and the local volatility ${\sigma }_D$ associated to the EuroSTOXX 50 implied volatility surface of 2009/04/02. Once function p, density of the couple $\left(S_t\mathrm{,}{y}_t\right)$ , is found, we compute the vanilla prices for different strikes and maturities using this density. To have a point of comparison, the same prices are also computed with the local volatility ${\sigma }_D$ , both of them are then compared to the targeted prices (column TP).

Let us then plot the error between the original vanillas and the ones obtained with the model. The calibration is quite efficient, the errors are equivalent to the ones of the local volatility model.

2.2.2. Cox-Ingersoll-Ross Process

We also focus on the calibration of a model inspired from the interest rates framework: the volatility is assumed to follow a CIR process.

$\frac{\text{d}{S}_t}{S_t}=r\left(t\right)\text{d}t+a\left(t,S_t\right)y_t\text{d}{W}_t^1$ (8)

$\text{d}{y}_t=\kappa \left(\alpha -y_t\right)\text{d}t+\gamma \sqrt{y_t}\text{d}{W}_t^2$ (9)

Detailled properties of this process are described by [8] . In particular, as long as $2\kappa \alpha >{\gamma }^2$ , $y_t$ is strictly positive a.e. Once again, Equation (4) is solved with this stochastic volatility.

3. Application to the “Local Correlation” Model

In this chapter, we are interested in the calibration of a market with n stocks and a basket on those stocks. The purpose is to define a diffusion model for those underlyings that is able to reproduce their implied volatility surface as well as the one of the basket.

The notations are the following, let ${\left(S_t^i\right)}_{1\le i\le n}$ denote the n stocks involved in our problem. The basket's value is given by

$B_t={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{w}_i{S}_t^i$ (10)

where the set ${\left(w_i\right)}_{1\le i\le n}$ stands for the weights of the different underlyings. They are assumed to be constant in the rest of our work. Let us also fix $n+1$ surfaces of vanillas ${\left(C_i\left(T\mathrm{,}K\right)\right)}_{1\le i\le n}$ and $C_B\left(T\mathrm{,}K\right)$ .

3.1. Inconsistencies between Stock and Basket Options

The naive approach to solve this problem is simply to consider n local volatility models

$\frac{\text{d}{S}_t^i}{S_t^i}=r\left(t\right)\text{d}t+{\sigma }_i\left(t,S_t^i\right)\text{d}{W}_t^i$ (11)

The functions ${\sigma }_i$ are easily determined to fit the surfaces ${\left(C_i\left(T\mathrm{,}K\right)\right)}_{1\le i\le n}$ with this diffusion. The correlation matrix $\rho ={\left({\rho }_{ij}\right)}_{1\le i,j\le n}$ associated to the standard brownian motions $W_t^i$ of each underlying can be estimated with historical data. The model is now entirely defined. By Equation (10) of $B_t$ , the vanilla prices for the basket are completely determined and are equal to

$\mathbb{E}\left[{\left({\displaystyle \sum }{w}_i{S}_T^i-K\right)}^{+}\right]$ . However, there is no particular reason for the surface

computed in this framework to be equal to $C_B\left(T\mathrm{,}K\right)$ . In fact, the skew of the basket is more pronounced on the market than in a model with constant correlations between the underlyings, [9] .

3.2. “Local Correlation” Model

In the manner of B. Dupire who decided to let the volatility depend on the level of the spot process, a degree of freedom is added to our model by distorting the matrix of correlation with a function of $B_t$ . This method appears for instance in [5] [10] where the basket level induces some feedback on the values of the different underlyings. In our context, the new correlation matrix is taken as a linear combination of $\rho$ and the constant matrix with only 1 as coefficients. This matrix $\tilde{\rho }$ is equal to

${\tilde{\rho }}_{ij}=\lambda +\left(1-\lambda \right){\rho }_{ij}={\rho }_{ij}+\lambda \left(1-{\rho }_{ij}\right)$ (12)

We shall see while writing the calibration equation that $\lambda$ has to be chosen as a function of the time and of $B_t$ . The matrix ${\tilde{\rho }}_{ij}$ can be seen as an analo- gous of Dupire's local volatility, a 'Local Correlation' so to speak. Assuming that $\rho$ is a proper correlation matrix (definite, positive...), and that the coefficients of the diffusion have a sufficient regularity, our diffusion model then possesses a density in the more general case of a matrix $\rho$ function of the couple $\mathrm{<mo\; stretchy="false">\; (\; </mo>}t\mathrm{,}B\mathrm{<mo\; stretchy="false">\; )\; </mo>}$ . It is also possible to write a condition for the vanillas of the model to be fitted.

Proposition 2. The diffusion model defined by (11) with a correlation function of the couple $\left(t\mathrm{,}B\right)$ , has a density $p\left(t\mathrm{,}{S}^1\mathrm{,}\cdots \mathrm{,}{S}^n\right)$ with respect to Lebesgue's measure. Furthermore, this model calibrates the surface $C_B\left(T\mathrm{,}K\right)$ of $B_t$ 's vanillas (represented by its local volatility surface ${\sigma }_B$ ) if and only if

$\begin{array}{l}{\sigma }_B^2\left(t,B\right)B^2{\displaystyle \int }\tilde{p}\left(t,B,S^2,\cdots ,S^n\right)\text{d}{S}^2\cdots \text{d}{S}^n\\ ={\displaystyle \underset{1\le i,j\le n}{\sum }}{w}_i{w}_j{\rho }_{ij}\left(t,B\right){\displaystyle \int }\left({\tilde{\beta }}_i{\tilde{\beta }}_j\tilde{p}\right)\left(t,B,S^2,\cdots ,S^n\right)\text{d}{S}^2\cdots \text{d}{S}^n\end{array}$ (13)

where

$\tilde{p}\left(t,B,S^2,\cdots ,S^n\right)=p\left(t,S^1,S^2,\cdots ,S^n\right),{\tilde{\beta }}_i\left(t,B,S^2,\cdots ,S^n\right)=S^i{\sigma }_i\left(t,S^i\right)$

with

$S^1=\frac{1}{w_1}\left(B-{\displaystyle \underset{i=2}{\overset{n}{\sum }}}{w}_i{S}^i\right)$ (14)

Proof. The existence of the transition density $p\left(t\mathrm{,}{S}^1\mathrm{,}\cdots \mathrm{,}{S}^n\right)$ stems from the assumptions made on the regularity of the coefficients, and on $\rho$ , for more details see [5] . We can now write the calibration problem for the vanillas of the basket $B_t$ . The density just defined satisfies Kolmogorov forward equation

$\frac{\partial p}{\partial t}-\frac{1}{2}{\displaystyle \underset{1\le i,j\le n}{\sum }}\frac{{\partial }^2}{\partial S^i\partial S^j}\left({\rho }_{ij}{S}_i{\sigma }_i{S}_j{\sigma }_{j}p\right)+{\displaystyle \underset{1\le i\le n}{\sum }}\frac{\partial }{\partial S^i}\left(r{S}^{i}p\right)+rp=0$

To ease the problem, it is useful to change the coordinates $\left(S^1\mathrm{,}\cdots \mathrm{,}{S}^n\right)$ into $\left(B\mathrm{,}{S}^2\mathrm{,}\cdots \mathrm{,}{S}^n\right)$ with $S^1$ defined by (14). After computations, the equation be- comes

$\begin{array}{l}\frac{\partial \tilde{p}}{\partial t}-\frac{1}{2}{\displaystyle \underset{1\le i,j\le n}{\sum }}{w}_i{w}_j\frac{{\partial }^2}{\partial B^2}\left({\rho }_{ij}{\tilde{\beta }}_i{\tilde{\beta }}_j\tilde{p}\right)\\ -\frac{1}{2}{\displaystyle \underset{2\le i,j\le n}{\sum }}\left(\frac{{\partial }^2}{\partial S^i\partial S^j}+w_i\frac{{\partial }^2}{\partial B\partial S^j}+w_j\frac{{\partial }^2}{\partial S^i\partial B}\right)\left({\rho }_{ij}{\tilde{\beta }}_i{\tilde{\beta }}_j\tilde{p}\right)\\ -{\displaystyle \underset{2\le i\le n}{\sum }}{w}_1\frac{{\partial }^2}{\partial B\partial S^j}\left({\rho }_{1i}{\tilde{\beta }}_1{\tilde{\beta }}_i\tilde{p}\right)+\frac{\partial }{\partial B}\left(rB\tilde{p}\right)+{\displaystyle \underset{2\le i\le n}{\sum }}\frac{\partial }{\partial S^i}\left(r{S}^i\tilde{p}\right)+r\tilde{p}=0\end{array}$

where $\tilde{p}\left(B\mathrm{,}{S}^2\mathrm{,}\cdots \mathrm{,}{S}^n\right)$ and ${\tilde{\beta }}_i\left(B\mathrm{,}{S}^2\mathrm{,}\cdots \mathrm{,}{S}^n\right)$ are defined above. Integrating the equation against the variables $\left(S^2\mathrm{,}\cdots \mathrm{,}{S}^n\right)$ , and writing $q=\frac{1}{w_1}{\displaystyle \int }\tilde{p}\text{d}{S}^2\cdots \text{d}{S}^n$ , the density of the marginal law of B satisfies

$\frac{\partial q}{\partial t}-\frac{1}{2}{\displaystyle \underset{1\le i,j\le n}{\sum }}{w}_i{w}_j\frac{{\partial }^2}{\partial B^2}\left({\rho }_{ij}{\displaystyle \int }{\tilde{\beta }}_i{\tilde{\beta }}_j\tilde{p}\text{d}{S}^2\cdots \text{d}{S}^n\right)+\frac{\partial }{\partial B}\left(rBq\right)+rq=0$

Comparing this equation to Dupire’s equation for the local volatility ${\sigma }_B$

$\frac{\partial q}{\partial t}-\frac{1}{2}\frac{{\partial }^2}{\partial B^2}\left({\sigma }_D^2{B}^{2}q\right)+\frac{\partial }{\partial B}\left(rBq\right)+rq=0$

if our model reproduces the vanillas ${\sigma }_D$ , then the following equality must be verified

${\sigma }_D^2{B}^2{\displaystyle \int }\tilde{p}\text{d}{S}^2\cdots \text{d}{S}^n={\displaystyle \underset{1\le i,j\le n}{\sum }}{w}_i{w}_j{\rho }_{ij}{\displaystyle \int }{\tilde{\beta }}_i{\tilde{\beta }}_j\tilde{p}\text{d}{S}^2\cdots \text{d}{S}^n$

Reciprocally, the condition we just wrote is clearly sufficient for the options to be calibrated

Remark 3.1. Let us note that this condition is written as an equality between two functions of the time and of B. The other variables are no longer repre- sented.

Assuming that condition (13) is not verified, the model defined by (11) does not fit the vanillas of the basket $B_t$ . It has to be enriched to solve the calibration problem. Our choice is to distort the correlation matrix. The new matrix $\tilde{\rho }$ is described by (12). Hence, let $\Theta$ denote the matrix ${\Theta }_{ij}=1$ for all $1\le i\mathrm{,}j\le n$ . We also notice that, the trace of $\rho$ being equal to n, its smallest eigenvalue $K_{\rho }$ is smaller than 1.

Lemma 1. The matrix $\tilde{\rho }=\left(1-\lambda \right)\rho +\lambda \Theta$ is also a correlation matrix as long as $\lambda$ is in $\left]r\mathrm{,1}\right[$ with

$r=-\min \left(\underset{i\ne j}{\max }\frac{1+{\rho }_{ij}}{1-{\rho }_{ij}},\frac{K_{\rho }}{n-K_{\rho }}\right)<0$ (15)

Proof. Clearly, for all $\left({\xi }_i\mathrm{,}{\xi }_j\right)\in {ℝ}^2$

${\displaystyle \underset{1\le i,j\le n}{\sum }}{\xi }_i\left(\left(1-\lambda \right){\rho }_{ij}+\lambda \right){\xi }_j\ge K_{\rho }\left(1-\lambda \right){\left|\xi \right|}^2+\lambda {\left({\displaystyle \underset{1\le i\le n}{\sum }}{\xi }_i\right)}^2$

If $\lambda$ is positive, since $K_{\rho }\left(1-\lambda \right)$ is stricty positive, the matrix remains definite positive. Now, if $\lambda <0$ , Cauchy-Schwarz gives

${\displaystyle \underset{1\le i,j\le n}{\sum }}{\xi }_i\left(\left(1-\lambda \right){\rho }_{ij}+\lambda \right){\xi }_j\ge \left(K_{\rho }\left(1-\lambda \right)+\lambda n\right){\left|\xi \right|}^2$

Since $K_{\rho }\le 1$ , $\lambda >-\frac{K_{\rho }}{n-K_{\rho }}$ is enough for $K_{\rho }\left(1-\lambda \right)+\lambda n$ to be stricly po- sitive.

The diagonal coefficients of $\tilde{\rho }$ are still 1. As for the other terms, thanks to the first term in relation (15), they still belong to the interval $\left]-\mathrm{1,1}\right[$ .

We introduce the new correlation matrix in condition (13), this gives

$\lambda \left(\tilde{p}\right)=\frac{{\sigma }_D^2{B}^2{\displaystyle \int }\tilde{p}\text{d}{S}^2\cdots \text{d}{S}^n-{\displaystyle \underset{1\le i,j\le n}{\sum }}{w}_i{w}_j{\rho }_{ij}{\displaystyle \int }{\tilde{\beta }}_i{\tilde{\beta }}_j\tilde{p}\text{d}{S}^2\cdots \text{d}{S}^n}{{\displaystyle \underset{1\le i,j\le n}{\sum }}{w}_i{w}_j\left(1-{\rho }_{ij}\right){\displaystyle \int }{\tilde{\beta }}_i{\tilde{\beta }}_j\tilde{p}\text{d}{S}^2\cdots \text{d}{S}^n}$ (16)

$\lambda \left(\tilde{p}\right)$ is a function of B and t. It is now possible to use this value to write a pide on the density $\tilde{p}$ . Any solution of the following equation is a density that calibrates the vanillas of the basket

$\frac{\partial \tilde{p}}{\partial t}+L_1\tilde{p}+L_2^{\lambda \left(\tilde{p}\right)}\left(\tilde{p}\right)=0$ (17)

where $L_1$ is linear and verifies

$\begin{array}{c}{L}_1\tilde{p}=-\frac{1}{2}{\displaystyle \underset{1\le i,j\le n}{\sum }}{w}_i{w}_j{\rho }_{ij}\frac{{\partial }^2}{\partial B^2}\left({\tilde{\beta }}_i{\tilde{\beta }}_j\tilde{p}\right)\\ -\frac{1}{2}{\displaystyle \underset{2\le i,j\le n}{\sum }}{\rho }_{ij}\left(\frac{{\partial }^2}{\partial S^i\partial S^j}+w_i\frac{{\partial }^2}{\partial B\partial S^j}+w_j\frac{{\partial }^2}{\partial S^i\partial B}\right)\left({\tilde{\beta }}_i{\tilde{\beta }}_j\tilde{p}\right)\\ -{\displaystyle \underset{2\le i\le n}{\sum }}{w}_1{\rho }_{1i}\frac{{\partial }^2}{\partial B\partial S^j}\left({\tilde{\beta }}_1{\tilde{\beta }}_i\tilde{p}\right)+\frac{\partial }{\partial B}\left(rB\tilde{p}\right)+{\displaystyle \underset{2\le i\le n}{\sum }}\frac{\partial }{\partial S^i}\left(r{S}^i\tilde{p}\right)+r\tilde{p}\end{array}$

and $L_2^{\lambda \left(\tilde{p}\right)}$ is the nonlinear part of the equation

$\begin{array}{c}{L}_2^{\lambda \left(\tilde{p}\right)}\left(\tilde{p}\right)=-\frac{1}{2}{\displaystyle \underset{1\le i,j\le n}{\sum }}{w}_i{w}_j\left(1-{\rho }_{ij}\right)\frac{{\partial }^2}{\partial B^2}\left(\lambda \left(\tilde{p}\right){\tilde{\beta }}_i{\tilde{\beta }}_j\tilde{p}\right)\\ -\frac{1}{2}{\displaystyle \underset{2\le i,j\le n}{\sum }}\left(1-{\rho }_{ij}\right)\left(\frac{{\partial }^2}{\partial S^i\partial S^j}+w_i\frac{{\partial }^2}{\partial B\partial S^j}+w_j\frac{{\partial }^2}{\partial S^i\partial B}\right)\left(\lambda \left(\tilde{p}\right){\tilde{\beta }}_i{\tilde{\beta }}_j\tilde{p}\right)\\ -{\displaystyle \underset{2\le i\le n}{\sum }}{w}_1\left(1-{\rho }_{1i}\right)\frac{{\partial }^2}{\partial B\partial S^j}\left(\lambda \left(\tilde{p}\right){\tilde{\beta }}_1{\tilde{\beta }}_i\tilde{p}\right)\end{array}$

Remark 3.2. The operator $L_1+L_2$ stems from a change of coordinates on a uniformly elliptic operator. It is also elliptic, uniformly on any domain where the ${\tilde{\beta }}_i$ are bounded away from 0 by a strictly positive constant.

Furthermore, the initial condition is

$\tilde{p}\left(0,B,S^2,\cdots ,S^n\right)=\delta \left({\displaystyle \sum }{w}_i{S}_0^i,S_0^2\cdots ,S_0^n\right)$

where $S_0^i$ is the market value at instant 0 of the i-th stock. Applying this initial condition to (16), the initial value of $\lambda$ is equal to

$\lambda \left(\tilde{p}\right)\left(0,B\right)=\frac{{\sigma }_D^2\left(0,B\right)B^2-{\displaystyle \underset{1\le i,j\le n}{\sum }}{w}_i{w}_j{\rho }_{ij}{S}_0^i{\sigma }_i\left(0,S_0^i\right)S_0^j{\sigma }_j\left(0,S_0^j\right)}{{\displaystyle \underset{1\le i,j\le n}{\sum }}{w}_i{w}_j\left(1-{\rho }_{ij}\right)S_0^i{\sigma }_i\left(0,S_0^i\right)S_0^j{\sigma }_j\left(0,S_0^j\right)}$ (18)

For a theoretical study of the calibration equation, we refer the reader to [5] .

3.3. Resolution of the Equation for the Calibration of a Basket

This subsection focuses on the results of the calibration for a two-underlyings basket. Let us consider two assets, both of them are assumed to generate the following implied volatility surface.

Using a Monte-Carlo simulation and the local volatilities stemming from those surfaces, the theoretical prices for the basket $B_t$ are computed, with weights $w_1=w_2=0.5$ and a correlation ${\rho }_{12}=-0.5$ .

Distorting this theoretical surface by a factor of 0.9, and thus making the prices of the basket inconsistent with the prices of the underlyings, the calibration algorithm is applied. Solving the partial integro-differential Equation (17) gives the following vanillas (quoted in implied volatility).

Here are some results for other tests. We calibrate a model with the following parameters ${\rho }_{12}=0$ , $w_1=0.3$ and $w_2=0.7$ , the targeted surface is the theo- retical one distorted with two factors: first 0.95 and second 1.05.

At last, a different surface for the second underlying is chosen, mutliplying the first one (described in 3.3) by 0.9, the correlation is this time taken as 0.5.

The results are rather satisfactory, especially at the money. To keep the computations to a reasonnable duration, a sparse initial surface was used. This explains why the calibration is not better far from the money, the fitting method is nontheless valid. Now follows an outlook of the values taken by the new correlation $\tilde{\rho }$ at different maturities when the theoretical surface is distorted by factors 0.95 and 1.05. The parameters are: ${\rho }_{12}=0.5$ , $w_1=0.7$ and $w_2=0.3$ .

As expected, the Local Correlation and the distorsion factor evolve in the same direction. The underlyings must be more correlated when the implied volatility of the basket is higher, and reciprocally. Furthermore, it appears that in the case of the 0.95 distorsion, the correlation has to violently decrease for high values of B: the two underlyings must be anti-correlated when they are both large.

As for the influence of the maturity, let us first state that in the computation of $\lambda$ , when the denominator is smaller than ${10}^{-6}$ , we chose not to change the correlation, to avoid numerical errors. It appears that, as long as B is in a zone where $\lambda$ was actually computed, the framework chosen to test the calibration actually generates a local correlation constant in time.

4. Application to Stochastic Interest Rates

This section is dedicated to hybrid local volatility models with stochastic rates. The interest rate is assumed to be stochastic and to follow a diffusion equation. The volatility depends on the level of the spot process exactly as in a local volatility model. The idea is to compute its exact value for the vanillas in this model to be calibrated.

4.1. Calibration of the Hybrid Local Volatility Model

The risk-neutral diffusion of the model is written as

$\frac{\text{d}{S}_t}{S_t}=r\left(t,y_t\right)\text{d}t+\sigma \left(t,S_t\right)\text{d}{W}_t^1$ (19)

$\text{d}{y}_t=\mu \left(t,y_t\right)\text{d}t+\lambda \left(t,y_t\right)\text{d}{W}_t^2$

The two brownian motions are correlated with a constant correlation denoted by $\rho$ . Classic regularity and ellipticity assumptions are made on the coefficients of the diffusion (described in [5] ) to get the

Proposition 3 The diffusion model defined above has a transition density with respect to Lebesgue’s measure. The value of $\sigma$ that fits its vanillas is given by

${\sigma }^2\left(t,S\right)={\sigma }_D^2\left(t,S\right)+2\frac{\bar{r}\left(t\right){\displaystyle {\int }_{ℝ}}{\displaystyle {\int }_S^{+\infty }}p\left(t,s,y\right)\text{d}s\text{d}y-{\displaystyle {\int }_{ℝ}}{\displaystyle {\int }_S^{+\infty }}r\left(t,y\right)p\left(t,s,y\right)dsdy}{S{\displaystyle {\int }_{ℝ}}p\left(t,S,y\right)\text{d}y}$ (20)

where $p\left(t\mathrm{,}S\mathrm{,}y\right)$ is the density of the couple $\left(S_t\mathrm{,}{y}_t\right)$ , and $\bar{r}\left(t\right)$ a determini- stic curve of rates used in the computation of Dupire’s local volatility ${\sigma }_D$ .

Proof. The existence of the density $p\mathrm{<mo\; stretchy="false">\; (\; </mo>}t\mathrm{,}S\mathrm{,}y\mathrm{<mo\; stretchy="false">\; )\; </mo>}$ stems from the assumptions on the coefficients. Let us prove formula (20). The function p solves the forward parabolic equation

$\begin{array}{l}\frac{\partial p}{\partial t}-\frac{{\partial }^2}{\partial S^2}\left(\frac{1}{2}{\sigma }^2{S}^{2}p\right)-\frac{{\partial }^2}{\partial S\partial y}\left(\rho \sigma \lambda Sp\right)-\frac{{\partial }^2}{\partial y^2}\left(\frac{1}{2}{\lambda }^{2}p\right)\\ +\frac{\partial }{\partial S}\left(rSp\right)+\frac{\partial }{\partial y}\left(\mu p\right)+rp=0\end{array}$

with the initial condition $p\left(0,S,y\right)={\delta }_{S_0,y_0}$ . As previously, the equation is inte- grated with respect to y, writing $q\left(t,S\right)={\displaystyle {\int }_{ℝ}}p\left(t,S,y\right)\text{d}y$

$\frac{\partial q}{\partial t}-\frac{{\partial }^2}{\partial S^2}\left(\frac{1}{2}{\sigma }^2{S}^{2}q\right)+\frac{\partial }{\partial S}\left(S{\displaystyle {\int }_{ℝ}}r\left(t,y\right)p\left(t,S,y\right)\text{d}y\right)+{\displaystyle {\int }_{ℝ}}r\left(t,y\right)p\left(t,S,y\right)\text{d}y=0$

This equation needs to be matched with

$\frac{\partial q_D}{\partial t}-\frac{{\partial }^2}{\partial S^2}\left(\frac{1}{2}{\sigma }_D^2{S}^2{q}_D\right)+\frac{\partial }{\partial S}\left(\bar{r}S{q}_D\right)+\bar{r}{q}_D=0$

$q_D\left(0,S\right)={\delta }_{S_0}$

Both of them can be written as

$\frac{1}{2}{\sigma }^2{S}^2=\frac{{\displaystyle {\int }_0^{+\infty }}{\left(s-S\right)}^{+}\left(\frac{\partial q}{\partial t}+\frac{\partial }{\partial s}\left(s{\displaystyle \int }rp\text{d}y\right)+{\displaystyle \int }rp\text{d}y\right)\text{d}s}{q}$

$\frac{1}{2}{\sigma }_D^2{S}^2=\frac{{\displaystyle {\int }_0^{+\infty }}{\left(s-S\right)}^{+}\left(\frac{\partial q_D}{\partial t}+\frac{\partial }{\partial s}\left(s\bar{r}{q}_D\right)+\bar{r}{q}_D\right)\text{d}s}{q_D}$

Computing

$\begin{array}{c}{\displaystyle {\int }_0^{+\infty }}{\left(s-S\right)}^{+}\frac{\partial }{\partial s}\left(s{\displaystyle \int }rp\text{d}y\right)\text{d}s={\displaystyle {\int }_0^{+\infty }}{\left(s-S\right)}^{+}\left(s\frac{\partial }{\partial s}\left({\displaystyle \int }rp\text{d}y\right)+{\displaystyle \int }rp\text{d}y\right)\text{d}s\\ =-{\displaystyle {\int }_0^{+\infty }}\frac{\partial }{\partial s}\left({\left(s-S\right)}^{+}\right)s{\displaystyle \int }rp\text{d}y\text{d}s\\ =-{\displaystyle {\int }_0^{+\infty }}s{1}_{s\ge S}{\displaystyle \int }rp\text{d}y\text{d}s\end{array}$

where the second line stems from a simple integration by parts. Reintroducing this into the previous equations, we get

$\frac{1}{2}{\sigma }^2{S}^2=\frac{{\displaystyle {\int }_0^{+\infty }}{\left(s-S\right)}^{+}\frac{\partial q}{\partial t}\text{d}s-S{\displaystyle {\int }_S^{+\infty }}{\displaystyle \int }rp\text{d}y\text{d}s}{q}$

$\frac{1}{2}{\sigma }_D^2{S}^2=\frac{{\displaystyle {\int }_0^{+\infty }}{\left(s-S\right)}^{+}\frac{\partial q_D}{\partial t}\text{d}s-\bar{r}S{\displaystyle {\int }_S^{+\infty }}{q}_D\text{d}s}{q_D}$

In order to calibrate the vanillas, all that remains to be done is match the marginal density q with $q_D$ , giving the necessary condition

$\frac{1}{2}{\sigma }^2{S}^2+\frac{S{\displaystyle {\int }_S^{+\infty }}{\displaystyle \int }rp\text{d}y\text{d}s}{q}=\frac{1}{2}{\sigma }_D^2{S}^2+\frac{\bar{r}S{\displaystyle {\int }_S^{+\infty }}q\text{d}s}{q}$

Replacing q by ${\displaystyle {\int }_{ℝ}}p\left(t\mathrm{,}S\mathrm{,}y\right)\text{d}y$ completes the proof.

The calibration equation for the vanillas of our hybrid model is thus

$\begin{array}{l}\frac{\partial p}{\partial t}-\frac{{\partial }^2}{\partial S^2}\left(\frac{1}{2}{\sigma }^2{S}^{2}p\right)-\frac{{\partial }^2}{\partial S\partial y}\left(\rho \sigma \lambda Sp\right)-\frac{{\partial }^2}{\partial y^2}\left(\frac{1}{2}{\lambda }^{2}p\right)\\ +\frac{\partial }{\partial S}\left(rSp\right)+\frac{\partial }{\partial y}\left(\mu p\right)+rp=0\end{array}$ (21)

with $\sigma$ given by Formula (20). Using similar technics to the other cases, an existence result can be obtained under certain assumptions, but this is not the scope of this paper. It is however noteworthy that one of the necessary hypothesis is the small variation of function r with respect to the deterministic curve $\bar{r}$ .

4.2. Numerical Calibration

In this section, the theoretical results above are applied to calibrate a given diffusion model. Assuming the instantaneous rate to obey a Vasicek model (or in other words an Ornstein-Uhlenbeck process), the diffusion equations become

$\frac{\text{d}{S}_t}{S_t}=r_t\text{d}t+\sigma \left(t,S_t\right)\text{d}{W}_t^1$

$\text{d}{r}_t=a\left(b-r_t\right)\text{d}t+\gamma \text{d}{W}_t^2$

The ADI algorithm described in section 5 is applied to Equation (21) with the coefficients associated to this diffusion. The initial condition is $p_0(S,r)={\delta }_{S_0,r_0}$ .

As in the two previous sections, this partial integro-differential equation is solved with a variable change for the spot process $x=\ln \left(S\right)$ . The grid chosen is $\left[-10\sigma \sqrt{t},10\sigma \sqrt{t}\right]\times \left[-0.1,0.2\right]$ with $\sigma =0.2$ and $r_0=0.04$ . We discretize it with 300 points in both the spot and the rate direction. The initial condition (Dirac mass at the point $\left(\ln \left(S_0\right),r_0\right)$ ) is approximated by a bivariate Gaussian centred at that point with a very small variance.

The following numerical values are taken for the diffusion

$a=0.5,b=0.7,\gamma =0.01,r_0=0.04$

These values generate the interest rate

To assess the quality of the calibration, call and put options on the spot process are computed by integration on the grid and compared to the targeted prices (target columns). Convergence is quite satisfactory. For instance, for 6 months and 1 year maturity vanillas, one finds

5. Algorithm for the Resolution of the Calibration Equation

In the previous sections, we wrote the calibration equations and gave graphs for their efficiency. In this one, we describe the algorithm used to solve them: a classic alternating direction implicit scheme. The nonlinear term is handled using a forward induction at first, and then a predictor-corrector method.

The strong feature of an ADI scheme is its convergence rate in time and space: $O\left(\delta x^2\right)+O\left(\delta t^2\right)$ . The nonlinearity of the equation challenges this assertion. It is however possible to prove that it remains true in this case too.

5.1. Alternating Direction Implicit Scheme

The calibration equation is a second-order parabolic equation. One of the most efficient method to solve such equations is a finite-difference approximation with alternating direction methods. For more informations on the subject, the reader can look into [4] , numerous articles have also been published, in particu- lar [11] [12] [13] [14] . Let us now consider the following equation

$\begin{array}{l}\frac{\partial p}{\partial t}-\frac{{\partial }^2}{\partial x^2}\left(\frac{1}{2}{f}^2{I}^2\left(p\right)p\right)-\frac{{\partial }^2}{\partial x\partial y}\left(\rho fgI\left(p\right)p\right)\\ -\frac{{\partial }^2}{\partial y^2}\left(\frac{1}{2}{g}^{2}p\right)+\frac{\partial }{\partial x}\left(\alpha p\right)+\frac{\partial }{\partial y}\left(\beta p\right)+\gamma p=0\end{array}$ (22)

$p\left(0,S,y\right)={\delta }_{S_0,y_0}$ (23)

where f, g, $\alpha$ , $\beta$ and $\gamma$ are functions of t, x and y, $\rho$ a constant and $I\left(p\right)$ the quotient of integrals

$I^2\left(p\right)\left(t,x\right)=\frac{{\displaystyle \int }{h}^2\left(y\right)p\left(t,x,y\right)\text{d}y}{{\displaystyle \int }p\left(t,x,y\right)\text{d}y}$ (24)

Remark 5.1 We restricted ourselves to two-dimensional equations since they cover all the concrete examples studied previously. But the computations that follow are true in the general case. The cost in time however becomes an issue in higher dimensions.

The domain for the numerical resolution is $\left]\mathrm{0,}T\right[\times \left]x_{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{x}^{\mathrm{<mo>\; *\; </mo>}}\right[\times \left]y_{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}^{\mathrm{<mo>\; *\; </mo>}}\right[$ . The first step is to take care of the initial condition. Instead of the Dirac, the initial condition $p_0$ is chosen as a gaussian distribution with very small variance. It obviously approximates our initial condition. It also verifies (on any bounded domain) the properties of regularity and strict positivity required in the theoretical study of the Equation (though in the present section, we are only interested in its numerical resolution). Let $p\left(t,x,y\right)=p_0\left(x,y\right)$ if $x=x_{*},x^{*}$ or $y=y_{*},y^{*}$ be the boundary conditions.

The algorithm is based upon a predictor-corrector approach. Let $\Delta x$ , $\Delta y$ and $\Delta t$ be increments of the variables x, y and t, where $\Delta x=\frac{x^{*}-x_{*}}{I}$ , $\Delta y=\frac{y^{*}-y_{*}}{J}$ and $\Delta t=\frac{T}{N}$ with I, J and N integers. The sets of points in the x,y,t-plane is given by $x_n=x_{*}+i\frac{\Delta x}{x^{*}-x_{*}}$ , $y_j=y_{*}+j\frac{\Delta y}{y^{*}-y_{*}}$ and $t_n=n\frac{\Delta t}{T}$ , for $0\le i\le I$ , $0\le j\le J$ and $0\le n\le N$ . We construct four sequences $p_n$ , $p_n^{\mathrm{<mo>\; *\; </mo>}}$ , $q_n$ and $q_n^{\mathrm{<mo>\; *\; </mo>}}$ of space-dependent functions with n between 0 and N.

A classic alternating direction implicit scheme functions as follows: let us define the initial functions

$p_0^{*}\left(i,j\right)=q_0\left(i,j\right)=q_0^{*}\left(i,j\right)=p_0\left(x_i,y_j\right)$

and then by induction

$\begin{array}{c}\frac{q_{n+1}^{*}-p_n}{\Delta t}={\delta }_x^2\left(I^2\left(p_n\right)\frac{f_{n+1}^2{q}_{n+1}^{*}+f_n^2{p}_n}{4}\right)\\ +{\delta }_{xy}^2\left(\rho f_n{g}_{n}I\left(p_n\right)p_n\right)+{\delta }_y^2\left(\frac{1}{2}{g}_n^2{p}_n\right)\\ -{\delta }_x\left(\frac{{\alpha }_{n+1}{q}_{n+1}^{*}+{\alpha }_n{p}_n}{2}\right)-{\delta }_y\left({\beta }_n{p}_n\right)-\frac{{\gamma }_{n+1}{q}_{n+1}^{*}+{\gamma }_n{p}_n}{2}\end{array}$ (25)

$\begin{array}{c}\frac{q_{n+1}-p_n}{\Delta t}={\delta }_x^2\left(I^2\left(p_n\right)\frac{f_{n+1}^2{q}_{n+1}^{*}+f_n^2{p}_n}{4}\right)+{\delta }_{xy}^2\left(\rho f_n{g}_{n}I\left(p_n\right)p_n\right)\\ +{\delta }_y^2\left(\frac{g_{n+1}^2{q}_{n+1}+g_n^2{p}_n}{4}\right)-{\delta }_x\left(\frac{{\alpha }_{n+1}{q}_{n+1}^{*}+{\alpha }_n{p}_n}{2}\right)\\ -{\delta }_y\left(\frac{{\beta }_{n+1}{q}_{n+1}+{\beta }_n{p}_n}{2}\right)-\frac{{\gamma }_{n+1}{q}_{n+1}+{\gamma }_n{p}_n}{2}\end{array}$ (26)

where $f_n\left(i\mathrm{,}j\right)$ designates $f\left(n\Delta t\mathrm{,}{x}_i\mathrm{,}{y}_j\right)$ (the same thing being true for the other coefficients of the equation). $\delta$ is a difference operator for the space derivatives. For instance, with $1\le i\le I-1$ and $1\le j\le J-1$

$\begin{array}{l}{\delta }_x^2{f}_n=\frac{f_n\left(i-1,j\right)-2f_n\left(i,j\right)+f_n\left(i+1,j\right)}{\Delta x^2}\\ {\delta }_{xy}^2{f}_n=\frac{f_n\left(i-1,j-1\right)-f_n\left(i+1,j\right)-f_n\left(i,j+1\right)+f_n\left(i+1,j+1\right)}{4\Delta x\Delta y}\\ {\delta }_x{f}_n=\frac{f_n\left(i+1,j\right)-f_n\left(i-1,j\right)}{2\Delta x}\\ ⋮\end{array}$

The Equations (25) and (26) form two tridiagonal systems that can be solved very efficiently. A recursion formula can be computed on the functions $q_n$ defined above

$\frac{q_{n+1}-q_{n+1}^{*}}{\Delta t}={\delta }_y^2\left(\frac{g_{n+1}^2{q}_{n+1}-g_n^2{p}_n}{4}\right)-{\delta }_y\left(\frac{{\beta }_{n+1}{q}_{n+1}-{\beta }_n{p}_n}{2}\right)-{\gamma }_{n+1}\frac{q_{n+1}-q_{n+1}^{*}}{2}$

Thus

$q_{n+1}^{*}=q_{n+1}-\frac{2\Delta t}{2+{\gamma }_{n+1}\Delta t}\left({\delta }_y^2\left(\frac{g_{n+1}^2{q}_{n+1}-g_n^2{p}_n}{4}\right)-{\delta }_y\left(\frac{{\beta }_{n+1}{q}_{n+1}-{\beta }_n{p}_n}{2}\right)\right)$

And eventually

$\begin{array}{l}\frac{q_{n+1}-p_n}{\Delta t}={\delta }_x^2\left(I^2\left(p_n\right)\frac{f_{n+1}^2{q}_{n+1}+f_n^2{p}_n}{4}\right)+{\delta }_{xy}^2\left(\rho f_n{g}_{n}I\left(p_n\right)p_n\right)\\ +{\delta }_y^2\left(\frac{g_{n+1}^2{q}_{n+1}+g_n^2{p}_n}{4}\right)-{\delta }_x\left(\frac{{\alpha }_{n+1}{q}_{n+1}+{\alpha }_n{p}_n}{2}\right)\\ -{\delta }_y\left(\frac{{\beta }_{n+1}{q}_{n+1}+{\beta }_n{p}_n}{2}\right)-\frac{{\gamma }_{n+1}{q}_{n+1}+{\gamma }_n{p}_n}{2}\\ -\frac{\Delta t}{4+2{\gamma }_{n+1}\Delta t}{\delta }_x^2\left(I^2\left(p_n\right)f_{n+1}^2\left({\delta }_y^2\left(\frac{g_{n+1}^2{q}_{n+1}-g_n^2{p}_n}{4}\right)-{\delta }_y\left(\frac{{\beta }_{n+1}{q}_{n+1}-{\beta }_n{p}_n}{2}\right)\right)\right)\\ +\frac{\Delta t}{2+{\gamma }_{n+1}\Delta t}{\delta }_x\left({\alpha }_{n+1}\left({\delta }_y^2\left(\frac{g_{n+1}^2{q}_{n+1}-g_n^2{p}_n}{4}\right)-{\delta }_y\left(\frac{{\beta }_{n+1}{q}_{n+1}-{\beta }_n{p}_n}{2}\right)\right)\right)\end{array}$ (27)

In the litterature concerning alternating direction implicit schemes, the functions $q_n$ are often called the predicted value of the solution p, a second “corrector” step usually follows. In our case, the equation being nonlinear, we start by chosing $p_{n+1}=q_{n+1}$ and study the finite-difference approximation that results.

Proposition 4. This last finite-difference equation is consistent with the partial differential Equation (22) on a bounded domain with smooth initial condition $p_0$ , the truncation error is

$O\left(\Delta t\right)+O\left(\Delta x^2\right)+O\left(\Delta y^2\right)+O\left(\frac{\Delta x^3}{\Delta y}\right)$

Proof. Let p be a classic solution of Equation (22) on a bounded domain, with ${p|}_{t=0}=p_0$ . We assume that p is strictly positive and sufficiently differentiable for all the quantities in the sequel to be properly defined. All the derivatives that appear are bounded as a consequence of the regularity assumptions on p. Writing $p_{i,j}^n=p^n\left(x_i,y_j\right)=p\left(t_n,x_i,y_j\right)$ , a simple Taylor expansion with re- mainder gives

$\begin{array}{c}\frac{p_{i,j}^{n+1}-p_{i,j}^n}{\Delta t}=\frac{\frac{\partial p}{\partial t}\left(t_n,x_i,y_j\right)+\frac{\partial p}{\partial t}\left(t_{n+1},x_i,y_j\right)}{2}\\ +\frac{\Delta t}{2}\left(\frac{{\partial }^{2}p}{\partial t^2}\left(t_n+\theta ,x_i,y_j\right)-\frac{{\partial }^{2}p}{\partial t^2}\left(t_n+{\theta }^{*},x_i,y_j\right)\right)\end{array}$

As for the space derivatives, clearly

${\delta }_x^2{p}^n=\frac{{\partial }^2{p}^n}{\partial x^2}\left(x_i,y_j\right)+\frac{\Delta x^2}{24}\left(\frac{{\partial }^4{p}^n}{\partial x^4}\left(x_i+{\theta }_1,y_j\right)+\frac{{\partial }^4{p}^n}{\partial x^4}\left(x_i-{\theta }_2,y_j\right)\right)$

$\begin{array}{c}{\delta }_{xy}^2{p}^n=\frac{{\partial }^2{p}^n}{\partial x\partial y}\left(x_i,y_j\right)+\frac{\Delta y^2}{12}\left(\frac{{\partial }^4{p}^n}{\partial x\partial y^3}\left(x_i,y_j+{\theta }_1\right)+\frac{{\partial }^4{p}^n}{\partial x\partial y^3}\left(x_i,y_j-{\theta }_2\right)\right)\\ +\frac{\Delta x^2}{24\Delta y}[\frac{{\partial }^3{p}^n}{\partial x^3}\left(x_i+{\theta }_3,y_{j-1}\right)-\frac{{\partial }^3{p}^n}{\partial x^3}\left(x_i+{\theta }_4,y_{j-1}\right)\\ +\frac{{\partial }^3{p}^n}{\partial x^3}\left(x_i-{\theta }_5,y_{j+1}\right)-\frac{{\partial }^3{p}^n}{\partial x^3}\left(x_i-{\theta }_6,y_{j+1}\right)]\end{array}$