function ImpRepo = bndfutimprepo(Price, FutPrice, FutSettle, ...
    Delivery, ConvFactor, CouponRate, Maturity,varargin)
%BNDFUTIMPREPO Implied repo rates for a bond future given price.
%
%   ImpRepo = bndfutimprepo(Price, FutPrice, FutSettle, ...
%                          Delivery, ConvFactor, CouponRate, Maturity)
%
%   ImpRepo = bndfutprice(Price, FutPrice, FutSettle, ...
%                          Delivery, ConvFactor, CouponRate, Maturity,...)
%
%   Optional Inputs: RepoBasis,ReinvestRate,ReinvestBasis,
%                    Period, Basis, EndMonthRule,
%                    IssueDate, FirstCouponDate, LastCouponDate, StartDate,
%                    Face
%
%   Description:
%       BNDFUTIMPREPO computes the implied repo rate for a bond future given
%       the price of a bond and the price of the future and bond properties, including the
%       bond's conversion factor.
%
%       Note that the default behavior is for the coupon reinvestment rate
%       to be identical to the repo rate -- however, the user can specify a
%       seprate reinvestment rate using optional inputs.
%
%   Inputs:
%       Price - numBonds X 1 vector of Bond Prices
%       FutPrice - numBonds X 1 vector of Future Prices
%       FutSettle - numBonds X 1 vector of Future Settle Dates
%       Delivery - numBonds X 1 vector of Future Delivery Dates
%       ConvFactor - numBonds X 1 vector of Bond Conversion Factors
%       CouponRate - numBonds X 1 vector of Coupon rates in decimal form.
%       Maturity - numBonds X 1 vector of Maturity dates.
%
%   Optional Inputs:
%       RepoBasis - Day count basis for repo rate (default is 2)
%       ReinvestRate - Rate for reinvesting intermediate coupons from the
%                           bond -- default is for this to be identical to
%                           the repo rate.
%       ReinvestBasis - Day count basis for resinvestment rate (default is
%                           identical to the repo basis)
%
%   Additionally, the following bond properties can also be set via
%   parameter value pair:
%
%            Period - Number of coupons payments per year.
%                     Possible values include:
%                     0, 1, 2 (default), 3, 4, 6, 12
%
%             Basis - Day-count basis.
%                     Possible values include:
%                     0 - actual/actual (default)
%                     1 - 30/360 SIA
%                     2 - actual/360
%                     3 - actual/365
%                     4 - 30/360 PSA
%                     5 - 30/360 ISDA
%                     6 - 30/360 European
%                     7 - actual/365 Japanese
%                     8 - actual/actual ISMA
%                     9 - actual/360 ISMA
%                    10 - actual/365 ISMA
%                    11 - 30/360 ISMA
%                    12 - actual/365 ISDA
%                    13 - bus/252
%
%      EndMonthRule - End-of-month rule; default is 1 (in effect)
%                     0 - Rule is NOT in effect for the bond(s)
%                     1 - (default) Rule is in effect for the bond(s) (meaning
%                         that a security that pays coupon interest on the last
%                         day of the month will always make payment on the last
%                         day of the month)
%
%         IssueDate - Bond issue date.
%
%   FirstCouponDate - Irregular or normal first coupon date.
%
%    LastCouponDate - Irregular or normal last coupon date.
%
%         StartDate - Forward starting date of payments.
%
%              Face - Face value of the bond; default is 100.  Note that
%                     Face has no impact on key rate duration -- this
%                     calling sequence is preserved for consistency.
%
%   Outputs:
%       ImpRepo - Implied repo rate, or the repo rate that would produce
%       the price input.
%
%   Example:
%
%   This example is based on example from Krgin:
%
%   bndfutimprepo(129,98,'9/21/2000','12/29/2000',1.3136,.0875,'8/15/2020')
%
%   See also BNDFUTPRICE, CONVFACTOR.

%   References:
%   [1] Krgin, Dragomir, Handbook of Global Fixed Income
%              Calculations, John Wiley & Sons, 2002

%   Copyright 1995-2017 The MathWorks, Inc.

% Parse input arguments
p = inputParser;

% Parse the bond inputs with no checking, do that in INSTARGBOND
p.addParameter('period',[]);
p.addParameter('basis',[]);
p.addParameter('endmonthrule',[]);
p.addParameter('issuedate',[]);
p.addParameter('firstcoupondate',[]);
p.addParameter('lastcoupondate',[]);
p.addParameter('startdate',[]);
p.addParameter('face',[]);

p.addParameter('repobasis',2,@(x) all(isvalidbasis(x)));
p.addParameter('reinvestrate',NaN,@isnumeric);
p.addParameter('reinvestbasis',[],@(x) all(isvalidbasis(x)));

% Check to make sure that Settle is a scalar
FutSettle = datenum(FutSettle);
Delivery = datenum(Delivery);

try
    p.parse(varargin{:});
catch ME
    newMsg = message('fininst:bndfutimprepo:optionalInputError');
    newME = MException(newMsg.Identifier,getString(newMsg));
    newME = addCause(newME,ME);
    throw(newME)
end

[CouponRate, FutSettle, Maturity, Period, Basis, EndMonthRule, IssueDate, ...
    FirstCouponDate, LastCouponDate, StartDate, Face] = ...
    instargbond(CouponRate,FutSettle,Maturity,...
    p.Results.period, p.Results.basis, p.Results.endmonthrule, ...
    p.Results.issuedate,p.Results.firstcoupondate,  ...
    p.Results.lastcoupondate, p.Results.startdate,  p.Results.face);

RepoBasis = p.Results.repobasis;
ReinvestRate = p.Results.reinvestrate;
ReinvestBasis = p.Results.reinvestbasis;

