// Determination de la volatilité implicite
// TP 1
// exec("TP1_smile.sci");

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)
  rho=sqrt(2*abs(mu));
  mu=log(S*exp(r*(T-t))/K);
  while (abs(h(mu,rho))>epsilon)
	rho=rho-h(mu,rho)/h_prime(mu,rho);
  end
  sigma=rho/sqrt(T-t);
endfunction
  
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];

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]);