Coverage for src / cvx / risk / sample / sample.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"""Risk models based on the sample covariance matrix. 

15 

16This module provides the SampleCovariance class, which implements a risk model 

17based on the Cholesky decomposition of the sample covariance matrix. This is 

18one of the most common approaches to portfolio risk estimation. 

19 

20Example: 

21 Create and use a sample covariance risk model: 

22 

23 >>> import numpy as np 

24 >>> from cvx.risk.sample import SampleCovariance 

25 >>> # Create risk model for up to 3 assets 

26 >>> model = SampleCovariance(num=3) 

27 >>> # Update with a covariance matrix 

28 >>> cov = np.array([[1.0, 0.5, 0.0], [0.5, 1.0, 0.5], [0.0, 0.5, 1.0]]) 

29 >>> model.update( 

30 ... cov=cov, 

31 ... lower_assets=np.zeros(3), 

32 ... upper_assets=np.ones(3) 

33 ... ) 

34 >>> # Estimate risk for a given portfolio 

35 >>> weights = np.array([0.4, 0.3, 0.3]) 

36 >>> risk = model.estimate(weights) 

37 >>> isinstance(risk, float) 

38 True 

39 

40""" 

41 

42from __future__ import annotations 

43 

44from dataclasses import dataclass 

45from typing import Any 

46 

47import clarabel 

48import numpy as np 

49from cvx.linalg import cholesky 

50from scipy import sparse 

51 

52from cvx.core import Bounds, Model, Parameter, Variable 

53 

54 

55@dataclass 

56class SampleCovariance(Model): 

57 """Risk model based on the Cholesky decomposition of the sample covariance matrix. 

58 

59 This model computes portfolio risk as the L2 norm of the product of the 

60 Cholesky factor and the weights vector. Mathematically, if R is the upper 

61 triangular Cholesky factor of the covariance matrix (R^T @ R = cov), then: 

62 

63 risk = ||R @ w||_2 = sqrt(w^T @ cov @ w) 

64 

65 This represents the portfolio standard deviation (volatility). 

66 

67 Attributes: 

68 num: Maximum number of assets the model can handle. The model can be 

69 updated with fewer assets, but not more. 

70 

71 Example: 

72 Basic usage: 

73 

74 >>> import numpy as np 

75 >>> from cvx.risk.sample import SampleCovariance 

76 >>> model = SampleCovariance(num=2) 

77 >>> model.update( 

78 ... cov=np.array([[1.0, 0.5], [0.5, 2.0]]), 

79 ... lower_assets=np.zeros(2), 

80 ... upper_assets=np.ones(2) 

81 ... ) 

82 >>> # Equal weight portfolio 

83 >>> weights = np.array([0.5, 0.5]) 

84 >>> risk = model.estimate(weights) 

85 >>> # Risk should be sqrt(0.5^2 * 1 + 0.5^2 * 2 + 2 * 0.5 * 0.5 * 0.5) 

86 >>> bool(np.isclose(risk, 1.0)) 

87 True 

88 

89 Using in optimization: 

90 

91 >>> from cvx.risk.portfolio import minrisk_problem 

92 >>> from cvx.core.variable import Variable 

93 >>> weights = Variable(2) 

94 >>> problem = minrisk_problem(model, weights) 

95 >>> problem.solve() 

96 >>> # Lower variance asset gets higher weight 

97 >>> bool(weights.value[0] > weights.value[1]) 

98 True 

99 

100 Mathematical verification - the risk estimate equals sqrt(w^T @ cov @ w): 

101 

102 >>> model = SampleCovariance(num=3) 

103 >>> cov = np.array([[0.04, 0.01, 0.02], 

104 ... [0.01, 0.09, 0.01], 

105 ... [0.02, 0.01, 0.16]]) 

106 >>> model.update( 

107 ... cov=cov, 

108 ... lower_assets=np.zeros(3), 

109 ... upper_assets=np.ones(3) 

110 ... ) 

111 >>> w = np.array([0.4, 0.35, 0.25]) 

112 >>> # Model estimate 

113 >>> model_risk = model.estimate(w) 

114 >>> # Manual calculation: sqrt(w^T @ cov @ w) 

115 >>> manual_risk = np.sqrt(w @ cov @ w) 

116 >>> bool(np.isclose(model_risk, manual_risk, rtol=1e-6)) 

117 True 

118 

119 """ 

