this is my first question so I hope I express myself clearly.

I'm trying to implement an Implicit and a Crank Nicolson algorithm for the generic PDE $\partial_\tau u(\tau,x)+a \partial_x^2 u(\tau,x) +d \partial_x u(\tau,x)+ c u(\tau,x)+d(\tau,x)=0$, where a,b,c are constants, x is the log-spot of an asset with black scholes dyanmics, that is $x=\log(S/S_0)$ and $\tau=T-t$ is the time to expiry. I solve the PDE forward from time 0 to expiry T.

Actually implementing both algorithms is quite straight forward, but my solution heavily depends on the choice of my spatial grid, namely left and right end points $x_\min$ and $x_\max$ appearing in

$x_\min < x_1 < \dots < x_N < x_\max $.

That means if I vary both values the solution variies, whereas if I fix both values and vary the discretization parameters in time $\Delta \tau$ and space $\Delta x$ the solution remains nearly unaffected.

Does someone know how to interpret this bevahiour? May it be possible that this is a direct consequence of the choice of my boundary conditions? They are probably not the best choice, but unfortunately the best I have. The solution of the PDE is the difference of a risk free OTC-derivative to that where counterparty credit risk, DVA and a funding cost adjustment is included.




Algorithm in full detail

Ok here are some details about the algorithm, but now in terms of Let us denote the values of the function at the grid points by $u_n^m= u(\tau_m,x_n)$. Then I approximate the spatial derivatives with central differences as usual and arrive at $ \partial_\tau u(\tau_m,x_n) = -a\frac{u_{n+1}^m-2u_n^m+u_{n-1}^m }{\delta x^2} - b \frac{u_{n+1}^m-u_{n-1}^m}{2\delta x}-cu_n^m -d_n^m\\[0.8em] =\left(-\frac{a}{\delta x^2}+\frac{b}{2\delta x}\right)u_{n-1}^m+\left(\frac{2a}{\delta x^2}-c\right)u_n^m+\left(-\frac{a}{\delta x^2}-\frac{b}{2 \delta x}\right) u_{n+1}^m-d_n^m $

or in more compact form

$\partial_\tau u(\tau_m,x_n)=\alpha_1 u_{n-1}^m+\alpha_2 u_n^m +\alpha_3 u_{n+1}^m -d_n^m :=g(u^m,d^m) $

The Crank Nicholson algorithm averages the forward time difference at point $(\tau_m,x_n)$ with the backward difference at $(\tau_{m+1},x_n)$

  • $\partial_\tau u(\tau_m,x_n) \approx \frac{u_n^{m+1}-u_n^m}{\delta t}$ forward difference at $(\tau_m,x_n)$
  • $\partial_t u(\tau_{m+1},x_n) \approx \frac{u_n^{m+1}-u_n^m}{\delta t}$ backward difference at $(\tau_{m+1},x_n)$

But values coincide and we average the spatial derivatives for both points

$\frac{u_n^{m+1}-u_n^m}{\delta t} = \frac{1}{2} g(u^m,d^m) +\frac{1}{2}g(u^{m+1},d^{m+1})$ or $u_n^{m+1} -\frac{\delta t }{2} g(u^{m+1},d^{m+1}) = u_n^m+\frac{\delta t }{2} g(u^m,d^m) $.

We can write the right hand side in matrix notation $\begin{pmatrix} u_1^m\\u_2^m \\ \vdots \\ \vdots \\ u^m_N \end{pmatrix} + \frac{\delta t }{2} \begin{pmatrix} \alpha_2 & \alpha_3 & 0 & \dots & 0 \\ \alpha_1 & \alpha_2 & \alpha_3 & & \vdots \\ 0 & \alpha_1 & \alpha_2 & \ddots & 0 \\ \vdots & 0 & \ddots & \ddots & \alpha_3 \\ 0 & 0 & & \alpha_1 & \alpha_2 \end{pmatrix} \begin{pmatrix} u_1^m\\u_2^m \\ \vdots \\ \vdots \\ u^m_N \end{pmatrix} + \frac{\delta t}{2} \begin{pmatrix} \alpha_1 u_0^m\\0 \\ \vdots \\ 0 \\ \alpha_3 u^m_{N+1} \end{pmatrix} - \frac{\delta t}{2} \begin{pmatrix} d_1^m\\d_2^m \\ \vdots \\ \vdots \\ d^m_N \end{pmatrix}$

