I have a problem numerically solving the following PDE with boundary conditions:
$$
u_t + \frac{x^2\sigma^2}2u_{xx} + rxu_x – ru = 0 \quad (x,t) \in (0,N) \times (0,T)
$$
with
$$
u(x,T) = \max\{0,x-K\}˛ \quad u(0,t) = 0, \quad u(N,t) = N – K.
$$
(This is the Black Scholes PDE to determine the fair price of an European call option.) Originally $x$ should vary over $\mathbb{R}^+$. However to numerically solve the problem, I should choose $N$ sufficiently large. The given parameters are $K = 12, r = 0.04, \sigma = 0.4$ and $T = 2$. I am supposed to use the explicit Euler method.
First, I discretise the underlying space using $m,n \in \mathbb{N}$, I define $x_i = h i$ and $t_j = \tau j$ with $h = \frac Nn, \tau = \frac Tm$ and $i \in \{0,\dots,n\}, j \in \{0,\dots,m\}$. I use the recursion, using $u_i^j := u(x_i,t_j)$,
$$
\frac1\tau(u_i^j – u_i^{j-1}) + \frac{x_i^2\sigma^2}2\frac1{h^2}(u_{i+1}^{j} – 2u_{i}^{j} + u_{i-1}^{j}) + rx_i \frac1{2h}(u_{i+1}^{j} – u_{i-1}^{j}) – r u_{i}^{j-1} = 0.
$$ (Notice that I use $w_i^{j-1}$ as the last term and not $w_i^j$). Rearranging, I get
the recursion
$$
w_i^{j-1} = \frac1{1+r\tau}\left(u_i^j + \frac{x_i^2\sigma^2}2\frac{\tau}{h^2}(u_{i+1}^{j} – 2u_{i}^{j} + u_{i-1}^{j}) + rx_i \frac{\tau}{2h}(u_{i+1}^{j} – u_{i-1}^{j})\right).
$$
Translating my boundary conditions yields
$$
w_0^j = 0, \quad w_n^j = N – K, \quad w_i^m = \max\{0, x_i – K\}.
$$
Now my Problem
When implementing this procedure (using R) (I use e.g. $N = 30$, $m = 5000$ and $n = 200$), I get a really unstable behaviour.
Here you can see 3 pictures, where you can see the (in that order) time steps $t_j$ for $j = 2962, 2959, 2956$
The solution seems to oscillate more and more and I am not exactly sure why. Maybe just due to numerical errors? This behaviour gets worse when I consider even smaller $t_j$, up until a point where it is unrecognisable.
For $n = 100$ the solutions looks fine, with no unstable behaviour at all. I would suspect that I get more errors when $n$ is too small, but here it is the other way around.
Unfortunately I have basically no experience for numerically solving a PDE and am not sure what I can try next and why I see this behaviour.
If you need more information (code, etc.) let me know.
Thanks in advance!
Best Answer
Edit:
Sorry, I just realise that the smooth graph below was obtained with a vol of $$ \sigma=0.2 $$ instead of your $\sigma=0.4\,.$ So it seems that the large vol is the reason for the instability, not the space step/time step configuration which seems totally OK. (end of edit)
I have experimented with a bit of python code myself and also consulted the classic book P. Willmot, S. Howison, J. Dewynne, The Mathematics of Financial Derivatives (1995). The authors write on p. 143 that an explicit scheme (yours) is stable if and only if (in your notation) $$ 0< \alpha:=\frac{\tau}{h^2}=\frac{T\,n^2}{N^2\,m}\le \frac{1}{2}\,. $$ In other words, when you have a stable configuration and double the number of space steps $n$ you have to quadruple the number $m$ of time steps to keep stability.
In your example we have $N=30,T=2,n=200,m=5000$ which gives $$ \alpha=0.01777 $$ which should give a stable result, and indeed I am getting a smooth graph which agrees very well with the closed form BS solution, unless the underlying is far ITM:
With $\sigma=0.2$ that I used the probability that the underlying is above $N=30$ in two years when it is at $13$ today is $$ \mathbb P(X_T>N)\approx 0.06\% $$ which is I think small enough. With your $\sigma=0.4$ this would be $2\%$ which I think is too large. Then also, my code fails.
I believe that the discrepancy of the graphs near $N=30$ shows that for large vols you have to take a far larger $N$ and ensure again that the stability parameter $\alpha$ is small.
Code: