%function vol_test = vol_test(S_today,Vola,Rate,Dividend,Kmin,Kmax,Knum,Tmax,Tmin,Tnum,num_plots)

function vol_test = vol_test()

S_today = 100;
Vola = 0.232;
Rate = 0.03;
Dividend = 0.01;
Kmin = 80;
Kmax = 130;
Knum = 20;
Tmin = 1/12;
Tmax = 12/12;
Tnum = 30;
num_plots = 3;
num_strikes = 10;

% vol_test(S_today,Vola,Rate,Dividend,Kmin,Kmax,Knum,Tmax,Tmin,Tnum,num_plots)
% This function plot the implied volatility and the local volatility (using DUPIRE) 
% for a set of call price generating from the Black Scholes models with a
% constant volatility. 
%
% DATA
% Underlying_asset_price = S_today;   % Data from the inputs
% Volatility = Vola;                   % Volatility data;
% Risk_free_interest_rate = Rate;     %Risk_free_interest_rate = 0.025;
% Dividend_rate = Dividend;           % Dividend Rate
% Exercise_price = [Kmin .......Kmax] with a number of colums = Knum
% Maturity = [Tmin .......Tmax] with a number of colums = Tnum
% Number_of_plots_T = num_plots 
%
% Internal Data:
% Newton Rapson method variables
% Maximum_number_of_iterations = 1000;
% Tolerance_of_convergence = 1e-6;
% Theta method variables.
% Theta_variable = 1;         % 0 <= Theta <= 1

% DATA
Underlying_asset_price = S_today;   % Data from the inputs
Volatility = Vola;                   % Volatility data;
Risk_free_interest_rate = Rate;     %Risk_free_interest_rate = 0.025;
Dividend_rate = Dividend;           % Dividend Rate

% Internal Data
% Newton Rapson method variables
Maximum_number_of_iterations = 1000;
Tolerance_of_convergence = 1e-6;
% Theta method variables.
Theta_variable = 0.5;       % 0 <= Theta <= 1

% Define K, T and Call_price
% Exercise_price = K
delta_K=(Kmax-Kmin)/(Knum-1);
for j=1:Knum   
Exercise_price(1,j) = Kmin+delta_K*(j-1);
end

% Maturity_in_years % months of the year, i.e. June = 6/12 
Month_of_today = 0/12;       % month of the year, i.e.  January = 1/12
delta_T=(Tmax-Tmin)/(Tnum-1);
for i=1:Tnum
Maturity_in_years(1,i) = Tmin+delta_T*(i-1);
end 

% Call_price
for i=1:Tnum
for j=1:Knum  
Call_price(i,j) = BLSPRICE(S_today,Exercise_price(1,j),Rate,Maturity_in_years(1,i),Vola,Dividend);
end
end

% NOTES IMPORTANTS
%
% A)This function use Black-Scholes put and call pricing
%  [CALL,PUT] = BLSPRICE(SO,K,R,T,VOL,Q) returns the value of
%  call and put options using the Black-Scholes pricing formula.  
%  Note: This function uses normcdf, the normal cumulative distribution function in the Statistics Toolbox.
%    SO is the current asset price, 
%    K is the exercise price, 
%    R is the risk-free interest rate, 
%    T is the time to maturity of the option in years, 
%    VOL is the standard deviation of the annualized continuously compounded rate of return of the asset (also known as volatility)
%    Q is the dividend rate of the asset.  The default Q is 0. 
%  
%  For example, the current price an asset is $100, the exercise price of the option is $95, the risk-free interest rate is 10
%    the time to maturity of the option is .25 years, and the standard deviation of the asset is 50.  Using this data,  
%    call = blsprice(100,95,.1,.25,.5,0) 
%    returns a call option price of $13.70.
% 
% B) Imply volatility
% This function solves for the implied volatility, using Newton's method and Black-Schooles equation.
% MAXITER is the maximum number of iterations used in solving for V.  By
% default, MAXITER = 50.
% By 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. 
%
% C) Local volatility
% This function solves for the local volatility, using dupire equation.
% For solving the PDE I use the theta method. 
% If Theta = 1   is called the explicit method.
% If Theta = 0.5 is called the Crank Nicholson method.

%VARIABLES IN OUR WORKSPACE

% Define variable names in our workspace
S0 = Underlying_asset_price;
K = Exercise_price;
C = Call_price;
T0 = Month_of_today;
T = Maturity_in_years;
R = Risk_free_interest_rate;
Q = Dividend_rate;
NK = Knum;  % = Num_of_strikes 
NT = Tnum;  % = Num_of_maturities
VOL = Volatility; % Volatility_from Data
warning off;

% Newton Rapson method variables in our workspace
MAXITER = Maximum_number_of_iterations;
TOL = Tolerance_of_convergence;

