function D = ZCB(p, R, N)
% Octave/MATLAB function to compute Zero Coupon
% Bond discounts from recombining binomial trees of
% risk neutral probabilities and riskless returns.
% INPUTS:                             (Example)
%   p =  NxN matrix of up probs.      (all 0.5)
%   R =  NxN matrix of riskless Rs  (all 1.002)
%   N = binomial tree height, must be >1    (9)
% OUTPUTS:
%   Z =  (N+1) rows of zero-coupon bond discounts
% EXAMPLE:
%   D = ZCB(0.5*ones(9,9), 1.002*ones(9,9), 9)
%
  D = zeros(N+1,N+1); % initialize output matrix
  for(j=0:N)         % level N
    D(N+1,j+1) = 1;
  end
  for(j=0:N-1)       % level N-1
    D(N,j+1) = 1/R(N,j+1);
  end
  for(n=N-2:-1:0)  % levels N-2,...,1,0
    for(j=0:n)  % Compute D by backward induction
      D(n+1,j+1) = ( p(n+1,j+1)*D(n+2,j+2) +
         (1-p(n+1,j+1))*D(n+2,j+1) )/R(n+1,j+1);
    end
  end
  return;  % D has rows 1:N+1
end