if isempty(ReinvestBasis)
    ReinvestBasis = RepoBasis;
end

[CouponRate, FutSettle, Maturity, Period, Basis, EndMonthRule, IssueDate, ...
    FirstCouponDate, LastCouponDate, StartDate, Face, Price,RepoBasis,...
    ReinvestBasis,Delivery,FutPrice,ConvFactor, ReinvestRate] = ...
    finargsz('all', CouponRate, FutSettle, Maturity, Period, Basis, ...
    EndMonthRule, IssueDate, FirstCouponDate,LastCouponDate, ...
    StartDate, Face, Price,RepoBasis,...
    ReinvestBasis,Delivery,FutPrice,ConvFactor, ReinvestRate);

AI_s = Face.*CouponRate./Period.*accrfrac(FutSettle, Maturity, Period, Basis,...
    EndMonthRule,IssueDate,FirstCouponDate,LastCouponDate,StartDate,Face);
AI_d = Face.*CouponRate./Period.*accrfrac(Delivery, Maturity, Period, Basis,...
    EndMonthRule,IssueDate,FirstCouponDate,LastCouponDate,StartDate,Face);

% Determine if there are any coupons to reinvest
[CFlowAmounts, CFlowDates] = cfamounts(CouponRate,...
    FutSettle, Maturity, Period, Basis, EndMonthRule, IssueDate,...
    FirstCouponDate, LastCouponDate,StartDate, Face);

nBonds = length(CouponRate);
ImpRepo = zeros(nBonds,1);

% The formulas used in the loop are from the Krgin reference above,
% Chapter 9, page 147
for bondidx=1:nBonds
    
    % Index of remaining coupons on bond between Settlement and Delivery
    intcoupidx = find(FutSettle(bondidx) < CFlowDates(bondidx,:) & ...
        CFlowDates(bondidx,:)<= Delivery(bondidx));
    nCoupons = length(intcoupidx);
    
    if isnan(ReinvestRate)
        % If reinvest is NaN, we assume that reinvest and repo rate are same
        switch nCoupons
            case 0
                ImpRepo(bondidx) = ((FutPrice(bondidx)*ConvFactor(bondidx) + AI_d(bondidx))./ ...
                    (Price(bondidx) + AI_s(bondidx)) - 1)./...
                    (yearfrac(FutSettle(bondidx),Delivery(bondidx),RepoBasis(bondidx)));
            case 1
                ImpRepo(bondidx) = (FutPrice(bondidx)*ConvFactor(bondidx) - ...
                    Price(bondidx) - AI_s(bondidx) + CFlowAmounts(bondidx,intcoupidx) ...
                    + AI_d(bondidx))./ ((Price(bondidx) + AI_s(bondidx))*...
                    yearfrac(FutSettle(bondidx),Delivery(bondidx),RepoBasis(bondidx)) - ...
                    CFlowAmounts(bondidx,intcoupidx)*yearfrac(CFlowDates(bondidx,intcoupidx),...
                    Delivery(bondidx),RepoBasis(bondidx)));
            otherwise
                dSD_D = yearfrac(FutSettle(bondidx),Delivery(bondidx),...
                                                        RepoBasis(bondidx));
                Days = yeardays(year(Delivery(bondidx)), RepoBasis(bondidx));
                b_D = yearfrac(CFlowDates(bondidx,intcoupidx(end)),...
                                    Delivery(bondidx),RepoBasis(bondidx));
                a_D = diff(CFlowDates(bondidx,intcoupidx))/Days;
                imprepofun = @(r) (Price(bondidx) + AI_s(bondidx))*(1 + r*dSD_D)...
                     - (FutPrice(bondidx)*ConvFactor(bondidx)) - AI_d(bondidx) - ...
                    (1 + r*b_D)*sum(CFlowAmounts(bondidx,intcoupidx).*(cumprod(1 + r*[1 a_D])));
                
                x0 = .05;
                [x,~,exitflag] = fzero(imprepofun,x0);
                if exitflag < 0
                    warning(message('fininst:bndfutimprepo:solutionConvergenceFailure', num2str( bondidx )))
                end
                ImpRepo(bondidx) = x;
                
        end
    else
        switch nCoupons
            case 0
                ImpRepo(bondidx) = ((FutPrice(bondidx)*ConvFactor(bondidx) + AI_d(bondidx))./ ...
                    (Price(bondidx) + AI_s(bondidx)) - 1)./...
                    (yearfrac(FutSettle(bondidx),Delivery(bondidx),RepoBasis(bondidx)));
            case 1
                ImpRepo(bondidx) = ((FutPrice(bondidx).*ConvFactor(bondidx) + ...
                    CFlowAmounts(bondidx,intcoupidx).*(1 + ReinvestRate(bondidx)*...
                    yearfrac(CFlowDates(bondidx,intcoupidx),Delivery(bondidx),ReinvestBasis(bondidx))) ...
                    + AI_d(bondidx))./(Price(bondidx) + AI_s(bondidx)) - 1)./ ...
                    yearfrac(FutSettle(bondidx),Delivery(bondidx),RepoBasis(bondidx));
            otherwise
                ImpRepo(bondidx) = ((FutPrice(bondidx).*ConvFactor(bondidx) + ...
                    sum(CFlowAmounts(bondidx,intcoupidx).*(cumprod(1 + ReinvestRate(bondidx)*...
                    yearfrac(CFlowDates(bondidx,intcoupidx),Delivery(bondidx),ReinvestBasis(bondidx))))) ...
                    + AI_d(bondidx))./(Price(bondidx) + AI_s(bondidx)) - 1)./ ...
                    yearfrac(FutSettle(bondidx),Delivery(bondidx),RepoBasis(bondidx));
        end 
    end
end