
function F = calibration(x,swaption_matrix)
%% xdata is the swaption_matrix which mainly used for specifying the ni, and nj,
% to represent 1*2, 2*2,3*3,... structure of the swap what we are
% calibrated to,
% for instance, 
% swaption_matrix =...
 %   [31.1  33    33.6  32.5  31.7  28.4  25.7;
 %    33.8  32.6  30.6  29.5  28.8  27    25.5;
 %    30.8  28.8  27.9  27.3  26.9  26    25.3;
 %    28.4  27.5  26.8  26.1  25.5  25.2  25;
 %    27.3  26    25.2  24.7  24.4  24.3  24.6;
 %    24.5  24.6  24.3  23.8  23.3  23.3  23.5;
 %    22.9  22.2  22    21.9  21.8  21.9 22]/100
% ydata is the swaptionprice_matrix which is the swaption price produced by
% the black76 formula from the market volatilty
% for instance,   
    % from Black model %swaptionprice_matrix =...
   % [ .00193674   .00467438    .00789733    .01102371  0.01413741;%   0.01883476   0.02483918;
   %  .00378584    .00795249    .01193239    .01599121  0.02006765 ;%  0.02691976   0.03609521;
   %  .00495465    .00975864    .01463299    .01947449  0.02413264 ;%  0.03260358   0.04431996;
   %  .00586024    .01149027    .01702368    .02215521  0.02695046 ;%  0.03655725   0.05026072;
   %  .00648863   0.01242197   0.01804322   0.02346658  0.02864702 ;]%  0.03894984   0.05427248;
   %  .00685729   0.01357967   0.01987754   0.02544701  0.03084793   0.04205001   0.05770927;
   %  .00692837   0.01331484   0.01971274   0.02565837  0.03157161   0.04267024   0.05817707]
   %% from model
   %swaptionprice_matrixmodel1 =...
    %[0.0092    0.0191    0.0290    0.0386    0.0474    0.0649    0.0881;
    %0.0114    0.0229    0.0342    0.0453    0.0562    0.0757    0.1039;
    %0.0116    0.0231    0.0341    0.0450    0.0547    0.0753    0.1018;
    %0.0120    0.0239    0.0350    0.0459    0.0567    0.0775    0.1045;
    %0.0117    0.0229    0.0335    0.0442    0.0545    0.0742    0.1011;
    %0.0111    0.0219    0.0323    0.0423    0.0520    0.0705    0.0960;
    %0.0102    0.0200    0.0295    0.0387    0.0476    0.0645    0.0878]
   
    
 % x=lsqcurvefit(@calibration1,[0.01 0.01],swaption_matrix,swaptionprice_matrix, [0.01 0.0001],[0.05 0.06])
 %%
 % dt = 0.2 swaption price matrix produced from calibration1
   %F = [ 0.0016    0.0037    0.0065    0.0101    0.0141    0.0202    0.0269
         %0.0024    0.0057    0.0100    0.0150    0.0203    0.0275    0.0350
         %0.0038    0.0088    0.0143    0.0204    0.0265    0.0337    0.0411
         %0.0054    0.0115    0.0180    0.0246    0.0311    0.0377    0.0452
         %0.0065    0.0133    0.0202    0.0271    0.0336    0.0405    0.0476
         %0.0075    0.0149    0.0220    0.0287    0.0358    0.0428    0.0489
         %0.0076    0.0153    0.0228    0.0301    0.0370    0.0431    0.0489]
  %% 
  %dt =0.02 swaption price matrix produced from calibration1
  %swaptionprice_matrixmodel2 =...
    %0.0034    0.0078    0.0133    0.0194    0.0257    0.0038    0.0551
    %0.0041    0.0091    0.0147    0.0206    0.0265    0.0379    0.0538
    %0.0047    0.0101    0.0157    0.0213    0.0268    0.0374    0.0522
    %0.0052    0.0105    0.0159    0.0212    0.0264    0.0361    0.0501
    %0.0052    0.0103    0.0155    0.0205    0.0253    0.0346    0.0476
    %0.0048    0.0095    0.0141    0.0184    0.0228    0.0311    0.0423
    %0.0040    0.0080    0.0118    0.0155    0.0191    0.0257    0.0348
  
  
%%
dr = 0.0001;  dt = 0.02;  rmax = 0.15; rmin = 0;rvec = rmin:dr:rmax; %theta = 0.01; 
N  = (rmax-rmin)/dr; % number of interest rate intervals
sig = x(1); alpha = x(2); theta=x(3)
 %sig = 0.03603; alpha = 0.09426; 
