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:
 Theta = 0.0: Euler's Explicit Method (EEM)
 Theta = 0.5: CrankNicolson Implicit Method (CNIM)
 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_{0} < x < x_{m+1},
t_{0} < t < t_{n};
where in the linear case,
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 BlackScholes 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 Ex, 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)^{2},
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
timetogo tg = Tt 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
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^{2});
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 timetogo 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 xdependent 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(x_{0},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}),
where
(Revised)
ThU_{i,j} = (1Theta)*U_{i,j} + Theta*U_{i,j+1} ,
while
DThU_{i,j}
and
DDThU_{i,j}, all constructed using ThU_{i,j}
as discrete values.
Use the Thomas Tridiagonal Solver Algorithm to solve the cases for which
Theta > 0.
Output Specifications:
 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.
 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 t_{j}=Ttg_{j}=Tj 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 Ueem_{ij}, Ucni_{ij} and
Ufim_{ij} are the
respective Euler's and CrankNicolson Implicit and Fully Implicit Method
results, respectively, for
U_{i,j}, while HybDifEem_{ij} is the
hybridrelativeabsolute difference of Ueem_{ij} compared to
Ucni_{ij}, i.e.,
HybdifEem_{ij} =
(Ueem_{ij}Ucni_{ij})/max(Ucni_{ij},eps),
where eps = 5.e1 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),

3D Scientific Visualization of Option Price versus (x,t):
Plot in 3dimensions 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.

2D Scientific Visualization of Initial Price, Final Price, Net Profit:
Plot in 2dimensions (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.

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.

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.

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:

Be Efficient in storage and floating point calculations.

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

Avoid calculating items like constant expressions in loops, but keep
code as general as practical.
 Use the Thomas Tridiagonal Algorithm for CNIM and FIM.
 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 predictorcorrector 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.