% Theta method variables in our workspace 
Theta = Theta_variable;

%CALCULATIONS
% Calculate Implied volatility with MATLAB
for i=1:NT
for j=1:NK
    if C(i,j)==0
    Imp_vol(i,j) = 0;
    else
    Imp_vol(i,j)=BLSIMPV(S0,K(1,j),R,T(i)-T0,C(i,j),MAXITER,Q,TOL);
    end
end
end

Imp_vol;

% Calculate de Local Volatility using dupire equation in terms of Strikes and Maturity

% define delta for maturities
for i=1:NT
     if i==NT
     delta_T(i)= 0; 
     else
     delta_T(i)=T(i+1)-T(i);
     end
end
% define delta for strikes
for j=1:NK
    if j==NK
    delta_K(j)= 0 ;
    else  
    delta_K(j)=K(j+1)-K(j);
    end
end
% parcial derivative of Call / respect MATURITY
for i=1:NT
for j=1:NK
    if i==NT | i==1    
    dp_C_r_T(i,j)=0;
    else    
    dp_C_r_T(i,j)=[C(i+1,j)-C(i-1,j)]/[2*delta_T(i)];
    end
end
end
% parcial derivative of Call / respect STRIKES
for i=1:NT
for j=1:NK
    if  i==1 | i==NT | j==1 | j==NK
    dp_C_r_K(i,j)=0;
    else  
    dp_C_r_K(i,j)=[Theta*(-C(i,j-1)+C(i,j+1))+(1-Theta)*(-C(i+1,j-1)+C(i+1,j+1))]/[2*delta_K(j)];
    end
end
end
% second parcial derivative of Call / respect STRIKES
for i=1:NT
for j=1:NK
    if i==1 | i==NT | j==1 | j==NK      
    sec_dp_C_r_K(i,j)=0;
    else
    sec_dp_C_r_K(i,j)=[Theta*(C(i,j-1)-2*C(i,j)+C(i,j+1))+(1-Theta)*(C(i+1,j-1)-2*C(i+1,j)+C(i+1,j+1))]/[delta_K(j)*delta_K(j)];
    end
end
end

%Using Dupire Equation in terms of Strikes and Maturity
% Loc_vol1 = dupire formula
% Loc_vol2 = Using Dupire Equation in terms call derivatives;
for i=1:NT
for j=1:NK
    if i==1 | i==NT | j==1 | j==NK
    Loc_vol1(i,j)=0;
    Loc_vol2(i,j)=0;
    else
    Loc_vol1(i,j)=sqrt([dp_C_r_T(i,j)]/[0.5*K(1,j)*K(1,j)*sec_dp_C_r_K(i,j)]);
    Loc_vol2(i,j)=sqrt([dp_C_r_T(i,j)+((R-Q)*K(1,j)*dp_C_r_K(i,j))+Q*C(i,j)]/[0.5*K(1,j)*K(1,j)*sec_dp_C_r_K(i,j)]);
    end
end
end
delta_K;
delta_T;
dp_C_r_T;
dp_C_r_K;
sec_dp_C_r_K;
Loc_vol1;
Loc_vol2;

% Calculate de Local Volatility using dupire equation in terms of imply volatility
% parcial derivative of imply volatility / respect MATURITY
for i=1:NT
for j=1:NK
    if i==1 | i==NT
    dp_IV_r_T(i,j)=0;
    else    
    dp_IV_r_T(i,j)=[Imp_vol(i+1,j)-Imp_vol(i-1,j)]/[2*delta_T(i)];
    end
end
end
% parcial derivative of imply volatility / respect STRIKES
for i=1:NT
for j=1:NK   
    if i==1 | i==NT | j==1 | j==NK      
    dp_IV_r_K(i,j)=0;  
    else   
    dp_IV_r_K(i,j)=[Theta*(-Imp_vol(i,j-1)+Imp_vol(i,j+1))+(1-Theta)*(-Imp_vol(i+1,j-1)+Imp_vol(i+1,j+1))]/[2*delta_K(j)]; 
    end
end
end
% second parcial derivative of imply volatility / respect STRIKES
for i=1:NT
for j=1:NK 
    if i==1 | i==NT | j==1 | j==NK     
    sec_dp_IV_r_K(i,j)=0; 
    else
    sec_dp_IV_r_K(i,j)=[Theta*(Imp_vol(i,j-1)-2*Imp_vol(i,j)+Imp_vol(i,j+1))+(1-Theta)*(Imp_vol(i+1,j-1)-2*Imp_vol(i+1,j)+Imp_vol(i+1,j+1))]/[delta_K(j)*delta_K(j)];
    end
