% drives the various finite difference methods from Chapter 8
% 
% 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( '../Chapter3' );   % <- in analytic solution to the black-scholes solution 

% initialize some parameters of our method:
% -- these are choosen to duplicate parts of Figure 8.6 (an European put) from Wilmott's 
%    book epage 157
% 
r     = 0.05;
sigma = 0.2; 
E     = 10.0; 
T     = 0.5;                % <- time to expiration (in years)
k     = r/(0.5*sigma^2); 
  
% 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 = 16.0;                % <- +\infty 
SLeft  = 1e-9;                % <- -\infty 
xLeft  = log(SLeft/E); 
xRight = log(SRight/E); 

Nx = 2000;
dx = (xRight-xLeft)/Nx; 
a  = 0.25; 
a  = 0.5;  
a  = 0.52; 
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);
fprintf('running alpha = %10.6f\n',a); 
[u,xgrid] = explicit_fd(@tran_payoff_put, @u_m_inf_put, @u_p_inf_put, r, sigma, xLeft, xRight, Nx, tau_Max, M );
%[u,xgrid] = explicit_fd(@tran_payoff_call, @u_m_inf_call, @u_p_inf_call, r, sigma, xLeft, xRight, Nx, tau_Max, M );

% solve some cash-or-nothing problems ... 
% 
if( 0 ) 
  B = 0.5*E; 
  global b; b = B/E; 
  %[u,xgrid] = explicit_fd(@tran_payoff_CON_put, @u_m_inf_CON_put, @u_p_inf_CON_put, r, sigma, xLeft, xRight, Nx, tau_Max, M );
  [u,xgrid] = explicit_fd(@tran_payoff_CON_call, @u_m_inf_CON_call, @u_p_inf_CON_call, r, sigma, xLeft, xRight, Nx, tau_Max, M );
end
  
% 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; ns=plot( S, V(end,:), 'x' ); grid on; hold on; xlabel( 'S' ); ylabel('C'); 

%return; 

[C,P] = blackscholes(E,T,S,r,sigma); 
%[Call,Put] = blsprice(S,E,r,T,sigma);    % <- use the MATLAB function ... as a check ... 
if( V(end,end-1)<1.0 )   % <- a heuristic as to which option we are pricing
  as=plot( S, P, '-og' ); 
  %plot( S, Put, '-dk' ); 
  title( ['The European Put Example from Figure 8.6.  \alpha=',num2str(a)] ); 
  fn=['fig_8_6_a_put_',num2str(a),'.eps']; 
else
  as=plot( S, C, '-og' ); 
  %plot( S, Call, '-dk' ); 
  title( ['The European Call Example.  \alpha=',num2str(a)] ); 
  fn=['fig_8_6_a_call_',num2str(a),'.eps']; 
end
legend( [ ns,as ], {'numerical solution', 'analytic solution'}, 'location', 'north' ); 
saveas( gcf, ['../../WriteUp/Graphics/Chapter8/',fn], 'eps' ); 

% evaluate for the S values's found in the table: 
Sgrid = [ 0.00, 2, 4, 6, 7, 8, 9, 10:16 ]; 
Vgrid = interp1( S, V(end,:), Sgrid, 'linear', 'extrap' ); 
Vgrid