% 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]