MATLAB Examples

Description

This is a MatLab (R2016a) code which builds a Hull-White Tree with equal time steps, constant rate of reversion, and constant volatility. It is based on Hull-White paper (HW-94) "Numerical Procedures for Implementing Term Structure Models I: Single-Factor Models." The paper can be downloaded from the Hull's web page http://www-2.rotman.utoronto.ca/~hull/downloadablepublications/Numerical%20procedures%20for%20implementing%20Term%20Structure%20Models%20I.pdf In our model the interest rate process r(t) satisfies the following equation $dr(t)=(\theta(t)-a*r(t))*dt+\sigma*dW(t)$, where $W(t)$ is a Wiener process.

Contents

Notations

Node $(i,j)$ is the node on the tree at time $t_{i}$ for which $x(i,j) = (i-j)*dx$. For any time step $i$, the highest node is $x(i,1) = (i-1)*dx$ and the lowest node is $x(i,3+2*(i-2)) = (1-i)*dx$. For n time steps we have $t_{1}=0 < t_{2}=dt < ... < t_{n+1}=n*dt$. The process starts from $0$, i.e. $x(1,1) = (1-1)*dx = 0$.

Defining HW tree parameters

Following HW paper we introduce $x(r(t),t)=r(t)-g(t)$ $dx(t)=-a*x(t)dt+\sigma*dW(t)$ In this step we generate trinomial tree for $x(t)$ and calculate branching probabilities. In order to reproduce numbers from the HW-94 paper (Exhibit 3 on p. 12), we set volatility parameter to $0.01$, mean reversion parameter to $0.1$, maturity to $3$ years, and number of steps to $3$ (i.e. the steps are annual). The user can override these parameter values.

sigma = 0.01;
a = 0.1;
% we are assuming that
% 1 year < maturity < 30 years
maturity = 3;
steps = 3;
dt = maturity/steps;
V = sigma^2*dt;
dx = sigma*sqrt(3*dt);
M = -a*dt;
jmax = ceil(.184/(a*dt));
if jmax < 2
    clear all
    input('The tree is too narrow. Pls, increase number of steps.')
    return
end

Initial Tree

In this step we will build initial tree (i.e. $\theta(t)$). This tree can be viewed in HW-94 in Exhibit 2 on p. 11.

for i = 1:(steps+1)
    nodes = min(3+2*(i-2),1+2*jmax);
    if 3+2*(i+1-2)<=1+2*jmax
        for j = 1:nodes
        x(i,j) = ((nodes+1)/2-j)*dx;
        Pu(i,j) = 1/6 + ((nodes+1)/2-j)^2*M^2/2 + ((nodes+1)/2-j)*M/2;
        Pm(i,j) = 2/3 - ((nodes+1)/2-j)^2*M^2;
        Pd(i,j) = 1/6 + ((nodes+1)/2-j)^2*M^2/2 - ((nodes+1)/2-j)*M/2;
        end
    else
        x(i,1) = ((nodes+1)/2-1)*dx;
        Pu(i,1) = 7/6 + ((nodes+1)/2-1)^2*M^2/2 + 3*((nodes+1)/2-1)*M/2;
        Pm(i,1) = -1/3 - ((nodes+1)/2-1)^2*M^2 - 2*((nodes+1)/2-1)*M;
        Pd(i,1) = 1/6 + ((nodes+1)/2-1)^2*M^2/2 + ((nodes+1)/2-1)*M/2;
        x(i,nodes) = ((nodes+1)/2-nodes)*dx;
        Pu(i,nodes) = 1/6 + ((nodes+1)/2-nodes)^2*M^2/2 - ((nodes+1)/2-nodes)*M/2;
        Pm(i,nodes) = -1/3 - ((nodes+1)/2-nodes)^2*M^2 + 2*((nodes+1)/2-nodes)*M;
        Pd(i,nodes) = 7/6 + ((nodes+1)/2-nodes)^2*M^2/2 - 3*((nodes+1)/2-nodes)*M/2;
        for j = 2:(nodes-1)
        x(i,j) = ((nodes+1)/2-j)*dx;
        Pu(i,j) = 1/6 + ((nodes+1)/2-j)^2*M^2/2 + ((nodes+1)/2-j)*M/2;
        Pm(i,j) = 2/3 - ((nodes+1)/2-j)^2*M^2;
        Pd(i,j) = 1/6 + ((nodes+1)/2-j)^2*M^2/2 - ((nodes+1)/2-j)*M/2;
        end
    end
