blsimpv(so,x,r,t,call,maxiter,q,tol)

function v = blsimpv(so,x,r,t,call,maxiter,q,tol) 
%BLSIMPV Black-Scholes implied volatility. 
%       V = BLSIMPV(SO,X,R,T,CALL,MAXITER,Q,TOL) returns the implied
%       volatility of an underlying asset given the current asset price SO,
%       the exercise price X, the risk free interest rate R, the time to
%       maturity in years T, and the call option value call.  This function
%       solves for the implied volatility V, using Newton's method.
%       MAXITER is the maximum number of iterations used in solving for V.  By
%       default, MAXITER = 50.  Q is the dividend rate for dividend-paying 
%       securities; default is 0.  TOL is the tolerance when the result of 
%       the Newton-Raphson method is considered to have converged; default 
%       is 1e-6.
%       
%              Note: This function uses normcdf and normpdf, the normal
%       cumulative distribution and normal probability density functions in
%       the Statistics Toolbox. 
% 
%       For example, an asset has a current price of $100, an 
%       exercise price of $95, the risk free interest rate is 
%       7.5% and the time to maturity of the option is .25 years. 
%       The call option has a value of $10.00.  
% 
%       Using blsimpv(100,95,.075,.25,10) returns an implied volatility 
%       of .31 or 31%. 
% 
%       See also BLSPRICE.
%
%       Reference: Chriss, Black-Scholes and Beyond: Option Pricing Models,
%                  Chapter 4 and 8. 
%
 
%       Author(s): C.F. Garvin, 2-23-95, P. N. Secakusuma, 12-15-2000 
%       Copyright (c) 1995-2000 The MathWorks, Inc. All Rights Reserved. 
%       $Revision: 1.9 $   $Date: 2000/12/15 19:01:13 $ 
 
if nargin < 5 
  error('Missing one of SO, X, R, T, and CALL.') 
end 

if ~exist('maxiter', 'var') | isempty(maxiter),
	maxiter = 50;
end

if ~exist('q', 'var') | isempty(q),
	q = 0;
end
if ~exist('tol', 'var') | isempty(tol),
	tol = 1e-6;
end

% Use the Newton-Raphson Method to find Implied Volatility
n = 1;
sig = 0.50;
c_sig = 0;
while abs(c_sig - call) > tol,
   psig = sig;
	
	c_sig = blsprice(so, x, r, t, psig, q);

    v_sig = blsvega(so, x, r, t, psig, q);
	
	sig = psig - ((c_sig - call) ./ v_sig);
	
	n = n + 1;
	if n > maxiter,
		error('No solution is found.  Increasing MAXITER might produce results.');
	end
end

v = sig;

return


%----- CODE BELOW IS NOT USED!!! -----%
%----- CODE BELOW IS NOT USED!!! -----%
%----- CODE BELOW IS NOT USED!!! -----%
% Comment out old version.....
if 0,

if nargin < 6 
  maxiter = 50;         % Maximum number of iterations 
end 
 
sig = .5*ones(size(so));  % Initialize volatility, delta, and iteration count 
f = 2*ones(size(so)); 
k = 1; 
 
% Using Newton's method to solve for implied volatility 
while 1 
  d1 = (log(so./x) + (r+sig.^2/2).*t)./(sig.*sqrt(t)); 
  d2 = d1 - (sig.*sqrt(t)); 
  d1p = sqrt(t)-(log(so./x)+(r+1/2.*sig.^2).*t)./sig.^2./sqrt(t); 
  d2p = d1p - sqrt(t); 
  % Black-Scholes equation set equal to zero. 
  f = so.*normcdf(d1) - x.*exp(-r.*t).*normcdf(d2)-call; 
  tmp = nan; 
  maxvol = 1000;               % Maximum volality is 1000% 
  litvol = blsprice(so,x,r,t,1e-6); 
  sind = find(call < litvol); 
  f(sind) = tmp(ones(size(sind))); 
  gind = find(call >= so); 
  f(gind) = tmp(ones(size(gind))); 
  if all(abs(f) < 1e-9) 
    v = sig; 
    return 
  end 
  % Derivative of above equation. 
  fp = so.*normpdf(d1).*d1p-x.*exp(-r.*t).*normpdf(d2).*d2p; 
 
  ind = find(abs(fp) < 1e-6); 
  fp(ind) = tmp(ones(size(ind)));  
 
  % Find new delta and repeat if necessary 
  del = -f./fp; 
  sig = sig+del; 
  k = k+1; 
  if k == maxiter+1 
    disp(char(7)) 
    fprintf(['Number of maximum iterations reached.\n',... 
                   'Please increase MAXITER.\n']) 
    sig(ind) = tmp(ones(size(ind))); 
    v = sig; 
    v(gind) = maxvol*ones(size(gind)); 
    return 
  end 
end 
 
v = sig; 
v(gind) = maxvol*ones(size(gind));

end
%--- End of OLD VERSION of BLSIMPV ---%
%--- End of OLD VERSION of BLSIMPV ---%
%--- End of OLD VERSION of BLSIMPV ---%