Coverage for src / cvxcla / first.py: 100%

1# Copyright 2023 Stanford University Convex Optimization Group 

2# 

3# Licensed under the Apache License, Version 2.0 (the "License"); 

4# you may not use this file except in compliance with the License. 

5# You may obtain a copy of the License at 

6# 

7# http://www.apache.org/licenses/LICENSE-2.0 

8# 

9# Unless required by applicable law or agreed to in writing, software 

10# distributed under the License is distributed on an "AS IS" BASIS, 

11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 

12# See the License for the specific language governing permissions and 

13# limitations under the License. 

14"""First turning point computation for the Critical Line Algorithm. 

15 

16This module provides functions to compute the first turning point on the efficient frontier, 

17which is the portfolio with the highest expected return that satisfies the constraints. 

18Two implementations are provided: a direct algorithm and a linear programming approach. 

19""" 

20 

21from __future__ import annotations 

22 

23import numpy as np 

24from numpy.typing import NDArray 

25 

26from .types import TurningPoint 

27 

28 

29# 

30def init_algo( 

31 mean: NDArray[np.float64], lower_bounds: NDArray[np.float64], upper_bounds: NDArray[np.float64] 

32) -> TurningPoint: 

33 """Compute the first turning point. 

34 

35 The key insight behind Markowitz's CLA is to find first the 

36 turning point associated with the highest expected return, and then 

37 compute the sequence of turning points, each with a lower expected 

38 return than the previous. That first turning point consists in the 

39 smallest subset of assets with highest return such that the sum of 

40 their upper boundaries equals or exceeds one. 

41 

42 We sort the expected returns in descending order. 

43 This gives us a sequence for searching for the 

44 first free asset. All weights are initially set to their lower bounds, 

45 and following the sequence from the previous step, we move those 

46 weights from the lower to the upper bound until the sum of weights 

47 exceeds one. If possible the last iterated weight is then reduced 

48 to comply with the constraint that the sum of weights equals one. 

49 This last weight is the first free asset, 

50 and the resulting vector of weights the first turning point. 

51 """ 

52 if np.any(lower_bounds > upper_bounds): 

53 msg = "Lower bounds must be less than or equal to upper bounds" 

54 raise ValueError(msg) 

55 

56 # Initialize weights to lower bounds 

57 weights = np.copy(lower_bounds) 

58 free = np.full_like(mean, False, dtype=np.bool_) 

59 

60 # Move weights from lower to upper bound 

61 # until sum of weights hits or exceeds 1 

62 for index in np.argsort(-mean): 

63 weights[index] += np.min([upper_bounds[index] - lower_bounds[index], 1.0 - np.sum(weights)]) 

64 if np.sum(weights) >= 1: 

65 free[index] = True 

66 break 

67 

68 # free = _free(weights, lower_bounds, upper_bounds) 

69 

70 if not np.any(free): 

71 # # We have not reached the sum of weights of 1... 

72 msg = "Could not construct a fully invested portfolio" 

73 raise ValueError(msg) 

74 

75 # Return first turning point, the point with the highest expected return. 

76 return TurningPoint(free=free, weights=weights) 

77 

78 

79def _free( 

80 w: NDArray[np.float64], lower_bounds: NDArray[np.float64], upper_bounds: NDArray[np.float64] 

81) -> NDArray[np.bool_]: 

82 """Determine which asset should be free in the turning point. 

83 

84 This helper function identifies the asset that should be marked as free 

85 in the turning point. It selects the asset that is furthest from its bounds, 

86 which helps ensure numerical stability in the algorithm. 

87 

88 Args: 

89 w: Vector of portfolio weights. 

90 lower_bounds: Vector of lower bounds for asset weights. 

91 upper_bounds: Vector of upper bounds for asset weights. 

92 

93 Returns: 

94 A boolean vector indicating which asset is free (True) and which are blocked (False). 

95 

96 """ 

97 # Calculate the distance from each weight to its nearest bound 

98 distance = np.min(np.array([np.abs(w - lower_bounds), np.abs(upper_bounds - w)]), axis=0) 

99 

100 # Find the index of the asset furthest from its bounds 

101 index = np.argmax(distance) 

102 

103 # Create a boolean vector with only that asset marked as free 

104 free = np.full_like(w, False, dtype=np.bool_) 

105 free[index] = True 

106 return free