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

u_t = F(x,t,u,u_x,u_xx),     x₀ < x < x_m+1,     t₀ < t < t_n;

F(x,t,u,u_x,u_xx) = a(x,t) u_xx(x,t) + b(x,t) u_x(x,t) + c(x,t) u(x,t) + d(x,t).

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

a(x,t) = - 0.5 (sigma x)²,

b(x,t) = -r x;

c(x,t) = r;

d(x,t) = 0,

The Boundary Conditions (BC) are

u(0,t) = E exp(-r(T-t));

u(q E,t) = 0,

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

u(x,T) = max[E-x,0];

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

sigma = 9.07%     (i.e., market volatility of sigma = 0.0907);

r = 3.75%     (i.e., riskless bond interest rate of 0.0375 and NOT the parabolic mesh ratio k/h²);

E = $50 per share     (i.e., exercise price of fifty dollars per share) with a contract for the option to sell 25 shares at the exercise time T;

q = 2.0     (i.e., an arbitrary infinite domain cutoff factor as multiple of exercise price E);

T = 0.5 years     (i.e., 6 months contract);

m = 100     (i.e., one hundred interior nodes in stock price mesh);

n = 100     (i.e., one hundred time-to-go subintervals, small enough to keep the parabolic mesh ratio sufficiently small);

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)² for the pure diffusion term, vanishes at x=0. The maximum here is max(rpmreff(x)) = 0.5k(sigma max(x)/h)², which can be very large of large max[x].

If your method needs corner values, i.e., u(x₀,0) and u(x_m+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

U_i,j+1 = U_i,j + k F(x_i,t_j+Theta;,ThU_i,j, DThU_i,j,DDThU_i,j),    

(Revised)
ThU_i,j = (1-Theta)*U_i,j + Theta*U_i,j+1 ,

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,

   j       t_j

   
   where t_j=T-tg_j=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,

   i     x_i     Ueem_ij     Ucni_ij     Ufim_ij     HybDifEem_ij     HybDifFIM_ij

   
   in 4 significant digit exponential format for the above floating point headings, where Ueem_ij, Ucni_ij and Ufim_ij are the respective Euler's and Crank-Nicolson Implicit and Fully Implicit Method results, respectively, for U_i,j, while HybDifEem_ij is the hybrid-relative-absolute difference of Ueem_ij compared to Ucni_ij, i.e.,

   HybdifEem_ij = (Ueem_ij-Ucni_ij)/max(Ucni_ij,eps),

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

   HybDifFim_ij = (Ufim_ij-Ucni_ij)/max(Ucni_ij,eps),

   
   

3. 3D Scientific Visualization of Option Price versus (x,t): Plot in 3-dimensions all the output data u=U_i,j versus Stock Price x=x_i and Forward Time t=t_j 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?).