Case study: compressing a factor model

The returns matrix

The returns are invented for illustration, not market data. They come from a factor model. Each month a few common moves hit every asset, each in its own proportion, and each asset adds noise of its own.

• A market move, with a monthly volatility of 4.5%, that every asset follows with a loading between 0.7 and 1.3.
• A move of 3% that lifts the first ten assets and pushes the last ten down.
• A third move of 3%, with loadings between $-1$ and 1.
• Each asset's own noise, with a volatility between 2% and 3.5% a month.

The cell draws them with np.random.default_rng(38). As in case study 3.2.5, the returns sit in a matrix $X$ with one row per month and one column per asset. Here $X$ is $120 \times 20$, in monthly decimals.

Centre first. Subtract each column's mean from that column, as case study 3.1.4 did before measuring angles. Call the result $X_c$. Every column now averages 0, and only the moves around the mean are left.

The sample covariance matrix of case study 3.2.5 is then a single product, $C = \frac{1}{n - 1}X_c^\top X_c$, with $n = 120$ months. Entry $(i, j)$ is the dot product of centred columns $i$ and $j$, over $n - 1$. Case study 3.2.5 wrote it $S$; in this module $S$ holds singular values, so it is $C$ here.

The SVD of the centred returns

Take the singular value decomposition $X_c = USV^\top$ (lesson 3.8.1). With full_matrices=False, $U$ is $120 \times 20$, there are 20 singular values, and $V^\top$ is $20 \times 20$.

Put it into $C$. Since $U^\top U = I$, $C = \frac{1}{n - 1}VS^\top U^\top USV^\top = V\,\frac{S^2}{n - 1}\,V^\top$. That is the spectral decomposition of $C$ (lesson 3.6.3), with $V$ in the place of $Q$.

So the SVD finds the principal components without forming $C$ at all. Here the three largest variances are 0.0376, 0.0091 and 0.0046 a month, out of a total of 0.0645.

It also gives each component's history: $X_c\mathbf{v}_i = \sigma_i\mathbf{u}_i$ (lesson 3.8.1) holds the 120 monthly returns of the "portfolio" with holdings $\mathbf{v}_i$.

The share of the total variance along $\mathbf{v}_i$ is $\lambda_i$ over the trace of $C$: the sum of the 20 variances, which equals the sum of the eigenvalues (lesson 3.6.1). The factor $n - 1$ cancels, so it is $\sigma_i^2 / \sum_j \sigma_j^2$: the share of lesson 3.8.2.

How many factors?

The first three components carry 58.33%, 14.10% and 7.06% of the variance, together 79.49%. After them comes a floor: 2.14%, 2.06%, 1.97% and so on, falling slowly to 0.44%.

The gap shows in the ratios of neighbouring singular values (lesson 3.8.2). $\sigma_3/\sigma_4 = 1.815$, while every later ratio lies between 1.004 and 1.103. So the data suggest $k = 3$, the number of factors that made them.

A target share would mislead here. Reaching 90% takes 9 components (90.39%), but components 4 to 9 carry only 1.43% to 2.14% each. They are the assets' own noise, which has no common direction; with noise this large, a gap is a better guide than a target.

What the components look like

• $\mathbf{v}_1$ has all 20 entries positive, from 0.16 to 0.30: the whole market moving together.
• $\mathbf{v}_2$ is positive for the first ten assets (0.12 to 0.30) and negative for the last ten ($-0.31$ to $-0.08$): the first ten against the rest.
• $\mathbf{v}_3$ mixes signs, with no simple pattern: the third move.

The model's three loading vectors are known here, because we made the data. Each lies almost inside the span of $\mathbf{v}_1$, $\mathbf{v}_2$ and $\mathbf{v}_3$ (lesson 3.5.1): projecting it there (lesson 3.5.2) keeps 99.86%, 99.32% and 99.08% of its length.

But the $\mathbf{v}_i$ are not the model's factors one by one. $\mathbf{v}_1$ lies 8.0° from the market's loading vector and $\mathbf{v}_2$ 10.1° from the split, but $\mathbf{v}_3$ is 44.5° from the third. The SVD forces its directions to be perpendicular, and the model's loading vectors are not: it finds the space of the factors, not their names.

A rank-3 copy of the returns

Keep three pieces (lesson 3.8.2): $X_3 = U_3S_3V_3^\top$. Each month's 20 centred returns become a mix of the three patterns $\mathbf{v}_1$, $\mathbf{v}_2$ and $\mathbf{v}_3$, with three amounts per month.

Storage. $X_c$ holds 2,400 numbers. $X_3$ needs $3 \times (120 + 20 + 1) = 423$, or 17.6%: three histories of 120 months, three patterns of 20 assets, and three singular values.

The residual. What is left, $E = X_c - X_3$, carries 20.51% of the squared size: 100% minus the 79.49% kept. It should look like the assets' own noise, and it does:

The common moves are gone: the residual correlations sit around 0, as in the figure. They lean slightly negative because $E$ is perpendicular to the three directions taken out, which pushes what is left a little apart.

Each asset's residual volatility tracks the noise the model gave it, from 2% to 3.5% a month: the two lists have a correlation of 0.89. The copy has split each asset's risk into a common part and its own part.

A compressed covariance matrix

The covariance matrix of the rank-3 copy keeps three terms of the spectral decomposition: $C_3 = V_3\Lambda_3V_3^\top$, with $\Lambda_3$ the diagonal of $\lambda_1$, $\lambda_2$ and $\lambda_3$. For a spread-out portfolio it is almost the whole story. With $\frac{1}{20}$ in each asset, $C$ gives a variance of 0.00182 a month, a volatility of 14.78% a year, and $C_3$ gives 0.001819.

The trap: a rank-3 matrix

