// Determination de la volatilité implicite
// TP 1

function [y]=N(x)
  y=cdfnor("PQ",x,0,1);
endfunction;

function [y]=f(mu,rho)
  y=N(mu/rho+rho/2)-exp(-mu)*N(mu/rho-rho/2);
endfunction;

function [y]=h(mu,rho)
  y=f(mu,rho)-P_obs/S;
endfunction;

function [y]=h_prime(mu,rho)
  y=(1/sqrt(2*%pi))*exp(-(mu/rho+rho/2)**2/2);
endfunction;

// Resolution par Newton
function [sigma]=Implied_vol(t,S,K,T,r,P_obs)
  //////////////////
  // A COMPLETER
  //////////////////
endfunction
  
// Parametres
K=[2650,2700,2750,2800,2850,2900,2950,3000];
S=2872;
t=0;
T=0.04;
r=0.05;
P_obs=[233,183,135,89,50,24,9,3];

// calcul et trace
for i=1:8
  sigma(i)=Implied_vol(t,S,K(i),T,r,P_obs(i));
end

plot2d(K,sigma,rect=[2600,0,3050,0.3]);