function [ITTree, OptionData] = itttree(StockSpec, RateSpec, ITTTimeSpec, StockOptSpec)
%ITTTREE Build an Implied Trinomial stock Tree.
%
% [ITTTREE] =  itttree(StockSpec, RateSpec, TimeSpec, StockOptSpec)
%                                  
% Inputs:
%   StockSpec     - Stock specification. Type "help stockspec" for 
%                   information on the stock specification. 
%   RateSpec      - Interest rate specification of the initial risk free rate curve.
%                   Type "help intenvset" for information on declaring an 
%                   interest rate variable. 
%   TimeSpec      - Tree time layout specification.  Type "help itttimespec"
%                   for information on the tree structure.  Defines the 
%                   observation dates of the implied trinomial tree.
%   StockOptSpec  - Option stock specification. Type "help stockoptspec" for 
%                   information on the option stock specification. 
%
% Outputs:
%   ITTTree  - Structure specifying stock and time information for an implied 
%              trinomial tree. 
%
% See also STOCKSPEC, ITTTIMESPEC, STOCKOPTSPEC, INTENVSET.

%   Author(s): M. Reyes-Kattar 15-March-2006
%   Copyright 1998-2006 The MathWorks, Inc.

%----------------------------------------------------------------
% Checking input arguments
%----------------------------------------------------------------
if nargin < 3
    error(message('fininst:itttree:InvalidInputs'));
end

if ~isafin(StockSpec,'StockSpec')
    error(message('fininst:itttree:InvalidStockSpec'));
end

if ~isafin(RateSpec,'RateSpec')
    error(message('fininst:itttree:InvalidRateSpec'));
end

if ~isafin(ITTTimeSpec,'ITTTimeSpec')
    error(message('fininst:itttree:InvalidTimeSpec'));
end

if nargin > 3 && ~isafin(StockOptSpec,'StockOptSpec')
    error(message('fininst:itttree:InvalidStockOptSpec'));
end

% Make sure the tree dates and the option dates are consistent
if(nargin > 3 && ITTTimeSpec.ValuationDate ~= StockOptSpec.Settle(1))
    error(message('fininst:itttree:InconsistentValuationSettleDates'));
end

% Set aside original RateSpec:
RateSpecOri = RateSpec;


% Stock Info
So          = StockSpec.AssetPrice;
Sigma0      = StockSpec.Sigma;
ExDates     = StockSpec.ExDividendDates;
DivAmounts   = StockSpec.DividendAmounts;
DivType     = StockSpec.DividendType;


if(length(So)>1)
    warning(message('fininst:binstocktree:MoreThanOneStock'));    
    So          = So (1);
    Sigma0      = Sigma0(1);
    DivType     = DivType{1};
    ExDates     = ExDates(1,:);
    DivAmounts  = DivAmounts(1,:);
end

if(iscell(DivType))
    DivType = DivType{1};
end

% Tree Info
ValuationDate   = ITTTimeSpec.ValuationDate;
Maturity        = ITTTimeSpec.Maturity;
NumPeriods      = ITTTimeSpec.NumPeriods;
Basis           = ITTTimeSpec.Basis;
EndMonthRule    = ITTTimeSpec.EndMonthRule;
TreeTimes       = ITTTimeSpec.tObs';
TreeDates       = ITTTimeSpec.dObs';

% Make sure we consider only dividend payment date that are to the right of
% the valuation date
if(~isempty(ExDates) && ~all(isnan(ExDates)))
    % Truncate payments before valuationDate
    TruncateIdx = ExDates<ValuationDate;

    if any(TruncateIdx)
        warning(message('fininst:itttree:IgnoreDividendsBeforeValuationDate'))
    end

    ExDates(TruncateIdx) = [];
    DivAmounts(TruncateIdx) = [];

    % Truncate payments after tree maturity
    TruncateIdx = ExDates>Maturity;

    if any(TruncateIdx)
        warning(message('fininst:itttree:DivAfterMaturity'))
    end

    ExDates(TruncateIdx) = [];
    DivAmounts(TruncateIdx) = [];

    % If there is anything left, find Times equivalents
    if(~isempty(ExDates))
        ExTimes = date2time(ValuationDate, ExDates, -1, Basis, EndMonthRule);

    else
        % Throw a warning and behave as if no dividends were specified
        warning(message('fininst:itttree:AllDividendsBeforeValuationDate'))
        DivAmounts = 0;
        DivType = [];
    end