Doing the same for the left hand side analogously and by defining By defining \begin{gather*} \boldsymbol{u}^m=\left(u_1^m,\dots,u_N^m\right)^T, \quad \boldsymbol{b}^m=\left(\alpha_1 u_0^m,0,\dots,0,\alpha_3 u_{N+1}^m\right)^T, \quad \boldsymbol{d}^m=\left(d_1^m,\dots,d_N^m\right)^T, \\ \eta= \frac{\delta t}{2} \end{gather*} we arrive at the Crank Nicholson algorithm. \begin{align*} \left(\mathbf{I}-\eta \mathbf{A}\right) \boldsymbol{u}^{m+1}- \eta \left(\boldsymbol{b}^{m+1}-\boldsymbol{d}^{m+1} \right)= \left(\boldsymbol{I}+\eta \mathbf{A}\right)\boldsymbol{u}^m + \lambda \left(\boldsymbol{b}^m-\boldsymbol{d}^m\right) \end{align*} Which can be written more compactly by \begin{align*} \mathbf{C}=&\mathbf{I}-\eta \mathbf{A} \qquad \text{and} \qquad \mathbf{D}=\mathbf{I}+\eta \mathbf{A},\\[0.6em] \boldsymbol{e}^m &=\eta` \left(\boldsymbol{b}^m+\boldsymbol{b}^{m+1}-\boldsymbol{d}^m-\boldsymbol{d}^{m+1}\right) \end{align*} the \textbf{Crank Nicolson Scheme} in matrix notation can be written in compact form \begin{align*} \mathbf{C} \boldsymbol{u}^{m+1} = \mathbf{D} \boldsymbol{u}^{m} +\boldsymbol{e}^m \end{align*}.

Boundary Conditions

The generic PDE is the result of a transformed PDE with coefficients depending on $S$.

I know by Feynman Kac theorem that the solution $U$ of the untransformed PDE satisfies

$U(t,S_t) =-c_1\int_t^T e^{-(r+\lambda_B+\lambda_C)(u-t)} \mathbb{E}_t \left[V^+(u,S(u))\right]du - c_2\int_t^T e^{-(r+\lambda_B+\lambda_C)(u-t)} \mathbb{E}_t \left[V^-(u,S(u))\right]du$

for some constants $\lambda_B,\lambda_C,c_1,c_2$. As I said $U$ is the difference from a risk free derivative $V$ and the same derivative but including counterparty credit risk. The derivative for which I want to determine the adjustment for the moment is a simple forward contract with value $V(t,S_t) = S_t - Ke^{-r(T-t)}$. I argue as follow:

We see that $V(t,S_t) \gg 0$ if $S_t \gg 0$. We assume then that \begin{align*} V(u,S_u) \gg 0 \quad \text{for all} \quad u\geq t \end{align*} Thus we can determine the upperboundary condition for $U$ by \begin{align*} U(t,S)&= -c_1\int_t^T D_{r+\lambda_B+\lambda_C}(t,u) \mathbb{E}_t \left[V(u,S(u)\right]du \\[0.8em] &= -c_1V(t,S_t)\int_t^T e^{-(\lambda_B+\lambda_C)(u-t)} du \\[0.8em] &= c_1\frac{1}{\lambda_B+\lambda_C} \left(1-e^{-(\lambda_B+\lambda_C)(T-t)}\right) \left(S_t - Ke^{-r(T-t)}\right). \end{align*} Analogously we can determine the lower boundary condition. We assume that if $S_t \to 0$ then $S_u \approx 0$ for all $u\geq t$. Therefore \begin{align*} V(u,S_u) = -Ke^{-r(T-t)} \quad \text{for all} \quad u\geq t \quad \text{as} \quad S_t \to 0 \end{align*} holds and therefore we can deduce \begin{align*} U(t,S)&= -c_2\int_t^T D_{r+\lambda_B+\lambda_C}(t,u) \mathbb{E}_t \left[V(u,S(u)\right]du \\[0.8em] &=- c_2V(t,S_t)\int_t^T e^{-(\lambda_B+\lambda_C)(u-t)} du \\[0.8em] &= c_2\frac{1 }{\lambda_B+\lambda_C} \left(1-e^{-(\lambda_B+\lambda_C)(T-t)}\right) Ke^{-r(T-t)} \end{align*} These considerations result in the boundary conditions. \begin{align*} u(\tau_m,x_{\text{max}} ) &=\frac{c_1}{\lambda_B+\lambda_C} \left(1-e^{-(\lambda_B+\lambda_C)\tau}\right) \left(x_0e^{x_{\text{max} } } - Ke^{-r\tau_m }\right) \\[0.8em] u(\tau_m,x_{\text{min}} )&= \frac{c_2}{\lambda_B+\lambda_C} \left(1-e^{-(\lambda_B+\lambda_C)\tau_m}\right) Ke^{-r\tau_m} \end{align*}

  • $\begingroup$ Could you provide a few more details about your numerical method? Can you write out the full finite difference discretization? I'm asking because the $\partial_x u$ can introduce some numerical instabilities if not done properly. $\endgroup$ – Tyler Olsen Dec 9 '15 at 16:19
  • $\begingroup$ Also, I believe you are seeing your boundary conditions. Can you write out what they are? You should not expect your solution to change much as you refine $\Delta x$ and $\Delta t$ (assuming they are sufficiently small to begin with). If you did see it changing, this would indicate that your numerical method is not stable, or that you have a bug in your implementation. $\endgroup$ – Tyler Olsen Dec 9 '15 at 16:23
  • $\begingroup$ push I hove that's allowed. $\endgroup$ – Simon Dec 11 '15 at 13:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.