n=12
T=1
delta_t=T/n
r=0.05
sigma=0.2
S0=140
K=S0
L=160
up=exp(sigma*sqrt(delta_t))
down=1/up
R=exp(r*delta_t)
p=(R-down)/(up-down)

// QUESTION 1

SS=zeros(n+1,n+1)
SS(1,1)=S0
for i=1:n
    SS(i+1,1)=SS(i,1)*down
    for j=1:i
        SS(i+1,j+1)=SS(i,j)*up
    end
end


function fi=phi(x)
    fi=max(0,x-K)
endfunction

function psu=psi(x)
    psu=max(0,K-x)
endfunction

CC=zeros(n+1,n+1)
for j=1:(n+1)
	CC(1+n,j)=phi(SS(1+n,j))
end
for i=n-1:-1:0
    for j=0:i
        CC(i+1,j+1)=(p*CC(i+2,j+2)+(1-p)*CC(i+2,j+1))/R;
	end
end

PP=zeros(n+1,n+1)
for j=1:(n+1)
	PP(1+n,j)=psi(SS(1+n,j))
end
for i=n-1:-1:0
    for j=0:i
        PP(i+1,j+1)=(p*PP(i+2,j+2)+(1-p)*PP(i+2,j+1))/R;
	end
end

// QUESTION 2

/* Code fourni par Francine Diener */
/*
function indi=sousL(i,j);
    if SS(i+1,j+1)<=L then indi=1;
        else indi=0;
    end;
endfunction;
DDIC=zeros(n+1,n+1,2);
for j=0:n
    for k=0:1
        if (k==1) then
            DDIC(n+1,j+1,k+1)=phi(SS(n+1,j+1));
        else
            DDIC(n+1,j+1,k+1)=0;
        end;
    end;
end;
for i=n-1:-1:0
    for j=0:i
        for k=0:1
            DDIC(i+1,j+1,k+1)=(p*DDIC(i+1+1,j+1+1,max(k,sousL(i+1,j+1))+1)+(1-p)*DDIC(i+1+1,j+1,max(k,sousL(i+1,j))+1))/R;
        end
    end   
end
*/

/* Une autre manière de coder */
function t=indic(x)
    if x<=L then
        t=1
    else
        t=0
    end
endfunction
DDIC=zeros(n+1,n+1)
for j=0:n
    DDIC(n+1,j+1)=phi(SS(n+1,j+1))*indic(SS(n+1,j+1))
end
for i=n-1:-1:0
    for j=0:i
        if SS(i+1,j+1)<=L then
            DDIC(i+1,j+1)=(p*CC(i+2,j+2)+(1-p)*CC(i+2,j+1))/R;
        else
            DDIC(i+1,j+1)=(p*DDIC(i+2,j+2)+(1-p)*DDIC(i+2,j+1))/R;
        end;
    end
end
disp(DDIC(1,1))


// QUESTION 5

function t=indicplus(x)
    if x>=L then
        t=1
    else
        t=0
    end
endfunction
DDOC=zeros(n+1,n+1)
for j=0:n
    DDOC(n+1,j+1)=phi(SS(n+1,j+1))*indicplus(SS(n+1,j+1))
end
for i=n-1:-1:0
    for j=0:i
        if SS(i+1,j+1)<=L then
            DDOC(i+1,j+1)=0;
        else
            DDOC(i+1,j+1)=(p*DDOC(i+2,j+2)+(1-p)*DDOC(i+2,j+1))/R;
        end;
    end
end
//disp(CC(3,:))
//disp((DDIC+DDOC)(3,:))

// QUESTION 6

DDIP=zeros(n+1,n+1)
for j=0:n
    DDIP(n+1,j+1)=psi(SS(n+1,j+1))*indic(SS(n+1,j+1))
end
for i=n-1:-1:0
    for j=0:i
        if SS(i+1,j+1)<=L then
            DDIP(i+1,j+1)=(p*PP(i+2,j+2)+(1-p)*PP(i+2,j+1))/R;
        else
            DDIP(i+1,j+1)=(p*DDIP(i+2,j+2)+(1-p)*DDIP(i+2,j+1))/R;
        end;
    end
end
disp(DDIP(1,1))

DDOP=zeros(n+1,n+1)
for j=0:n
    DDOP(n+1,j+1)=psi(SS(n+1,j+1))*indicplus(SS(n+1,j+1))
end
for i=n-1:-1:0
    for j=0:i
        if SS(i+1,j+1)<=L then
            DDOP(i+1,j+1)=0;
        else
            DDOP(i+1,j+1)=(p*DDOP(i+2,j+2)+(1-p)*DDOP(i+2,j+1))/R;
        end;
    end
end
disp(PP(3,:))
disp((DDIP+DDOP)(3,:))
