% matlab code for solving                                                                 %%
% v_tau = v_{yy} + c v_y - d v         (d = 2 r/sigma^2,   c = d-1 for link to options)   %% 
% range:  ymin < y < ymax  .....  take -ymin = ymax = ln 3                                %%
%         0 < tau < 1/2 sigma^2 T                                                         %%
% initial cond v(y,0) = (1-e^y)_+                                                         %%
% bdry conds   v(ymin,tau) = e^{-d tau} - exp(ymin)                                       %%
%              v(ymax,tau) = 0                                                            %%
% and exploring the link with European put option prices                                  %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r     = 0.05;
sigma = 0.20;
T     = 1;
K     = 10;
d     = 2*r/sigma^2;               % can take any d here
c     = d-1;                       % can take any c here

ymin  = -log(3);
ymax  = +log(3);
N     = 200;                       
h     = (ymax-ymin)/N;
k     = h^2/3;                     % initial choice of k - well within stability
M     = ceil(sigma^2*T/2/k);       % "ceil" to make this an integer
k     = sigma^2*T/2/M;             % final choice of k 

y     = ymin:h:ymax;
tau   = 0:k:sigma^2*T/2; 
v     = zeros(N+1,M+1);

v(1,:) = exp(-d*tau)-exp(ymin);     % can take any a(tau) here - explain this choice 
v(:,1) = (1 - exp(y)).*(y<0);       % can take any f(y) here 

% solving by the Euler method 
for j=1:M
for i=2:N
    v(i,j+1) = v(i,j) + k/h^2*(v(i+1,j)-2*v(i,j)+v(i-1,j)) + ...
        c*k/2/h*(v(i+1,j)-v(i-1,j)) - d*k*v(i,j); 
end
end 

figure(1)    % three-d plot of solution 
Y   = y'*ones(1,M+1);
Tau = ones(N+1,1)*tau;
mesh(Y,Tau,v) 

figure(2)    % three-d plot of associated (European) option prices 
u = K * v;
t = T - 2*Tau/sigma^2; 
x = K * exp(Y);
mesh(x,t,u)

figure(3)    % 2-d plot of option prices for 9 values of tau 
hold on 
plot(x(:,1),u(:,1))
plot(x(:,ceil(M/8)),u(:,ceil(M/8)))
plot(x(:,ceil(2*M/8)),u(:,ceil(2*M/8)))
plot(x(:,ceil(3*M/8)),u(:,ceil(3*M/8)))
plot(x(:,ceil(4*M/8)),u(:,ceil(4*M/8)))
plot(x(:,ceil(5*M/8)),u(:,ceil(5*M/8)))
plot(x(:,ceil(6*M/8)),u(:,ceil(6*M/8)))
plot(x(:,ceil(7*M/8)),u(:,ceil(7*M/8)))
plot(x(:,M),u(:,M))
% note how options prices can be BELOW (K-x)_+ (not the case for American) 

% compare with the BS formula 
% [Call,Put] = blsprice(Price, Strike, Rate, Time, Volatility)
% note their time is "time remaining till expiration" = T-t in our approach 

[call put] = blsprice( x, K, r, T-t, sigma );

figure(4) 
mesh(x,t,put) 

% look at errors 
qq = abs(put - u);

figure(5)
mesh(x,t,qq) 

figure(6) 
hold on
plot(x(:,end),qq(:,end),'b')
plot(x(:,2),qq(:,2),'r')
% explain the source of error and look at how this goes down as N increases 