end

% Define tree and time variables
NumLevels = NumPeriods + 1;
TMaturity = date2time(ValuationDate, Maturity, -1, Basis, EndMonthRule);

% Those ExDates that coincide with treedates have their
% corresponding times snap with the TreeTimes:
[MemberMask, Loc] = ismember(ExDates, TreeDates);
if any(MemberMask)
    ExTimes(MemberMask) = TreeTimes(Loc(Loc ~= 0));
end

% If the ValuationDate of the tree is between two maturities of the
% RateSpec, we need to truncate the information in the RateSpec so that we
% can set the ValuationDate of the RateSpec equal to the valuationDate of
% the tree:
RSValuationDate = intenvget(RateSpec, 'ValuationDate');
if(RSValuationDate < ValuationDate)
    RSStartDate = intenvget(RateSpec, 'StartDate');
    RSEndDate = intenvget(RateSpec, 'EndDate');
    RSRates = intenvget(RateSpec, 'Rates');

    % Those segments that are completely to the left of the ValuationDate
    % can be eliminated (emptied out)
    outIdx = find(RSStartDate < ValuationDate & RSEndDate <= ValuationDate);
    RSStartDate(outIdx) = [];
    RSEndDate(outIdx) = [];
    RSRates(outIdx)=[];

    % Segments that start before ValuatinDate but that end after
    % ValuationDate can be truncated to start on ValuationDate.
    RSStartDate(RSStartDate < ValuationDate) = ValuationDate;
    RateSpec = intenvset(RateSpec, 'ValuationDate', ValuationDate, 'StartDate', RSStartDate, 'EndDate', RSEndDate, 'Rates', RSRates);
end

% Throw a warning if RateSpec is not continuously compounded in case the
% users are not aware that the rates will change when compounding is set to
% -1
Compound  = intenvget(RateSpec, 'Compounding');
if Compound ~= -1
    warning(message('fininst:itttree:NonContinuousCompounding'));
end

% Obtain the risk-free rates based on the given risk-free curve
RateSpec = intenvset(RateSpec, 'Basis', Basis, 'EndMonthRule', EndMonthRule, ...
    'Compounding', -1, 'ValuationDate', ValuationDate);
RateSpec = intenvset(RateSpec, 'StartTimes', TreeTimes(1:end-1), 'EndTimes', TreeTimes(2:end));
RiskFreeRates = intenvget(RateSpec, 'Rates');
RiskFreeDisc = intenvget(RateSpec, 'Disc');

% -------------------------------------------------------------
% Calculate model parameters
% -------------------------------------------------------------

% Delta Time
dT = TMaturity/NumPeriods;

% up and down factors
u = exp(Sigma0*sqrt(2*dT));
u = u * ones(NumPeriods,1);

% -------------------------------------------------------------
% Define outputs
% -------------------------------------------------------------

% Build actual tree
STree =  mktrintree(NumLevels  ,1,1:2:2*NumLevels-1);
PTree =  mktrintree(NumLevels-1,3,1:2:2*NumLevels-3);

if strcmpi(DivType, 'continuous')
    
    BRateSpec = intenvset(RateSpec, 'Rates', RiskFreeRates-DivAmounts);
    BDisc = intenvget(BRateSpec, 'Disc');

    % The continuous model affects only the calculation of forward prices. 
    % For the rest of the algorithm, behave as if no dividends were specified:
    DivAmounts = 0;
    DivType = [];
    
else
    BDisc = RiskFreeDisc;
end


% Populate the tree
if isempty(DivType) || strcmpi(DivType, 'constant')
    % Calculate Dividend Factors
    if isempty(DivType)
        DivFactors = ones(1, NumLevels);
    else
        % DivFactors = (1 - Delta(i))
        DivFactors = getdivfactors(TreeTimes, ExTimes, DivAmounts);
    end

    STree{1} = So;
    for iLevel = 2:NumLevels

        limit=floor((2*iLevel-1)/2);
        idx = limit:-1:-limit;
        STree{iLevel}(:) = DivFactors(iLevel) * STree{1} * (u(iLevel-1) .^ idx);
    end
