IV_calculator

% Implied Volatility Calculator
%            Daijun Xu
%            Project 3
%%
% The function will estimate the market implied volatility of put and call
% options based on Black Scholers model. I got CL and MSFT option data from
% Yahoo!Finance for vectorization test and plot. And you can also use your
% own data(as scalars) to calculate.
%%
function []=IV_calculator()

%% Set the Black Scholers methods to calculate options price
        flag=1; % Default value of option type is call (call when flag=1, put when flag=0)
        d1=@(S,K,r,vol,T)(log(S/K)+(r+vol^2/2)*T)/(vol*sqrt(T));
        d2=@(S,K,r,vol,T)(log(S/K)+(r-vol^2/2)*T)/(vol*sqrt(T));
        bscallpx=@(S,K,r,vol,T)S*cdf('norm',d1(S,K,r,vol,T),0,1)-K*exp(-r*T)*cdf('norm',d2(S,K,r,vol,T),0,1);
        bsputpx=@(S,K,r,vol,T)K*exp(-r*T)*cdf('norm',-d2(S,K,r,vol,T),0,1)-S*cdf('norm',-d1(S,K,r,vol,T),0,1);
        bspx=@(S,K,r,vol,T,flag)flag*bscallpx(S,K,r,vol,T)+(1-flag)*bsputpx(S,K,r,vol,T);
        bsvega=@(S,K,r,vol,T)S*pdf('norm',d1(S,K,r,vol,T),0,1)*sqrt(T);
    quit=0;   % Default value of quitting the program is 0  
   
   
    
%% Get parameters from user and check the validity of parameters
      while(quit==0)

            flag=input('Please choose the type of option \n (call is 1, put is 0):  ');

            T=input('Please input Tenor(in year) of the option:  ');
            K=input('Please input Strike price(in US dollar) of the option:  ');
            r=input('Please input riskless rate(in decimal digits) in the market:  ');
            S=input('Please input current price(in US dollar) of the underlying stock:  ');
            optpx=input('Please input current price(in US dollar) of the option:    ');
            vol_int=input('Please input a guess of volatility(in decimal digits)for newton methods:  ');
            vol=vol_int;
                if T>0&&K>0&&r>-1&&r<1&&S>0&&optpx>0&&vol_int>0&&vol_int<3      % Check the validity of parameters
                    quit=2;
                 fprintf('\nThe results of bisection methods are: \n')
                 [ EstimatedPrice,Interations,ImpliedVolatility]=bisection_iv(optpx,bspx,S,K,r,vol,T,flag)
                 fprintf('\nThe results of newton methods are: \n')
                 [ EstimatedPrice,Interations,ImpliedVolatility]=newton_iv(optpx,bspx,S,K,r,vol,T,flag,bsvega,vol_int)
                 
                else
                    quit=input('Parameters Error. To quit, press 1. To continue inputing parameters,press 0 \n');   
                    if quit==1
                       break
                    end    
 %% To ask if the user want to compute other options
                quit=input('Do you want compute another option? To quit press 1. To continue,press 0 \n');  
                end
      end
end
%% Subfunction bisection_iv

function [ EstimatedPrice,Interations,ImpliedVolatility ] = bisection_iv(optpx,bspx,S,K,r,vol,T,flag)

%   Using bisection methods to find the implied volatility

    Threshold1=0.000001;    %   Set the accuracy
    vol_lb=0;    %   define a lowerbound to guess implied volatility
    vol_ub=2;    %   define a upbound to guess implied volatility

    vol=(vol_ub+vol_lb)./2; % First guess
    i=1;        %   set the counter for times of guess
 
       while(abs(optpx-bspx(S,K,r,vol,T,flag))/optpx>Threshold1) 
       %  Controling the accuracy by the difference between estimated
       %  option price under estimated volatility and real price 
               
         if optpx>bspx(S,K,r,vol,T,flag)
            vol_lb=vol;             % Change the lowerbound to the guess point
            vol=(vol+vol_ub)/2;     % Change the guess point
         else
            vol_ub=vol;             % Change the upperbound to the guess point
            vol=(vol+vol_lb)/2;     % Change the guess point
         end
          i=i+1;

       end
        EstimatedPrice=bspx(S,K,r,vol,T,flag); % Price under estimated volatility
        Interations=i;                         % Interations of guess
        ImpliedVolatility=vol;                 % Estimated volatility
end

function [ EstimatedPrice,Interations,ImpliedVolatility ] = newton_iv(optpx,bspx,S,K,r,vol,T,flag,bsvega,vol_int)

%   Using Newton's methoad to find the implied volatility
  
        Threshold2=0.001;   %   Set the accuracy
        vol=vol_int;        %   First guess
        i=1;                %   Set a counter for times of guess

                while (abs(optpx-bspx(S,K,r,vol,T,flag))/optpx>Threshold2)
                 %  Controling the accuracy by the difference between estimated
                 %  option price under estimated volatility and real price 
                    
                    vol=vol-(bspx(S,K,r,vol,T,flag)-optpx)./bsvega(S,K,r,vol,T);    %   The change of implied vol is made
                    i=i+1;
                end
                           
        EstimatedPrice=bspx(S,K,r,vol,T,flag);
        Interations=i;
        ImpliedVolatility=vol;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%