Using 'fmincon' for Portfolio optimization with nonlinear constraint

Hi, I have get stuck with a problem where I try to maximize the Sharpe Ratio for a portfolio, but with the linear constraint that each asset maximum weight is 10%, but also that the sum of weights over 5% can't exceed 40% (UTICS III). This is the reason I haven't used quadprog() or similar solver (that only deal with linear constraints).

Right now I am trying to write a function that estimates optimal level of lambda/riskaversion using fminbnd(), while also using fmincon() to find optimal level of weights for a vector pwgt (which should have length N). I haven't included any more constraints than lower- and upper bound in the code below yet.

My questions are:

1. How do I get Matlab to interpret "pwgt" as a vector of length N and not a scalar?

2. Any suggestion of how I should write the 5/40% constraint? I am thinking of something like:

function [c,ceq] = nonlincon(pwgt, obj)
c=[];
r=[];
Q=[];
r=obj.AssetMean;
Q=obj.AssetCovar;
c=[];
%Nonlinear inequality constraints
for i=1:26
    if pwgt(i) > 0.05
    c(i)=pwgt(i);
    else
    c(i)=0;
    end
end

% Nonlinear equality constraints
ceq = [];
c=sum(c(1,1:length(r)))-0.40;

The code for the whole function:

function [lambdaopt]=maxSharpeUtics(obj)

fhandle=@(lambda) local_objective(lambda, obj);

options = optimset('fminbnd');
options = optimset(options, 'Display', 'Iter', 'TolX', 1.0e-10);
[lambdaopt, ~, exitflag] = fminbnd(fhandle, 0, 3, options);

[lambdaopt]=lambdaopt;

function sratio = local_objective(lambda, obj)
r=obj.AssetMean;
Q=obj.AssetCovar;
N=length(r);
lb=zeros(N,1);
ub(1:N,1)=0.10;
x0=1/obj.NumAssets;
pwgt=[];
f = r'*pwgt-lambda*(pwgt'*Q*pwgt);

fhand = @(pwgt) f;
[pwgt]=fmincon(fhand,x0,lb,ub);

ret=r'*pwgt*12;
s=sqrt((pwgt'*Q*pwgt*12));

if ~isempty(obj.RiskFreeRate)
  sratio = -(ret - obj.RiskFreeRate)/s;
else
  sratio = -ret/s;
end

"obj" is a Portfolio object which contains the means and covariances needed to solve the problem.

I am using R2014b.

Anything else that I have done wrong in writing the functions?

Many thanks in advance.

Br, Sten

Answer 1

[Matt J]

1. How do I get Matlab to interpret "pwgt" as a vector of length N and not a scalar?

If your initial guess of pwgt is of length N, fmincon (or any other Optimization Toolbox solver) will know that it is length N and not a scalar.

2. Any suggestion of how I should write the 5/40% constraint? I am thinking of something like:

You can't do it that way, because it's a highly non-smooth constraint. I think you will have to do a combinatoric search, i.e., loop through all subsets of pwgt(i) and run the solver on a different sub-problem for each subset. In each sub-problem, you will impose a lower bound of 5% on the fixed subset, an upper bound of 5% on pwgt(i), not in the subset, and an inequality constraint on the sum of the weights in the subset. Decomposing the problem this way, I guess you might be able use quadprog...

Answer 2

[Johan Löfberg]

Solving this using fmincon sounds like the wrong tool.

First, the Sharpe ratio can be rewritten as a convex quadratic by suitable redefinition of the weights. The original nonlinear form is nonconvex so it should be avoided as far as possible. Using nonlinear solvers for this simple problem is a common mistake, so I added some info about this in the portfolio example in the YALMIP Wiki (a toolbox for optimization)

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Examples.Portfolio

With that reformulation, you can use quadprog instead.

You constraint on maximal weight will be z_i <= sum(z)*.1, while the other set of constraints most likely are best to implement using standard integer programming ideas. If you are lazy and want to implement it using YALMIP, post a question on the YALMIP Google groups and I might cook something up for you

[Sten]

Thanks for your input Johan. Yes, I have found that fmincon() is the wrong tool in this case. One solution I have found is using intlinprog() to estimate the optimal portfolio for a given risk aversion/lambda in this way:

http://se.mathworks.com/help/optim/examples/mixed-integer-quadratic-programming-portfolio-optimization.html

With this method I can implement my UTICS 5/40% constraint and it works very well. To find the maximum Sharpe-ratio I, above this, tried to use fminbnd() in this way:

function [pwgt, lambdaopt, fval] = maximizeSharpeUtics(obj)

% minimize the local objective function to obtain a value on lambda that maximizes the Sharpe ratio

fhandle = @(lambda) local_objective(lambda, obj);

options = optimset('fminbnd');
options = optimset(options, 'Display', 'Iter', 'TolX', 1.0e-2); 
[lambdaopt, fval, exitflag] = fminbnd(fhandle, 0, 100, options);

[fval]=fval;
[lamdaopt]=lambdaopt;
[pwgt]=EstimateFrontierByLambda(obj, lambdaopt);

end

function sratio = local_objective(lambda, obj)

r=obj.AssetMean;
Q=obj.AssetCovar;

[fwgt]=EstimateFrontierByLambda(obj, lambda);

ret=r'*fwgt*12;
s=sqrt((fwgt'*Q*fwgt*12)); 

if ~isempty(obj.RiskFreeRate)
  sratio = -(ret - obj.RiskFreeRate)/s;
else
  sratio = -ret/s;
end
end

However, the problem here is (as you say) that the Sharpe-ratio in this form is nonconvex and the fminbnd()-optimizer fails to find the maximum Sharpe-ratio.

The "EstimateFrontierByLambda(obj, lambda);" is my intlinprog-optimizer with all constraints included. "obj" is a Portfolio-object from the Financial Toolbox.

Could you please explain the reformulation of the Sharpe-ratio as you mention? I don't really get the intuition. Would the reformulation let me use fminbnd to find a value at lambda that maximizes the Sharpe-ratio?

Many thanks & br,

Sten

[Johan Löfberg]

Why are you still talking about general nonlinear solvers (fminbnd) when I've told you it is possible to formulate as a convex QP (possibly with mixed-integer extensions leading to MIQP). Also, solving a MIQP using intlinprog is really a band-aid solution, motivated entirely from the fact that mathworks do not ship a MIQP solver with the optimization toolbox.

The intuition is simple. sum(f)=1 is entirely arbitrary. You simply want to say that f should not be zero, and all your constraints are simply related to the relation of the individual elements. You can find a solid derivation in e.g. www.ie.bilkent.edu.tr/~mustafap/courses/OIF.pdf

When you start adding combinatorial 5/40 constraints to the convex QP Sharpe model, you will need a MIQP solver. If you are academia, you have all the main players available for free, such as gurobi, cplex, mosek. If not, scip probably works ok.

To summarize, if you go for a MIQP appraoch (which is the only reasonable approach in my opinion) and model it using YALMIP, it will be a pretty trivial problem to setup (post a question on the YALMIP forum and I'll help you set up the logic model)