Fix MaximumSharpeRatioPortfolioOptimizer to maximize the Sharpe ratio

Algorithm.Framework/Portfolio/MaximumSharpeRatioPortfolioOptimizer.cs

@@ -14,6 +14,7 @@
 */
 
 using System.Collections.Generic;
+using System.Linq;
 using Accord.Math;
 using Accord.Math.Optimization;
 using Accord.Statistics;
@@ -44,41 +45,36 @@ public MaximumSharpeRatioPortfolioOptimizer(double lower = -1, double upper = 1,
             _riskFreeRate = riskFreeRate;
         }
 
-        /// <summary>
-        /// Sum of all weight is one: 1^T w = 1 / Σw = 1
-        /// </summary>
-        /// <param name="size">number of variables</param>
-        /// <returns>linear constraint object</returns>
-        protected LinearConstraint GetBudgetConstraint(int size)
-        {
-            return new LinearConstraint(size)
-            {
-                CombinedAs = Vector.Create(size, 1.0),
-                ShouldBe = ConstraintType.EqualTo,
-                Value = 1.0
-            };
-        }
-
         /// <summary>
         /// Boundary constraints on weights: lw ≤ w ≤ up
         /// </summary>
+        /// <remarks>
+        /// Expressed in the substituted variable y = κw (κ = 1ᵀy &gt; 0), the per-weight bounds
+        /// become linear: yᵢ − up·(1ᵀy) ≤ 0 and yᵢ − lw·(1ᵀy) ≥ 0.
+        /// </remarks>
         /// <param name="size">number of variables</param>
         /// <returns>enumeration of linear constraint objects</returns>
         protected IEnumerable<LinearConstraint> GetBoundaryConditions(int size)
         {
             for (int i = 0; i < size; i++)
             {
-                yield return new LinearConstraint(1)
+                // yᵢ − up·(1ᵀy) ≤ 0
+                var upper = Vector.Create(size, -_upper);
+                upper[i] += 1.0;
+                yield return new LinearConstraint(size)
                 {
-                    VariablesAtIndices = new int[] { i },
-                    ShouldBe = ConstraintType.GreaterThanOrEqualTo,
-                    Value = _lower
+                    CombinedAs = upper,
+                    ShouldBe = ConstraintType.LesserThanOrEqualTo,
+                    Value = 0.0
                 };
-                yield return new LinearConstraint(1)
+                // yᵢ − lw·(1ᵀy) ≥ 0
+                var lower = Vector.Create(size, -_lower);
+                lower[i] += 1.0;
+                yield return new LinearConstraint(size)
                 {
-                    VariablesAtIndices = new int[] { i },
-                    ShouldBe = ConstraintType.LesserThanOrEqualTo,
-                    Value = _upper
+                    CombinedAs = lower,
+                    ShouldBe = ConstraintType.GreaterThanOrEqualTo,
+                    Value = 0.0
                 };
             }
         }
@@ -96,36 +92,49 @@ public double[] Optimize(double[,] historicalReturns, double[] expectedReturns =
             var returns = (expectedReturns ?? historicalReturns.Mean(0)).Subtract(_riskFreeRate);
 
             var size = covariance.GetLength(0);
-            var x0 = Vector.Create(size, 1.0 / size);
-            var k = returns.Dot(x0);
+            var equalWeights = Vector.Create(size, 1.0 / size);
+
+            // The Charnes-Cooper substitution needs a portfolio with positive expected excess
+            // return to exist, otherwise the Sharpe ratio cannot be maximized.
+            var feasible = _lower >= 0 ? returns.Any(x => x > 0) : returns.Any(x => x != 0);
+            if (!feasible)
+            {
+                return equalWeights;
+            }
 
+            // Charnes-Cooper substitution y = κw (κ = 1ᵀy): maximizing the Sharpe ratio
+            // (µ − r_f)ᵀw / √(wᵀΣw) becomes minimizing wᵀΣw subject to (µ − r_f)ᵀy = 1,
+            // recovering the weights afterwards as w = y / (1ᵀy).
+            // https://quant.stackexchange.com/questions/18521/sharpe-maximization-under-quadratic-constraints
             var constraints = new List<LinearConstraint>
             {
-                // Sharpe Maximization under Quadratic Constraints
-                // https://quant.stackexchange.com/questions/18521/sharpe-maximization-under-quadratic-constraints
-                // (µ − r_f)^T w = k
+                // (µ − r_f)ᵀy = 1
                 new LinearConstraint(size)
                 {
                     CombinedAs = returns,
                     ShouldBe = ConstraintType.EqualTo,
-                    Value = k
+                    Value = 1.0
                 }
             };
 
-            // Σw = 1
-            constraints.Add(GetBudgetConstraint(size));
-
             // lw ≤ w ≤ up
             constraints.AddRange(GetBoundaryConditions(size));
 
-            // Setup solver
+            // Setup solver: minimize yᵀΣy
             var optfunc = new QuadraticObjectiveFunction(covariance, Vector.Create(size, 0.0));
             var solver = new GoldfarbIdnani(optfunc, constraints);
 
             // Solve problem
-            var success = solver.Minimize(Vector.Copy(x0));
-            var sharpeRatio = returns.Dot(solver.Solution) / solver.Value;
-            return success ? solver.Solution : x0;
+            var success = solver.Minimize(Vector.Copy(equalWeights));
+            if (!success)
+            {
+                return equalWeights;
+            }
+
+            // Recover the portfolio weights: w = y / (1ᵀy)
+            var y = solver.Solution;
+            var sum = y.Sum();
+            return sum > 0 ? y.Divide(sum) : equalWeights;
         }
     }
 }