else

    % Find the present value of all future cash flows at each
    % tree node time where there is a dividend payment
    CashFlows = TreeCashFlows(RateSpec, ExTimes, DivAmounts, TreeTimes);

    STree{1} = So - CashFlows(1);
    for iLevel = 2:NumLevels

        limit=floor((2*iLevel-1)/2);
        idx = limit:-1:-limit;
        STree{iLevel}(:) = STree{1} * (u(iLevel-1) .^ idx);
    end

    % Add the present value of future dividends to each node:
    for iLevel = 2:NumLevels

        STree{iLevel}(:) = STree{iLevel}(:) + CashFlows(iLevel);

        if any(STree{iLevel}(:) < 0)
           error(message('fininst:itttree:InvalidDividends', sprintf( '%f', TreeTimes( iLevel ) )));
        end
    end

    % Finally, replace value at node zero
    STree{1} = So;
end

% -------------------------------------------------------------------------
% If requested, return the ideal list of options that we need to build 
% the tree without having to interpolate/extrapolate along the volatility
% surface
if nargout>1
    OptionData = StockOptRequirements(ITTTimeSpec, STree);
end


% if StockOptSpec was not specified, we are done
if(nargin<4)
    ITTree=[];
    return;
end

% Extract interpolation preference
InterpMethod = StockOptSpec.InterpMethod;

% -------------------------------------------------------------------------
% Obtain the implied volatility from the stock and option data. This yields
% data point for the volatility smile (do this only if we are using vol
% interpolation):
if strcmp(InterpMethod, 'vol')
    [Sigma_i] = GetImpliedVolatility(StockSpec, RateSpecOri, StockOptSpec);
else
    Sigma_i = [];
end

%--------------------------------------------------------------------------
% Now build the AD Tree and the probability tree. To do this we will need
% the option prices at the tree nodes. We obtain those from interpolating
% from the input data
% -------------------------------------------------------------------------
[CallPrices, PutPrices] = InterpOptPrices(STree, StockSpec, RateSpecOri, ...
    ITTTimeSpec, StockOptSpec, TreeDates, TreeTimes, Sigma_i, InterpMethod);

% Initialize A-D Tree
Q =  mktrintree(NumLevels,1,1:2:2*NumLevels-1);
Q{1} = 1;

jMax = 2*(1:NumLevels)-1;
jMid = ceil(jMax/2);

