Incorrect negative portfolio weights for constrained minimum-variance portfolios

Hi Community, 

I implemented the Jagannathan & Ma (Journal of Finance, 2003) algorithm, which reduces the error in the sample covariance matrix to improve the out-of-sample performance of portfolios (at the end of this post is the excerpt from the paper with the theoretical derivation). The resulting portfolio weights are supposed to be nonnegative and in particular it should work for large portfolios where the number of assets exceed the number of observations (N>T). 

However, the results from the Gurobi optimizer show negative weights. And in particular if I increase N to N=100, the results get even worse, i.e. more negative. 

Could you please check the implementation as it appears to be an issue with the Lagrangian multipliers, especially if N>T?

Below is the code I am using as well as the set up of the Gurobi function. 

Thank you very much in advance.

% This is a simple code to reproduce the issue

% Generate test data
N=30;
data = randn(60,N).*repmat(rand(1,N)/10,60,1);

% Construct sample covariance matrix
S = cov(data); 

N = size(data, 2);
ConstrL = 0; % non-negativity constraint

% Run optimization using Gurobi function 
[x,fval,exitflag,output,lambda]=gurobiquadprog(S,zeros(1,N),[],[],ones(1,N),1, repmat(ConstrL,1,N));

% Jagannathan and Ma correction (New JF-Article, no impact on treatment of shortsale constraints)
J = -(lambda.upper*ones(1,N)+ones(N,1)*lambda.upper')...
+(lambda.lower*ones(1,N)+ones(N,1)*lambda.lower');

% Shrink the covariance matrix
SS = S-J;

% Compute portfolio weights
w = (SS\ones(N,1))/(ones(N,1)'/SS*ones(N,1))

% This is the gurobi optimization function
function [x,fval,exitflag,output,lambda] = gurobiquadprog(H,f,A,b,Aeq,beq,lb,ub,x0,options)



% Initialize missing arguments
if nargin < 10
    options = struct();
    if nargin == 9
        if isa(x0,'struct') || isa(x0,'optim.options.SolverOptions')
            options = x0; % x0 was omitted and options were passed instead
            x0 = [];
        end
    else
        x0 = [];
        if nargin < 8
            ub = [];
            if nargin < 7
                lb = [];
                if nargin < 6
                    beq = [];
                    if nargin < 5
                        Aeq = [];
                    end
                end
            end
        end
    end
end

% Warn user if x0 argument is ignored
if ~isempty(x0)
    warning('LINPROG will ignore non-empty starting point X0');
end

% Build Gurobi model
model.obj = f;
model.Q = sparse(H);
model.A = [sparse(A); sparse(Aeq)]; % A must be sparse
model.sense = [repmat('<',size(A,1),1); repmat('=',size(Aeq,1),1)];
model.rhs = full([b(:); beq(:)]); % rhs must be dense
if ~isempty(lb)
    model.lb = lb;
else
    model.lb = -inf(size(model.A,2),1); % default lb for MATLAB is -inf
end
if ~isempty(ub)
    model.ub = ub;
end

% Extract relevant Gurobi parameters from (subset of) options
params = struct();

if isfield(options,'Display') || isa(options,'optim.options.SolverOptions')
    if any(strcmp(options.Display,{'off','none'}))
        params.OutputFlag = 0;
    end
end

if isfield(options,'MaxTime') || isa(options,'optim.options.SolverOptions')
    params.TimeLimit = options.MaxTime;
end

params.outputflag = 0;
params.NonConvex = 2; % avoid issue for non-convex cases but leads to
% matlab freeze
params.PSDTol = inf;

% Solve model with Gurobi
result = gurobi(model,params);

% Resolve model if status is INF_OR_UNBD
if strcmp(result.status,'INF_OR_UNBD')
    params.DualReductions = 0;
    warning('Infeasible or unbounded, resolve without dual reductions to determine...');
    result = gurobi(model,params);
end

% Collect results
x = [];
output.message = result.status;
output.constrviolation = [];

if isfield(result,'x')
    x = result.x;
    if nargout > 3
        slack = model.A*x-model.rhs;
        violA = slack(1:size(A,1));
        violAeq = norm(slack((size(A,1)+1):end),inf);
        viollb = model.lb(:)-x;
        violub = 0;
        if isfield(model,'ub')
            violub = x-model.ub(:);
        end
        output.constrviolation = max([0; violA; violAeq; viollb; violub]);
    end
end

fval = [];

if isfield(result,'objval')
    fval = result.objval;
end

if strcmp(result.status,'OPTIMAL')
    exitflag = 1; % converged to a solution
elseif strcmp(result.status,'UNBOUNDED')
    exitflag = -3; % problem is unbounded
elseif strcmp(result.status,'ITERATION_LIMIT')
    exitflag = 0; % maximum number of iterations reached
else
    exitflag = -2; % no feasible point found
end

lambda.lower = [];
lambda.upper = [];
lambda.ineqlin = [];
lambda.eqlin = [];

if nargout > 4
    if isfield(result,'rc')
        lambda.lower = max(0,result.rc);
        lambda.upper = -min(0,result.rc);
    end
    if isfield(result,'pi')
        if ~isempty(A)
            lambda.ineqlin = -result.pi(1:size(A,1));
        end
        if ~isempty(Aeq)
            lambda.eqlin = -result.pi((size(A,1)+1):end);
        end
    end
end

% Local Functions =========================================================

function [f,A,b,Aeq,beq,lb,ub,x0,options] = probstruct2args(s)
%PROBSTRUCT2ARGS Get problem structure fields ([] is returned when missing)

f = getstructfield(s,'f');
A = getstructfield(s,'Aineq');
b = getstructfield(s,'bineq');
Aeq = getstructfield(s,'Aeq');
beq = getstructfield(s,'beq');
lb = getstructfield(s,'lb');
ub = getstructfield(s,'ub');
x0 = getstructfield(s,'x0');
options = getstructfield(s,'options');

function f = getstructfield(s,field)
%GETSTRUCTFIELD Get structure field ([] is returned when missing)

if isfield(s,field)
f = getfield(s,field);
else
f = [];
end