function [S,Q1,q,up,down,p0,N1] = IBT123(S0,Ks,Cs,rho) 
% Octave/Matlab function to compute the implied binomial
% tree from spot prices and market European-style Call
% option ask prices, by Rubinstein's "1,2,3" method. 
% INPUTS:                                      (Example)
%   S0= spot stock price (BAC 2021/12/11)       (44.52)
%   Ks= equispaced, increasing strike prices    (40:46) 
%   Cs= Call asks   ([5.6,4.8,4.1,3.45,2.87,2.29,1.81])
%   rho= riskless return to 2022/3/20  (exp(r*109/365))
% OUTPUT:
%   S =  binomial tree of stock prices at all (n,j)
%   Q1 = time T risk-neutral probabilities.
%   q = risk-neutral probs along one path to (n,j)
%   up = upfactors at (n,j)
%   down = downfactors at (n,j)
%   p0 = risk-neutral up probabilities at (n,j)
%   N1 = depth of the pruned tree
% NOTE: since Octave array indices start at 1, output
%   array value for (n,j) is stored at index (n+1,j+1).
% EXAMPLE:
%   S0=44.22; Ks=40:46; % BAC spot price and strikes
%   Cs=[5.6,4.8,4.1,3.45,2.87,2.29,1.81]; % Call asks
%   r=0.06/100; rho=exp(r*109/365); % 13-wk T-bill APR
%   [S,Q1,q,up,down,p0,N1] = IBT123(S0, Ks, Cs, rho)
%
m = length(Ks);        % expect valid Ks(1),...,Ks(m)
% Check that Ks and Cs have the same length
if(length(Cs) != m ) % expect valid Cs(1),...,Cs(m)
  error ("length(Cs)!=length(Ks)");
end
N=m+1;             % initial tree depth, before pruning
delta=Ks(2)-Ks(1); % expect delta=Ks(j)-Ks(j-1), all j
if(!(delta>0))     % expect Ks(1)<Ks(2)<...
  error("must have Ks(1)<Ks(2)")
end
for j=3:m          % expect equispaced Ks
 if( Ks(j)!= Ks(j-1)+delta )
   error("must have equispaced Ks")
 end
end
%% Generate time-T risk-neutral probabilities
Q=1:(N+1);   % so Q(j+1)=Qj for j=0,...,N=m+1
for j=2:m-1    % from Eq.12.24, p.162:
   Q(j+1) = rho*(Cs(j-1)-2*Cs(j)+Cs(j+1))/delta;
end
Q(2)=0; Q(m+1)=0;   % Q1 and Qm are special
Q(N+1) = rho*(Cs(m-1)-Cs(m))/delta; % QN=Qm+1 is special
Q(1)=1-sum(Q(2:(N+1)));  % so sum(Q)=1, as prob.
%% Generate the initial bottom row of stock prices
Sb=1:(N+1);  % so Sb(j+1)=S(N,j) for the initial tree
b=(Cs(2)-Cs(1))/delta;   % temporary "slope" variable
Sb(1)=rho*(S0-Cs(1)+Ks(1)*b)/(1+rho*b); % S(N,0) special
for j=1:m
  Sb(j+1)= Ks(j); % S(N,1),...,S(N,m) are Ks(1),...,Ks(m)
end
Sb(N+1) = rho*Cs(m)/Q(m+2) +Sb(m+1); % S(N,N) is special
%% Prune the bottom row of S to use only positive Qs:
Q1=Q(Q>0);          %  subvector of the positive Q(j)
Sb1 = Sb(Q>0);      % ...and the corresponding S(N,j)
N1 = length(Q1)-1;  % depth of the pruned tree
R = rho^(1/N1);     % new one-step risk-free return
%% Initialize and assign the output binomial trees
% Initialze:
S = zeros(N1+1,N1+1);    % array of stock prices.
q = zeros(N1+1,N1+1);    % array of one-path probs.
p0 = zeros(N1+1,N1+1);   % array of one-step up probs.
up = zeros(N1+1,N1+1);   % up factors
down = zeros(N1+1,N1+1); % down factors
%% Generate the pruned implied binomial trees
% Last row of the pruned stock price tree
for j=0:N1
 S(N1+1,j+1) = Sb1(j+1);
end
% Last row of single-path risk neutral probabilities
for j=0:N1
 q(N1+1,j+1) = Q1(j+1)/nchoosek(N1,j); % path-indep.
end
% Backward induction to fill the trees
for n=N1-1:-1:0
  for j=0:n
   q(n+1,j+1) = q(n+2,j+2)+q(n+2,j+1);  
   p0(n+1,j+1) = q(n+2,j+2)/q(n+1,j+1); 
   S(n+1,j+1) = (p0(n+1,j+1)*S(n+2,j+2)
                  + (1-p0(n+1,j+1))*S(n+2,j+1))/R;
   up(n+1,j+1) = S(n+2,j+2)/S(n+1,j+1);   % upfactor
   down(n+1,j+1) = S(n+2,j+1)/S(n+1,j+1); % downfactor
  end
end
return;
end
