I'm trying to calibrate Local stochastic volatility model using finite difference method, and I'm mainly following this referece: Tian (2015).
I met a problem when calibrating leverage function - the discritized fokker planck equation with ADI Douglas scheme will generate negative density, which leads to negative value in the integrals of leverage function and this will then break the square root calculation. However, there seems not much discussion in the references about this issue.
So my question is:
- Is ADI supposed to preserve the positivity of the density function? (If so, this indicates my current implementation is subject to some error.)
- If not, what should I do to avoid this?
-----------------------------------------
Updates (0401):
I found that applying zero-flux boundary condition (described in QuantLib Document) helps a lot to stabilize the calibration process - no more negative integrals in square root function. The Monte Carlo roundtrip of five-year option has an error of 0.5 (compared to BS price) with ATM strike.
However, It is only stable when the time step is small (dt=1/300). If I try to use a non-uniform temporal grid which has larger time step for longer maturities, the density will lose its structure and calibration then fails. Another concern is I do see the density getting smaller as the calibration proceed in time and currently I'm fixing this by rescaling the density to one in each step.
So my current question is:
I implemented the zero-flux boundary condition only on boundary of v and applied a simpler boundary condition over S - assume diminishing second order term in PDE. I'm wondering if this is the correct way to do it. Should I apply zero flux on boundary of S as well? if not, what should I use?
How should I deal with the loss of density in calibration?
How can I improve the calibration stability for large time step?
Thanks a lot!
-----------------------------------------
Updates (0405):
The calibration process is still not stable. Truncation at zero only helps for small time step. I think there should be something to do with the way I implement ADI scheme. Here is the details about how I deal with Fokker Planck in ADI.
I start from the Fokker Planck equation as follow: \begin{align*} \partial_{t}f = &\frac{1}{2}v\partial_{x}^2L^2f \\ &+ (-r+q)\partial_{x}f + \frac{v}{2}\partial_{x}L^2f\\ &+ \frac{\sigma^2}{2}\frac{1}{v}\partial_{z}^2f \\ &+ ((-\kappa\theta-\frac{\sigma^2}{2})\frac{1}{v} + \kappa)\partial_{z}f + \kappa\theta\frac{1}{v}f\\ &+\rho\sigma\partial_{x}\partial_{z}Lf \end{align*} where $x=log(S/S_0)$ and $z=log(v/v0)$.
Then I rewrite this into matrix form: $$\partial_{t}f = F_0+F_1+F_2 $$ $$F_0 = \sigma*\rho*dx@f@dz$$
\begin{align*} F_1=f @ &[(-\kappa*\theta-0.5*\sigma^2)*dz@Diag(1/v\_vec)\\ &+ \kappa*dz\\ &+ 0.5*\sigma^2*dz2@Diag(1/v\_vec)\\ &+\kappa*\theta*Diag(1/v\_vec)] \end{align*}
for $ j = 1,2,...,Nz:$ \begin{align*} F_2[:,j]=&[-r*dx\\ &+ 0.5*v\_vec[j]*dx@Diag(L^2)\\ &+ 0.5*v\_vec[j]*dx2@Diag(L^2)]@f[:.j] \end{align*}
where:
- $f$ is an $m\times n$ matrix. $L$ is a vector with length $m$. $v\_vec$ is a vector with length $n$.
- $dx,dx2,dz,dz2$ are all differential operators (tri-diagonal matrix).
- $Diag(A)$ represents a square matrix with $A$ being diagonal.
- $'*'$ represents normal multiplication and $'@'$ represents matrix multiplication.
Then this form can directly fit into the ADI structure described in Tian(2015).
Could anyone see any problem in these ADI steps? Thanks a lot!
-----------------------------------------
Updates (0415):
Problem solved!
The key factor is the range of the variables - especially the variable $z=log(v/v0)$. It needs to be negative enough so $v$ can be close enough to zero.