$C_3$ has rank 3, so it is singular (lesson 3.4.1): there are portfolios it calls riskless. Take a pair trade, long asset 1 and short asset 2, and adjust it slightly so that it has no exposure to $\mathbf{v}_1$, $\mathbf{v}_2$ or $\mathbf{v}_3$. Its holdings are $(0.988, -0.983, \ldots)$, with the other 18 between $-0.07$ and 0.05.

$C_3$ gives this portfolio a variance of 0, up to rounding. $C$ gives it 0.001212, a volatility of 3.48% a month. An optimiser handed $C_3$ would load up on such trades, believing them free of risk.

The fix: put back each asset's own risk

Add back, on the diagonal, the variance that $C_3$ misses for each asset: $D = \operatorname{diag}(C) - \operatorname{diag}(C_3)$. These are the residual variances.

$C_f$ gives the pair trade 3.39% a month, close to the 3.48% from $C$, and the spread-out portfolio 14.91% a year. Its smallest eigenvalue is 0.000416, above 0, so no portfolio is called riskless.

It needs $20 \times 3 + 3 + 20 = 83$ numbers, against 210 for $C$. With 500 assets the gap is wider: 2,003 against 125,250.

Closer to the truth. We made the data, so the model's own covariance matrix is known. The Frobenius distance (lesson 3.8.2) from $C$ to it is 0.0074; from $C_f$ it is 0.00685. Replacing the 17 small components by a diagonal threw away mostly noise, so the compressed matrix is the better estimate.

On the desk

Case study 3.6.5 found the principal components of yield changes from their covariance matrix with eigh. Run the steps of this case study on those changes, centred, and the right singular vectors are the level, slope and curvature it found, with no covariance matrix formed. Working on $X_c$ itself is also the more accurate route, for the reason lesson 3.5.3 gave for lstsq.

import numpy as np


def _returns():
    """The invented returns of case study 3.8.4: 120 months of 20 assets, seed 38."""
    rng = np.random.default_rng(38)
    n, m = 120, 20
    loadings = np.column_stack([
        rng.uniform(0.7, 1.3, m),
        np.r_[rng.uniform(0.4, 1.0, 10), -rng.uniform(0.4, 1.0, 10)],
        rng.uniform(-1.0, 1.0, m),
    ])
    factors = rng.standard_normal((n, 3)) * [0.045, 0.03, 0.03]
    specific = rng.uniform(0.02, 0.035, m)
    means = rng.uniform(0.004, 0.012, m)
    return means + factors @ loadings.T + rng.standard_normal((n, m)) * specific


X = _returns()
C = np.cov(X, rowvar=False)


def _reference(k):
    lam, Q = np.linalg.eigh(C)
    lam, Q = lam[::-1][:k], Q[:, ::-1][:, :k]
    common = Q @ np.diag(lam) @ Q.T
    return common + np.diag(np.diag(C) - np.diag(common))


def _raises(f, what):
    try:
        f()
    except ValueError:
        return
    raise AssertionError(f"expected ValueError for {what}")


def test_share_of_three():
    """Three components carry 79.49% of the variance of the case study's returns"""
    expect_close(variance_share(X, 3), 0.794933, what="variance_share(X, 3)", rel=1e-5)


def test_share_of_all():
    """All 20 components carry all of it"""
    expect_close(variance_share(X, 20), 1.0, what="variance_share(X, 20)", rel=1e-9)


def test_matches_eigh():
    """factor_covariance(X, 3) matches the same matrix built from eigh of np.cov"""
    expect_close(factor_covariance(X, 3), _reference(3), what="factor_covariance(X, 3)", rel=1e-7, abs_tol=1e-12)


def test_diagonal_of_cov():
    """Its diagonal is the diagonal of np.cov: each asset keeps its own variance"""
    expect_close(np.diag(factor_covariance(X, 3)), np.diag(C), what="diagonal of factor_covariance(X, 3)", rel=1e-9)


def test_positive_definite():
    """With 3 components it is symmetric and positive definite (smallest eigenvalue above 0)"""
    Cf = np.asarray(factor_covariance(X, 3))
    expect_equal(bool(np.allclose(Cf, Cf.T)), True, what="factor_covariance(X, 3) is symmetric")
    expect_equal(bool(np.linalg.eigvalsh(Cf).min() > 1e-5), True, what="smallest eigenvalue above 1e-5")


def test_all_components_give_cov():
    """With all 20 components it is np.cov itself"""
    expect_close(factor_covariance(X, 20), C, what="factor_covariance(X, 20)", rel=1e-7, abs_tol=1e-12)


def test_gross_returns():
    """Gross returns 1 + r give the same share and the same matrix: the returns are centred"""
    expect_close(variance_share(1 + X, 3), variance_share(X, 3), what="variance_share(1 + X, 3)", rel=1e-9)
    expect_close(factor_covariance(1 + X, 3), _reference(3), what="factor_covariance(1 + X, 3)", rel=1e-6, abs_tol=1e-12)


def test_bad_k():
    """k = 0, k = 21 for 20 assets and k = 2.5 raise ValueError"""
    for k in (0, 21, 2.5):
        _raises(lambda: factor_covariance(X, k), f"k = {k}")
        _raises(lambda: variance_share(X, k), f"variance_share with k = {k}")


def test_bad_returns():
    """A 1-D array, a single row and a missing value raise ValueError"""
    _raises(lambda: factor_covariance(X[:, 0], 1), "a 1-D array")
    _raises(lambda: factor_covariance(X[:1], 1), "a single row")
    gap = X.copy()
    gap[5, 3] = np.nan
    _raises(lambda: factor_covariance(gap, 3), "a nan in the returns")
    _raises(lambda: variance_share(gap, 3), "variance_share with a nan")