for iLevel=1:NumLevels-1
    Fi=STree{iLevel}/BDisc(iLevel);
    for jState = jMax(iLevel):-1:1
        if(jState >= jMid(iLevel))
            % Below and including the central node
            kIdx=jState+1:jMax(iLevel);
            Sum=sum(Q{iLevel}(kIdx) .* (STree{iLevel+1}(jState+1) - STree{iLevel}(kIdx)/BDisc(iLevel)));

            PTree{iLevel}(3,jState) = (PutPrices{iLevel+1}(jState+1)/RiskFreeDisc(iLevel) - Sum)/(Q{iLevel}(jState)*(STree{iLevel+1}(jState+1)-STree{iLevel+1}(jState+2)));
            PTree{iLevel}(1,jState) = (Fi(jState) + PTree{iLevel}(3,jState)*(STree{iLevel+1}(jState+1)-STree{iLevel+1}(jState+2))-STree{iLevel+1}(jState+1)) / ...
                (STree{iLevel+1}(jState) - STree{iLevel+1}(jState+1));
        else
            % Above and excluding the central node
            kIdx = 1:jState-1;
            Sum=sum(Q{iLevel}(kIdx) .* (STree{iLevel}(kIdx)/BDisc(iLevel) - STree{iLevel+1}(jState+1)));

            PTree{iLevel}(1,jState) = (CallPrices{iLevel+1}(jState+1)/RiskFreeDisc(iLevel) - Sum)/(Q{iLevel}(jState)*(STree{iLevel+1}(jState)-STree{iLevel+1}(jState+1)));
            PTree{iLevel}(3,jState) = (Fi(jState) - PTree{iLevel}(1,jState)*(STree{iLevel+1}(jState)-STree{iLevel+1}(jState+1))-STree{iLevel+1}(jState+1)) / ...
                (STree{iLevel+1}(jState+2)-STree{iLevel+1}(jState+1));
        end
        
        % Calculate middle node
        PTree{iLevel}(2,jState) = 1 - PTree{iLevel}(1,jState) - PTree{iLevel}(3,jState);
        
        % Correct faulty probabilities
        if(any(PTree{iLevel}(:,jState)<0) || any(PTree{iLevel}(:,jState)>1))

            if(Fi(jState)>STree{iLevel+1}(jState+1) && Fi(jState) < STree{iLevel+1}(jState))

                PTree{iLevel}(1,jState) = 1/2 * ((Fi(jState)-STree{iLevel+1}(jState+1))/(STree{iLevel+1}(jState)-STree{iLevel+1}(jState+1)) + ...
                    (Fi(jState)-STree{iLevel+1}(jState+2))/(STree{iLevel+1}(jState)-STree{iLevel+1}(jState+2)));
                PTree{iLevel}(3,jState) = 1/2 * (STree{iLevel+1}(jState)-Fi(jState))/(STree{iLevel+1}(jState)-STree{iLevel+1}(jState+2));

            elseif(Fi(jState)>STree{iLevel+1}(jState+2) && Fi(jState)<STree{iLevel+1}(jState+1))

                PTree{iLevel}(1,jState) = 1/2*(Fi(jState)-STree{iLevel+1}(jState+2)) / (STree{iLevel+1}(jState)-STree{iLevel+1}(jState+2));
                PTree{iLevel}(3,jState) = 1/2*((STree{iLevel+1}(jState)-Fi(jState))/(STree{iLevel+1}(jState)-STree{iLevel+1}(jState+2)) + ...
                    (STree{iLevel+1}(jState+1)-Fi(jState))/(STree{iLevel+1}(jState+1)-STree{iLevel+1}(jState+2)));

            end
        
        PTree{iLevel}(2,jState) = 1 - PTree{iLevel}(1,jState) - PTree{iLevel}(3,jState);
        end
                
        
        if(any(any(PTree{iLevel}<0)))
           error(message('fininst:itttree:NegativeProbabilities'));

        end
        
        % Arrow-Debreu Prices:
        if(iLevel==1)
            % First level
            Q{iLevel+1}(:) = PTree{iLevel}(:,1)*RiskFreeDisc(iLevel);
        else
            if(jState==jMax(iLevel))  % Last Node
                Q{iLevel+1}(end) = Q{iLevel}(jState)*PTree{iLevel}(3,jState)*RiskFreeDisc(iLevel);
            elseif(jState==1) % First Node
                
                Q{iLevel+1}(jState+2) = (PTree{iLevel}(1,jState+2)*Q{iLevel}(jState+2) + ...
                                         PTree{iLevel}(2,jState+1)*Q{iLevel}(jState+1) + ...
                                         PTree{iLevel}(3,jState)*Q{iLevel}(jState))*RiskFreeDisc(iLevel);
                                     
                Q{iLevel+1}(jState+1) = (PTree{iLevel}(1,jState+1)*Q{iLevel}(jState+1) + ...
                                         PTree{iLevel}(2,jState)*Q{iLevel}(jState))*RiskFreeDisc(iLevel);

                Q{iLevel+1}(jState) =   PTree{iLevel}(1,jState)*Q{iLevel}(jState)*RiskFreeDisc(iLevel);
                                     
            elseif(jState==jMax(iLevel)-1) % Previous-to-last node
                
                Q{iLevel+1}(jState+2) = (PTree{iLevel}(2,jState+1)*Q{iLevel}(jState+1) + ...
                                         PTree{iLevel}(3,jState)*Q{iLevel}(jState))*RiskFreeDisc(iLevel);

            else % Everyone else
                
                Q{iLevel+1}(jState+2) = (PTree{iLevel}(1,jState+2)*Q{iLevel}(jState+2) + ...
                                         PTree{iLevel}(2,jState+1)*Q{iLevel}(jState+1) + ...
                                         PTree{iLevel}(3,jState)*Q{iLevel}(jState))*RiskFreeDisc(iLevel);
            end
        end % if iLevel=1
    end % jState
    if(any(Q{iLevel+1}<0))
        % We should NEVER get here.
        error(message('fininst:itttree:NegativeArrowDebreuPrices'));
    end