Algorithm.Framework/Portfolio/MaximumSharpeRatioPortfolioOptimizer.py

@@ -55,24 +55,23 @@ def optimize(self, historical_returns, expected_returns = None, covariance = Non
 
         size = covariance.columns.size   # K x 1
         x0 = np.array(size * [1. / size])
-        k = expected_returns.dot(x0)
 
-        # Sharpe Maximization under Quadratic Constraints
+        # SLSQP maximizes the Sharpe ratio (µ − r_f)^T w / √(w^T Σ w) directly, so the fractional
+        # objective is optimized in place without any substitution. The budget constraint Σw = 1 and
+        # the per-weight bounds lw ≤ w ≤ up are applied as-is. The previous implementation instead
+        # fixed (µ − r_f)^T w to the equal-weight return, which collapsed the optimizer to minimum
+        # variance. The C# implementation uses the Charnes-Cooper QP substitution because its solver
+        # only handles quadratic objectives.
         # https://quant.stackexchange.com/questions/18521/sharpe-maximization-under-quadratic-constraints
-        # (µ − r_f)^T w = k
         constraints = [
-            {'type': 'eq', 'fun': lambda weights: expected_returns.dot(weights) - k}]
+            # Σw = 1
+            {'type': 'eq', 'fun': lambda weights: self.get_budget_constraint(weights)}]
 
-        # Σw = 1
-        constraints.append(
-            {'type': 'eq', 'fun': lambda weights: self.get_budget_constraint(weights)})
-
-        opt = minimize(lambda weights: self.portfolio_variance(weights, covariance),   # Objective function
+        opt = minimize(lambda weights: -expected_returns.dot(weights) / np.sqrt(self.portfolio_variance(weights, covariance)),   # Objective function: −Sharpe ratio
                        x0,                                                        # Initial guess
                        bounds = self.get_boundary_conditions(size),               # Bounds for variables: lw ≤ w ≤ up
                        constraints = constraints,                                 # Constraints definition
                        method='SLSQP')        # Optimization method:  Sequential Least SQuares Programming
-        sharpe_ratio = expected_returns.dot(opt['x']) / opt.fun
 
         return opt['x'] if opt['success'] else x0
 

Tests/Algorithm/Framework/Portfolio/MaximumSharpeRatioPortfolioOptimizerTests.cs

@@ -65,14 +65,14 @@ public void SetUp()
 
             ExpectedResults = new List<double[]>
             {
-                new double[] { -0.562396, 0.608942, 0.953453 },
-                new double[] { 0.686025, -0.269589, 0.583023 },
-                new double[] { 0.26394, -0.043374, 0.779434 },
-                new double[] { -0.223905, 0.401036, 1, 0.065329, -0.24246 },
-                new double[] { 0.5, 0.5 },
                 new double[] { -0.5, 0.5, 1 },
-                new double[] { -0.242647, 1, 0.242647 },
-                new double[] { -1, 0.922902, 0.364512, 0.712585 },
+                new double[] { 0, 0, 1 },
+                new double[] { -0.404692, 0.404692, 1 },
+                new double[] { -0.418338, 0.023261, 1, 0.040668, 0.35441 },
+                new double[] { 0.5, 0.5 },
+                new double[] { -0.670213, 0.670213, 1 },
+                new double[] { -1, 1, 1 },
+                new double[] { -1, 0.315476, 0.684524, 1 },
             };
         }
 
@@ -98,21 +98,22 @@ public override void OptimizeWeightings(int testCaseNumber)
         public void OptimizeWeightingsSpecifyingLowerBoundAndRiskFreeRate(int testCaseNumber)
         {
             var testOptimizer = new MaximumSharpeRatioPortfolioOptimizer(lower: 0, riskFreeRate: 0.04);
-            var expectedResult = new double[] { 0, 0.44898, 0.55102 };
+            var expectedResult = new double[] { 0, 0, 1 };
 
             var result = testOptimizer.Optimize(HistoricalReturns[testCaseNumber]);
 
             Assert.AreEqual(expectedResult, result.Select(x => Math.Round(x, 6)));
         }
 
         [Test]
-        public void SingleSecurityPortfolioReturnsNaN()
+        public void SingleSecurityPortfolioReturnsFullWeight()
         {
             var testOptimizer = new MaximumSharpeRatioPortfolioOptimizer();
             var historicalReturns = new double[,] { { -0.1 } };
             var expectedReturns = new double[] { -0.1 };
 
-            var expectedResult = new double[] { double.NaN };
+            // With a single security the budget constraint Σw = 1 leaves it fully invested
+            var expectedResult = new double[] { 1 };
 
             var result = testOptimizer.Optimize(historicalReturns, expectedReturns);