Coverage for src / cvx / risk / factor / factor.py: 100%
131 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-13 06:46 +0000
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"""Factor risk model.
16This module provides the FactorModel class, which implements a factor-based
17risk model for portfolio optimization. Factor models decompose portfolio risk
18into systematic (factor) risk and idiosyncratic (residual) risk.
20Example:
21 Create a factor model and estimate portfolio risk:
23 >>> import numpy as np
24 >>> from cvx.risk.factor import FactorModel
25 >>> # Create factor model with 10 assets and 3 factors
26 >>> model = FactorModel(assets=10, k=3)
27 >>> # Set up factor exposure and covariance
28 >>> np.random.seed(42)
29 >>> exposure = np.random.randn(3, 10) # 3 factors x 10 assets
30 >>> factor_cov = np.eye(3) # Factor covariance matrix
31 >>> idio_risk = np.abs(np.random.randn(10)) # Idiosyncratic risk
32 >>> model.update(
33 ... exposure=exposure,
34 ... cov=factor_cov,
35 ... idiosyncratic_risk=idio_risk,
36 ... lower_assets=np.zeros(10),
37 ... upper_assets=np.ones(10),
38 ... lower_factors=-0.1 * np.ones(3),
39 ... upper_factors=0.1 * np.ones(3)
40 ... )
41 >>> # Model is ready for optimization
42 >>> w = np.zeros(10)
43 >>> w[:5] = 0.2
44 >>> risk = model.estimate(w)
45 >>> isinstance(risk, float)
46 True
48"""
50from __future__ import annotations
52from dataclasses import dataclass
53from typing import Any
55import clarabel
56import numpy as np
57from cvx.linalg import cholesky
58from scipy import sparse
60from cvx.core import Bounds, Model, Parameter, Variable
63@dataclass
64class FactorModel(Model):
65 """Factor risk model for portfolio optimization.
67 Factor models decompose portfolio risk into systematic risk (from factor
68 exposures) and idiosyncratic risk (residual risk). The total portfolio
69 variance is:
71 Var(w) = w' @ exposure' @ cov @ exposure @ w + sum((idio_risk * w)^2)
73 This implementation uses the Cholesky decomposition of the factor covariance
74 matrix for efficient risk computation.
76 Attributes:
77 assets: Maximum number of assets in the portfolio.
78 k: Maximum number of factors in the model.
80 Example:
81 Create and use a factor model:
83 >>> import numpy as np
84 >>> from cvx.risk.factor import FactorModel
85 >>> # Create model
86 >>> model = FactorModel(assets=5, k=2)
87 >>> # Factor exposure: 2 factors x 5 assets
88 >>> exposure = np.array([[1.0, 0.8, 0.6, 0.4, 0.2],
89 ... [0.2, 0.4, 0.6, 0.8, 1.0]])
90 >>> # Factor covariance
91 >>> factor_cov = np.array([[1.0, 0.3], [0.3, 1.0]])
92 >>> # Idiosyncratic risk per asset
93 >>> idio_risk = np.array([0.1, 0.1, 0.1, 0.1, 0.1])
94 >>> model.update(
95 ... exposure=exposure,
96 ... cov=factor_cov,
97 ... idiosyncratic_risk=idio_risk,
98 ... lower_assets=np.zeros(5),
99 ... upper_assets=np.ones(5),
100 ... lower_factors=-0.5 * np.ones(2),
101 ... upper_factors=0.5 * np.ones(2)
102 ... )
103 >>> w = np.array([0.2, 0.2, 0.2, 0.2, 0.2])
104 >>> risk = model.estimate(w)
105 >>> isinstance(risk, float)
106 True
108 Mathematical verification of risk decomposition:
110 >>> model = FactorModel(assets=3, k=2)
111 >>> # Factor exposure: how much each asset is exposed to each factor
112 >>> exposure = np.array([[1.0, 0.5, 0.0], # Market factor
113 ... [0.0, 0.5, 1.0]]) # Sector factor
114 >>> # Factor covariance (diagonal = uncorrelated factors)
115 >>> factor_cov = np.array([[0.04, 0.0], # Market vol = 20%
116 ... [0.0, 0.0225]]) # Sector vol = 15%
117 >>> # Idiosyncratic risk per asset
118 >>> idio = np.array([0.10, 0.12, 0.08])
119 >>> model.update(
120 ... exposure=exposure,
121 ... cov=factor_cov,
122 ... idiosyncratic_risk=idio,
123 ... lower_assets=np.zeros(3),
124 ... upper_assets=np.ones(3),
125 ... lower_factors=-np.ones(2),
126 ... upper_factors=np.ones(2)
127 ... )
128 >>> # Equal weight portfolio
129 >>> w = np.array([1/3, 1/3, 1/3])
130 >>> model_risk = model.estimate(w)
131 >>> # Manual: total_var = y^T @ cov @ y + sum((idio * w)^2)
132 >>> y = exposure @ w # Factor exposures
133 >>> systematic_var = y @ factor_cov @ y
134 >>> idio_var = np.sum((idio * w)**2)
135 >>> manual_risk = np.sqrt(systematic_var + idio_var)
136 >>> bool(np.isclose(model_risk, manual_risk, rtol=1e-5))
137 True
139 Error handling for dimension violations:
141 >>> model = FactorModel(assets=3, k=2)
142 >>> try:
143 ... model.update(
144 ... exposure=np.random.randn(5, 3), # 5 factors > k=2
145 ... cov=np.eye(5),
146 ... idiosyncratic_risk=np.ones(3),
147 ... lower_assets=np.zeros(3),
148 ... upper_assets=np.ones(3),
149 ... lower_factors=-np.ones(5),
150 ... upper_factors=np.ones(5)
151 ... )
152 ... except ValueError as e:
153 ... print("Caught:", str(e))
154 Caught: Too many factors
156 """
158 assets: int = 0
159 """Maximum number of assets in the portfolio."""
161 k: int = 0
162 """Maximum number of factors in the model."""
164 def __post_init__(self) -> None:
165 """Initialize the parameters after the class is instantiated.
167 Creates parameters for factor exposure, idiosyncratic risk, and the Cholesky
168 decomposition of the factor covariance matrix. Also initializes bounds for
169 both assets and factors.
171 Example:
172 >>> from cvx.risk.factor import FactorModel
173 >>> model = FactorModel(assets=10, k=3)
174 >>> # Parameters are automatically created
175 >>> model.parameter["exposure"].shape
176 (3, 10)
177 >>> model.parameter["idiosyncratic_risk"].shape
178 10
179 >>> model.parameter["chol"].shape
180 (3, 3)
182 """
183 self.parameter["exposure"] = Parameter(
184 shape=(self.k, self.assets),
185 name="exposure",
186 )
188 self.parameter["idiosyncratic_risk"] = Parameter(shape=self.assets, name="idiosyncratic risk")
190 self.parameter["chol"] = Parameter(
191 shape=(self.k, self.k),
192 name="cholesky of covariance",
193 )
195 self.bounds_assets = Bounds(m=self.assets, name="assets")
196 self.bounds_factors = Bounds(m=self.k, name="factors")
198 def estimate(self, weights: np.ndarray, **kwargs: Any) -> float:
199 """Compute the total portfolio risk using the factor model.
201 Combines systematic risk (from factor exposures) and idiosyncratic risk
202 to calculate the total portfolio risk. The formula is:
204 risk = sqrt(||chol @ y||^2 + ||idio_risk * w||^2)
206 where y = exposure @ weights (factor exposures).
208 Args:
209 weights: Numpy array representing portfolio weights.
210 **kwargs: Additional keyword arguments, may include:
211 - y: Factor exposures as a numpy array. If not provided, calculated
212 as exposure @ weights.
214 Returns:
215 Float representing the total portfolio risk.
217 Example:
218 >>> import numpy as np
219 >>> from cvx.risk.factor import FactorModel
220 >>> model = FactorModel(assets=3, k=2)
221 >>> model.update(
222 ... exposure=np.array([[1.0, 0.5, 0.0], [0.0, 0.5, 1.0]]),
223 ... cov=np.eye(2),
224 ... idiosyncratic_risk=np.array([0.1, 0.1, 0.1]),
225 ... lower_assets=np.zeros(3),
226 ... upper_assets=np.ones(3),
227 ... lower_factors=-np.ones(2),
228 ... upper_factors=np.ones(2)
229 ... )
230 >>> w = np.array([0.4, 0.3, 0.3])
231 >>> risk = model.estimate(w)
232 >>> isinstance(risk, float)
233 True
235 """
236 w = np.asarray(weights)
237 y = np.asarray(kwargs.get("y", self.parameter["exposure"].value @ w))
239 var_systematic = np.linalg.norm(self.parameter["chol"].value @ y)
240 var_residual = np.linalg.norm(self.parameter["idiosyncratic_risk"].value * w)
242 return float(np.sqrt(var_systematic**2 + var_residual**2))
244 def update(self, **kwargs: Any) -> None:
245 """Update the factor model parameters.
247 Updates the factor exposure matrix, idiosyncratic risk vector, and
248 factor covariance Cholesky decomposition. The input dimensions can
249 be smaller than the maximum dimensions.
251 Args:
252 **kwargs: Keyword arguments containing:
253 - exposure: Factor exposure matrix (k x assets).
254 - idiosyncratic_risk: Vector of idiosyncratic risks.
255 - cov: Factor covariance matrix.
256 - lower_assets: Array of lower bounds for asset weights.
257 - upper_assets: Array of upper bounds for asset weights.
258 - lower_factors: Array of lower bounds for factor exposures.
259 - upper_factors: Array of upper bounds for factor exposures.
261 Raises:
262 ValueError: If number of factors or assets exceeds maximum.
264 Example:
265 >>> import numpy as np
266 >>> from cvx.risk.factor import FactorModel
267 >>> model = FactorModel(assets=5, k=3)
268 >>> # Update with 2 factors and 4 assets
269 >>> model.update(
270 ... exposure=np.random.randn(2, 4),
271 ... cov=np.eye(2),
272 ... idiosyncratic_risk=np.abs(np.random.randn(4)),
273 ... lower_assets=np.zeros(4),
274 ... upper_assets=np.ones(4),
275 ... lower_factors=-np.ones(2),
276 ... upper_factors=np.ones(2)
277 ... )
279 """
280 self.parameter["exposure"].value = np.zeros((self.k, self.assets))
281 self.parameter["chol"].value = np.zeros((self.k, self.k))
282 self.parameter["idiosyncratic_risk"].value = np.zeros(self.assets)
284 # get the exposure
285 exposure = kwargs["exposure"]
287 # extract dimensions
288 k, assets = exposure.shape
289 if k > self.k:
290 msg = "Too many factors"
291 raise ValueError(msg)
292 if assets > self.assets:
293 msg = "Too many assets"
294 raise ValueError(msg)
296 self.parameter["exposure"].value[:k, :assets] = kwargs["exposure"]
297 self.parameter["idiosyncratic_risk"].value[:assets] = kwargs["idiosyncratic_risk"]
298 self.parameter["chol"].value[:k, :k] = cholesky(kwargs["cov"])
299 self.bounds_assets.update(**kwargs)
300 self.bounds_factors.update(**kwargs)
302 def solve_minrisk(
303 self,
304 weights: Variable,
305 base: np.ndarray,
306 extra_constraints: list[tuple[np.ndarray, float | None, float | None]],
307 y_var: Variable | None = None,
308 ) -> tuple[float | None, float | None, str]:
309 """Build and solve the Clarabel SOC problem for this model."""
310 n = weights.n
311 k = self.k
313 chol = self.parameter["chol"].value
314 exposure = self.parameter["exposure"].value
315 idio = self.parameter["idiosyncratic_risk"].value
317 lb_w, ub_w = self.bounds_assets.get_bounds()
318 lb_y, ub_y = self.bounds_factors.get_bounds()
320 n_vars = 1 + n + k
321 P = sparse.csc_matrix((n_vars, n_vars)) # noqa: N806
322 q = np.zeros(n_vars)
323 q[0] = 1.0
325 A_rows: list[np.ndarray] = [] # noqa: N806
326 b_rows: list[np.ndarray] = []
327 cones: list[Any] = []
329 soc_size = 1 + k + n
330 A_soc = np.zeros((soc_size, n_vars)) # noqa: N806
331 A_soc[0, 0] = -1.0
332 A_soc[1 : 1 + k, 1 + n :] = -chol
333 A_soc[1 + k :, 1 : 1 + n] = -np.diag(idio)
334 b_soc = np.zeros(soc_size)
335 b_soc[1 + k :] = -idio * base
336 A_rows.append(A_soc)
337 b_rows.append(b_soc)
338 cones.append(clarabel.SecondOrderConeT(soc_size))
340 A_eq_sum = np.zeros((1, n_vars)) # noqa: N806
341 A_eq_sum[0, 1 : 1 + n] = 1.0
342 A_rows.append(A_eq_sum)
343 b_rows.append(np.array([1.0]))
344 cones.append(clarabel.ZeroConeT(1))
346 A_eq_exp = np.zeros((k, n_vars)) # noqa: N806
347 A_eq_exp[:, 1 : 1 + n] = -exposure
348 A_eq_exp[:, 1 + n :] = np.eye(k)
349 A_rows.append(A_eq_exp)
350 b_rows.append(np.zeros(k))
351 cones.append(clarabel.ZeroConeT(k))
353 A_wlb = np.zeros((n, n_vars)) # noqa: N806
354 A_wlb[:, 1 : 1 + n] = -np.eye(n)
355 A_rows.append(A_wlb)
356 b_rows.append(-lb_w)
357 cones.append(clarabel.NonnegativeConeT(n))
359 A_wub = np.zeros((n, n_vars)) # noqa: N806
360 A_wub[:, 1 : 1 + n] = np.eye(n)
361 A_rows.append(A_wub)
362 b_rows.append(ub_w)
363 cones.append(clarabel.NonnegativeConeT(n))
365 A_ylb = np.zeros((k, n_vars)) # noqa: N806
366 A_ylb[:, 1 + n :] = -np.eye(k)
367 A_rows.append(A_ylb)
368 b_rows.append(-lb_y)
369 cones.append(clarabel.NonnegativeConeT(k))
371 A_yub = np.zeros((k, n_vars)) # noqa: N806
372 A_yub[:, 1 + n :] = np.eye(k)
373 A_rows.append(A_yub)
374 b_rows.append(ub_y)
375 cones.append(clarabel.NonnegativeConeT(k))
377 for a, lb_val, ub_val in extra_constraints:
378 a = np.asarray(a)
379 if lb_val is not None and ub_val is not None and lb_val == ub_val:
380 A_extra = np.zeros((1, n_vars)) # noqa: N806
381 A_extra[0, 1 : 1 + n] = a
382 A_rows.append(A_extra)
383 b_rows.append(np.array([lb_val]))
384 cones.append(clarabel.ZeroConeT(1))
385 else:
386 if lb_val is not None:
387 A_extra = np.zeros((1, n_vars)) # noqa: N806
388 A_extra[0, 1 : 1 + n] = -a
389 A_rows.append(A_extra)
390 b_rows.append(np.array([-lb_val]))
391 cones.append(clarabel.NonnegativeConeT(1))
392 if ub_val is not None:
393 A_extra = np.zeros((1, n_vars)) # noqa: N806
394 A_extra[0, 1 : 1 + n] = a
395 A_rows.append(A_extra)
396 b_rows.append(np.array([ub_val]))
397 cones.append(clarabel.NonnegativeConeT(1))
399 A = sparse.csc_matrix(np.vstack(A_rows)) # noqa: N806
400 b = np.concatenate(b_rows)
402 settings = clarabel.DefaultSettings()
403 settings.verbose = False
404 sol = clarabel.DefaultSolver(P, q, A, b, cones, settings).solve()
405 status = str(sol.status)
407 if "Solved" in status:
408 weights.value = np.array(sol.x[1 : 1 + n])
409 if y_var is not None:
410 y_var.value = np.array(sol.x[1 + n : 1 + n + k])
411 return float(sol.obj_val), float(sol.x[0]), status
412 return None, None, status