end % iLevel


% Build Tree Structure
ITTree = classfin('ITStockTree');
ITTree.StockSpec = StockSpec;
ITTree.StockOptSpec = StockOptSpec;
ITTree.TimeSpec = ITTTimeSpec;
ITTree.RateSpec = RateSpecOri;
ITTree.tObs = TreeTimes';
ITTree.dObs = TreeDates';
ITTree.STree = STree;
ITTree.Probs = PTree;

return;


% =============================================================
% AUX Function Block
% =============================================================
function CashFlows = TreeCashFlows(RateSpec, ExTimes, DivAmounts, TreeTimes)

% Ignore cash flows occurring before and at tree ValuationDate
BadCFIdx = ExTimes <= TreeTimes(1);
if any(BadCFIdx)
    warning(message('fininst:itttree:DivBeforeValuationDate'));
    ExTimes(BadCFIdx) = [];
    DivAmounts(BadCFIdx) = [];
end

% Pull out div at ValuationDate. They don't need to be discounted:
ValDatePayment = 0;
if ExTimes(1) == TreeTimes(1)
    ValDatePayment = DivAmounts(1);
    DivAmounts(1)=[];
    ExTimes(1)=[];
end


% Union all times but retrieve location of tree times and ExTimes
AllTimes = union(TreeTimes, ExTimes);
[dummy, TTIdx] = ismember(TreeTimes, AllTimes);%#ok
[dummy, ETIdx] = ismember(ExTimes, AllTimes);%#ok


% Find applicable discounts to bring CFs to present value
RateSpecNew = intenvset(RateSpec, 'StartTimes', AllTimes(1:end-1), ...
    'EndTimes', AllTimes(2:end));

Disc = intenvget(RateSpecNew, 'Disc');

% Build a discount matrix to find the present value of future
% dividend payments on each node

% Find the index of elements in the lower right triangle
% of matrix, excluding the anti-diagonal
MaskIdx = find(~flipud(triu(ones(length(Disc)))));