120 

121 num: int = 0 

122 """Maximum number of assets the model can handle.""" 

123 

124 def __post_init__(self) -> None: 

125 """Initialize the parameters after the class is instantiated. 

126 

127 Creates the Cholesky decomposition parameter and initializes the bounds. 

128 The Cholesky parameter is a square matrix of size (num, num), and bounds 

129 are created for asset weights. 

130 

131 Example: 

132 >>> from cvx.risk.sample import SampleCovariance 

133 >>> model = SampleCovariance(num=5) 

134 >>> # Parameters are automatically created 

135 >>> model.parameter["chol"].shape 

136 (5, 5) 

137 

138 """ 

139 self.parameter["chol"] = Parameter( 

140 shape=(self.num, self.num), 

141 name="cholesky of covariance", 

142 ) 

143 self.bounds = Bounds(m=self.num, name="assets") 

144 

145 def estimate(self, weights: np.ndarray, **kwargs: Any) -> float: 

146 """Estimate the portfolio risk using the Cholesky decomposition. 

147 

148 Computes the L2 norm of the product of the Cholesky factor and the 

149 weights vector. This is equivalent to the square root of the portfolio 

150 variance (i.e., portfolio volatility). 

151 

152 Args: 

153 weights: Numpy array representing portfolio weights. 

154 **kwargs: Additional keyword arguments (not used). 

155 

156 Returns: 

157 Float representing the portfolio risk (standard deviation). 

158 

159 Example: 

160 >>> import numpy as np 

161 >>> from cvx.risk.sample import SampleCovariance 

162 >>> model = SampleCovariance(num=2) 

163 >>> # Identity covariance (uncorrelated assets with unit variance) 

164 >>> model.update( 

165 ... cov=np.eye(2), 

166 ... lower_assets=np.zeros(2), 

167 ... upper_assets=np.ones(2) 

168 ... ) 

169 >>> risk = model.estimate(np.array([0.5, 0.5])) 

170 >>> isinstance(risk, float) 

171 True 

172 

173 """ 

174 return float(np.linalg.norm(self.parameter["chol"].value @ np.asarray(weights))) 

175 

176 def update(self, **kwargs: Any) -> None: 

177 """Update the Cholesky decomposition parameter and bounds. 

178 

179 Computes the Cholesky decomposition of the provided covariance matrix 

180 and updates the model parameters. The covariance matrix can be smaller 

181 than num x num. 

182 

183 Args: 

184 **kwargs: Keyword arguments containing: 

185 - cov: Covariance matrix (numpy.ndarray). Must be positive definite. 

186 - lower_assets: Array of lower bounds for asset weights. 

187 - upper_assets: Array of upper bounds for asset weights. 

188 

189 Example: 

190 >>> import numpy as np 

191 >>> from cvx.risk.sample import SampleCovariance 

192 >>> model = SampleCovariance(num=5) 

193 >>> # Update with a 3x3 covariance (smaller than max) 

194 >>> cov = np.array([[1.0, 0.3, 0.1], 

195 ... [0.3, 1.0, 0.2], 

196 ... [0.1, 0.2, 1.0]]) 

197 >>> model.update( 

198 ... cov=cov, 

199 ... lower_assets=np.zeros(3), 

200 ... upper_assets=np.ones(3) 

201 ... ) 

202 >>> # Cholesky factor is updated 

203 >>> model.parameter["chol"].value[:3, :3].shape 

204 (3, 3) 

205 

206 """ 

207 cov = kwargs["cov"] 

208 n = cov.shape[0] 

209 

210 chol = np.zeros((self.num, self.num)) 

