function [S,Sx,Ipv] = CRRmD(T,S0,Di,ri,v,N) 
% Octave/MATLAB function to model an asset price
% as a risky ex-dividend price Sx plus a riskless
% dividend cash flow present value Ipv, using the
% modified Cox-Ross-Rubinstein (CRR) binomial
% pricing model with time-varying rates ri.
% INPUTS:                               (Example)
%   T  =  expiration time                (1 year)
%   S0 =  spot stock price                  (100)
%   Di =  2:N+2 dividends in a vector (see below)
%   ri =  1:N+1 annual yield curve    (see below)
%   v  =  volatility; must be >0           (0.15)
%   N  =  height of the binomial tree        (12)
% OUTPUTS:
%   S =  cum-dividend binomial tree at all (n,j).
%   Sx =  ex-dividend binomial tree at all (n,j).
%   Ipv =  dividend present value at all n.
% EXAMPLE:
%   Di=[0 1 0 0 1 0 0 2 0 0 2 0 0 2]; % Di(k+1)=D_k
%   ri=[1 1 1 1 1 2 2 2 2 2 3 3 3]/100; % ri(k+1)=r_k
%   [S,Sx,Ipv]=CRRD(1,100,Di, ri, 0.15, 12);
%
  [pui,up,Ri] = CRRparams(T,ri,v,N); % Modified CRR
  Sx=zeros(N+1,N+1); S=zeros(N+1,N+1);
  Ipv=zeros(1,N+1); % ...allocate output matrices
  Ipv(N+1)=Di(N+2)/Ri(N+1); % ignore d_n if n>N+1
  for n=(N-1):(-1):0  % textbook time index n is..
    Ipv(n+1)=(Ipv(n+2)+Di(n+2))/Ri(n+1); % at n+1
  end  % ...backward recursion for Ipv.
  Sx(1,1)=S0-Ipv(1); S(1,1)=S0; % initial values
  for n=1:N     % textbook time indices
    for j=0:n   % states j={0,1,...,n} at time n
      Sx(n+1,j+1)=Sx(1,1)*up^(2*j-n);
      S(n+1,j+1)=Sx(n+1,j+1)+Ipv(n+1);
    end
  end  % ...forward recursion for Sx and S.
  return; % Prices are in S, Sx, and Ipv.
end