end
x
dim = min(3+2*(i-2),1+2*jmax);
Pu=Pu(1:steps,1:dim)
Pm=Pm(1:steps,1:dim)
Pd=Pd(1:steps,1:dim)
x =

         0         0         0         0         0
    0.0173         0   -0.0173         0         0
    0.0346    0.0173         0   -0.0173   -0.0346
    0.0346    0.0173         0   -0.0173   -0.0346


Pu =

    0.1667         0         0         0         0
    0.1217    0.1667    0.2217         0         0
    0.8867    0.1217    0.1667    0.2217    0.0867


Pm =

    0.6667         0         0         0         0
    0.6567    0.6667    0.6567         0         0
    0.0267    0.6567    0.6667    0.6567    0.0267


Pd =

    0.1667         0         0         0         0
    0.2217    0.1667    0.1217         0         0
    0.0867    0.2217    0.1667    0.1217    0.8867

Final Tree

In this step we are adjusting the tree by fitting the term structure of the interest rates. The term structure is defined as in HW-94 paper.

for i = 1:steps+1
    t(i) = i*dt;
    L(i) = .08 - .05*exp(-.18*i*dt);
end
% The code can be used for any other term structure. For example,
% the term structure can be pulled from a csv file as shown below.
% Loading spot interest rates using a csv file.
%    load irusd.csv;
%    T = irusd(1,:)/365;
%    L = irusd(2,:);
%    t = [dt:dt:maturity+dt];
%    L = interp1(T, L, t, 'linear', 'extrap');
% We will now calculate root Arrow-Debreu prices following
% using formula (5) on p. 10 in HW-94
% for step 1 (t_{o})
Q(1,1) = 1;
C(1,1) = 1;
g(1) = L(1);
% for step 2 (t_{1})
Q(2,1) = Q(1,1)*Pu(1,1)*exp(-L(1)*dt);
Q(2,2) = Q(1,1)*Pm(1,1)*exp(-L(1)*dt);
Q(2,3) = Q(1,1)*Pd(1,1)*exp(-L(1)*dt);
C(2,1) = exp(-dt*x(2,1));
C(2,2) = exp(-dt*x(2,2));
C(2,3) = exp(-dt*x(2,3));
g(2) = (2*dt*L(2)+log(Q(2,:)*transpose(C(2,:))))/dt;
% for the other steps
% we will split them into 5 categories
for i = 3:(steps+1)
    nodes = min(3+2*(i-2),1+2*jmax);
    if 3+2*(i-2)<=1+2*jmax
        % Q(i,j) on the upper bound
        Q(i,1) = Q(i-1,1)*Pu(i-1,1)*exp(-dt*(x(i-1,1)+g(i-1)));
        C(i,1) = exp(-dt*x(i,1));
        % Q(i,j) one level below the upper bound
        Q(i,2) = Q(i-1,1)*Pm(i-1,1)*exp(-dt*(x(i-1,1)+g(i-1)))+ ...
                 Q(i-1,2)*Pu(i-1,2)*exp(-dt*(x(i-1,2)+g(i-1)));
        C(i,2) = exp(-dt*x(i,2));
        % Q(i,j) on the lower bound
        Q(i,nodes) = Q(i-1,nodes-2)*Pd(i-1,nodes-2)*exp(-dt*(x(i-1,nodes-2)+g(i-1)));
        C(i,nodes) = exp(-dt*x(i,nodes));
        % Q(i,j) one level above the lower bound
        Q(i,nodes-1) = Q(i-1,nodes-2)*Pm(i-1,nodes-2)*exp(-dt*(x(i-1,nodes-2)+g(i-1)))+ ...
                       Q(i-1,nodes-3)*Pd(i-1,nodes-3)*exp(-dt*(x(i-1,nodes-3)+g(i-1)));
        C(i,nodes-1) = exp(-dt*x(i,nodes-1));
        % remaining Q(i,j)
        for j = 3:(nodes-2)
            Q(i,j) = Q(i-1,j-2)*Pd(i-1,j-2)*exp(-dt*(x(i-1,j-2)+g(i-1)))+ ...
                     Q(i-1,j-1)*Pm(i-1,j-1)*exp(-dt*(x(i-1,j-1)+g(i-1)))+ ...
                     Q(i-1,j)*Pu(i-1,j)*exp(-dt*(x(i-1,j)+g(i-1)));
            C(i,j) = exp(-dt*x(i,j));
        end
        g(i) = (i*dt*L(i)+log(Q(i,:)*transpose(C(i,:))))/dt;
    else
        % Q(i,j) on the upper bound
        Q(i,1) = Q(i-1,1)*Pu(i-1,1)*exp(-dt*(x(i-1,1)+g(i-1)))+ ...
                 Q(i-1,2)*Pu(i-1,2)*exp(-dt*(x(i-1,2)+g(i-1)));
        C(i,1) = exp(-dt*x(i,1));
        % Q(i,j) one level below the upper bound
        Q(i,2) = Q(i-1,1)*Pm(i-1,1)*exp(-dt*(x(i-1,1)+g(i-1)))+ ...
                 Q(i-1,2)*Pm(i-1,2)*exp(-dt*(x(i-1,2)+g(i-1)))+ ...
                 Q(i-1,3)*Pu(i-1,3)*exp(-dt*(x(i-1,3)+g(i-1)));
        C(i,2) = exp(-dt*x(i,2));
        % Q(i,j) two levels below the upper bound
        Q(i,3) = Q(i-1,1)*Pd(i-1,1)*exp(-dt*(x(i-1,1)+g(i-1)))+ ...
                 Q(i-1,2)*Pd(i-1,2)*exp(-dt*(x(i-1,2)+g(i-1)))+ ...
                 Q(i-1,3)*Pm(i-1,3)*exp(-dt*(x(i-1,3)+g(i-1)))+ ...
                 Q(i-1,4)*Pu(i-1,4)*exp(-dt*(x(i-1,4)+g(i-1)));
        C(i,3) = exp(-dt*x(i,3));
        % Q(i,j) on the lower bound
        Q(i,nodes) = Q(i-1,nodes)*Pd(i-1,nodes)*exp(-dt*(x(i-1,nodes)+g(i-1)))+ ...
                     Q(i-1,nodes-1)*Pd(i-1,nodes-1)*exp(-dt*(x(i-1,nodes-1)+g(i-1)));
        C(i,nodes) = exp(-dt*x(i,nodes));
        % Q(i,j) one level above the lower bound
        Q(i,nodes-1) = Q(i-1,nodes)*Pm(i-1,nodes)*exp(-dt*(x(i-1,nodes)+g(i-1)))+ ...
                       Q(i-1,nodes-1)*Pm(i-1,nodes-1)*exp(-dt*(x(i-1,nodes-1)+g(i-1)))+ ...
                       Q(i-1,nodes-2)*Pd(i-1,nodes-2)*exp(-dt*(x(i-1,nodes-2)+g(i-1)));
        C(i,nodes-1) = exp(-dt*x(i,nodes-1));
        % Q(i,j) two levels above the lower bound
        Q(i,nodes-2) = Q(i-1,nodes-2)*Pm(i-1,nodes-2)*exp(-dt*(x(i-1,nodes-2)+g(i-1)))+ ...
                       Q(i-1,nodes-3)*Pd(i-1,nodes-3)*exp(-dt*(x(i-1,nodes-3)+g(i-1)))+ ...
                       Q(i-1,nodes-1)*Pu(i-1,nodes-1)*exp(-dt*(x(i-1,nodes-1)+g(i-1)))+ ...
                       Q(i-1,nodes)*Pu(i-1,nodes)*exp(-dt*(x(i-1,nodes)+g(i-1)));
        C(i,nodes-2) = exp(-dt*x(i,nodes-2));
        % remaining Q(i,j)
        if jmax > 2
            for j = 4:(nodes-3)
                Q(i,j) = Q(i-1,j-1)*Pd(i-1,j-1)*exp(-dt*(x(i-1,j-1)+g(i-1)))+ ...
                         Q(i-1,j)*Pm(i-1,j)*exp(-dt*(x(i-1,j)+g(i-1)))+ ...
                         Q(i-1,j+1)*Pu(i-1,j+1)*exp(-dt*(x(i-1,j+1)+g(i-1)));
                C(i,j) = exp(-dt*x(i,j));
            end
            g(i) = (i*dt*L(i)+log(Q(i,:)*transpose(C(i,:))))/dt;
        else
            g(i) = (i*dt*L(i)+log(Q(i,:)*transpose(C(i,:))))/dt;
        end
    end