211 chol[:n, :n] = cholesky(cov) 

212 self.parameter["chol"].value = chol 

213 self.bounds.update(**kwargs) 

214 

215 def solve_minrisk( 

216 self, 

217 weights: Variable, 

218 base: np.ndarray, 

219 extra_constraints: list[tuple[np.ndarray, float | None, float | None]], 

220 y_var: Variable | None = None, 

221 ) -> tuple[float | None, float | None, str]: 

222 """Build and solve the Clarabel SOC problem for this model.""" 

223 n = weights.n 

224 chol = self.parameter["chol"].value 

225 lb, ub = self.bounds.get_bounds() 

226 

227 n_vars = 1 + n 

228 P = sparse.csc_matrix((n_vars, n_vars)) # noqa: N806 

229 q = np.zeros(n_vars) 

230 q[0] = 1.0 

231 

232 A_rows: list[np.ndarray] = [] # noqa: N806 

233 b_rows: list[np.ndarray] = [] 

234 cones: list[Any] = [] 

235 

236 A_soc = np.zeros((n + 1, n_vars)) # noqa: N806 

237 A_soc[0, 0] = -1.0 

238 A_soc[1:, 1:] = -chol 

239 b_soc = np.zeros(n + 1) 

240 b_soc[1:] = -chol @ base 

241 A_rows.append(A_soc) 

242 b_rows.append(b_soc) 

243 cones.append(clarabel.SecondOrderConeT(n + 1)) 

244 

245 A_eq = np.zeros((1, n_vars)) # noqa: N806 

246 A_eq[0, 1:] = 1.0 

247 A_rows.append(A_eq) 

248 b_rows.append(np.array([1.0])) 

249 cones.append(clarabel.ZeroConeT(1)) 

250 

251 A_lb = np.zeros((n, n_vars)) # noqa: N806 

252 A_lb[:, 1:] = -np.eye(n) 

253 A_rows.append(A_lb) 

254 b_rows.append(-lb) 

255 cones.append(clarabel.NonnegativeConeT(n)) 

256 

257 A_ub = np.zeros((n, n_vars)) # noqa: N806 

258 A_ub[:, 1:] = np.eye(n) 

259 A_rows.append(A_ub) 

260 b_rows.append(ub) 

261 cones.append(clarabel.NonnegativeConeT(n)) 

262 

263 for a, lb_val, ub_val in extra_constraints: 

264 a = np.asarray(a) 

265 if lb_val is not None and ub_val is not None and lb_val == ub_val: 

266 A_extra = np.zeros((1, n_vars)) # noqa: N806 

267 A_extra[0, 1:] = a 

268 A_rows.append(A_extra) 

269 b_rows.append(np.array([lb_val])) 

270 cones.append(clarabel.ZeroConeT(1)) 

271 else: 

272 if lb_val is not None: 

273 A_extra = np.zeros((1, n_vars)) # noqa: N806 

274 A_extra[0, 1:] = -a 

275 A_rows.append(A_extra) 

276 b_rows.append(np.array([-lb_val])) 

277 cones.append(clarabel.NonnegativeConeT(1)) 

278 if ub_val is not None: 

279 A_extra = np.zeros((1, n_vars)) # noqa: N806 

280 A_extra[0, 1:] = a 

281 A_rows.append(A_extra) 

282 b_rows.append(np.array([ub_val])) 

283 cones.append(clarabel.NonnegativeConeT(1)) 

284 

285 A = sparse.csc_matrix(np.vstack(A_rows)) # noqa: N806 

286 b = np.concatenate(b_rows) 

287 

288 settings = clarabel.DefaultSettings() 

289 settings.verbose = False 

290 sol = clarabel.DefaultSolver(P, q, A, b, cones, settings).solve() 

291 status = str(sol.status) 

292 

293 if "Solved" in status: 

294 weights.value = np.array(sol.x[1 : 1 + n]) 

295 return float(sol.obj_val), float(sol.x[0]), status 

296 return None, None, status