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

// Parametres
t=0;
S=95;
K=100;
T=1;
r=0.1;
P_obs=5.;

// Precision calcul
epsilon=1.e-10;

mu=log(S*exp(r*(T-t))/K);

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;

//////////////////////////
// Iterations de Newton
//////////////////////////

rho=sqrt(2*abs(mu));
count=0;
while (abs(h(mu,rho))>epsilon)
  //////////////////
  // A COMPLETER
  //////////////////
  count=count+1;
end
count
sigma=rho/sqrt(T-t)

//////////////////////////
// Verification: prix du call
//////////////////////////

d1=(1/(sigma*sqrt(T-t)))*(log(S/K)+(r+sigma**2/2)*(T-t));
d2=d1-sigma*sqrt(T-t);
S*N(d1)-K*exp(-r*(T-t))*N(d2)

// $$$ rho=sigma*sqrt(T-t);
// $$$ S*f(mu,rho)

//////////////////////////
// Methode de point fixe
//////////////////////////

rho=sqrt(2*abs(mu));
count=0;
while (abs(h(mu,rho))>epsilon)
  //////////////////
  // A COMPLETER
  //////////////////
  count=count+1;
end
count
sigma=rho/sqrt(T-t)

//////////////////////////
// Methode de bissection
//////////////////////////

rho_inf=0;
rho_sup=sqrt(T-t)*0.5;
rho=(rho_sup+rho_inf)/2;
count=0;
while (abs(h(mu,rho))>epsilon)
  //////////////////
  // A COMPLETER
  //////////////////
  count=count+1;
end
count
sigma=rho/sqrt(T-t)