end
Q
g
r=bsxfun(@plus,x,g')
% Now we plot the tree for x(t)
figure(1);
hold on
for i = 1:steps
    t = [i,i+1];
    nodes = min(3+2*(i-2),1+2*jmax);
    if 3+2*(i-1)<=1+2*jmax
        for j = 1:nodes
            plot(t,[x(i,j),x(i+1,j)])
            plot(t,[x(i,j),x(i+1,j+1)])
            plot(t,[x(i,j),x(i+1,j+2)])
        end
    else
            plot(t,[x(i,1),x(i+1,1)])
            plot(t,[x(i,1),x(i+1,2)])
            plot(t,[x(i,1),x(i+1,3)])
            plot(t,[x(i,nodes),x(i+1,nodes)])
            plot(t,[x(i,nodes),x(i+1,nodes-1)])
            plot(t,[x(i,nodes),x(i+1,nodes-2)])
        for j = 2:(nodes-1)
            plot(t,[x(i,j),x(i+1,j-1)])
            plot(t,[x(i,j),x(i+1,j)])
            plot(t,[x(i,j),x(i+1,j+1)])
        end
    end
    axis tight
    title('Branching Process for x(r,t)');
end
% Now we plot the tree for r(t)
figure(2);
hold on
for i = 1:steps
    t = [i,i+1];
    plot(t,[L(i),L(i+1)],'r')
    plot(t,[g(i),g(i+1)],'g')
    nodes = min(3+2*(i-2),1+2*jmax);
    if 3+2*(i-1)<=1+2*jmax
        for j = 1:nodes
            plot(t,[x(i,j)+g(i),x(i+1,j)+g(i+1)])
            plot(t,[x(i,j)+g(i),x(i+1,j+1)+g(i+1)])
            plot(t,[x(i,j)+g(i),x(i+1,j+2)+g(i+1)])
        end
    else
            plot(t,[x(i,1)+g(i),x(i+1,1)+g(i+1)])
            plot(t,[x(i,1)+g(i),x(i+1,2)+g(i+1)])
            plot(t,[x(i,1)+g(i),x(i+1,3)+g(i+1)])
            plot(t,[x(i,nodes)+g(i),x(i+1,nodes)+g(i+1)])
            plot(t,[x(i,nodes)+g(i),x(i+1,nodes-1)+g(i+1)])
            plot(t,[x(i,nodes)+g(i),x(i+1,nodes-2)+g(i+1)])
        for j = 2:(nodes-1)
            plot(t,[x(i,j)+g(i),x(i+1,j-1)+g(i+1)])
            plot(t,[x(i,j)+g(i),x(i+1,j)+g(i+1)])
            plot(t,[x(i,j)+g(i),x(i+1,j+1)+g(i+1)])
        end
    end
    axis tight;
    title('Branching Process for r(t)');
end
% Now let's check that calculated g_{i}and Q_{ij}
% satisfy equation (6)
for i = 1:steps+1
    test(i) = 0;
    nodes = min(3+2*(i-2),1+2*jmax);
    for j = 1:nodes
        test(i) = test(i)+Q(i,j)*exp(-dt*(x(i,j)+g(i)));
    end
    test(i) = -log(test(i))/(i*dt);
end
t = [dt:dt:(maturity+dt)];
figure(3);
hold on
subplot(3,1,1);
plot(t,L,'--r+')
title('Actual Libor')
subplot(3,1,2);
plot(t,test,'--go')
title('Fitted Libor')
subplot(3,1,3);
axis tight
error = test - L;
plot(t,error,'-')
title('Error')
Q =

    1.0000         0         0         0         0
    0.1604    0.6417    0.1604         0         0
    0.0182    0.1998    0.4736    0.2033    0.0189
    0.0371    0.1957    0.3821    0.2022    0.0399


g =

    0.0382    0.0520    0.0625    0.0688


r =

    0.0382    0.0382    0.0382    0.0382    0.0382
    0.0694    0.0520    0.0347    0.0520    0.0520
    0.0972    0.0799    0.0625    0.0452    0.0279
    0.1034    0.0861    0.0688    0.0514    0.0341