% drives the projected SOR algorithm for American options problems
% 
% our parameters are choosen to duplicate the Figure 9.4 (an American vs European put) 
% from Wilmott's book epage 194
%
% Written by:
% -- 
% John L. Weatherwax                2008-06-11
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; drawnow; 
clear all; 
clear functions; 

addpath( '../Chapter8' ); 

sigma = 0.4;
r     = 0.1; 
E     = 10.0; 
k     = r/(0.5*sigma^2); 
T     = 3/12;                % <- time to expiration (in years)
T     = 6/12;                % <- time to expiration (in years)
t     = 0;

% the change of variables from financial variables to dimensionless variables is given by 
% 
% x   = log(S/E)
% tau = 0.5*sigma^2 * (T-t) 
% v   = P/E
% u   = v * exp( -alpha * x - beta * tau )
% 

% initialize the "x=log(S/E)" grid ... we will convert to financial variables after solving 
%
SRight = 20.0;                % <- +\infty 
SLeft  = 1e-9;                % <- -\infty 
xLeft  = log(SLeft/E); 
xRight = log(SRight/E); 

Nx = 1000;
dx = (xRight-xLeft)/Nx; 

a  = 1.0; 

dt = a*dx^2;
% we integrate the diffusion equation from tau=0 (t=T expiration) to tau=0.5*sigma^2 (t=0 now)
tau_Max = (0.5*sigma^2)*T; 
M       = ceil(tau_Max/dt);

% 1) Solve the European put problem: 
%
%
[u,xgrid] = implicit_fd_LU(@tran_payoff_put, @u_m_inf_put, @u_p_inf_put, r, sigma, xLeft, xRight, Nx, tau_Max, M );

% the change of variables from dimensionless variables to financial variables is given by 
%
S   = E*exp( xgrid ); 
t   = 0;                       % <- the financial time when we want to evaluate our options value 
tau = 0.5*(sigma^2)*(T-t);     % <- translates into this scaled time 

Spow = (S.^(0.5*(1-k))); 
Smat = repmat( Spow(:).', [M+1, 1] ); 
V    = (E^(0.5*(1+k))) * Smat * exp( -(1/4)*((k+1)^2)*tau ).*u; 

figure; es=plot( S, V(end,:), '-x' ); grid on; hold on; xlabel( 'S' ); ylabel('P'); 

% evaluate for the S values's found in the table: 
Sgrid = [ 0.0+1e-6, 2, 4, 6, 8, 10:2:16 ]; 
Vgrid_eu = interp1( S, V(end,:), Sgrid, 'linear', 'extrap' ); 

% 2) Solve the American put problem  
%
%
[u,xgrid] = crank_fd_PSOR(@tran_payoff_put, @u_m_inf_put, @u_p_inf_put, r, sigma, xLeft, xRight, Nx, tau_Max, M );

% the change of variables from dimensionless variables to financial variables is given by 
%
S   = E*exp( xgrid ); 
t   = 0;                       % <- the financial time when we want to evaluate our options value 
tau = 0.5*(sigma^2)*(T-t);     % <- translates into this scaled time 

Spow = (S.^(0.5*(1-k))); 
Smat = repmat( Spow(:).', [M+1, 1] ); 
V    = (E^(0.5*(1+k))) * Smat * exp( -(1/4)*((k+1)^2)*tau ).*u; 

as=plot( S, V(end,:), '-or' ); grid on; hold on; xlabel( 'S' ); ylabel('P'); 

title( ['The Americna/European Put Example from Figure 9.4.  \alpha=',num2str(a),'T=',num2str(T)] ); 
legend( [ es,as ], {'European Put', 'American Put'}, 'location', 'north' ); 

fn=['fig_9_4.eps'];
saveas( gcf, ['../../WriteUp/Graphics/Chapter9/',fn], 'epsc' ); 

% evaluate for the S values's found in the table: 
Sgrid = [ 0.0+1.e-6, 2, 4, 6, 8, 10:2:16 ]; 
Vgrid_am = interp1( S, V(end,:), Sgrid, 'linear', 'extrap' ); 

fprintf('the option value for various asset prices.\n');
fprintf('First row are the prices S\n');
fprintf('Second row are the European option prices\n');
fprintf('Third row are the American option prices\n'); 
[Sgrid; Vgrid_eu; Vgrid_am]