Coverage for src/basanos/math/_factor_model.py: 100%
75 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-08-04 07:53 +0000
1r"""Factor risk model decomposition (Section 4.1 of basanos.pdf).
3This private module provides the `FactorModel` frozen dataclass, which
4encapsulates the three-component factor model
6$$
7\\bm{\\Sigma} = \\mathbf{B}\\mathbf{F}\\mathbf{B}^\\top + \\mathbf{D}
8$$
10and a class method for fitting the model from a return matrix via the
11Singular Value Decomposition (Section 4.2).
12"""
14from __future__ import annotations
16import dataclasses
17from typing import cast
19import numpy as np
21# cvx-linalg exposes this constant as a public top-level name since >= 0.7;
22# the project floor is cvx-linalg>=0.9.0 (see pyproject.toml), so the public
23# import always resolves and no legacy-layout fallback is required.
24from cvx.linalg import DEFAULT_COND_THRESHOLD as _DEFAULT_COND_THRESHOLD
25from cvx.linalg import DimensionMismatchError, SingularMatrixError
26from cvx.linalg import check_and_warn_condition as _check_and_warn_condition
27from cvx.linalg import inv as _inv
28from cvx.linalg import solve as _solve
30from basanos.exceptions import FactorModelError
33@dataclasses.dataclass(frozen=True)
34class FactorModel:
35 r"""Frozen dataclass for a factor risk model decomposition (Section 4.1).
37 Encapsulates the three components of the factor model
39 $$
40 \bm{\Sigma} = \mathbf{B}\mathbf{F}\mathbf{B}^\top + \mathbf{D}
41 $$
43 where
45 - $\mathbf{B} \in \mathbb{R}^{n \times k}$ is the *factor loading
46 matrix*: column $j$ gives the sensitivity of each asset to
47 factor $j$.
48 - $\mathbf{F} \in \mathbb{R}^{k \times k}$ is the *factor covariance
49 matrix* (positive definite), capturing how the $k$ factors
50 co-vary.
51 - $\mathbf{D} = \operatorname{diag}(d_1, \dots, d_n)$ with
52 $d_i > 0$ is the *idiosyncratic variance* diagonal, capturing
53 the asset-specific variance unexplained by the common factors.
55 The central assumption is $k \ll n$: the dominant systematic sources
56 of risk are captured by a handful of factors while the idiosyncratic
57 component is, by construction, uncorrelated across assets.
59 Attributes:
60 factor_loadings: Factor loading matrix $\mathbf{B}$,
61 shape ``(n, k)``.
62 factor_covariance: Factor covariance matrix $\mathbf{F}$,
63 shape ``(k, k)``.
64 idiosyncratic_var: Idiosyncratic variance vector
65 $(d_1, \dots, d_n)$, shape ``(n,)``. All entries must be
66 strictly positive.
68 Examples:
69 >>> import numpy as np
70 >>> loadings = np.eye(3, 2)
71 >>> cov = np.eye(2) * 0.5
72 >>> idio = np.array([0.5, 0.5, 1.0])
73 >>> fm = FactorModel(factor_loadings=loadings, factor_covariance=cov, idiosyncratic_var=idio)
74 >>> fm.n_assets
75 3
76 >>> fm.n_factors
77 2
78 >>> fm.covariance.shape
79 (3, 3)
80 """
82 factor_loadings: np.ndarray
83 factor_covariance: np.ndarray
84 idiosyncratic_var: np.ndarray
86 def __post_init__(self) -> None:
87 """Validate shape consistency and strict positivity after initialization.
89 Raises:
90 FactorModelError: If ``factor_loadings`` is not 2-D.
91 FactorModelError: If ``factor_covariance`` shape does not
92 match the number of factors inferred from ``factor_loadings``.
93 FactorModelError: If ``idiosyncratic_var`` length does
94 not match the number of assets inferred from ``factor_loadings``.
95 FactorModelError: If any element of
96 ``idiosyncratic_var`` is not strictly positive.
97 """
98 if self.factor_loadings.ndim != 2:
99 raise FactorModelError(f"factor_loadings must be 2-D, got ndim={self.factor_loadings.ndim}.") # noqa: TRY003
100 n, k = self.factor_loadings.shape
101 if self.factor_covariance.shape != (k, k):
102 raise FactorModelError( # noqa: TRY003
103 f"factor_covariance must have shape ({k}, {k}) to match "
104 f"factor_loadings columns, got {self.factor_covariance.shape}."
105 )
106 if self.idiosyncratic_var.shape != (n,):
107 raise FactorModelError( # noqa: TRY003
108 f"idiosyncratic_var must have shape ({n},) to match factor_loadings rows, "
109 f"got {self.idiosyncratic_var.shape}."
110 )
111 if not np.all(self.idiosyncratic_var > 0):
112 raise FactorModelError("All entries of idiosyncratic_var must be strictly positive.") # noqa: TRY003
114 @property
115 def n_assets(self) -> int:
116 """Number of assets *n* (rows of ``factor_loadings``)."""
117 return int(self.factor_loadings.shape[0])
119 @property
120 def n_factors(self) -> int:
121 """Number of factors *k* (columns of ``factor_loadings``)."""
122 return int(self.factor_loadings.shape[1])
124 @property
125 def covariance(self) -> np.ndarray:
126 r"""Reconstruct the full $n \times n$ covariance matrix.
128 Computes $\bm{\Sigma} = \mathbf{B}\mathbf{F}\mathbf{B}^\top +
129 \mathbf{D}$ by combining the low-rank systematic component with the
130 diagonal idiosyncratic component.
132 Returns:
133 np.ndarray: Shape ``(n, n)`` symmetric covariance matrix.
135 Examples:
136 >>> import numpy as np
137 >>> loadings = np.array([[1.0, 0.0], [0.0, 1.0], [0.0, 0.0]])
138 >>> cov = np.eye(2)
139 >>> idio = np.ones(3)
140 >>> fm = FactorModel(factor_loadings=loadings, factor_covariance=cov, idiosyncratic_var=idio)
141 >>> fm.covariance.diagonal().tolist()
142 [2.0, 2.0, 1.0]
143 """
144 return cast(
145 "np.ndarray",
146 self.factor_loadings @ self.factor_covariance @ self.factor_loadings.T + np.diag(self.idiosyncratic_var),
147 )
149 @property
150 def woodbury_condition_number(self) -> float:
151 r"""Condition number of the inner $k \times k$ Woodbury matrix.
153 Returns the condition number of the matrix
155 $$
156 \mathbf{M} = \mathbf{F}^{-1} + \mathbf{B}^\top\mathbf{D}^{-1}\mathbf{B}
157 $$
159 which is the matrix actually inverted during `solve`. A large
160 value (above ``_DEFAULT_COND_THRESHOLD`` ≈ 1e12) indicates that the
161 Woodbury solve is numerically unreliable.
163 This property gives callers a way to inspect the numerical health of
164 the model without performing a full solve. Unlike the condition number
165 of the full $n \times n$ covariance matrix, this measure is
166 specific to the $k \times k$ inner system solved inside the
167 Woodbury identity.
169 Returns:
170 float: Condition number $\kappa(\mathbf{M})$. Returns
171 ``inf`` when $\mathbf{F}$ is not positive-definite (e.g.
172 singular or indefinite), as the Cholesky decomposition used to
173 form $\mathbf{F}^{-1}$ fails in that case.
175 Examples:
176 >>> import numpy as np
177 >>> loadings = np.eye(3, 1)
178 >>> cov = np.eye(1)
179 >>> idio = np.ones(3)
180 >>> fm = FactorModel(factor_loadings=loadings, factor_covariance=cov, idiosyncratic_var=idio)
181 >>> fm.woodbury_condition_number > 0
182 True
183 """
184 d_inv = 1.0 / self.idiosyncratic_var # (n,)
185 d_inv_b_mat = d_inv[:, None] * self.factor_loadings # D^{-1} B, shape (n, k)
186 try:
187 f_inv = _inv(self.factor_covariance)
188 mid = f_inv + self.factor_loadings.T @ d_inv_b_mat # (k, k)
189 except (np.linalg.LinAlgError, SingularMatrixError):
190 return float("inf")
191 return float(np.linalg.cond(mid))
193 def solve(
194 self,
195 rhs: np.ndarray,
196 cond_threshold: float = _DEFAULT_COND_THRESHOLD,
197 ) -> np.ndarray:
198 r"""Solve $\bm{\Sigma}\,\mathbf{x} = \mathbf{b}$ via the Woodbury identity.
200 Applies the Sherman--Morrison--Woodbury formula (Section 4.3 of
201 basanos.pdf) to avoid forming or factorising the full
202 $n \times n$ covariance matrix:
204 $$
205 (\mathbf{D} + \mathbf{B}\mathbf{F}\mathbf{B}^\top)^{-1}
206 = \mathbf{D}^{-1}
207 - \mathbf{D}^{-1}\mathbf{B}
208 \bigl(\mathbf{F}^{-1} + \mathbf{B}^\top\mathbf{D}^{-1}\mathbf{B}\bigr)^{-1}
209 \mathbf{B}^\top\mathbf{D}^{-1}.
210 $$
212 Because $\mathbf{D}$ is diagonal, $\mathbf{D}^{-1}$ is
213 free. The inner matrix is $k \times k$ with cost
214 $O(k^3)$, and the surrounding multiplications cost
215 $O(kn)$. Total cost is $O(k^3 + kn)$ rather than
216 $O(n^3)$.
218 Args:
219 rhs: Right-hand side vector $\mathbf{b}$, shape ``(n,)``.
220 cond_threshold: Condition-number threshold above which an
221 `IllConditionedMatrixWarning` is
222 emitted. The check is applied to both ``factor_covariance``
223 ($\mathbf{F}$) and to the inner $k \times k$
224 Woodbury matrix $\mathbf{F}^{-1} + \mathbf{B}^\top
225 \mathbf{D}^{-1}\mathbf{B}$. Defaults to ``1e12``.
227 Returns:
228 np.ndarray: Solution vector $\mathbf{x}$, shape ``(n,)``.
230 Raises:
231 DimensionMismatchError: If ``rhs`` length does not match
232 ``n_assets``.
233 SingularMatrixError: If the inner $k \times k$ matrix is
234 singular.
236 Examples:
237 >>> import numpy as np
238 >>> loadings = np.eye(3, 1)
239 >>> cov = np.eye(1)
240 >>> idio = np.ones(3)
241 >>> fm = FactorModel(factor_loadings=loadings, factor_covariance=cov, idiosyncratic_var=idio)
242 >>> rhs = np.array([1.0, 2.0, 3.0])
243 >>> x = fm.solve(rhs)
244 >>> np.allclose(fm.covariance @ x, rhs)
245 True
246 """
247 n = self.n_assets
248 if rhs.shape != (n,):
249 raise DimensionMismatchError(rhs.size, n)
251 # D^{-1} is free because D is diagonal
252 d_inv = 1.0 / self.idiosyncratic_var # (n,)
253 d_inv_rhs = d_inv * rhs # D^{-1} b, shape (n,)
254 d_inv_b_mat = d_inv[:, None] * self.factor_loadings # D^{-1} B, shape (n, k)
256 # Solve mid * w = B^T D^{-1} b, where mid = F^{-1} + B^T D^{-1} B.
257 # F^{-1} is obtained via a Cholesky solve rather than an explicit
258 # inversion, consistent with the Cholesky-first discipline in _linalg.py.
259 # A condition-number check on factor_covariance is applied first so
260 # that ill-conditioned F is flagged before its inverse enters mid.
261 rhs_k = self.factor_loadings.T @ d_inv_rhs # (k,)
262 try:
263 _check_and_warn_condition(self.factor_covariance, cond_threshold)
264 mid = _inv(self.factor_covariance, cond_threshold) + self.factor_loadings.T @ d_inv_b_mat # (k, k)
265 _check_and_warn_condition(mid, cond_threshold)
266 w = _solve(mid, rhs_k, cond_threshold) # (k,)
267 except np.linalg.LinAlgError as exc:
268 raise SingularMatrixError(str(exc)) from exc
270 # x = D^{-1} b - D^{-1} B w
271 return cast("np.ndarray", d_inv_rhs - d_inv_b_mat @ w)
273 @classmethod
274 def from_returns(cls, returns: np.ndarray, k: int) -> FactorModel:
275 r"""Fit a rank-*k* factor model from a return matrix via truncated SVD.
277 Extracts latent factors from the return matrix
278 $\mathbf{R} \in \mathbb{R}^{T \times n}$ using the Singular
279 Value Decomposition (SVD). The top-*k* singular triplets define the
280 factor model components:
282 $$
283 \mathbf{B} = \mathbf{V}_k, \quad
284 \mathbf{F} = \bm{\Sigma}_k^2 / T, \quad
285 \hat{d}_i = 1 - \bigl(\mathbf{B}\mathbf{F}\mathbf{B}^\top\bigr)_{ii}
286 $$
288 where $\mathbf{V}_k$ and $\bm{\Sigma}_k$ are the top-*k*
289 right singular vectors and singular values of $\mathbf{R}$
290 respectively. When *returns* contains unit-variance columns (as
291 produced by `vol_adj`), the sample
292 covariance has unit diagonal; the idiosyncratic term
293 $\hat{d}_i = 1 - (\mathbf{B}\mathbf{F}\mathbf{B}^\top)_{ii}$
294 absorbs the residual so the full covariance $\hat{\mathbf{C}}^{(k)}$
295 also has unit diagonal. Each $\hat{d}_i$ is clamped from below
296 at machine epsilon to guarantee strict positivity.
298 Args:
299 returns: Return matrix of shape ``(T, n)``, typically
300 volatility-adjusted log returns with rows as timestamps and
301 columns as assets.
302 k: Number of factors to retain. Must satisfy
303 ``1 <= k <= min(T, n)``.
305 Returns:
306 FactorModel: Fitted factor model with ``n_assets = n`` and
307 ``n_factors = k``.
309 Raises:
310 FactorModelError: If *returns* is not 2-D.
311 FactorModelError: If *k* is outside the range ``[1, min(T, n)]``.
313 Examples:
314 >>> import numpy as np
315 >>> rng = np.random.default_rng(0)
316 >>> ret = rng.standard_normal((50, 5))
317 >>> fm = FactorModel.from_returns(ret, k=2)
318 >>> fm.n_factors
319 2
320 >>> fm.n_assets
321 5
322 >>> fm.covariance.shape
323 (5, 5)
324 """
325 if returns.ndim != 2:
326 raise FactorModelError(f"Return matrix must be 2-D, got ndim={returns.ndim}.") # noqa: TRY003
327 t_len, n = returns.shape
328 if not (1 <= k <= min(t_len, n)):
329 raise FactorModelError(f"k must satisfy 1 <= k <= min(T, n) = {min(t_len, n)}, got k={k}.") # noqa: TRY003
331 _, s, vt = np.linalg.svd(returns, full_matrices=False)
333 # Top-k right singular vectors as columns: shape (n, k)
334 v_k = vt[:k].T
335 s_k = s[:k]
337 # Factor covariance: diagonal matrix with entries s_j**2 / T
338 factor_cov = np.diag(s_k**2 / t_len)
340 # Diagonal of B*F*B^T = sum_j (s_j**2/T) * B[:,j]**2
341 factor_diag = (v_k**2) @ (s_k**2 / t_len)
343 # Idiosyncratic variance: target diagonal is 1.0 (unit-variance columns
344 # assumed); residual = 1.0 - systematic contribution, clamped to (0, inf)
345 _unit_variance = 1.0
346 d = np.maximum(_unit_variance - factor_diag, np.finfo(float).eps)
348 return cls(factor_loadings=v_k, factor_covariance=factor_cov, idiosyncratic_var=d)