classdef DieboldLi < IRCurve
    %DIEBOLDLI Diebold Li model of yield curve
    %   Implementation of the Diebold Li model from the following paper:
    %   Francis X. Diebold, Canlin Li, Forecasting the term structure of
    %   government bond yields, Journal of Econometrics, Volume 130,
    %   Issue 2, February 2006, Pages 337-364
    %
    %   DieboldLi has the following properties:
    %
    %   Type: Type of Curve, either Zero or ParYield
    %   Settle: Settle date for curve (scalar)
    %   Beta1:Beta1 parameter
    %   Beta2:Beta2 parameter
    %   Beta3:Beta3 parameter
    %   Lambda: Time parameter
    %   Compounding: Compounding for object
    %   Basis: Day-count basis.
    %
    %   DieboldLi has the following methods, in addition to its constructor:
    %
    %   GETFORWARDRATES
    %   GETZERORATES
    %   GETDISCOUNTFACTORS
    %   GETPARYIELDS
    %   TORATESPEC
    %
    %   And the static method:
    %
    %   FITBETASFROMYIELDS
    
    %   Copyright 2011 The MathWorks, Inc.
    
    properties (SetAccess = protected, GetAccess = public)
        Beta1
        Beta2
        Beta3
        Lambda
    end
    
    methods
        function obj = DieboldLi(Type,Settle,Params,varargin)
            %DIEBOLDLI Create a DieboldLi object
            %   DLI = DIEBOLDLI(TYPE,SETTLE, PARAMS) creates a curve with specified
            %   Type and Settle date and Params vector.
            %
            %   PARAMS is a 1 x 4 vector of Beta1, Beta2, Beta3 and Lambda
            %   parameters of the Diebold Li model.
            %
            %   Optional inputs include BASIS and COMPOUNDING, specified in parameter
            %   value pairs
            %
            %   Example:  Create a DieboldLi object
            %
            %           c = DieboldLi('zero',today,[5 2 1 2]);
            %
            %
            
            if all(~strcmpi(Type,{'zero','paryield'}))
                error(message('fininst:DieboldLi:invalidType'))
            end
            
            obj.Type = Type;
            
            try
                obj.Settle = datenum(Settle);
            catch E
                error(message('fininst:DieboldLi:invalidSettle'))
            end
            
            if length(obj.Settle) > 1
                error(message('fininst:DieboldLi:invalidSettleSize'))
            end
            
            if length(Params) ~= 4
                error(message('fininst:DieboldLi:invalidParamSize'))
            end
            
            obj.Beta1 = Params(1);
            obj.Beta2 = Params(2);
            obj.Beta3 = Params(3);
            obj.Lambda = Params(4);
            
            % Parse inputs
            p = inputParser;
            
            p.addParamValue('compounding',2,@(x) ismember(x,[-1 1 2 3 4 6 12]));
            p.addParamValue('basis',0,@isvalidbasis);
            
            try
                p.parse(varargin{:});
            catch ME
                newMsg = message('fininst:DieboldLi:optionalInputError');
                newME = MException(newMsg.Identifier,getString(newMsg));
                newME = addCause(newME,ME);
                throw(newME)
            end
            
            obj.Compounding = p.Results.compounding;
            obj.Basis = p.Results.basis;
        end
        
        function outForwardRates = getForwardRates(obj,inpDates)
            Params = [obj.Beta1 obj.Beta2 obj.Beta3 obj.Lambda];
            switch lower(obj.Type)
                case 'zero'
                    ZeroRates = locDieLi(Params,obj.Settle,inpDates,obj.Basis);
                case 'paryield'
                    ParRates = locDieLi(Params,obj.Settle,inpDates,obj.Basis);
                    ZeroRates = pyld2zero(ParRates, inpDates, obj.Settle, obj.Compounding, obj.Basis);
            end
            outForwardRates = zero2fwd(ZeroRates, inpDates, obj.Settle);
        end
        
        function outZeroRates = getZeroRates(obj,inpDates)
            Params = [obj.Beta1 obj.Beta2 obj.Beta3 obj.Lambda];
            switch lower(obj.Type)
                case 'zero'
                    outZeroRates = locDieLi(Params,obj.Settle,inpDates,obj.Basis);
                case 'paryield'
                    ParRates = locDieLi(Params,obj.Settle,inpDates,obj.Basis);
                    outZeroRates = pyld2zero(ParRates, inpDates, obj.Settle, obj.Compounding, obj.Basis);
            end
        end
        
        function outDiscountFactors = getDiscountFactors(obj,inpDates)
            Params = [obj.Beta1 obj.Beta2 obj.Beta3 obj.Lambda];
            switch lower(obj.Type)
                case 'zero'
                    ZeroRates = locDieLi(Params,obj.Settle,inpDates,obj.Basis);
                case 'paryield'
                    ParRates = locDieLi(Params,obj.Settle,inpDates,obj.Basis);
                    ZeroRates = pyld2zero(ParRates, inpDates, obj.Settle, obj.Compounding, obj.Basis);
            end
            outDiscountFactors = zero2disc(ZeroRates, inpDates, obj.Settle);
        end
        function outParYields = getParYields(obj,inpDates)
            Params = [obj.Beta1 obj.Beta2 obj.Beta3 obj.Lambda];
            switch lower(obj.Type)
                case 'zero'
                    ZeroRates = locDieLi(Params,obj.Settle,inpDates,obj.Basis);
                    outParYields = zero2pyld(ZeroRates, inpDates, obj.Settle, obj.Compounding, obj.Basis);
                case 'paryield'
                    outParYields = locDieLi(Params,obj.Settle,inpDates,obj.Basis);
            end
        end
        function outRateSpec = toRateSpec(obj,inpDates)
            Params = [obj.Beta1 obj.Beta2 obj.Beta3 obj.Lambda];
            switch lower(obj.Type)
                case 'zero'
                    ZeroRates = locDieLi(Params,obj.Settle,inpDates,obj.Basis);
                case 'paryield'
                    ParRates = locDieLi(Params,obj.Settle,inpDates,obj.Basis);
                    ZeroRates = pyld2zero(ParRates, inpDates, obj.Settle, obj.Compounding, obj.Basis);
            end
            outRateSpec = intenvset('Rates', ZeroRates(:), 'EndDates',...
                inpDates(:),'StartDates',obj.Settle,...
                'Basis',obj.Basis,'Compounding',obj.Compounding);
        end
    end
    methods(Static = true)
        function obj = fitBetasFromYields(inSettle,lambda_t,inDates,inYields)
            %DIEBOLDLI/FITBETASFROMYIELDS Fits Diebold Li model from data and  lambda
            %   IRFC = DIEBOLDLI.FITBETASFROMYIELDS(SETTLE,LAMBDA,DATES,YIELDS)
            %   fits a Diebold Li model to market data given a lambda time
            %   parameter and returns a DieboldLi object.
            %
            %   SETTLE: Settle date for curve (scalar)
            %   LAMBDA: Time parameter from DieboldLi model -- recommended
            %   to be .0609 by the paper.
            %   DATES: Dates of input Yields
            %   YIELDS: Yields to be fit
            %
            TimeToMat = yearfrac(inSettle,inDates,1);
            X = [ones(size(TimeToMat)) (1 - exp(-lambda_t*TimeToMat))./(lambda_t*TimeToMat) ...
                ((1 - exp(-lambda_t*TimeToMat))./(lambda_t*TimeToMat) - exp(-lambda_t*TimeToMat))];
            Beta = regress(inYields,X);
            obj = DieboldLi('ParYield',inSettle,[Beta' lambda_t]);
        end
    end
end

function outZeroRates = locDieLi(Params,Settle,inpDates,Basis)
    inTimes = yearfrac(Settle,inpDates,Basis);
    Beta1 = Params(1);
    Beta2 = Params(2);
    Beta3 = Params(3);
    Lambda = Params(4);
    outZeroRates = Beta1 + Beta2*(1 - exp(-Lambda*inTimes))./(Lambda*inTimes) + Beta3*((1 - exp(-Lambda*inTimes))./(Lambda*inTimes) - exp(-Lambda*inTimes));
end