# Black Scholes heat equation form) Crank Nicolson

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 any one 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:

Code:
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
Gold Member
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.

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##.

Last edited:
##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.

Last edited:
hilbert2
Gold Member
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.

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:
##

A= \left[ \begin{matrix}
\psi & \gamma & & 0\\
\alpha & \ddots & \ddots & \\
& \ddots & \ddots & \gamma \\
0 & & \alpha & \psi \end{matrix} \right]
,##
##

B= \left[ \begin{matrix}
-\xi & -\gamma & & 0\\
-\alpha & \ddots & \ddots & \\
& \ddots & \ddots & - \gamma \\
0 & &- \alpha & -\xi \end{matrix} \right]
,##

##
\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

Last edited: