% This MATLAB script drives the projected SOR algorithm for pricing American options
% 
% These parameters are choosen to duplicate the Figure 9.5 (an American vs European call)
% 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; 

%close all; drawnow; 
clear all; 

addpath( '../Chapter8' ); addpath( '../Chapter3' ); 

sigma = 0.8;
r     = 0.25; 
E     = 10.0; 
D0    = 0.2; 
%D0    = 0; 
k     = (r-D0)/(0.5*sigma^2); 
T     = 1/3;                 % <- 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   = C/E
% u   = v * exp( -alpha * x - beta * tau )
% 

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

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

a  = 1.0; 

dtau = 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/dtau);

% 1) Solve the European call problem (with dividends): 
%    -- note I'm not sure that this is correct ... 
%       the European call with dividends maybe slightly more complicated that this ... 
%       see the section in the book 6.2.2 -> A Constant Dividend Yield 
%
if( 0 ) 
  [u,xgrid] = implicit_fd_LU(@tran_payoff_call, @u_m_inf_call, @u_p_inf_call, r-D0, 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 ... this is also called tau_Max 

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

  % This is the function V_1(S,t) option prices ... discussed in the text 
  % -- modify them to incorporate the dividend yield D0
  % 
  tArray = 0:( tau_Max/M ): tau_Max; tArray = tArray.'; tArray = flipud( tArray ); 
  expPt  = repmat( exp( -D0*( T - tArray ) ), [1,Nx+1] ); 
  V      = expPt .* V_1; 
  fh=figure; es=plot( S, V(end,:), '-x' ); grid on; hold on; xlabel( 'S' ); ylabel('C'); 

  % 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' ); 

  clear V u; 
end

% 2) Solve the American call problem  
%
[u,xgrid] = crank_fd_PSOR(@tran_payoff_call, @u_m_inf_call, @u_p_inf_call, r-D0, 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_1  = (E^(0.5*(1+k))) * Smat * exp( -(1/4)*((k+1)^2)*tau ).*u; 

% This is the function V_1(S,t) option prices ... discussed in the text 
% -- modify them to incorporate the dividend yield D0
% 
tArray = 0:( tau_Max/M ): tau_Max; tArray = tArray.'; tArray = flipud( tArray ); 
expPt  = repmat( exp( -D0*( T - tArray ) ), [1,Nx+1] ); 
V      = expPt .* V_1; 

fh=gcf; 
figure(fh); as=plot( S, V(end,:), '-or' ); grid on; hold on; xlabel( 'S' ); ylabel('C'); 

[C,P] = blackscholes(E,0,S,r-D0,sigma); 
C = C .* expPt(end,:); 
%[C,P] = blsprice(S,E,r,T,sigma,D0);    % <- use the MATLAB function ... as a check ... 


figure(fh); bss=plot( S, C, '-k', 'LineWidth', 2 ); 

title( ['The American/European Call Example from Figure 9.5.  \alpha=',num2str(a),' T=',num2str(T)] ); 
%legend( [ es,as,bss ], {'European Call', 'American Call', 'Black-Scholes Analytic Solution'}, 'location', 'northwest' ); 
legend( [ as,bss ], {'American Call', 'Payoff at Expiration'}, 'location', 'northwest' ); 
axis( [0,30,0,22] ); 

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

% evaluate for the S values's found in the table: 
Sgrid = [ 0.0+1e-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 stock prices S\n');
%fprintf('Second row are the European option prices\n');
fprintf('The next row are the American option prices\n'); 
[Sgrid; Vgrid_am]