%% Inputs
dr        = 0.001;
dt        = 0.02;
rmax      = 0.15;
rmin      = 0;      
theta     = 0.00072836%0.025;

a         = 0.058115%0.023 % 0.05--0.1
sig       = 0.0065852%0.044 %0.015--0.02

T         = 4% inception of swap
t         = [ 4 5 6 7 8 9 10] ; % tenor of the swap 6 years swap
K         = 3.18547 /100% the strike of the swaption;
fswaprate = K;
M         = T/dt;           % number of time intervals
N         = (rmax-rmin)/dr; % number of interest rate intervals
rvec      = rmin:dr:rmax
%% 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))

annuity = sum(discountfactors(T+1:t(end)))
%% specifies the Initial Boundary at Maturity
f=zeros(N+1,M+1)
f(:,1) = max( K-rvec',0)*annuity
f(N+1,2:M+1) = 1

 %% Defining Coefficients
for j = 1 : N
pu(j)  =  -0.25*dt*((sig/dr)*(sig/dr) +  theta/dr         -  a*(j-1));
pd(j)  =  -0.25*dt*((sig/dr)*(sig/dr) -  theta/dr         +  a*(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-a*rmin)/dr

B(1, 1) = -(1 + dt*rmin/2-2)
B(1, 2) =0.25*dt*(theta-a*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 Price of The Swaption
swaptionprice = f(:,end)
swaptionprice = interp1(rvec,f(:,end),K)
mesh(0:dt:T,rvec,f(:,3:end))
title('recieveer Swaption Price under Hull - White Model')
xlabel('Maturity T')
xlim([0 T])
ylabel('Interest Rate r')
ylim([rmin rmax])
legend('alpha = 0.05 sigma = 0.1 Theta = 0.025','location','northwest')

