Black Scholes heat equation form) Crank Nicolson

[tomatoyeungHK]

Dear all:

I am studying the finite difference method for black scholes equation(heat equation form), and my Prof. asked me to use the Crank Nicolson and Compact Scheme to do it. But I cannot reach second order for the space, i.e., x . The error only reach 1st order and sometimes not stable, it is totally contradict with the unconditionally stable property of Crank Nicolson Method. Is there anyone can give e some help in solving this problem? Deadline is coming and I cannot find out any problems. Thanks a lot. My Matlab code is here:

clear

%the length of the interval , set: x_min = -1.4, x_max = 5;
%the length of the interval is 6.4
%Parameters
K = 100; sigma = 0.25; r = 0.05;
beta = sigma*sigma/2;

%Discretization of x and t
step_x = 640;
step_t = 640;

dx = 6.4 /(step_x);
%the T =0.25
dt = 0.25/(step_t);

%Creating Vectors
x = linspace(-1.4,5,(step_x-1));
t = linspace(0,0.25,step_t );
U = zeros(step_x -1, step_t);

%Initial value and Boundary Value

u0_C= K*(max(exp(x) -1,0));

%Setting the left hand matrix, i.e. coeff. of k+1.th layer vectors
%upper diagonal
upp = (-(r-beta)/(4*dx) - beta/(2*dx*dx));
%main diagonal
main = (beta/(dx*dx) + r/2 + 1/dt);      
%lower diagonal
low = ((r-beta)/(4*dx) - beta/(2*dx*dx)); 

%main diagonal for right hand side,i.e. coeff. of k.th layer vectors
main_R = (beta/(dx*dx) + r/2 - 1/dt);

%Set up left hand matrix
l1 = [main, low, zeros(1,step_x -3)];
l2 = [main, upp, zeros(1,step_x -3)];
L_Mat = toeplitz(l1,l2);

%Set up right hand matrix
r1 = [-main_R,  -low, zeros(1,step_x -3)];
r2 = [-main_R,  -upp, zeros(1,step_x -3)];
R_Mat = toeplitz(r1,r2);

v = u0_C;
v = v';
for i = 1:step_t

      %boundary conditions for call option
      %Left points
      uL0 = K*(max(exp(-1.4)-exp(-r*(i-1)*dt),0));
      uL1 = K*(max(exp(-1.4)-exp(-r*(i*dt)),0));

      %Right points
      uR0 = K*(max(exp(5)-exp(-r*(i-1)*dt),0));
      uR1 = K*(max(exp(5)-exp(-r*(i*dt)),0));

      rhs = R_Mat*v;

      %Updating the boundary conditions
      rhs(1) = rhs(1) - low*uL0 - upp*uL1;
      rhs(end) = rhs(end)-upp*uR0 - low*uR1;

      v = L_Mat\rhs;
      U(:,i)=v;
  end

%U(:,i)
S=K*exp(x);

%U
[Call,Put] = blsprice(S,K,0.05,0.25,0.25);

plot(S,Call);
hold on
sol = U(:,end);
plot(S,sol,'r');
%finding the error at x=0,i.e. the error at exercise price K,
[x0,I0]=min(abs(x));

error_0=abs(sol(I0)-Call(I0));
error_0

