1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Black Scholes heat equation form) Crank Nicolson

  1. Apr 22, 2017 #1
    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 (Text):
    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
     
     
  2. jcsd
  3. Apr 22, 2017 #2

    hilbert2

    User Avatar
    Science Advisor
    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.
     
  4. Apr 23, 2017 #3
    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: Apr 23, 2017
  5. Apr 23, 2017 #4
    ##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: Apr 23, 2017
  6. Apr 23, 2017 #5

    hilbert2

    User Avatar
    Science Advisor
    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.
     
  7. Apr 23, 2017 #6
    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:blushing:
     
    Last edited: Apr 23, 2017
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Black Scholes heat equation form) Crank Nicolson
  1. Heat Equation (Replies: 1)

  2. Heat equation (Replies: 2)

  3. The heat equation (Replies: 17)

Loading...