MCS 571 Computer Project 1 -- Numerical PPDE Methods -- Spring 2004

Theta Algorithm Applied to Financial RisK Elimination in Option Pricing:
MATLAB or C/Fortran Programming Computation


MATLAB or C/Fortran Programming Code and Output are Due Monday 05 April 2004 in Class .


General Instructions:

Write a single, efficient program, that compares the performance of the three (3) values of Theta for the Theta Algorithm:

  1. Theta = 0.0: Euler's Explicit Method (EEM)

  2. Theta = 0.5: Crank-Nicolson Implicit Method (CNIM)

  3. Theta = 1.0: Fully Implicit Method (FIM)

Apply these methods to the inhomogeneous, Parabolic Partial Differential Equation (PPDE):


where in the linear case,


Special Test Coefficients:

Consider the special test example to test your general code. The application is to the Black-Scholes financial options pricing model for the case of European put options. The put option contract gives the contractee the optional right to sell a specified number of shares of stock at a price E at exercise time T to the option contractor. A rational put option contractee bets that the market stock price will fall below price E at T so that stock at the new price x can be sold at the exercise price E to the contractor at a profit of E-x, otherwise the option would not be exercised by a rational contractee.

Unlike the usual forward PPDE, the signs on the terms are different since the PDE is a backward PPDE that starts at the final or exercise time T and is integrated backwards to the initial time to find the best starting price.

The coefficients (COEFF) are





on 0 < x < q E and 0 < x < T.

The Boundary Conditions (BC) are



on 0 < t < T.

Unlike PPDEs discussed in class this problem does NOT have an initial condition, but a Final Condition (FC) that is the Put Option Payoff:


on 0 < x < q E. Due to (FC), this is a Backward Diffusion Problem and needs to be solved backward, starting from t = T and ending at t = 0. Hence, it is suggested that you change the time variable from t to the time-to-go tg = T-t which is equivalent to reversing the signs in the spatial coefficients (a,b,c,d), starting from tg = 0 and ending tg = T, making a proper forward diffusion equation with positive diffusion coefficient a(x,t).


Test Parameter Specifications:

Your code should be written for general parameter variables, but the code should be tested with the model parameters that have been chosen to be








Note: The problem is called a Singular Diffusion, since the diffusion coefficient a(x,t) vanishes as x --> 0, so that the x-dependent effective parabolic mesh ratio, here rpmreff(x) = 0.5k(sigma x/h)2 for the pure diffusion term, vanishes at x=0. The maximum here is max(rpmreff(x)) = 0.5k(sigma max(x)/h)2, which can be very large of large max[x].

If your method needs corner values, i.e., u(x0,0) and u(xm+1,0), use the average values of the boundary and initial conditions, since that is what the numerical approximations will converge to eventually. However, corner values should not be necessary.

Warning: Do not write your code only for numerical values given here, but write in general variables and not just for this test case.


Method Specifications:

The Theta Algorithm is


where
while DThUi,j and DDThUi,j, all constructed using ThUi,j as discrete values.

Use the Thomas Tridiagonal Solver Algorithm to solve the cases for which Theta > 0.


Output Specifications:

  1. Input Parameters: Print out the values of the input parameters including the model parameters and the numerical parameters. Also, print out the diffusion term max(rpmreff(x)) as an indicator of stability.

  2. Output Table: Tabulate, with good labels, for j = 0:4 (i.e., the first few points for debugging and accuracy) and j=25:25:n (i.e., with step 25 in this variable time sample of option prices) using time headings,
    where tj=T-tgj=T-j dt and dt is the time step, while for i = 0 : 10 : m (i.e., in steps of 10) for each sampled time form the table,
    in 4 significant digit exponential format for the above floating point headings, where Ueemij, Ucniij and Ufimij are the respective Euler's and Crank-Nicolson Implicit and Fully Implicit Method results, respectively, for Ui,j, while HybDifEemij is the hybrid-relative-absolute difference of Ueemij compared to Ucniij, i.e.,

      HybdifEemij = (Ueemij-Ucniij)/max(Ucniij,eps),


    where eps = 5.e-1 is used to avoid Ucniij = 0 or bad properties of the relative error for very small values, by keeping the denominator no smaller than 50 cents. Similarly,

      HybDifFimij = (Ufimij-Ucniij)/max(Ucniij,eps),


  3. 3D Scientific Visualization of Option Price versus (x,t): Plot in 3-dimensions all the output data u=Ui,j versus Stock Price x=xi and Forward Time t=tj for i=0:m+1 and j=0:n, separately for EEM, CNIM and FIM. If you use MATLAB, see the function mesh(t,x,u) along with view(45,30) and axis tight;, say. Be sure to label your plot with informative title, xlabel and ylabels. If you use C, C++ or Fortran language programming, you can export your output to a file and then input it into MATLAB just for plotting or use some other comparable plotting program.

  4. 2D Scientific Visualization of Initial Price, Final Price, Net Profit: Plot in 2-dimensions (MATLAB plot function) the numerical approximations to the initial option price u(x,0), the final option price u(x,T), the net option profit per share u(x,T)-u(x,0) and the net option profit for the contracted 25 shares 25*(u(x,T)-u(x,0)) on a single graphs with labels and legends, with one graph for each of the three methods.

  5. 2D Scientific Visualization of Hybrid Differences: Also, plot HybDifEem and HybDifFim on the same plot at the initial time t = 0 on the same plot versus the stock price x with labels and a legend.

  6. Timings: In addition, report the timings in seconds (see the the clock function with etime elapsed time function, or other timing functions with higher resolution) for each method, counting common initializations of terms like coefficient arrays for {a, b, c, d} and output in your timing loop.

  7. Discussion of Results: Discuss the significance of your results, e.g., stability and the comparison of the results for the three Theta methods.


Some Proper Efficiency Suggestions:

  1. Be Efficient in storage and floating point calculations.

  2. Compute once and store the spatially dependent coefficients in an array to use in all time steps.

  3. Avoid calculating items like constant expressions in loops, but keep code as general as practical.

  4. Use the Thomas Tridiagonal Algorithm for CNIM and FIM.

  5. In C or Fortran programming, you must use double floating point precision (MATLAB uses double precision automatically). Float or single precision (real) is inadequate for scientific, statistical or engineering computations.

Added Note: These methods can also be used for nonlinear PPDEs, directly for EEM, but with a predictor-corrector modification for CNIM and FIM. Upwinding does not seem necessary here (Why?).


Web Source: http://www.math.uic.edu/~hanson/mcs571/pcp571S04.html

MCS 571 HomePage: http://www.math.uic.edu/~hanson/mcs571/

Email Comments or Questions to Professor Hanson.