%Norm_2 for the average error
error_2 = sqrt(6.4/step_x)*norm(sol'-Call,2);
error_2

 

[hilbert2]

Could you explicitly write the equation that you are solving, preferably using LaTeX code? If it's a linear PDE for a function f(x,t) of some two variables, x and t, you could test your code against some exactly solvable case, like a x-independent initial state f(x,0) = A or a linearly varying f(x,0) = Bx.

 

[tomatoyeungHK]

First of all, thanks for your reply.
I am so sorry for that I should state the problem much clear.
Here are some illustrations.

1st.
The Black-Scholes equation is a PDE in calculating the options in stock markets.
This is the original form of BS equation:
$\frac{\partial V}{\partial t}+ \frac{1}{2}\sigma^2S^2\frac{\partial^2 V}{\partial S^2}+ rS\frac{\partial V}{\partial S}-rV = 0$
Obviously, it is a backward equation and non-constant coefficient.

2nd
Here is the initial and boundary condition for call option:
$E$ is the exercise price and expiry date $T$ in the following.
$C(S,T) = max(S-E,0),\: for\: all\: S\geq0,$

We use the asset price boundary condition at zero asset price, $S =0, as\, S\rightarrow \infty.$ The PDE shows that if $S$ is ever $0$, $dS$ is zero also, it means $S$ can never change. As if $S = 0$ at expiry, the payoff will be 0, so the Call option is worthless on $S = 0$, even if there is a long time to expiry date. Then, on $S = 0$

$C(0,t) = 0,\: for\: all\: t\geq 0.$

At the end, as $S \rightarrow \infty$, it is more possible that the option will be exercised and the magnitude of the exercise price becomes less and less important. Therefore, as $S\rightarrow \infty$, the value of the option becomes: $asset - exercise\,price$, we need to pay to exchange for the asset.\\
For all $t>0$,

$C(S,t) \sim S - Ee^{-r(T-t)},\: as \: S\rightarrow \infty$

$V(S,t) = S$ itself is a solution to the Black-Scholes PDE through direct calculation.
Actually, $V(S,t) \equiv S$ for all $S$ and $t$ and then,

when $S = 0$, $V(0,t) = 0$ for all $t$;

when $S\rightarrow \infty$, $V(S,t) = S$ for all $t$,

and for $t = T$,$V(S,T) = S(T)$ for all $S$.

Unfortunately, it doesn't mean that the stock price can be obtained through solving the Black-Scholes PDE, because we do not know the exact value of the boundary conditions. For example, we don't know $S(T)$ at time $T$, when we are at time $t=0$.

 

[tomatoyeungHK]

$C(S,t) = C(S(t), E, T-t, r) = SN(d_{1}) - Ee^{-r(T-t)}N(d_{2}),$

where,

$N(d) = \frac{1}{\sqrt{2\pi}}\int^{d}_{-\infty}e^{{-\frac{1}{2}}s^{2}}ds,$

is the cumulative distribution function for the standard normal distribution,

$d_{1} = \frac{\ln\frac{S}{E} + (r + \frac{1}{2}\sigma^{2})(T-t)}{\sigma \sqrt{T-t}}$

$d_{2} = \frac{\ln\frac{S}{E} + (r - \frac{1}{2}\sigma^{2})(T-t)}{\sigma \sqrt{T-t}}$ For the boundary condition above tells us that,
$d_{1}, d_{2} \rightarrow -\infty \,as\, S \rightarrow 0,\, and\, N(-\infty) = 0.$
So, $C(0,t) = 0N(-\infty) - Ee^{-r(T-t)}N(-\infty) = 0.$

For the boundary condition
$C(S,t) \sim S - Ee^{-r(T-t)},\: as \: S\rightarrow \infty,$
$d_{1},\, d_{2} \rightarrow \infty\, as\, S \rightarrow \infty,$
and $N(\infty) = 1$. then,
$C(S,t) \rightarrow SN(\infty) - Ee^{-r(T-t)}N(\infty) \sim S - Ee^{-r(T-t)} \qquad as\:S\rightarrow\infty.$

3rd
We change the equation to a heat equation

For a Call option, the Black-Scholes PDE become:

$\left.\Biggl\{\begin{array}{ll} \frac{\partial C}{\partial t}+ \frac{1}{2}\sigma^2S^2\frac{\partial^2 C}{\partial S^2}+ rS\frac{\partial C}{\partial S}-rC = 0,\\ C(0,t) = 0,\quad C(S,t)\sim S-Ee^{-r(T-t)},\: as\: S\rightarrow\infty, \\ C(S,T) = max(S-E,0) \end{array}\right.$

To make the variables dimensionless and reverse in time, we can let
$S = Ee^{x}, \quad t = T - \tau\frac{1}{2}\sigma^{2},\quad C(S,t) = Ev(x,\tau).$

Substitute the variables into the equation, and we get

$\frac{\partial v}{\partial \tau} = \frac{\partial (C/E)}{\partial t}\frac{\partial t}{\partial \tau} = -\frac{1}{E}\frac{\partial C}{\partial t} (\frac{1}{\frac{1}{2}\sigma^2}),$
$\frac{\partial v}{\partial x} = \frac{\partial (C/E)}{\partial S}\frac{\partial S}{\partial x} = \frac{1}{E}\frac{\partial C}{\partial S}Ee^{x}= \frac{\partial C}{\partial S}e^{x},$
$\frac{\partial^2 v}{\partial x^2} = e^{x}\frac{\partial C}{\partial S} + e^{x}\frac{\partial^2 S}{\partial x^2}\frac{\partial S}{\partial x} = e^{x}\frac{\partial C}{\partial S} + e^{x}\frac{\partial^2 C}{\partial S^2}Ee^{x},$

Finally, the PDE becomes
$\frac{\partial v}{\partial \tau} = \frac{\partial^2 v}{\partial x^2}+(k-1)\frac{\partial v}{\partial x} - kv,$
where,$k = r/\frac{1}{2}\sigma^2.$The equation is already transformed into a forward parabolic equation.}
The initial condition becomes:
$v(x,0) = max(e^{x} - 1,0).$
The boundary conditions after transformation:

$\left\{\begin{array}{ll} v(x,\tau) = 0, & \qquad x\rightarrow -\infty \\ v(x, \tau)\sim(e^{x} - e^{-r\tau}), & \qquad x\rightarrow \infty \end{array} \right.$


The above are the whole derivation of BS equation to heat equation.

 

[hilbert2]

So, I guess you are solving the equation on some rectangular domain in $(S,t)$ plane: $S_{min} \leq S \leq S_{max}$, $t_{min} \leq t \leq t_{max}$.

If you take the original equation
$\frac{\partial V}{\partial t}+ \frac{1}{2}\sigma^2S^2\frac{\partial^2 V}{\partial S^2}+ rS\frac{\partial V}{\partial S}-rV = 0$
and try to find a stationary solution by setting $\frac{\partial V}{\partial t}=0$, you get an ordinary 2nd order DE that can be solved with Mathematica and you'll find out that it's solutions are of the form $V(S)=C_1 S^{a} + C_2 S^{b}$, where $a,b$ are exponents that depend on $\sigma , r$.

If you set $\frac{\partial V}{\partial S}=0$ everywhere, you get a solution that doesn't depend on $S$ but grows exponentially as a function of $t$, with a doubling time that depends on $r$.

This kind of solvable results can be compared to what your code predicts. It's best if you describe the problem as a purely mathematical question without putting any particular real-world meaning on the parameters, as most people on PF don't necessarily know much about economics.

 

[tomatoyeungHK]

In fact, my prof. he asked me to do the heat equation form, but I am afraid of the no one would understand what I am saying,
that's why I provide the some basic ideas on the previous post.
This is the equation form I am using actually:
$\frac{\partial v}{\partial \tau} = \frac{\sigma^2}{2}\frac{\partial^2 v}{\partial x^2}+(r-\frac{\sigma^2}{2})\frac{\partial v}{\partial x} - rv$,
and the initial condition:
$v(x,0) = max(Ke^{x} - K,0),$
boundary condition:
$v(x,\tau) \approx \left\{\begin{array}{ll} 0, & \qquad x\rightarrow -\infty \\ Ke^{x} - Ke^{-r\tau}, & \qquad x\rightarrow \infty \end{array} \right.$
Now, the equation can be computed from time = 0 to the terminal time = 0.25.

And my professor said that I am just a bachelor student who is doing a final year paper, no need to make the things too complicated.
Therefore, according to my understanding, I only use the difference to replace the differential term and obtain the coefficient matrix is enough.
Just like this:
$\frac{v^{k+1}_{i} - v^{k}_{i}}{\Delta t} =\frac{1}{2}[ \frac{\sigma ^2}{2}(\frac{v^{k}_{i-1} - 2v^{k}_{i} + v^{k}_{i+1}}{\Delta x^2} + \frac{v^{k+1}_{i-1} - 2v^{k+1}_{i} + v^{k+1}_{i+1}}{\Delta x^2})] +\frac{1}{2}[(r-\frac{\sigma ^2}{2})(\frac{v^{k}_{i+1} - v^{k}_{i-1}}{2\Delta x} + \frac{v^{k+1}_{i+1} - v^{k+1}_{i-1}}{2\Delta x}) ] - \frac{r}{2}(v^{k+1}_{i} + v^{k}_{i})$
** in the equation, I used $\beta$ to replace the $\frac{\sigma^2}{2}$ for easier calculation.
and then we obtain the equation:
$[\frac{r-\beta}{4(\Delta x)} - \frac{\beta}{2\Delta x^2}]v^{k+1}_{i-1} + [\frac{\beta}{\Delta x^2} + \frac{r}{2} + \frac{1}{\Delta t}]v^{k+1}_{i} + [\frac{-(r-\beta)}{4(\Delta x)} - \frac{\beta}{2\Delta x^2}]v^{k+1}_{i+1}$
=
$-[\frac{r-\beta}{4(\Delta x)} - \frac{\beta}{2\Delta x^2}]v^{k}_{i-1} + -[\frac{\beta}{\Delta x^2} + \frac{r}{2} - \frac{1}{\Delta t}]v^{k}_{i} -[\frac{-(r-\beta)}{4(\Delta x)} - \frac{\beta}{2\Delta x^2}]v^{k}_{i+1}$

then:

Let $\alpha = [\frac{r-\beta}{4(\Delta x)} - \frac{\beta}{2\Delta x^2}], \psi = [\frac{\beta}{\Delta x^2} + \frac{r}{2} + \frac{1}{\Delta t}], \gamma = [\frac{-(r-\beta)}{4(\Delta x)} - \frac{\beta}{2\Delta x^2}], \xi = [\frac{\beta}{\Delta x^2} + \frac{r}{2} - \frac{1}{\Delta t}]$

the left hand matrix:
$\begin{equation} A= \left[ \begin{matrix} \psi & \gamma & & 0\\ \alpha & \ddots & \ddots & \\ & \ddots & \ddots & \gamma \\ 0 & & \alpha & \psi \end{matrix} \right] \end{equation},$
$\begin{equation} B= \left[ \begin{matrix} -\xi & -\gamma & & 0\\ -\alpha & \ddots & \ddots & \\ & \ddots & \ddots & - \gamma \\ 0 & &- \alpha & -\xi \end{matrix} \right] \end{equation},$

$\begin{align} A\cdot \begin{bmatrix} v_{1}^{k+1} \\ v_{2}^{k+1} \\ \vdots \\ v_{m-1}^{k+1} \end{bmatrix} \nonumber = B\cdot \begin{bmatrix} v_{1}^{k} \\ v_{2}^{k} \\ \vdots \\ v_{m-1}^{k} \end{bmatrix} \nonumber + \begin{bmatrix} \gamma( v_{0}^{k}) + \alpha(v_{0}^{k+1}) \\ 0 \\ \vdots \\ \alpha(v_{m-1}^{k}) + \gamma(v_{m-1}^{k+1}) \end{bmatrix} \nonumber \end{align}$

almost like this,
but I am just think that the zeros in the last column vector of right hand side,do I need to add something.
because I just read one of my book again which writes the last column matrix in its derivation in the following way:
(of course the coefficients are not r/2 in this case)
$\begin{bmatrix} \frac{r}{2}(u_{0}^{k} + u_{0}^{k+1}) + \tau f(x_{1},t_{k+1/2}) \\ \tau f(x_{2},t_{k+1/2}) \\ \vdots \\ \tau f(x_{m-2},t_{k+1/2})\\ \frac{r}{2}(u_{m}^{k} + u_{m}^{k+1}) + \tau f(x_{m-1},t_{k+1/2}) \end{bmatrix} \nonumber$

I hope that you understand what I am saying