end
end
%Using Dupire Equation in terms of implied volatility
% Loc_vol3 = using Dupire Equation in terms of implied volatility;
for i=1:NT
for j=1:NK
    if i==1 | i==NT | j==1 | j==NK
    Loc_vol3(i,j)=0;
    else
    d1(i,j)=[log(S0/K(1,j))+(R-Q+0.5*Imp_vol(i,j)^2*(T(i)-T0))]/[Imp_vol(i,j)*(sqrt(T(i)-T0))];   
    d2(i,j)=[1+K(1,j)*d1(i,j)*dp_IV_r_K(i,j)*sqrt(T(i)-T0)]^2+[Imp_vol(i,j)^2*K(1,j)^2*(T(i)-T0)*(sec_dp_IV_r_K(i,j)-d1(i,j)*dp_IV_r_K(i,j)^2*(sqrt(T(i)-T0)))];
    Loc_vol3(i,j)=sqrt([Imp_vol(i,j)^2+2*Imp_vol(i,j)*(T(i)-T0)*(dp_IV_r_T(i,j)+(R-Q)*K(1,j)*dp_IV_r_K(i,j))]/d2(i,j));
    end
end
end

dp_IV_r_T;
dp_IV_r_K;
sec_dp_IV_r_K;
Loc_vol3;

% Diference betwwen Local and Impled Volatilit
for i=1:NT
for j=1:NK
    if i==1 | i==NT | j==1 | j==NK
    Dif_imp(i,j)=0;
    Dif_vol1(i,j)=0;
    Dif_vol2(i,j)=0;
    Dif_vol3(i,j)=0;
    else
    Dif_imp(i,j)=VOL-Imp_vol(i,j);
    Dif_vol1(i,j)=VOL-Loc_vol1(i,j);
    Dif_vol2(i,j)=VOL-Loc_vol2(i,j);
    Dif_vol3(i,j)=VOL-Loc_vol3(i,j);
    end
end
end

% Plots and worksapce results
% Loc_vol1 = dupire formula
% Loc_vol2 = Using Dupire Equation in terms call derivatives;
% Loc_vol3 = Using Dupire Equation in terms of implied volatility;
% Dif_vol1 = Difference between volatility given and volatility 1
% Dif_vol2 = Difference between volatility given and volatility 2
% Dif_vol3 = Difference between volatility given and volatility 3

if num_plots > NT-1  
num_plots=NT-1;
end

for i=1:num_plots
    
    if num_plots == 1       deltaplot(1)=NT-1; 
    elseif num_plots == 2      
        deltaplot = [NT/2,NT-1];
    else
        if i==1                 deltaplot(1)=round(NT/num_plots);
        elseif i==num_plots     deltaplot(num_plots)=NT-1;
        else
        deltaplot(i)=round(deltaplot(i-1)+(NT-1)/(num_plots));
        end
    end
end

if num_strikes > NK-1  
    num_strikes=NK-1;
end
for j=1:num_strikes
   
    if num_strikes == 1       deltastrike(1)=NK-1; 
    elseif num_strikes == 2      
        deltastrike = [NK/2,NK-1];
    else
        if j==1                   deltastrike(1)=round(NK/num_strikes);
        elseif j==num_strikes     deltastrike(num_strikes)=NK-1;
        else
        deltastrike(j)=round(deltastrike(j-1)+(NK-1)/(num_strikes));
        end
    end
end

%define worksapce

for i=1:num_plots
for j=1:num_strikes
        Maturity_in_years(i) = T(deltaplot(i)); 
        Exercise_price(j) = K(deltastrike(j));
        Call_price(i,j) = C(deltaplot(i),deltastrike(j));
        Imp_vol_X(i,j) = Imp_vol(deltaplot(i),deltastrike(j));
        Loc_vol_formula(i,j) = Loc_vol1(deltaplot(i),deltastrike(j));
        Loc_vol_in_call_terms(i,j) = Loc_vol2(deltaplot(i),deltastrike(j));
        Loc_vol_in_imp_vol_terms(i,j) = Loc_vol3(deltaplot(i),deltastrike(j));
        Dif_imp_X(i,j) = Dif_imp(deltaplot(i),deltastrike(j));
        Dif_vol_formula(i,j) = Dif_vol1(deltaplot(i),deltastrike(j));
        Dif_vol_in_call_terms(i,j) = Dif_vol2(deltaplot(i),deltastrike(j));
        Dif_vol_in_imp_vol_terms(i,j) = Dif_vol3(deltaplot(i),deltastrike(j));
end
end

%worksapce display

Exercise_price;
Maturity_in_years;
Call_price;
Imp_vol;
Loc_vol1;
Loc_vol2;
Loc_vol3;
Dif_imp;
Dif_vol1;
Dif_vol2;
Dif_vol3;

