function pu = RiskNeut(S, R, N)
% Octave/MATLAB function to compute risk neutral up
% probabilities from recombining binomial trees S,R
% of underlying asset prices and riskless returns.
% INPUTS:                                  (Example)
%   S = asset price matrix               (N+1 x N+1)
%   R = riskless returns matrix              (N x N)
%   N = binomial tree height, must be >1         (3)
% OUTPUT:
%   pu =  N row matrix of risk neutral up probs.
% EXAMPLE:
%  S=[100,0,0,0;87,115,0,0;75,100,133,0;65,87,115,152]
%  R=[1.025,0,0;1.03,1.02,0;1.035,1.025,1.01]
%  pu=RiskNeut(S,R,3)
%
  pu = zeros(N,N); % allocate the output matrix
  for(n=0:N-1)
    for(j=0:n)
      up = S(n+2,j+2)/S(n+1,j+1); % up factor
      dn = S(n+2,j+1)/S(n+1,j+1); % down factor
      pu(n+1,j+1) = (R(n+1,j+1)-dn)/(up-dn);
    end
  end
  return
end
