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

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 raise ValueError("Lower bounds must be less than or equal to upper bounds") 

54 

55 # Initialize weights to lower bounds 

56 weights = np.copy(lower_bounds) 

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

58 

59 # Move weights from lower to upper bound 

60 # until sum of weights hits or exceeds 1 

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

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

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

64 free[index] = True 

65 break 

66 

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

68 

69 if not np.any(free): 

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

71 raise ValueError("Could not construct a fully invested portfolio") 

72 

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

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

75 

76 

77def _free( 

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

79) -> NDArray[np.bool_]: 

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

81 

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

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

84 which helps ensure numerical stability in the algorithm. 

85 

86 Args: 

87 w: Vector of portfolio weights. 

88 lower_bounds: Vector of lower bounds for asset weights. 

89 upper_bounds: Vector of upper bounds for asset weights. 

90 

91 Returns: 

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

93 

94 """ 

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

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

97 

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

99 index = np.argmax(distance) 

100 

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

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

103 free[index] = True 

104 return free