cds_price(initDate, endDate, initSpread, curSpreads, discountrate, recovery, freq, useAcc, dcc)

function [MTM RPV01 UFF] = cds_price(initDate, endDate, initSpread, curSpreads, discountrate, recovery, freq, useAcc, dcc)
% 
% Version 0.2
% Author: Rogier Swierstra
% [e-mail]
%
% function MTM =  cds_price(initDate, endDate, initSpread, curSpreads,
%                           discountrate, recovery, freq, useAcc, dcc)
%
% CDS_PRICE calculates the MTM price of a CDS deal, given change in spreads
% It is based on the pricing model of Dominic O'Kane and Stuart Turnbull,
% "Valuation of Credit Default Swaps" (April, 2003) available at:
% http://www.nuclearphynance.com/User%20Files/256/cds.pdf
%
% The input variables are as follows:
% - dates are in Matlab format 'DD-MMM-YY'
% - spreads are in BPs (basepoints)
% - recovery is the recovery assumption fraction (decimal)
% - discountrate is either a fixed discount rate, or a term structure
% - freq is the frequency of coupon payments, per year
% - useAcc is a boolean, whether to compute with accrued payments
% - dcc is the choice of day count convention, default ACT/ACT
% 
% The result of the calculation is the change in value for the protection
% buyer, i.e. an increase in spreads gives a positive MTM. 
%
% RPV01 is the risky present value of 1bp ending at a default 
% UFF is the up-front-fee corresponding to 500bp running spread
%
% example: cds_price('01-Jul-2008','01-Jul-2013',600,2000,.04,6,4,1,0)
%
% This function requires financial toolbox for date manipulations.

% Define primary variables
if nargin < 9, dcc = 0; end        % default: use ACT/ACT
if nargin < 8, useAcc = 1; end     % default: use accrued payments
if nargin < 7, freq = 4; end       % default: quarterly installments
if nargin < 6, recovery = 0.4; end % default: 40% recovery

% Date determinations
% cfdates, daysdif and yearfrac are FT-functions
curDate     = datenum(date());     % or specify the current date
start       = datenum(initDate);
firstcoupon = addtodate(start, 12/freq, 'month');
coupons     = [start cfdates(initDate,endDate,freq,dcc,0,initDate,firstcoupon)];
thisperiod  = find(coupons<curDate, 1,'last');
periods     = yearfrac(coupons(1:end-1), coupons(2:end),dcc);
accrperiod  = accrfrac(coupons(thisperiod),curDate,freq,dcc);

% Determine the discount vector
if isscalar(discountrate) % fixed rate
    discount = exp(-discountrate*periods);
elseif isvector(discountrate) && length(discountrate) == length(periods) % supplied term structure
    discount = cumprod(1./(1+discountrate).^periods);
else error('Discount vector of incorrect size')
end % if

% Compute survival probabilities.
% Note the conversion from basis points to decimal notation!
curTerms = terms(curSpreads, length(periods), freq); % see functions below
q = hazardrate(curTerms/10000, discount, recovery, periods, useAcc); % see functions below

% Risky present value of a basis point, and net mark-to-market value of a
% long position (protection buyer)
RPV01 = discount.*periods*(q(2:end) - useAcc*diff(q)/2)'/10000;
marketSpread = curTerms(thisperiod);
MTM = (marketSpread - initSpread)*RPV01 + useAcc*accrperiod*initSpread/10000;
UFF = MTM + (initSpread - 500)*RPV01;
end % function cds_price


% Given spreads, determine the survival probabilities during each period.
% Strictly speaking, this is not the hazard rate.
%
% This calculation cannot be explained in a comment line, but see the
% referenced article, section 9, and solve inductively for 
%   q(t) = exp(-lambda tau)
% Note that if accr=1 we do assume accrual payments, at halfway the
% timestep. If accr=0, we assume no accrual payments.
function q = hazardrate(spreads, disc, recov, d, accr)
    n = length(disc);   % time periods
    lgd = 1-recov;      % loss given default
    q = zeros(1, n+1);  % survival probabilities
    DL = zeros(1, n);   % Default Leg
    ppl = zeros(1, n);  % pre-spread premium Leg
    PL = zeros(1, n);   % Premium Leg, including accrual terms

    q(1) = 1; s = spreads(1);
    alpha = lgd - accr*s*d(1)/2;
    q(2) = alpha/(alpha + s*d(1));
    
    a = (lgd - (accr/2)*spreads.*d).*disc;
    b = (lgd + (1-accr/2)*spreads.*d).*disc;
    
    for i = 2:n
        DL(i) = DL(i-1) + lgd*disc(i)*(q(i-1)-q(i));
        ppl(i)= ppl(i-1) + disc(i)*(q(i)*d(i) + accr*d(i)*((q(i-1)-q(i))/2)); 
        PL(i) = ppl(i)*spreads(i);
        % a = lgd - (accr/2)*spreads(i)*d(i) ;
        % b = lgd + (1-accr/2)*spreads(i)*d(i);
        q(i+1) = (DL(i)-PL(i) + q(i)*a(i))/b(i);    
    end % for
end % function hazardrate


% Given spreads (5Y/1-3-5-7-10Y/term structure) in bps, convert to a vector
% of suitable length 
function spreads = terms(spr, n, f)
    if isscalar(spr) % Only CDS5Y
        spreads = repmat(spr, 1, n);
    elseif length(spr) == 5 % spreads for 1Y 3Y 5Y 7Y 10Y
        temp = [repmat(spr(1),1,f) repmat(spr(2),1,2*f) repmat(spr(3),1,2*f) repmat(spr(4),1,2*f) repmat(spr(5),1,3*f)];
        spreads = temp(1:n);
    elseif length(spr) == 6 % spreads for 1Y 2Y 3Y 5Y 7Y 10Y
        temp = [repmat(spr(1),1,f) repmat(spr(2),1,f) repmat(spr(3),1,f) repmat(spr(4),1,2*f) repmat(spr(5),1,2*f) repmat(spr(6),1,3*f)];
        spreads = temp(1:n);
    elseif length(spr) == n % if a complete spread structure is available
        spreads = spr;
    end % if
end % function terms

% NOTES
%
% * The function "terms" still works best with flat spreads. No 
% interpolation is done between given data points.
%
% * Fractional maturities may still cause problems.
%

% Author: Rogier Swierstra
% [e-mail]