omega(returns,targets)

function S = omega(returns,targets)
% OMEGA Computes portfolio omega 
%
%   S = OMEGA(RETURNS) gives a structure S containing data so that:
%       S.returnLevel = input return levels
%       S.omegaValues = omega value at that return level
%   RETURNS must be a vector.
%
%   S = OMEGA(RETURNS,LEVELS) returns the same structure, but only provides
%   values at the return levels specified in LEVELS.
%
%   For example, to reproduce one of the images in the Keating and Shadwick
%   paper (page 9):
%
%       means = [6 7]; stds = [4 3]; npts = 5000;
%       C = {}; Lgnds = {};
%       for i = 1:length(means)
%           returns = randn(npts,1)*stds(i) + means(i);
%           S = omega(returns);
%           C{end+1} = S.returnLevel; C{end+1} = S.omegaValues;
%           Lgnds{end+1} = ['\mu = ' num2str(means(i))  ', \sigma = ' num2str(stds(i))];
%       end
%       plot(C{:}); legend(Lgnds{:}); xlim([3 6])

%   Version 1.2
%   Based on the omega value described in Keating and Shadwick,
%   "An Introduction to Omega", from the Finance Development Centre
%   Limited.   
%
%   Nabeel Azar
%   http://www.nabeelazar.com
%   [e-mail]

%   Parse inputs
if ~isvector(returns)
    error(['RETURNS input must be a vector']);
end
returns = returns(:);

%   Need to compute a cdf function for these returns.
%   First, get a sorted list of unique return values.
sortedReturns = unique(returns);

%   Get an estimate of the CDF.
%   The adjustmentVector accounts for duplicate returns
%   Only apply if needed
adjustmentVector = 0;
if (length(sortedReturns)~=length(returns))
    adjustmentVector = cumsum(hist(returns,sortedReturns) - 1);
end
numerator = (1:length(sortedReturns)) + adjustmentVector;
cdfValues = numerator.'./length(returns);

%   Get estimates of the areas relevant to computing omega
areaBelowAndToLeft = cumtrapz(sortedReturns,cdfValues);
areaAboveAndToRight = trapz(sortedReturns,1-cdfValues) - ... 
    cumtrapz(sortedReturns,1-cdfValues);

%   Now compute an omega value for each return
warningState = warning;
warning('off','MATLAB:divideByZero');
omegaValues = areaAboveAndToRight./areaBelowAndToLeft;
warning(warningState)

%   Put together the appropriate outputs.
switch nargin
    case 1
        S.returnLevel = sortedReturns;
        S.omegaValues = omegaValues;
    case 2
        %   Linearly interpolate to get the values at the 
        %   requested target points.
        samples = zeros(size(targets));
        samples(:) = interp1(sortedReturns,omegaValues,targets(:),'linear');
        S.returnLevel = targets;
        S.omegaValues = samples;
end