function [S] = IVT(S0, isigma, r, T, N) 
%  Octave/Matlab function to compute the Implied Volatility Tree
%
%  Inputs:
%    S0 =  spot stock price
%    isigma = implied volatility percentage function handle
%    r  = riskless return per time step
%    T = time to expiry
%    N = number of time steps
%
%  Outputs:
%    S = binomial tree
% 
%  NOTE: since Octave array indices start at 1, the output array
%        value for (n,j) is stored at index (n+1,j+1).
%
%
% Example 11.3 from the textbook, p.141:
%
%   S0=100; r=0.05; T=4; N=4;
%   isigma = @(K) ifelse(K<110, 20-K/20, 25.5-K/10);
%   [S] = IVT(S0, isigma, r, T, N)
%
% M.V.Wickerhauser, 2020-11-24
%
%%%
% Allocate matrices for the output and intermediate trees:
S = zeros(N+1,N+1);              % Modeled asset prices
lambda = zeros(N+1,N+1);         % Arrow-Debreu prices
Vput = zeros(N+1,N+1);           % One-step a.t.m. Put premium
% Minor time step
dt = T/N;
% Riskless return per minor time step
R = exp(r*dt);

% Base case
S(1,1) = S0;
[CE,PE]=BlackScholes(1,S0,S0,r,isigma(S0)/100);
lambda(1,1) = 1;
Vput(1,1) = PE/lambda(1,1);
up = (S(1,1)+Vput(1,1))/(S(1,1)/R-Vput(1,1));
down = 1/up;
S(2,2) = up*S(1,1);
S(2,1) = down*S(1,1);
pi_up = (R*S(1,1)-S(2,1))/(S(2,2)-S(2,1));
pi_down = 1-pi_up;
lambda(2,2) = pi_up/R;
lambda(2,1) = pi_down/R;

% Forward induction
for n = 1:N
  K = S(n,n);
  [CEup,PEup]=BlackScholes(n,K,K,r,isigma(K)/100);
  for j = 0:n


  end
end

endfunction
