ShapleyRiskDecomposition(C,w,createplot,names)

function PCTR = ShapleyRiskDecomposition(C,w,createplot,names)
% PCTR = ShapleyRiskDecomposition(C,w,createplot,names);
%
% Calculates the Shapley percentage contibution to risk (=volatility) 
% of a set N number of assets
% Warning    : Takes too long to end if N > 20  
%
% Inputs     :
% ------------
% C          : NxN covariance matrix
% w          : Nx1 weights vector
% createplot : (=1) if a plot of risk contributions is needed 0 otherwise
% names      : Nx1 cell array of the names of assets 
%
%
% Outputs:
% ------------
% PCTR       : Nx2 matrix of percentage contribution to risk of w
%              The first column is Shapley PCTR, 
%              the second column is Euler PCTR. 
%
%
% Author: Hakan Kaya
% Date  : February 2nd, 2013

% Demos
%
% Sensitivity to correlation. What is the difference between Shapley and
% Euler decompositions if correlation (\rho) varies between two assets.
% PCTR = ShapleyRiskDecomposition; 
%
% Create random covariance and calculate risk contrib of equal weight.
% N = 10; % number of assets
% U = randn(100,N); C = U'*U/2; w = ones(N,1)/N;
% PCTR = ShapleyRiskDecomposition(C,w,1,[]);

if nargin < 1
    x = [90 45 0 -45 -90]/100;
    w = (0:1:100)/100;
    lt = cell(length(x),1);
    Y = zeros(length(w),length(x));
    X = Y;
    for i = 1:length(x)
        lt(i,1) = {['\rho = ',num2str(x(i))]};
        C = [1 x(i);x(i) 1];
        for j = 1:length(w)
            PCTR = ShapleyRiskDecomposition(C,[w(j) 1-w(j)],0,[]);
            Y(j,i) = PCTR(1,1);
            X(j,i) = PCTR(1,2);
        end
    end
    figure;
    m = min([min(min(Y)) min(min(X))]);
    M = max([max(max(Y)) max(max(X))]);
    plot(w,Y,'linewidth',3); grid on
    
    xlabel('Weight of Asset 1')
    ylabel('PCTR of Asset 1')
    hold on;
    plot(w,X,'linewidth',1); grid on
    title('Euler(Light) vs Shapley(Thick) PCT Volatility','Fontsize',15)
    
    xlabel('Weight of Asset 1')
    ylabel('PCTR of Asset 1')
    ylim([m M])
    legend(lt,'location','nw')
    drawnow;
    return;
end

if exist('createplot','var')==0
    createplot = 0;
    if exist('names','var')==0
        names = [];
    end
end

w = w(:);
n = size(C,1);
S = 1:n;
CTR = zeros(n,1);
[YY,NB] = nchoose(setdiff(S,1));
for k = 1:n
    c = 0;
    if k > 1
        YY = replaceelement(YY,k,k-1);
    end
    for s = 0:(n-1)      
        Y = YY(NB==s);
        for j = 1:size(Y,1)
            D = DeltaF(k,Y(j),w,C);
            c = c + (1 / (nchoosek(n,s)*(n-s))) * D;
        end
    end
    CTR(k,1) = c;
