hullwhite.m

% Implementation of Hull-White interest rate tree assuming symmetrical
% branching
% Author : Sandeep Gore
% Date : 21st June 2007 V 1.0

sigma = 0.01;a = 0.1;dt = 1;
% M = -a*dt;
% V = sigma^2*dt;
M = exp(-a*dt)-1;
V = sigma^2*(1-exp(-2*a*dt))/(2*a);
dr = sqrt(3*V);
% dr = 0.0173;
jmax = ceil(-0.184/M); jmin = -jmax;
% jmax = 2; jmin = -jmax;
maxnodes = length(jmax:-1:jmin);
numlevels = 3;
n = 0:numlevels-1;
numnodes = 2*n+1;

% Calculate the no. of nodes at each time step
numnodes(numnodes>maxnodes) = maxnodes;

% The probabilities of upward, straight and downward transition from initial
% time step to the next time step 
prob{1}(1,1) = 1/6;
prob{1}(2,1) = 2/3;
prob{1}(3,1) = 1/6;

R{1} = 0;
nextnodes{1} = (1:3)';

for i = 2:numlevels
    j = floor(numnodes(i)/2);
    jrange = -j:1:j;
    
    R{i} = jrange.*dr;
    nodes = 1:numnodes(i);
%   Probabilities of transition in case of standard nodes
    prob{i}(1, nodes) = 1/6+(1/2)*(jrange.^2*M^2 + jrange*M);    % Up probability
    prob{i}(2, nodes) = 2/3-jrange.^2*M^2;                       % Middle probability
    prob{i}(3, nodes) = fliplr(prob{i}(1,:));                    % Down probability
    
    % Map the nodes in next stage that are connected by nodes in this
    % stage. The elements of each column represent the nodes in next step
    % connected by each of the nodes in current time step, starting from
    % bottom, going up. Hence the first column of nextnodes{i} indicates
    % the nodes connected by first node in current step starting from
    % bottom of the tree to the nodes in next step. 
    
    nextnodes{i} = [nodes; nodes+1; nodes+2];
    
    % Remember the number of nodes in a time step where term structure
    % introduces non-standard branching symmetrically
    if j == jmax
        flatnumnodes = 1:numnodes(i-1);
    end

    % Calculate probabilities at modified(non-standard) branches. 
    % The modified branches at top and bottom edges are mirror images of each other

    if j >= jmax
        prob{i}(1,1) = 1/6+(1/2)*((-j).^2*M^2 - (-j)*M);
        prob{i}(2,[1 end]) = -1/3-(-j)^2*M^2 + 2*(-j)*M;
        prob{i}(3,1) = 7/6 + (1/2)*((-j).^2*M^2 - 3*(-j)*M);
        prob{i}(1,end) = prob{i}(3,1);
        prob{i}(3,end) = prob{i}(1,1);
        
        % Remember the nodes in next step that are connected by nodes in
        % current step in case of non-standard branching.
       
        nextnodes{i} = [flatnumnodes; flatnumnodes+1; flatnumnodes+2];
        temp  = [nextnodes{i}(:, 1) nextnodes{i}(:, end)];
        nextnodes{i} = [temp(:, 1) nextnodes{i} temp(:, 2)];
    end
    
end

% Define term structure and price of bond
% Note : dT and the time spacing in term structure must be same. Else, use
% interpolation
t = 0:5;
termrate = 0.08 - 0.05*exp(-0.18*t);
price = exp(-termrate(2:end).*t(2:end));

% Determine the term structure that matches with the initial term
% structure.

% Initial shift(alpha0) and interest rate. 
shift = zeros(1, numlevels + 1);
shift(1) = termrate(2);                     
R{1} = R{1} + shift(1);
Q = R;
Q{1} = 1;
temp = prob;

for i = 2:numlevels 
    % Calculate value of Q at each node of the term structure
    temp{i-1} = flipud(temp{i-1});
    for j  = 1:numnodes(i)
        % Find out which nodes from previous time step are connected to the
        % node referred by j in current time step. 
        [r,c] = find(nextnodes{i-1} == j);
       
        % Map the connecting nodes from previous time step to corresponding
        % probabilities of connection and calculate the value of Q
        Q{i}(1,j) = sum(temp{i-1}(nextnodes{i-1} == j)'.*exp(-R{i-1}(c)*dt).*Q{i-1}(c));
    end  
    
    shift(i) = (log(cumsum(Q{i}*exp(-R{i}(:)*dt)))-log(price(i)))/dt;
    R{i} = R{i} + shift(i);
end