% T = 4;% inception of swap
% t = [4  5  6  7  8  9 10 ] % tenor of the swap 6 years swap
 K = [1.61265  1.8469   2.07475  2.28107  2.44275;%  2.67431 2.86094;
      2.07998  2.30581  2.5043   2.65086  2.76404;%  2.917   3.03967;
      2.53415  2.71882  2.84327  2.93659  3.00746;%  3.09263 3.17047;
      2.90349  2.99783  3.07037  3.12554  3.16202; %  3.19868 3.255;
      3.09217  3.15346  3.19936  3.22638  3.24149]/100%   3.26924 3.30201
      %3.29153  3.29929  3.30014  3.29484  3.31562   3.33412  3.32819
      %3.30123  3.33894  3.35249  3.35969  3.36241   3.34015  3.32047]/100 ;
  % the strike of the swaption;
fswaprate = K
%fswaprate = 0.02864;
%% discount factors
zerorates=...
    [1       0.012458;
    2        0.014123;
    3        0.016354;
    4        0.018689;
    5        0.020749;
    6        0.022442;
    7        0.023830;
    8        0.024962;
    9        0.025863;
   10        0.026581;
   11        0.027144;
   12        0.027713;
   13        0.028181;
   14        0.028584;
   15        0.028930;
   16        0.029176;
   17        0.029388;
   18        0.029576;
   19        0.029746;
   20        0.029893;];
discountfactors = exp(-zerorates(:,1).*zerorates(:,2));

 %% 
 tao = 1; % make the calibration process easier, I make the tao = 1
        % which means the swap payment occurs once a year
        swaptionprice1=zeros(size(swaption_matrix(:,1:end)))
        
 for ni=1:size(swaption_matrix,2)    
     % = [1 2 3 4 ] column of swaptionprice_matirx option maturity (T)
      M  = ni/dt     
      % number of time intervals
    for nj=1:size(swaption_matrix,1) 
     % = [1 2 3 4 ] row of swaptionprice_matrix  swap maturity   (t)
      % annuity factor
      annuity(ni,nj) = sum(discountfactors(ni+1:ni+nj))     
      
%% specifies the Initial Boundary at Maturity
f = zeros(N+1,M+1)
f(:,1) = max(K(ni,nj) - rvec',0 )
f(N+1,2:M+1) =1
% r= max 
for i = 2 : M+1
    f(N,i) = 0; 
end
%% Calculating Theta
%for i = 1: N% length(forwardrates)
%theta(i) = ft(i) + a*forwardrates(i) +0.5*(1-exp(-2*a*dt))*(sig*sig)/a 
%end 
 %% Defining Coefficients
for j = 1 : N
pu(j)  =  -0.25*dt*((sig/dr)*(sig/dr) +  theta/dr         -  alpha*(j-1));   
pd(j)  =  -0.25*dt*((sig/dr)*(sig/dr) -  theta/dr         +  alpha*(j-1));   
pm(j)  =  0.5*dt*(sig/dr)*(sig/dr)    +  0.5*(j-1)*dr*dt  +  1;                 
pmm(j) =  pm(j )-2;
end 
%% Forming Matrix A and B
% form matrix A; B*U(t) = A*U(t+1)
for j = 2 : N+1
A(j, j - 1) = pd(j - 1);
A(j, j) = pm(j - 1);
end;
for j = 2 : N 
A(j, j + 1) = pu(j - 1);
end;
% form matrix A; B*U(t) = A*U(t+1)
for j = 2 : N+1
B(j, j - 1) = -pd(j - 1);
B(j, j) = -pmm(j -1);
end;
for j = 2 : N 
B(j, j + 1) = -pu(j - 1);
end;

%% Boundary Condition
%A(1, 1) = 1 + dt*rmin/2;  A(1, 2) = -0.25*dt*(theta-alpha*rmin)/dr;
%B(1, 1) = 1 + dt*rmin/2;  B(1, 2) = -0.25*dt*(theta-alpha*rmin)/dr;

A(1, 1) =1 + dt*rmin/2;
A(1, 2) =-0.25*dt*(theta-alpha*rmin)/dr;

B(1, 1) = -(1 + dt*rmin/2-2);
B(1, 2) =-0.25*dt*(theta-alpha*rmin)/dr;

%% Solving System of Equations
[L,U] = lu(A)
for i = 1 : 1 : M + 1
f(:, i + 1) = U\(L\(B*f(:, i)));
end;
f=[rvec',f];
%% Interpolating The Pricing of The Swaption
swaptionprice = f(:,end)
% swaptionprice = interp1(rvec,f(:,end),K)
swaptionprice = interp1(rvec,f(:,end),fswaprate(ni,nj));
% the rvec needs to be the forward swap rate which calculates from the curve

swaptionprice1(ni,nj)= swaptionprice*annuity(ni,nj)

  end
     
 end
 F = swaptionprice1
end
 


