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 , where is a Wiener process.

## Notations

Node is the node on the tree at time for which . For any time step , the highest node is and the lowest node is . For n time steps we have . The process starts from , i.e. .

## Defining HW tree parameters

Following HW paper we introduce In this step we generate trinomial tree for and calculate branching probabilities. In order to reproduce numbers from the HW-94 paper (Exhibit 3 on p. 12), we set volatility parameter to , mean reversion parameter to , maturity to years, and number of steps to (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. ). 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