end
PCTR = CTR / GetF(C,w);
PCTR = [PCTR,diag(w)*C*w/ ( w'*C*w )];


if createplot == 1
   if isempty(names)
       names = 1:size(PCTR,1);
   end
   figure;
   subplot(3,1,1)
   bar(PCTR); grid on;
   xlim([0 size(PCTR,1)+1])
   title('Euler and Shapley Risk Decomposition','FontSize',15);
   legend('Shapley','Euler')
   set(gca,'xtick',1:size(PCTR,1),'xticklabel',names)
   subplot(3,1,2)
   [sd,Cr] = cov2corr(C);
   bar(sd); grid on;
   xlim([0 size(PCTR,1)+1])
   title('Volatility','FontSize',15);
   set(gca,'xtick',1:size(PCTR,1),'xticklabel',names)
   subplot(3,1,3)  
   imagesc(Cr); grid on; colorbar('Location','SouthOutside');
   set(gca,'xtick',1:size(PCTR,1),'xticklabel',names)
   set(gca,'ytick',1:size(PCTR,1),'yticklabel',names)
   title('Correlation Heatmap','FontSize',15);   
end


function YY = replaceelement(YY,k,l)
for i = 1:length(YY)
    YY{i}(YY{i}==k) = l;
end

function x = DeltaF(k,I,w,C)
J = [I{1},k];
if sum(I{1})>0
    x = GetF(C(J,J),w(J)) - GetF(C(I{1},I{1}),w(I{1}));
else
    x = GetF(C(J,J),w(J));
end


function x = GetF(C,w)
if isempty(w)==0
    x = sqrt(w'*C*w);
else
    x = 0;
end


function [W,C] = nchoose(S)
% NCHOOSE - all combinations of the elements of a set
%   W = nchoose(S) returns all possible combinations of one or more
%   elements of the set S. In total there are 2^N-1 combinations, where N
%   is the number of elements in S. W is a cell array: each cell holds one
%   of these combination (as a row vector). S can be a cell array, and each
%   cell of W will then contain a cell array.  
%
%   Example:
%      nchoose([2 4 6 8])
%        % -> { [2] ; 
%        %      [4] ; 
%        %      [2 4] ; 
%        %      [6] ;
%        %      ...
%        %      [2 6 8] ;
%        %      [4 6 8] ;
%        %      [2 4 6 8]} ; % in total (2^4)-1 = 15 different combinations
% 
%   Notes: 
%   - For sets containing more than 18 elements a warning is given, as this
%     can take some time. On my PC, a set of 20 elements took 20 seconds.
%     Hit Ctrl-C to intterupt calculations.
%   - If S contain non-unique elements (e.g. S = [1 1 2]), NCHOOSE will
%     return non-unique cells. In other words, NCHOOSE treats all elements
%     of S as being unique. One could use NCHOOSE(UNIQUE(S)) to avoid that.
%   - Loosely speaking, NCHOOSE(S) collects all output of multiple calls to
%     NCHOOSEK(S,K) where K is looping from 1 to the number of elements of
%     S. The implementation of NCHOOSE, however, does rely of a different
%     method and is much faster than such a loop.
%   - By adding an empty cell to the output, the power set of a cell array
%     is formed:
%       S = {1 ; 'hello' ; { 1 2 3}}
%       PowerS = nchoose(S) 
%       PowerS(end+1) = {[]}
%     See: http://en.wikipedia.org/wiki/Power_set
%   
%   See also NCHOOSEK, PERMS
%            COMBN, ALLCOMB on the File Exchange


N = numel(S) ; 

% How many unique combinations of 1 to N elements are there
M = (2^N)-1 ;
if N > 18,
    warning('Nchoose:LargeOutput', ...
        'There are %d unique combinations. Please be patient ...',M) ;
end

S = S(:).' ;   % make the set a row vector, for uniform output

% The selection of elements is based on the binary representation of all
% numbers between 1 and M. The binary representation of numbers stored in X
% can be calculated using the formula (see DEC2BIN):
%    B1 = rem(floor(X(:) * pow2(1-N:0)),2) 
% (* means matrix multiplication) in which row B(k,:) represents X(k)
% Based on the suggestion by US, this formula can be rewritten as:
%    B2 = bitget(X2 * (2.^(N-1:-1:0)),N) > 0
% for a specific number X2. Although the bits of B2 are in reversed order
% (compared to B1), this is, however, exactly what we want, as we can now
% select the elements of S directly. We do not need to flip the order
% afterwards. 

% Looping over the numbers is faster than creating a whole matrix at once
W = cell(M+1,1) ;    % Pre-allocation of output
p2=2.^(N-1:-1:0) ; % This part of the formula can be taken out of the loop
C = zeros(M+1,1);
for i=1:M,
    % calculate the (reversed) binary representation of i
    % select the elements of the set based on this representation
    W{i} = S(bitget(i*p2,N) > 0) ; 
    C(i,1) = length(W{i});
end
W{M+1} = [];

% Two other options, which are less efficient
% % 1) looping over nchoosek (slow)
% for i=N:-1:1,
%     W{i} = nchoosek(S,i) ;
% end
% % and now W still has to be split into distinct cells ...
% 
% % 2) in one big step (memory inefficient)
% Q = rem(floor([1:M].'*pow2(1-N:0)),2).' ;
% [R,C] = find(Q) ;
% W = mat2cell(S(R),1,sum(Q)).' ;