# Black Scholes heat equation form) Crank Nicolson

• tomatoyeungHK
In summary, the conversation is about using the Crank Nicolson and Compact Scheme to solve the finite difference method for the Black Scholes equation, but the user is facing issues with achieving a second order for the space variable. They also mention using Matlab code and testing it against an exactly solvable case. The equation being solved is the Black Scholes equation for a call option with initial and boundary conditions stated, and it is explained that the stock price cannot be obtained through solving this equation due to unknown boundary conditions.
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:
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

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.

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

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

## What is the Black-Scholes heat equation in Crank-Nicolson form?

The Black-Scholes heat equation is a partial differential equation used in mathematical finance to model the behavior of financial options. In Crank-Nicolson form, it is a finite difference approximation of the original equation, which allows for efficient numerical computation.

## What is the significance of the Black-Scholes heat equation in financial modeling?

The Black-Scholes heat equation is significant because it provides a theoretical framework for pricing financial options. It takes into account factors such as stock price, time to expiration, and volatility, and allows for the calculation of a fair price for an option.

## What is the Crank-Nicolson method and how is it related to the Black-Scholes heat equation?

The Crank-Nicolson method is a numerical technique used to solve partial differential equations, including the Black-Scholes heat equation. It uses a combination of forward and backward time steps to approximate the solution, resulting in a more accurate and stable result compared to other methods.

## What are the limitations of the Black-Scholes heat equation in Crank-Nicolson form?

One limitation of the Black-Scholes heat equation in Crank-Nicolson form is that it assumes a constant volatility and interest rate, which may not always hold true in real-world financial markets. It also does not take into account factors such as transaction costs, dividends, and market crashes.

## How is the Black-Scholes heat equation in Crank-Nicolson form used in practice?

The Black-Scholes heat equation in Crank-Nicolson form is often used in finance to price options and calculate risk measures. It is also used in the calibration of option pricing models and in the development of trading strategies. Additionally, it is a key component in the valuation of other financial instruments, such as convertible bonds and credit derivatives.

• Differential Equations
Replies
23
Views
4K
• Differential Equations
Replies
2
Views
2K
• Programming and Computer Science
Replies
4
Views
4K
• MATLAB, Maple, Mathematica, LaTeX
Replies
3
Views
3K
• MATLAB, Maple, Mathematica, LaTeX
Replies
41
Views
8K
• Differential Equations
Replies
12
Views
8K
• Engineering and Comp Sci Homework Help
Replies
6
Views
2K
• Differential Equations
Replies
7
Views
3K
• MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
4K
• Programming and Computer Science
Replies
2
Views
4K