deltaplot;
deltastrike;
Exercise_price;
Maturity_in_years;
Call_price;
Imp_vol_X;
Loc_vol_formula;
Loc_vol_in_call_terms;
Loc_vol_in_imp_vol_terms;
Dif_imp_X
Dif_vol_formula;
Dif_vol_in_call_terms
Dif_vol_in_imp_vol_terms

close all
warning off

figure
for i=1:num_plots
    subplot(num_plots,1,i)
plot(K,Dif_vol2(deltaplot(i),:),K,Dif_vol3(deltaplot(i),:));
%axis([K(1,2),K(1,NK-1),-1/(NK*10),1/(NK*10)]);
YLABEL(+ Maturity_in_years(deltaplot(i)));
GRID ON
if i==1    
TITLE('Error = volatility given minus Local Volatility                    ');
legend('Local Volatility in terms of calls','Local Volatility in terms of implied volatility');
end
if i==num_plots    
XLABEL('Maturity vs Strikes');
end
end
figure
for i=1:num_plots
    subplot(num_plots,1,i)
plot(K,Imp_vol(deltaplot(i),:),K,Loc_vol1(deltaplot(i),:),K,Loc_vol2(deltaplot(i),:),K,Loc_vol3(deltaplot(i),:));
axis([K(1,2),K(1,NK-1),VOL-1/NK,VOL+1/NK]);
YLABEL(+ Maturity_in_years(deltaplot(i)));
GRID ON
if i==1 
TITLE('Implied volatility and Local Volatilities                   '); 
legend('Implied volatility','Local vol(dupire formula)','Local vol(call prices)','Local vol(implied volatility)');
end
if i==num_plots    
XLABEL('Maturity vs Strikes');
end
end

% Plot 3D

for i=2:NT-1
for j=2:NK-1
        T_plot(i-1) = real(T(i)); 
        K_plot(j-1) = real(K(j));
        Call_price_plot(i-1,j-1) = real(C(i,j));
        Imp_vol_plot(i-1,j-1) = real(Imp_vol(i,j));
        Loc_vol1_plot(i-1,j-1) = real(Loc_vol1(i,j));
        Loc_vol2_plot(i-1,j-1) = real(Loc_vol2(i,j));
        Loc_vol3_plot(i-1,j-1) = real(Loc_vol3(i,j));
        Dif_imp_plot(i-1,j-1) = real(Dif_imp(i,j));
        Dif_vol1_plot(i-1,j-1) = real(Dif_vol1(i,j));
        Dif_vol2_plot(i-1,j-1) = real(Dif_vol2(i,j));
        Dif_vol3_plot(i-1,j-1) = real(Dif_vol3(i,j));
end
end

figure
subplot(3,2,2);
surf(K,T,Dif_imp);
TITLE('                            0.232 - Imp vol');
ZLABEL(' Error ');
colormap cool;
axis([Kmin,Kmax,Tmin,Tmax,-2*10^-14,2*10^-14]);
subplot(3,2,4);
surf(K_plot,T_plot,Dif_vol1_plot);
TITLE('            0.232 - Loc vol (Call prices)');
ZLABEL(' Error ')
colormap cool;
axis([Kmin,Kmax,Tmin,Tmax,-0.2,0.2]);
subplot(3,2,6); 
surf(K_plot,T_plot,Dif_vol3_plot);
TITLE('                          0.232 - Loc vol (Imp)');
XLABEL('Strikes');YLABEL('Maturity');ZLABEL(' Error ');
colormap cool;
axis([Kmin,Kmax,Tmin,Tmax,-2*10^-13,2*10^-13]);

subplot(3,2,1);
surf(K,T,Imp_vol);
TITLE('Implied volatility');
colormap cool; axis([Kmin,Kmax,Tmin,Tmax,VOL-0.001,VOL+0.001]);
subplot(3,2,3);
surf(K_plot,T_plot,Loc_vol1_plot);
TITLE('Local volatility (Call prices)');ZLABEL('Volatility')
colormap cool; axis([Kmin,Kmax,Tmin,Tmax,0.1,0.3]);
subplot(3,2,5);
surf(K_plot,T_plot,Loc_vol3_plot);
TITLE('Local volatility (Implied volatility)');
XLABEL('Strikes');YLABEL('Maturity');
colormap cool;axis([Kmin,Kmax,Tmin,Tmax,VOL-0.001,VOL+0.001]);

K;
T;
Call_price;
Imp_vol;
Loc_vol1;
Loc_vol2;
Loc_vol3;
Dif_imp;
Dif_vol1;
Dif_vol2;
Dif_vol3;

T_plot;
K_plot;
Call_price_plot;
Imp_vol_plot;
Loc_vol1_plot;
Loc_vol2_plot;
Loc_vol3_plot;
Dif_imp_plot;
Dif_vol1_plot;
Dif_vol2_plot;
Dif_vol3_plot;