% Daniel Sevcovic, Beata Stehlikova, Karol Mikula
% Analyticke a numericke metody ocenovania financnych derivatov
% Nakladatelstvo STU, Bratislava, 2009. 200 stran
% ISBN 978-80-229-3014-3

clear all; close all;

sigma=0.5; r=0.04; q=0.2;
T=0.5; e=80;

alfa=(r-q)/(sigma^2)-0.5;
beta=0.5*(r+q)+(sigma^2)/8 + ((r-q)^2)/(2*sigma^2);

NN=100; n=2*NN-1;
m=100;

L=4;
k=T/m;
h=L/NN;

gam=k*sigma^2/(2*h^2);

A=zeros(n,n);
for i=1:n
    A(i,i)=1+2*gam;
end;
for i=1:n-1
    A(i,i+1)=-gam;
    A(i+1,i)=-gam;
end;


index=1;
for i=-NN+1:NN-1
    u0(index)=exp(alfa*i*h)*max(exp(i*h)-1,0);
    index=index+1;
end;
u0=u0';

% TREBA UPRAVIT OKRAJOVE PODMIENKY - NEMOZU KLESNUT POD TRANSFORMOVANY PAYOFF
for j=0:m-1
    phi(j+1)=0;
    psi(j+1)=exp((alfa+1)*NN*h+(beta-q)*j*k);
end;

ustare=u0;

for j=1:m
    b=zeros(n,1); b(1)=gam*phi(j); b(n)=gam*psi(j);
       
    % gvec bude vektor, s ktorym budeme porovnavat iteracie
    % kazda zlozka vektora u musi byt >= prislusna zlozka vektora gvec
    gvec=zeros(1,n);
    index=1;
    for i=-NN+1:NN-1
        gvec(index)= % ... VYTVORTE VEKTOR gvec
        index=index+1;
    end;
    
    omega=1.7;
    
    uvec=ustare;
    fvec=ustare+b;
    
    % POROVNAJTE ZACIATOCNU ITERACIU S gvec 
    % AK JE NIEKDE MENSIA, PREPISTE JU PRISLUSNOU ZLOZKOU gvec
    
    % .....
    
    
    for kiter=1:10
        uvec(1)=(1-omega)*uvec(1)+(omega/A(1,1))*(fvec(1)-A(1,2)*uvec(2));
        uvec(1)= % DOPLNTE POROVNANIE S PRISLUSNOU ZLOZKOU gvec
    
        for i=2:n-2
            uvec(i)=(1-omega)*uvec(i)+(omega/A(i,i))*( fvec(i)-A(i,i-1)*uvec(i-1)-A(i,i+1)*uvec(i+1));
            uvec(i)= % DOPLNTE POROVNANIE S PRISLUSNOU ZLOZKOU gvec
        end;
        uvec(n-1)=(1-omega)*uvec(n-1)+(omega/A(n-1,n-1))*( fvec(n-1)-A(n-1,n-2)*uvec(n-2));
        uvec(n-1)= % DOPLNTE POROVNANIE S PRISLUSNOU ZLOZKOU gvec
    end;

    unove=uvec;        
    ustare=unove;
end;

callPayoff=inline('max(0,s-e)','s','e');
call=inline('s.*exp(-d*tau).*normcdf(((r-d+0.5*sigma^2)*tau+log(s/e))/(sigma*sqrt(tau)))-e.*exp(-r*tau).*normcdf(((r-d-0.5*sigma^2)*tau+log(s/e))/(sigma*sqrt(tau)))','s','e','tau','sigma','r','d');

j=m-1;
ind=1;
for i=-NN/2+1:NN-1
    Vnove(ind,1)=e*exp(i*h);
    Vnove(ind,2)=e*exp(-alfa*i*h-beta*j*k)*unove(i+NN);
    Vcall(ind)=call(e*exp(i*h),e,T,sigma,r,q);
    ind=ind+1;
end;

figure(1);
subplot(1,2,1);
plot(Vnove(:,1),Vnove(:,2),'blue'); xlabel('S'); ylabel('V'); hold on; % (S,V numericke pre US call)
sVect=50:0.01:120;
plot(sVect, callPayoff(sVect,e),'red'); hold on; % (S, payoff)
plot(Vnove(:,1),Vcall,'green');  % (S, europsky call)
axis([50 120 -1 40]);

subplot(1,2,2);
plot(Vnove(:,1),Vnove(:,2),'blue'); xlabel('S'); ylabel('V'); % (S,V numericke pre US call)
hold on;
sVect=50:0.01:120;
plot(sVect, callPayoff(sVect,e),'red'); hold on;   % (S, payoff)
plot(Vnove(:,1),Vpresne,'green');   % (S, europsky call)
axis([90 120 -1 40]);

