%METODO IMPLICITO PARA LA RESOLUCION DE LA ECUACION DE BLACK-SCHOLES MODIFICADA (CON LA ESTRUCTURA DE COSTOS DE TRANSACCION DE LELAND)
%ECUACION:  (dV/dt)+1/2*sigma^2*S^2*(d^2V/dS^2)-kapa*sigma*S^2*sqrt(2/(pi*deltat))*|d^2V/dS^2|+r*S*(dV/dS)-r*V=0
%Este modelo contempla tambien el caso sin costos, poniendo kapa=0 (y deltat distinto de 0 por como es el algoritmo)
function V=implicito(sigma,r,kapa,deltat,dt,dS,M,N)
% INPUT: sigma es la volatilidad, r la tasa de interes libre de riesgo, kapa la constante asociada al costo, delta t es el intervalo 
% entre rehedging, dt*M=T (donde T es el stike price), dt es el intervalo de tiempo en la grilla y M es la cantidad de pasos, dS es
% el intervalo de la grilla y N es la cantidad de pasos en esa direccion.
% OUTPUT: V(n) es V(n*dS,0)=V(S,0) 
for n=0:N
    if dS*n>55
        V(n+1)=10;
    elseif dS*n<45
        V(n+1)=0;
    else
        V(n+1)=dS*n-45;
    end
    %V(n+1)=pay_off(n*dS);
    %El Payoff es una funcion de S, hay que definirla
end
%Me construyo la matriz del sistema y el vector solucion Y
A=zeros(N-1,N-1);
Y=zeros(N-1);
if (V(1)-2*V(2)+V(3)>=0)
    a1=(-1/2*sigma^2*1^2+kapa*sigma*1^2*sqrt(2/(pi*deltat))+(r*1)/2)*dt;
    b=(1+(sigma^2*1^2-2*kapa*sigma^2*1^2*sqrt(2/(pi*deltat))+r)*dt);
    c=(-1/2*sigma^2*1^2+kapa*sigma*1^2*sqrt(2/(pi*deltat))-(r*1)/2)*dt;
else
    a1=(-1/2*sigma^2*1^2-kapa*sigma*1^2*sqrt(2/(pi*deltat))+(r*1)/2)*dt; 
    b=(1+(sigma^2*1^2+2*kapa*sigma^2*1^2*sqrt(2/(pi*deltat))+r)*dt);
    c=(-1/2*sigma^2*1^2-kapa*sigma*1^2*sqrt(2/(pi*deltat))-(r*1)/2)*dt;
end
A(1,1)=b;
A(1,2)=c;
Y(1)=V(2)-a1*V(1);
if (V(N-1)-2*V(N)+V(N+1)>=0)
    a=(-1/2*sigma^2*(N-1)^2+kapa*sigma*(N-1)^2*sqrt(2/(pi*deltat))+(r*(N-1))/2)*dt;
    b=(1+(sigma^2*(N-1)^2-2*kapa*sigma^2*(N-1)^2*sqrt(2/(pi*deltat))+r)*dt);
    cN=(-1/2*sigma^2*(N-1)^2+kapa*sigma*(N-1)^2*sqrt(2/(pi*deltat))-(r*(N-1))/2)*dt;
else
    a=(-1/2*sigma^2*(N-1)^2-kapa*sigma*(N-1)^2*sqrt(2/(pi*deltat))+(r*(N-1))/2)*dt;
    b=(1+(sigma^2*1^2+2*kapa*sigma^2*1^2*sqrt(2/(pi*deltat))+r)*dt);
    cN=(-1/2*sigma^2*(N-1)^2-kapa*sigma*(N-1)^2*sqrt(2/(pi*deltat))-(r*(N-1))/2)*dt;
end
A(N-1,N-1)=b;
A(N-1,N-2)=a;
Y(N-1)=V(N)-V(N+1)*c;
for n=2:N-2
    if (V(n)-2*V(n+1)+V(n+2)>=0)
        a=(-1/2*sigma^2*n^2+kapa*sigma*n^2*sqrt(2/(pi*deltat))+(r*n)/2)*dt;
        b=(1+(sigma^2*n^2-2*kapa*sigma^2*n^2*sqrt(2/(pi*deltat))+r)*dt);
        c=(-1/2*sigma^2*n^2+kapa*sigma*n^2*sqrt(2/(pi*deltat))-(r*n)/2)*dt;
            
    else
        a=(-1/2*sigma^2*n^2-kapa*sigma*n^2*sqrt(2/(pi*deltat))+(r*n)/2)*dt;
        b=(1+(sigma^2*n^2+2*kapa*sigma^2*n^2*sqrt(2/(pi*deltat))+r)*dt);
        c=(-1/2*sigma^2*n^2-kapa*sigma*n^2*sqrt(2/(pi*deltat))-(r*n)/2)*dt;
    end
    A(n,n-1)=a;
    A(n,n)=b;
    A(n,n+1)=c;
    Y(n)=V(n+1);
end
%Usa la descomposicion LU, se puede mejorar el algoritmos aprovechando que la matriz es tridiagonal o que es esparsa
[L,U,P]=lu(A);
X=zeros(N-1);
%Itera
for m=M:-1:1
    %Resolucion del sistema AX=Y
    X=L\(P*Y);
    X=U\X;
    V(1)=0;
    V(N+1)=10;
    %V(1)=V_0(0,m*dt);
    %V(N+1)=V_inf(N*dS,m*dt);
    %Estas dos funciones son las condiciones de contorno, hay que definirlas
    for n=2:N
        V(n)=X(n-1);
    end
    Y(1)=V(2)-a1*V(1);
    Y(N-1)=V(N)-V(N+1)*cN;
    for n=2:N-2
        Y(n)=V(n+1);
    end
end

for n=2:N
    V(n)=X(n-1);
end

%Grafica S versus V el precio del activo contra el valor del portfolio
for n=0:N
    S(n+1)=dS*n;
end
plot(S,V) 

%Realizado por Manuel Maurette y Javier Kreiner para Introduccion a Finanzas, dictada el primer cuatrimestre de 2004 por Pablo Amster.