% Calculate cummulative discounts at each node
CummDisc = flipud(Disc(:)) * ones(size(Disc'));
CummDisc(MaskIdx) = 1;
CummDisc = cumprod(CummDisc);
CummDisc(MaskIdx) = 0;
CummDisc = flipud(CummDisc);

% Calculate cummulative discounted future dividends
AllAmounts  = zeros(length(AllTimes)-1,1);
AllAmounts(ETIdx-1) = DivAmounts;

% Last node always has zero since there are no future dividends to consider.
% Hence, pad the result vector of discounted dividends with a zero.
CashFlows = [CummDisc * AllAmounts; 0];

% Extract Cash Flows on tree times
CashFlows = CashFlows(TTIdx);

% Add any payments made at ValDate
CashFlows(1) = CashFlows(1);

return;

%----------------------------
function DivFactors = getdivfactors(TreeTimes, ExTimes, DivAmounts)

NDiv = length(DivAmounts);
NumLevels = length(TreeTimes);

TreeTimes = (TreeTimes(:))';
TreeTimesEx = TreeTimes(ones(NDiv,1), :);


ExTimes = ExTimes(:);
ExTimesEx = ExTimes(:, ones(1, NumLevels));

DivAmounts = 1-DivAmounts(:);
DivAmountsEx = DivAmounts(:, ones(1, NumLevels));

% Put ones only where the ExDate snaps into the TreeDat
Mask = diff([false(NDiv,1) TreeTimesEx >= ExTimesEx], 1, 2);

DivFactors = Mask .* DivAmountsEx;
DivFactors(~Mask) = 1;

if(size(DivFactors,1)>1)
    DivFactors = prod(DivFactors);
end

return;

% -------------------------------------------------------------------------
% Display the list of the required range of option prices.
function OptionData = StockOptRequirements(ITTTimeSpec, STree)
TreeDates = datestr(ITTTimeSpec.dObs,2);
NumLevels = length(STree);
NumStates = 2*(1:NumLevels)-1;
idx=1;
for iLevel=2:NumLevels
    for jState=1:iLevel-1
        % calls
        OptionData(idx)=struct('Type', 'call', 'Settle', TreeDates(1,:), 'Maturity', TreeDates(iLevel,:), 'Strike', STree{iLevel}(jState)); %#ok
        idx=idx+1;
    end
    for jState=iLevel:NumStates(iLevel)
        % puts
        OptionData(idx)=struct('Type', 'put', 'Settle', TreeDates(1,:), 'Maturity', TreeDates(iLevel,:), 'Strike', STree{iLevel}(jState));
        idx=idx+1;
    end
end
return

% -------------------------------------------------------------------------
% Obtain the implied volatility based on the input data
function [Sigma_i, Strike, Maturity] = GetImpliedVolatility(StockSpec, RateSpecOri, StockOptSpec)

% Save verified data
So = StockSpec.AssetPrice;
OptPrice = StockOptSpec.OptPrice;
Strike = StockOptSpec.Strike;
Settle = StockOptSpec.Settle;
Maturity = StockOptSpec.Maturity;
MaturityTimes = date2time(Settle, Maturity, -1);
OptSpec = StockOptSpec.OptSpec;

% Find the Risk-free rates for the maturity of all options:
RS = intenvset(RateSpecOri, 'StartDates', Settle, 'EndDates', Maturity);

RiskFreeRates = intenvget(RS, 'Rates');
Sigma_i = blsimpv(So, Strike, RiskFreeRates, MaturityTimes, OptPrice, [], [], [], OptSpec);
return

% -------------------------------------------------------------------------
% Main interpolation function for option prices. The interpolation is
% actually done in the volatility domain because it is less prone to errors
% due to extrapolation.
% -------------------------------------------------------------------------
function [CallPrices, PutPrices] = InterpOptPrices(STree, StockSpec, RateSpecOri, ...
    TimeSpec, StockOptSpec, TreeDates, TreeTimes, Sigma_i, InterpMethod)

InOptPrices = StockOptSpec.OptPrice;
Strike = StockOptSpec.Strike;
Settle = StockOptSpec.Settle;
Maturity = date2time(Settle, StockOptSpec.Maturity, -1, TimeSpec.Basis, TimeSpec.EndMonthRule);
OptSpec = StockOptSpec.OptSpec;
So = STree{1};

% Initialize OptPrices Tree
NumLevels = length(STree);
CallPrices =  mktrintree(NumLevels,1,1:2:2*NumLevels-1, NaN);
PutPrices = CallPrices;

% Separate puts from calls
callMask = strcmpi(OptSpec, 'call');
putMask = ~callMask;

% Use flags to indicate the type of interpolation
bPriceInterp = strcmp(InterpMethod, 'price');
bVolInterp = ~bPriceInterp;

% If the volatility will be needed, eliminate cases where the it could not 
% be calculated.
callSigma_i = [];
putSigma_i = [];
callPrice = [];
putPrice = [];

if bVolInterp
    nanVolMask = isnan(Sigma_i);
    
    % Some of the vols may be NaN  - consequence of using blsimpv. Ignore
    % these data points for the purposes of interpolation.
    callMask = callMask & ~nanVolMask;
    putMask = putMask & ~nanVolMask;

callSigma_i = Sigma_i(callMask);
putSigma_i = Sigma_i(putMask);
end

if bPriceInterp
callPrice = InOptPrices(callMask);
putPrice = InOptPrices(putMask);
end

% Separate the rest of the call and put data
callStrike = Strike(callMask);
putStrike = Strike(putMask);
callMaturity = Maturity(callMask);
putMaturity = Maturity(putMask);

% Make sure compounding is set to continuous, and that the rates are zero
% rates.
if(intenvget(RateSpecOri, 'Compounding') ~= -1)
    RateSpecOri = intenvset(RateSpecOri,'Compounding', -1);
end
RS = intenvset(RateSpecOri, 'StartDates', TreeDates(1), 'EndDates', TreeDates(2:end));
Rates = intenvget(RS, 'Rates');

% -------------------------------------------------------------------------
% We will first try regular interpolation, hoping that the option price
% data passed in by the user covers the range needed by the tree. If that
% is not the case, we need to extrapolate, and this may not a good thing -
% especially if the dataset if coarse. If this happens, we may end up with
% physically unfeasible values, such as negative prices.
% -------------------------------------------------------------------------
SnappedToZero=false;
ExtrapOpts = [];
ExtrapOptsIdx = 1;
for iLevel=2:NumLevels
    CallTreeMask=false(size(STree{iLevel}(:)));
    CallTreeMask(1:iLevel-1)=true;
    PutTreeMask=~CallTreeMask;
    
    % We interpolate based on the interpolation method specified by the
    % user:
    if(bPriceInterp)
        inCallData = callPrice;
        inPutData = putPrice;
    else
        inCallData = callSigma_i;
        inPutData = putSigma_i;
    end
    
    callInterpData = griddata(callStrike,callMaturity,inCallData,STree{iLevel}(CallTreeMask),TreeTimes(iLevel), 'cubic');
    if(any(isnan(callInterpData)))        
        % Save the option information so that we can list it at the end
        NanOpts = find(isnan(callInterpData));
        TStrikes  = STree{iLevel}(CallTreeMask);
        for NanIdx = NanOpts
            ExtrapOpts(ExtrapOptsIdx).OptSpec = 'call';
            ExtrapOpts(ExtrapOptsIdx).Strike = TStrikes(NanIdx); %#ok
            ExtrapOpts(ExtrapOptsIdx).Maturity = TreeDates(iLevel);
            ExtrapOptsIdx = ExtrapOptsIdx + 1;
        end
        callInterpData = griddata(callStrike,callMaturity,inCallData,STree{iLevel}(CallTreeMask),TreeTimes(iLevel), 'v4');        
    end

    putInterpData = griddata(putStrike,putMaturity,inPutData,STree{iLevel}(PutTreeMask)',TreeTimes(iLevel), 'cubic');
    if(any(isnan(putInterpData)))        
         % Save the option information so that we can list it at the end
        NanOpts = find(isnan(putInterpData));
        TStrikes  = STree{iLevel}(PutTreeMask);
        for NanIdx = NanOpts
            ExtrapOpts(ExtrapOptsIdx).OptSpec = 'put';
            ExtrapOpts(ExtrapOptsIdx).Strike = TStrikes(NanIdx); %#ok
            ExtrapOpts(ExtrapOptsIdx).Maturity = TreeDates(iLevel);
            ExtrapOptsIdx = ExtrapOptsIdx + 1;
        end
        putInterpData = griddata(putStrike,putMaturity,inPutData,STree{iLevel}(PutTreeMask)',TreeTimes(iLevel), 'v4');        
    end
    
    % Snap to zero if interp leads to neg values - this may happen if not
    % enough data points (coarse interpolation)
    if(any(callInterpData<0) || any(putInterpData<0))
        % We should never get negative vols or prices. Still, interpolation
        % may lead to very small negative values when extrapolating. Snap
        % them to zero if this happens.
        callInterpData(callInterpData<0)=0;
        putInterpData(putInterpData<0)=0;
    end

    if(bVolInterp)
    % -----------------------------------------------------------------------
    % Now go back to the price domain using blsprice:
        Call = blsprice(So, STree{iLevel}(CallTreeMask)', Rates(iLevel-1), TreeTimes(iLevel), callInterpData(:));
        [Dummy, Put] = blsprice(So, STree{iLevel}(PutTreeMask)', Rates(iLevel-1), TreeTimes(iLevel), putInterpData(:)); %#ok
    elseif(bPriceInterp)
        Call = callInterpData;
        Put = putInterpData;
    end    
    CallPrices{iLevel}(CallTreeMask) = Call; 
    PutPrices{iLevel}(PutTreeMask) = Put;           
end

% Warn the user about possible problems with the interpolation process, and
% show the options that caused extrapolation
if(length(ExtrapOpts)>0)
    msg = '';
    for optIdx = 1:length(ExtrapOpts)
        msg = [msg sprintf('Option Type: ''%s''   Maturity: %s  Strike=%.4f\n', ExtrapOpts(optIdx).OptSpec, ...
           datestr(ExtrapOpts(optIdx).Maturity), ExtrapOpts(optIdx).Strike)];
    end            
    warning(message('fininst:itttree:Extrapolation', msg));
end

return
