Covariance Matrix Estimation for Portfolio Construction

The portfolio optimizers in Portfolio Management Fundamentals treated expected returns and covariance as supplied inputs. Style and Factor Exposures estimated exposures from historical returns. We now examine the risk estimates themselves: how much structure should we impose on a noisy covariance matrix, and how should we evaluate the resulting allocations?

We move from hand-checkable covariance calculations to structured estimators, time-varying volatility, and a common-sample portfolio comparison. The objective is to understand what each estimate assumes and how estimation errors affect an investment decision. There is no universally best estimator.

Run the examples in order. Historical observations are monthly decimal returns. Covariance is reported in squared return units unless explicitly rescaled; multiplying monthly covariance by 12 is a reporting convention under zero serial cross-covariances.

import json
import numpy as np
import pandas as pd
import PortfolioOptimizationKit as pok

From Returns to Portfolio Risk

The inputs below describe returns over a specified observation interval. Their covariance supports linear portfolio-return calculations with beginning-of-period weights. For an option portfolio, a fixed delta exposure is only a local approximation: Deriving and Using Black–Scholes shows that delta changes with price and time. Large scenario moves or nonlinear tail-risk calculations can require repricing each contract before measuring portfolio losses. Historical simple-return volatility is also not exactly a GBM diffusion parameter at a finite observation interval.

For an $N$-asset return vector $\mathbf{R}$ over one common horizon, define:

$$\mathit{\mu }=\mathbb{E}[\mathbf{R}],\mathrm{\Sigma }=\mathbb{E}[(\mathbf{R}-\mathit{\mu })(\mathbf{R}-\mathit{\mu }{)}^T].$$

Its diagonal entries are variances; off-diagonal entries are covariances. For assets with positive volatility, ${\mathrm{\Sigma }}_{ij}={\rho }_{ij}{\sigma }_i{\sigma }_j$. Fixed beginning-period portfolio weights give:

$${\sigma }_p^2={\mathbf{w}}^T\mathrm{\Sigma }\mathbf{w},{\sigma }_p=\sqrt{{\mathbf{w}}^T\mathrm{\Sigma }\mathbf{w}}.$$

A covariance matrix must be symmetric and positive semidefinite (PSD): no portfolio can have negative variance. Positive definiteness is stronger: it rules out nonzero zero-variance combinations and permits an ordinary matrix inverse. A singular PSD matrix can still define a valid constrained minimum-variance problem.

Sample Covariance: A Worked Calculation

With $T>1$ common observations, estimate the mean of each asset and calculate:

$$S=\frac{1}{T-1}\sum _{t=1}^T({\mathbf{r}}_t-\overline{\mathbf{r}})({\mathbf{r}}_t-\overline{\mathbf{r}}{)}^T=\frac{X^{T}X}{T-1},$$

where row $t$ of $X$ contains the demeaned returns. The $T-1$ denominator accounts for estimating the mean. Under iid sampling with finite second moments, this is unbiased for the population covariance. Historical observations still need to be representative of the risk we want to estimate: sample covariance is not free of modeling assumptions.

Take three monthly returns for Asset A, $(-2\mathrm{\%},0\mathrm{\%},2\mathrm{\%})$, and Asset B, $(1\mathrm{\%},-1\mathrm{\%},0\mathrm{\%})$. Both means are zero. Then:

$$s_A^2=\frac{{0.02}^2+0+{0.02}^2}{2}=0.0004,s_B^2=\frac{{0.01}^2+{0.01}^2+0}{2}=0.0001,$$$$s_{AB}=\frac{(-0.02)(0.01)+0(-0.01)+(0.02)0}{2}=-0.0001.$$

Correlation is $-0.0001/(0.02\times 0.01)=-0.5$. A 60/40 portfolio has monthly variance ${0.6}^2(0.0004)+{0.4}^2(0.0001)+2(0.6)(0.4)(-0.0001)=0.000112$, hence monthly volatility about 1.058%.

worked_returns = pd.DataFrame({"A": [-0.02, 0.0, 0.02], "B": [0.01, -0.01, 0.0]})
centered = worked_returns - worked_returns.mean()
worked_covariance = centered.T @ centered / (len(worked_returns) - 1)
worked_weights = np.array([0.6, 0.4])
worked_variance = worked_weights @ worked_covariance @ worked_weights
print("Monthly sample covariance:")
print(pok.sample_cov(worked_returns).to_string(float_format=lambda x: f"{x:.6f}"))
print(f"60/40 monthly variance: {worked_variance:.6f}")
print(f"Monthly volatility: {np.sqrt(worked_variance):.3%}")
print(f"Annualized volatility: {np.sqrt(12 * worked_variance):.3%}")

The toolkit requires a complete common sample, with unique observation and asset labels. Computing each covariance from a different subset of available dates can produce an indefinite matrix. Decide how to handle missing observations before estimation, rather than letting pairwise deletion silently change the sample for each asset pair.

Dimensionality and Unstable Allocations

A symmetric $N\times N$ covariance matrix has $N(N+1)/2$ distinct entries: $N$ variances and $N(N-1)/2$ covariances. Expected returns add another $N$ parameters if the allocation model needs them. For 500 assets, covariance alone has 125,250 entries.

The number of entries is not a literal minimum number of dates needed for estimation. A sharper algebraic limitation follows from the demeaned data matrix:

$$\mathrm{rank}(S)\le min(N,T-1).$$

If $N\ge T$, sample covariance is necessarily singular. Even when invertible, very small estimated eigenvalues can make unconstrained allocations highly sensitive to estimation error. For example, the unconstrained fully invested GMV solution is proportional to $S^{-1}\mathbf{1}$; inversion amplifies directions with small eigenvalues. Long-only constraints limit some extreme allocations but do not remove uncertainty in the risk estimate.

rng = np.random.default_rng(42)
short_history = pd.DataFrame(rng.normal(0, 0.02, size=(6, 12)))
short_covariance = pok.sample_cov(short_history)
short_shrinkage = pok.shrinkage_cov(short_history, delta=0.5)
print(f"Assets: {short_history.shape[1]}; observed periods: {len(short_history)}")
print(f"Distinct covariance entries: {12 * 13 // 2}")
print(f"Sample covariance rank: {np.linalg.matrix_rank(short_covariance)} (at most 5)")
print(f"Smallest sample eigenvalue: {np.linalg.eigvalsh(short_covariance).min():.3e}")
print(f"Smallest 50% shrinkage eigenvalue: {np.linalg.eigvalsh(short_shrinkage).min():.3e}")

Tiny negative sample eigenvalues at floating-point precision can be roundoff, rather than economically negative variance. Here shrinkage restores positive definiteness because its target is positive definite and receives positive weight. This result is not automatic for every target, especially with zero-variance assets or perfect dependence.

Increasing the historical window adds observations, but also includes older regimes. A shorter window may adapt more quickly while increasing sampling noise. That trade-off motivates structured covariance estimates.

Constant Correlation and Shrinkage

A Constant-Correlation Target

Retain each sample volatility $s_i$ and replace the distinct off-diagonal correlations with their average:

$$\overline{\rho }=\frac{2}{N(N-1)}\sum _{i<j}{\hat{\rho }}_{ij},F_{ij}=\{\begin{array}{ll}{s}_i^2,& i=j,\\ \overline{\rho }{s}_i{s}_j,& i\ne j.\end{array}$$

The target reduces the covariance description to $N$ variances and one common correlation. It does not assume identical volatilities. Its correlation matrix has eigenvalues $1-\overline{\rho }$ with multiplicity $N-1$ and $1+(N-1)\overline{\rho }$ once, so it is PSD when $-1/(N-1)\le \overline{\rho }\le 1$.

With complete common data, averaging the valid sample correlations respects that range. Constant-return assets have undefined correlations; the toolkit gives them zero covariance rows and columns and averages correlations only among variable assets. For a single asset, the target reduces to its variance. For two variable assets, there is only one pair correlation, so this target equals sample covariance: shrinking toward it cannot change the estimate.

Combining Sample Detail with Structure

Define a shrinkage estimate:

$${\hat{\mathrm{\Sigma }}}_{\delta }=(1-\delta )S+\delta F,0\le \delta \le 1.$$

At $\delta =0$, retain the sample matrix; at $\delta =1$, use the target. A convex combination of PSD matrices remains PSD. The target can reduce sampling variability but introduce structural bias if true correlations differ substantially. Shrinkage exchanges some sample detail for stability; it does not establish that the target is the true covariance.

The helper shrinkage_cov() uses a user-specified intensity, defaulting to 0.5. It does not estimate an optimal Ledoit–Wolf intensity. Select any tuned intensity using information available before the evaluation period, or a properly nested validation procedure.

Hand check. Consider three assets with annual volatilities 20%, 30%, and 25%, and pair correlations 0.75, 0.25, and 0.50. Their average is 0.50. For the first pair, sample covariance is $0.75(0.2)(0.3)=0.045$, target covariance is $0.50(0.2)(0.3)=0.030$, and 50% shrinkage gives 0.0375. The third asset makes the target correlation meaningfully different from the first pair's sample correlation.

annual_volatilities = np.array([0.20, 0.30, 0.25])
correlation_matrix = np.array([[1.0, 0.75, 0.25], [0.75, 1.0, 0.50], [0.25, 0.50, 1.0]])
mean_correlation = correlation_matrix[np.triu_indices(3, k=1)].mean()
sample_matrix = correlation_matrix * np.outer(annual_volatilities, annual_volatilities)
target_matrix = np.full((3, 3), mean_correlation) * np.outer(annual_volatilities, annual_volatilities)
np.fill_diagonal(target_matrix, annual_volatilities ** 2)
delta = 0.5
shrunk_matrix = (1 - delta) * sample_matrix + delta * target_matrix
print(f"Average pair correlation: {mean_correlation:.2f}")
print(f"First-pair covariance: sample {sample_matrix[0, 1]:.4f}; target {target_matrix[0, 1]:.4f}; shrinkage {shrunk_matrix[0, 1]:.4f}")
print("Retained annual variances:", np.diag(shrunk_matrix))

Try it: set delta to 0 and 1. The off-diagonal covariance moves between the two endpoints, while each diagonal variance remains unchanged. The historical examples below use the helper to estimate both matrices from observed monthly returns.

Factor-Based Covariance

The factor model from the previous chapter also supplies a covariance structure. For demeaned returns, write:

$$\mathbf{R}-\mathit{\mu }=B\mathbf{f}+\mathit{\epsilon }.$$

Let $B$ be an $N\times K$ loading matrix, $\mathrm{\Omega }=\mathrm{Cov}(\mathbf{f})$, and $D=\mathrm{Cov}(\mathit{\epsilon })$. If factors and residuals are uncorrelated:

$$\mathrm{\Sigma }=B\mathrm{\Omega }{B}^T+D.$$

The two terms represent common-factor and residual risk. Assuming $D$ is diagonal reduces the number of residual covariance parameters, but is a substantive restriction: omitted industry or firm connections may leave correlated residuals.

For a specified factor basis, the covariance construction uses $NK$ loadings, $K(K+1)/2$ factor covariance entries, and $N$ residual variances. With two potentially correlated factors, that is $3N+3$, excluding means and intercepts. If factor covariances are constrained to zero, it is $3N+2$. These counts depend on what is estimated, fixed, or normalized; a latent-factor model also has identification restrictions.

For a concrete one-year model, use factor variances 0.04 and 0.01, uncorrelated factors, and the loadings below. Asset A's variance is $1^2(0.04)+{0.2}^2(0.01)+0.01=0.0504$. The A–B covariance is $(1)(0.8)(0.04)+(0.2)(-0.3)(0.01)=0.0314$.

loadings = np.array([[1.0, 0.2], [0.8, -0.3], [1.2, 0.5]])
factor_covariance = np.diag([0.04, 0.01])
residual_covariance = np.diag([0.01, 0.02, 0.015])
factor_model_covariance = loadings @ factor_covariance @ loadings.T + residual_covariance
print(pd.DataFrame(factor_model_covariance, index=["A", "B", "C"], columns=["A", "B", "C"]).to_string(float_format=lambda x: f"{x:.5f}"))
print(f"Asset A annual volatility: {np.sqrt(factor_model_covariance[0, 0]):.2%}")

Positive semidefinite $\mathrm{\Omega }$ and $D$ make the constructed matrix PSD. The calculation shows the structure, rather than estimating a factor model from this three-asset illustration. Estimation still requires common dates, identified exposures, and an appropriate residual model.

Time-Varying Volatility

A Historical Market Proxy

Use the same bundled 30-industry proxy as Part 1. It combines each month's value-weighted industry returns using the preceding month's implied capitalization fractions. This is an educational market proxy, not an official market index. The data do not establish when all capitalization metadata became public.

market_returns = pok.get_total_market_index_returns()
periods_per_year = 12
full_sample_volatility = pok.annualize_vol(market_returns, periods_per_year)
print(f"Proxy: {market_returns.index[0]} to {market_returns.index[-1]}, {len(market_returns)} monthly returns")
print(f"Full-sample annualized sample volatility: {full_sample_volatility:.2%}")
print("\n<ECHARTS_DATA>" + json.dumps({
    "dates": market_returns.index.astype(str).tolist(),
    "marketProxyReturns": {
        "title": "Lagged Industry Market Proxy: Monthly Returns",
        "yAxisName": "Monthly return (%)",
        "series": {"Market proxy": (market_returns * 100).tolist()},
    },
}, allow_nan=False))

A full-sample standard deviation summarizes one entire history. It is not the volatility estimate that would have been available at every earlier date.

Rolling Windows and Annualization

For a trailing window of $L$ monthly observations, compute its sample standard deviation and multiply by $\sqrt{12}$, regardless of whether $L$ is 12, 36, or 60:

$${\hat{\sigma }}_{t,\mathrm{ann}}^{(L)}=\sqrt{12}s(r_{t-L+1},\dots ,r_t).$$

The window length determines how much history is estimated; the observation frequency determines annualization. Using $\sqrt{60}$ for 60 monthly observations incorrectly changes the reporting horizon. As before, square-root scaling is a convention for sums under zero serial covariance, not an exact formula for compounded annual return volatility.

rolling_volatility = pd.DataFrame({
    "12-month window": market_returns.rolling(12).std(ddof=1) * np.sqrt(periods_per_year),
    "60-month window": market_returns.rolling(60).std(ddof=1) * np.sqrt(periods_per_year),
})
print("Mean of available rolling estimates (%):")
print((rolling_volatility.mean() * 100).to_string(float_format=lambda x: f"{x:.2f}"))
print("\n<ECHARTS_DATA>" + json.dumps({
    "dates": rolling_volatility.index.astype(str).tolist(),
    "rollingVolatility": {
        "title": "Trailing Volatility: Same Annual Units, Different Windows",
        "yAxisName": "Annualized sample volatility (%)",
        "series": {name: [None if pd.isna(value) else float(value * 100) for value in rolling_volatility[name]]
                   for name in rolling_volatility},
    },
}, allow_nan=False))

Each value uses returns through the labeled month. To use it before the next month's return, shift the estimate forward by one observation. The longer window generally changes more slowly but can retain information from a regime that has ended. The mean of rolling volatilities need not equal full-sample volatility: the windows overlap, cover different dates, and take square roots before averaging.

Volatility Cones

A volatility cone summarizes the distribution of rolling estimates for several window lengths. Here the minimum, maximum, median, and central historical quantiles describe overlapping estimates from this one history. They are not confidence limits for next year's volatility.

window_years = [1, 2, 3, 5, 10, 20]
cone_rows = []
for years in window_years:
    estimates = market_returns.rolling(12 * years).std(ddof=1).dropna() * np.sqrt(12)
    cone_rows.append([estimates.min(), estimates.quantile(0.1), estimates.median(),
                      estimates.quantile(0.9), estimates.max(), len(estimates)])
volatility_cone = pd.DataFrame(cone_rows, index=window_years,
                              columns=["Minimum", "10th percentile", "Median", "90th percentile", "Maximum", "Windows"])
print(volatility_cone.to_string(float_format=lambda x: f"{x:.4f}"))
print("\n<ECHARTS_DATA>" + json.dumps({
    "volatilityCone": {
        "title": "Historical Volatility Cone",
        "xAxis": {"type": "value", "name": "Trailing window (years)", "min": 1},
        "yAxisName": "Annualized sample volatility (%)",
        "series": {name: [[years, float(volatility_cone.loc[years, name] * 100)] for years in window_years]
                   for name in volatility_cone.columns[:-1]},
    },
}, allow_nan=False))

Longer windows often give a narrower observed range, but nesting or monotonicity is not guaranteed. There are also fewer long-window estimates. Even iid constant-volatility returns produce variation in estimated volatility; a cone alone does not prove heteroskedasticity or an ARCH mechanism.

Exponential Weighting and Conditional Variance

Finite-History Exponential Weights

Equal weighting assigns $1/T$ to each observation. With a decay factor $0<\lambda \le 1$, normalized finite-history weights are:

$$a_t=\frac{{\lambda }^{T-t}}{\sum _{j=1}^T{\lambda }^{T-j}},\sum _t{a}_t=1.$$

Smaller $\lambda$ discounts older observations faster. For $\lambda <1$, the relative-weight half-life is $\log (0.5)/\log (\lambda )$ observation periods. A decay convention such as 0.94 historically used for daily RiskMetrics updates should not be transferred to monthly data without considering its very different calendar half-life.

observations = 24
decay_values = [1.0, 0.97, 0.94, 0.8]
exponential_weights = pd.DataFrame(index=np.arange(1, observations + 1))
for decay in decay_values:
    weights = decay ** np.arange(observations - 1, -1, -1)
    exponential_weights[f"Decay {decay:g}"] = weights / weights.sum()
print("Effective sample sizes 1/sum(weight²):")
print((1 / (exponential_weights ** 2).sum()).to_string(float_format=lambda x: f"{x:.2f}"))
print("\n<ECHARTS_DATA>" + json.dumps({
    "exponentialWeights": {
        "title": "Weights on 24 Monthly Observations",
        "xAxis": {"type": "value", "name": "Observation: oldest to newest", "min": 1, "max": observations},
        "yAxisName": "Observation weight (%)",
        "series": {name: [[int(t), float(value * 100)] for t, value in exponential_weights[name].items()]
                   for name in exponential_weights},
    },
}, allow_nan=False))

Effective sample size is a weight-concentration diagnostic; serial dependence can reduce information further. With $\lambda =1$, all weights are equal and the value is 24.

What the Toolkit's EWA Helper Measures

The existing helper annualize_vol_ewa() measures normalized weighted squared deviations around the ordinary arithmetic mean of the supplied history:

$${\hat{v}}_{\lambda }=\sum _t{a}_t(r_t-\overline{r}{)}^2,{\hat{\sigma }}_{\lambda ,\mathrm{ann}}=\sqrt{12{\hat{v}}_{\lambda }}.$$

There is no degrees-of-freedom correction. At $\lambda =1$, it equals standard deviation with ddof=0, annualized at the same frequency. It is neither variance around a weighted mean nor the recursive conditional-variance model below. In particular, choosing a center and estimating a forecast are separate modeling decisions.

Hand check. For returns $(-2\mathrm{\%},0\mathrm{\%},2\mathrm{\%})$ and $\lambda =0.5$, the ordinary mean is zero and normalized weights are $(1,2,4)/7$. Weighted variance is $(1/7+4/7)(0.02{)}^2=0.000285714$, giving annualized volatility about 5.855%.

ewa_example = pd.Series([-0.02, 0.0, 0.02], index=pd.period_range("2024-01", periods=3, freq="M"), name="Example")
worked_ewa = pok.annualize_vol_ewa(ewa_example, decay=0.5, periods_per_year=12)
equal_weight_volatility = pok.annualize_vol_ewa(ewa_example, decay=1.0, periods_per_year=12)
print(f"Decay 0.5 annualized dispersion: {worked_ewa:.3%}")
print(f"Decay 1 annualized dispersion: {equal_weight_volatility:.3%}")
print(f"Ordinary ddof=0 annualized volatility: {pok.annualize_vol(ewa_example, 12, ddof=0):.3%}")

Now compare endpoint estimates from the same market-proxy history. We reuse the loaded series; these numbers summarize information through December 2018, not estimates available throughout the displayed past.

endpoint_ewa = pd.Series({decay: pok.annualize_vol_ewa(market_returns, decay=decay, periods_per_year=12)
                         for decay in decay_values}, name="Annualized dispersion (%)") * 100
print(endpoint_ewa.to_string(float_format=lambda x: f"{x:.3f}"))

Try it: replace the input with market_returns.loc[:"2008-12"]. The estimated center and endpoint change because the available history changes. A lower decay can raise or lower estimated dispersion; it depends on the observations receiving more weight.

ARCH and GARCH: Explicit Conditional Dynamics

Let ${\epsilon }_t=R_t-\mathbb{E}[R_t\mid {\mathcal{F}}_{t-1}]$ be a return innovation, and let $h_t=\mathrm{Var}({\epsilon }_t\mid {\mathcal{F}}_{t-1})$. An ARCH($p$) model specifies:

$$h_t=\omega +\sum _{j=1}^p{\alpha }_j{\epsilon }_{t-j}^2.$$

Nonnegative ${\alpha }_j$ and positive $\omega$ ensure positive conditional variance. Under the usual standardized-innovation assumptions, $\sum _j{\alpha }_j<1$ gives a finite unconditional variance $\omega /(1-\sum _j{\alpha }_j)$. Writing $\omega =(1-\sum _j{\alpha }_j)V_L$ is a reparameterization in terms of that long-run variance, not a general requirement that a separately chosen set of coefficients sum to one.

Using $p$ for shock lags and $q$ for variance lags, a GARCH($p,q$) model is:

$$h_t=\omega +\sum _{j=1}^p{\alpha }_j{\epsilon }_{t-j}^2+\sum _{j=1}^q{\beta }_j{h}_{t-j}.$$

For GARCH(1,1), positive $\omega$, nonnegative $\alpha ,\beta$, and $\alpha +\beta <1$ give a finite long-run variance:

$$V_L=\frac{\omega }{1-\alpha -\beta },h_{t+1}=\omega +\alpha {\epsilon }_t^2+\beta h_t.$$

The response to a new squared innovation is controlled by $\alpha$; persistence of expected future variance depends on $\alpha +\beta$. Estimation also requires a conditional mean, a distributional or quasi-likelihood specification, initialization, and parameter uncertainty. Fitting univariate GARCH models alone does not estimate cross-asset correlations.

Hand check. Set $V_L=0.0004$, $\alpha =0.10$, and $\beta =0.85$, so $\omega =0.00002$. If the previous variance is 0.0004 and the observed innovation is 8%, the next variance is $0.00002+0.1(0.08{)}^2+0.85(0.0004)=0.001$, giving next-month volatility 3.162%.

long_run_variance = 0.02 ** 2
arch_coefficient, garch_coefficient = 0.10, 0.85
omega = (1 - arch_coefficient - garch_coefficient) * long_run_variance
observed_shock = 0.08
next_variance = omega + arch_coefficient * observed_shock ** 2 + garch_coefficient * long_run_variance
print(f"Omega: {omega:.6f}")
print(f"Next monthly variance: {next_variance:.6f}; volatility: {np.sqrt(next_variance):.3%}")

Recursive EWMA instead uses $h_{t+1}=\lambda h_t+(1-\lambda ){\epsilon }_t^2$, with a specified initial $h_0$. This corresponds to the zero-intercept boundary $\alpha +\beta =1$ and does not imply mean reversion to a fixed positive long-run variance. It also differs from the finite normalized weights and sample-mean centering in annualize_vol_ewa().

shocks = np.tile([0.02, -0.02], 12)
shocks[4] = 0.08
recursive_decay = 0.90
garch_variances = [long_run_variance]
ewma_variances = [long_run_variance]
for shock in shocks:
    garch_variances.append(omega + arch_coefficient * shock ** 2 + garch_coefficient * garch_variances[-1])
    ewma_variances.append(recursive_decay * ewma_variances[-1] + (1 - recursive_decay) * shock ** 2)
print("\n<ECHARTS_DATA>" + json.dumps({
    "conditionalVariance": {
        "title": "Specified Shock Path: GARCH and Recursive EWMA",
        "xAxis": {"type": "value", "name": "Number of observed monthly shocks", "min": 0},
        "yAxisName": "Next-month volatility (%)",
        "series": {
            "GARCH(1,1)": [[i, float(np.sqrt(v) * 100)] for i, v in enumerate(garch_variances)],
            "Recursive EWMA": [[i, float(np.sqrt(v) * 100)] for i, v in enumerate(ewma_variances)],
        },
    },
}, allow_nan=False))

Point zero is the initial variance forecast; point five incorporates the fifth observed shock. These are illustrative parameters applied to a prescribed shock sequence, not fitted market forecasts. A particular realized path need not converge to the unconditional variance.

Factor GARCH

Factor GARCH combines factor structure with time-varying risk, for example:

$${\mathrm{\Sigma }}_t=B{\mathrm{\Omega }}_t{B}^T+D_t.$$

Conditional factor covariance ${\mathrm{\Omega }}_t$ and residual variances in $D_t$ evolve over time; loadings may be fixed or modeled separately. The dimension reduction comes from modeling $K$ common factors for $N$ assets, usually with $K\ll N$. It does not mean replacing the $q$ variance lags of a univariate GARCH model with $K<q$ factors. Dynamic correlations and residual dependence remain substantive choices.

Portfolio Construction with Different Covariance Estimates

One Universe, One Evaluation Calendar

Use the 49 value-weighted within-industry return portfolios from 1974–2018. We allocate across these same industry portfolios using equal weights, lagged capitalization weights, or long-only minimum variance. This avoids confusing equal weighting within industries with equal allocations across industries.

Use 48 monthly observations for each estimation window. That leaves 492 evaluated months, January 1978–December 2018. The first window ends in December 1977, before the first evaluated return. With 49 assets and 48 observations, each sample covariance has rank at most 47: an intentional example of high-dimensional estimation.

industry_returns = pok.get_ind_file(filetype="rets", nind=49, ew=False).loc["1974":"2018"]
industry_counts = pok.get_ind_file(filetype="nfirms", nind=49).loc[industry_returns.index]
industry_sizes = pok.get_ind_file(filetype="size", nind=49).loc[industry_returns.index]
industry_capitalization = industry_counts * industry_sizes
industry_cap_weights = industry_capitalization.div(industry_capitalization.sum(axis=1), axis=0)
cash_returns = pok.get_fff_returns()["RF"]
estimation_months = 48
first_window = industry_returns.iloc[:estimation_months]
print(f"Universe: {industry_returns.shape[1]} industries; history: {len(industry_returns)} months")
print(f"First training window: {first_window.index[0]} to {first_window.index[-1]}")
print(f"First evaluation month: {industry_returns.index[estimation_months]}")

Capitalization is firm count times average firm size; select the common period before calculating these fractions because the earliest 49-industry metadata contain missing observations. The following snapshot has complete inputs. Lagging capitalization avoids using the evaluated month's end-of-month weights, but does not establish point-in-time publication availability.

Compare the Initial Risk Estimates

Inspect the first training window before looking at strategy performance. All three matrices estimate monthly covariance on the same data. They share diagonal sample variances; they differ in the cross-asset structure.

initial_covariances = {
    "Sample": pok.sample_cov(first_window),
    "Constant correlation": pok.cc_cov(first_window),
    "50% shrinkage": pok.shrinkage_cov(first_window, delta=0.5),
}
eigenvalues = pd.DataFrame({name: np.linalg.eigvalsh(cov) for name, cov in initial_covariances.items()},
                         index=np.arange(1, industry_returns.shape[1] + 1))
print(pd.DataFrame({name: [np.linalg.matrix_rank(cov), np.linalg.eigvalsh(cov).min(), np.linalg.cond(cov)]
                    for name, cov in initial_covariances.items()},
                   index=["Numerical rank", "Minimum eigenvalue", "Condition number"]).T.to_string())
print("\n<ECHARTS_DATA>" + json.dumps({
    "covarianceSpectrum": {
        "title": "First Training Window: Covariance Eigenvalues",
        "xAxis": {"type": "value", "name": "Eigenvalue rank (ascending)", "min": 1},
        "yAxisName": "Monthly variance (squared percentage points)",
        "series": {name: [[int(i), float(value * 10000)] for i, value in eigenvalues[name].items()]
                   for name in eigenvalues},
    },
}, allow_nan=False))

Variance in decimal-return units is multiplied by ${100}^2=\mathrm{10,000}$ for squared percentage points. The condition number of a singular matrix is theoretically infinite; a finite but enormous numerical value is another indication of near-zero computed eigenvalues. Better conditioning alone does not guarantee lower realized portfolio risk.

Baselines and Minimum-Variance Allocations

The long-only minimum-variance problem is:

$$\underset{\mathbf{w}}{min}{\mathbf{w}}^T\hat{\mathrm{\Sigma }}\mathbf{w},{\mathbf{1}}^T\mathbf{w}=1,w_i\ge 0.$$

No expected-return forecast is needed. Multiplying covariance by a common positive annualization factor changes reported variance, not the minimizing weights. The helper solves the convex quadratic objective, returns labeled weights, and checks the solver result. It does not invert a singular sample covariance matrix.

First compute EW and capitalization baselines. Both use the same warm-up and evaluation months as the covariance-based strategies. Capitalization weights come from the last training month, rather than the beginning of the 48-month window.

portfolio_returns = pd.DataFrame({
    "EW": pok.backtest_weight_scheme(industry_returns, window=estimation_months, weight_scheme=pok.weight_ew),
    "CW": pok.backtest_weight_scheme(industry_returns, window=estimation_months,
                                    weight_scheme=pok.weight_cw, cap_ws=industry_cap_weights),
})
print(portfolio_returns.iloc[estimation_months:estimation_months+3].to_string())

Each of the next three cells estimates its covariance from months $t-48$ through $t-1$, chooses weights, then earns the industry returns in month $t$. Keep the same universe and constraints to isolate the estimator change.

portfolio_returns["GMV sample"] = pok.backtest_weight_scheme(
    industry_returns, window=estimation_months, weight_scheme=pok.weight_minvar,
    cov_estimator=pok.sample_cov, periods_per_year=12,
)
print(portfolio_returns["GMV sample"].dropna().head().to_string())

portfolio_returns["GMV constant correlation"] = pok.backtest_weight_scheme(
    industry_returns, window=estimation_months, weight_scheme=pok.weight_minvar,
    cov_estimator=pok.cc_cov, periods_per_year=12,
)
print(portfolio_returns["GMV constant correlation"].dropna().head().to_string())

shrinkage_intensity = 0.5
portfolio_returns["GMV shrinkage"] = pok.backtest_weight_scheme(
    industry_returns, window=estimation_months, weight_scheme=pok.weight_minvar,
    cov_estimator=pok.shrinkage_cov, periods_per_year=12, delta=shrinkage_intensity,
)
print(portfolio_returns["GMV shrinkage"].dropna().head().to_string())

Evaluate Realized Outcomes on the Common Sample

Remove the explicitly identified warm-up rather than allowing each strategy to choose a different evaluation sample. Use the observed monthly cash series for all Sharpe calculations. Report compounded growth separately from the arithmetic excess-return measure used by Sharpe.

evaluation_returns = portfolio_returns.iloc[estimation_months:]
portfolio_excess = evaluation_returns.sub(cash_returns.loc[evaluation_returns.index], axis=0)
portfolio_statistics = pd.DataFrame({
    "CAGR (%)": pok.annualize_rets(evaluation_returns, 12) * 100,
    "Annualized volatility (%)": pok.annualize_vol(evaluation_returns, 12) * 100,
    "Sharpe": np.sqrt(12) * portfolio_excess.mean() / portfolio_excess.std(ddof=1),
    "Maximum drawdown (%)": evaluation_returns.apply(lambda r: -pok.drawdown(r)["Drawdown"].min() * 100),
})
print(f"Evaluation: {evaluation_returns.index[0]} to {evaluation_returns.index[-1]}, {len(evaluation_returns)} months")
print(portfolio_statistics.to_string(float_format=lambda x: f"{x:.3f}"))
portfolio_wealth = 100 * (1 + evaluation_returns).cumprod()
portfolio_wealth.loc[evaluation_returns.index[0] - 1] = 100.0
portfolio_wealth = portfolio_wealth.sort_index()
print("\n<ECHARTS_DATA>" + json.dumps({
    "dates": portfolio_wealth.index.astype(str).tolist(),
    "covariancePortfolioWealth": {
        "title": "Common Evaluation Sample: Covariance-Based Allocations",
        "yAxis": {"type": "log", "name": "Wealth per 100 initial units"},
        "series": {name: portfolio_wealth[name].tolist() for name in portfolio_wealth},
    },
}, allow_nan=False))

All strategies rebalance monthly and omit transaction costs. A retrospective comparison can reveal realized differences, but the best ending wealth does not establish the best covariance forecast. GMV targets variance, not return; realized risk can differ from estimated risk. Selection of the universe, window, and shrinkage intensity also affects the result.

With the default sample, annualized realized volatility is approximately 16.22% for EW, 15.18% for CW, 11.59% for sample GMV, 11.90% for constant-correlation GMV, and 11.66% for shrinkage GMV. Shrinkage improves the initial matrix's conditioning but does not deliver the lowest realized volatility in this comparison. That distinction is why an estimator's numerical properties and its investment outcomes must be assessed separately.

Try it: change shrinkage_intensity to 0 or 1 and rerun the shrinkage and evaluation cells. These endpoints should reproduce the sample or constant-correlation strategy, within numerical tolerance. Choosing the best intensity after examining this same evaluation period would turn the exercise into in-sample tuning.

Inspect Weights at a Stated Decision Date

To understand how the estimators allocate risk, compare the portfolios chosen for January 2009, using only returns through December 2008. This is a single decision snapshot, not a fit using the entire post-2009 sample.

decision_month = pd.Period("2009-01", freq="M")
decision_window = industry_returns.loc[:decision_month-1].iloc[-estimation_months:]
decision_weights = pd.DataFrame({
    "EW": pok.weight_ew(decision_window),
    "CW": pok.weight_cw(decision_window, cap_ws=industry_cap_weights),
    "GMV sample": pok.weight_minvar(decision_window, cov_estimator=pok.sample_cov),
    "GMV constant correlation": pok.weight_minvar(decision_window, cov_estimator=pok.cc_cov),
    "GMV shrinkage": pok.weight_minvar(decision_window, cov_estimator=pok.shrinkage_cov, delta=shrinkage_intensity),
})
print(f"Training: {decision_window.index[0]} to {decision_window.index[-1]}; decision: {decision_month}")
print((decision_weights * 100).to_string(float_format=lambda x: f"{x:.2f}"))
print("Effective number of holdings, 1/sum(weight²):")
print((1 / (decision_weights ** 2).sum()).to_string(float_format=lambda x: f"{x:.2f}"))
largest_industries = decision_weights.max(axis=1).nlargest(8).index
print("\n<ECHARTS_DATA>" + json.dumps({
    "decisionWeights": {
        "type": "bar", "title": "January 2009 Decision: Eight Largest Industry Exposures",
        "xAxis": {"type": "category", "data": largest_industries.tolist()},
        "yAxisName": "Portfolio weight (%)",
        "series": {name: (decision_weights.loc[largest_industries, name] * 100).tolist() for name in decision_weights},
    },
}, allow_nan=False))

The table contains all industries; the chart selects the eight with the largest maximum weight across strategies. The effective number of holdings describes weight concentration, not independent risk sources. A highly concentrated GMV portfolio can still have low estimated variance, while a widely spread portfolio can retain substantial common-factor exposure.

For this decision date, sample GMV has about 5.74 effective holdings, constant-correlation GMV 3.48, and shrinkage GMV 4.45, versus 49 for equal weighting. Imposing a simpler covariance structure does not necessarily produce a more evenly spread allocation.

Additional Risk Diagnostics: VaR and Expected Shortfall

Covariance summarizes second moments. It does not determine tail risk without further distributional assumptions. Revisit the historical loss measures from Understanding Returns and Assessing Risks with Value at Risk using the same market proxy: VaR is a loss threshold, and ES averages the worst tail. Monthly returns produce monthly risk measures.

tail_levels = [0.10, 0.05, 0.01]
tail_summary = pd.DataFrame({
    "Confidence": [f"{1-level:.0%}" for level in tail_levels],
    "Historical VaR (%)": [pok.var_historic(market_returns, level=level) * 100 for level in tail_levels],
    "Historical ES (%)": [pok.cvar_historic(market_returns, level=level) * 100 for level in tail_levels],
    "Tail mass (observations)": [len(market_returns) * level for level in tail_levels],
}).set_index("Confidence")
print(tail_summary.to_string(float_format=lambda x: f"{x:.3f}"))

Now estimate 99% VaR and ES from each previous 36-month window, then compare the next observed loss with those estimates. The return sign matters: define loss $L_t=-R_t$ and plot positive loss thresholds against losses. Mirroring VaR and ES above and below a return series would incorrectly suggest symmetric uncertainty bands.

tail_window, tail_probability = 36, 0.01
lagged_var = market_returns.rolling(tail_window).apply(
    lambda r: pok.var_historic(r, level=tail_probability), raw=False).shift(1)
lagged_es = market_returns.rolling(tail_window).apply(
    lambda r: pok.cvar_historic(r, level=tail_probability), raw=False).shift(1)
tail_diagnostics = pd.DataFrame({"Observed loss": -market_returns, "Prior-window VaR": lagged_var,
                                "Prior-window ES": lagged_es}).loc["1990":"2018"]
exceedances = tail_diagnostics["Observed loss"] > tail_diagnostics["Prior-window VaR"]
print(f"Tail mass per estimation window: {tail_window * tail_probability:.2f} observations")
print(f"Realized VaR exceedances: {exceedances.sum()} of {len(exceedances)} months ({exceedances.mean():.2%})")
print("\n<ECHARTS_DATA>" + json.dumps({
    "dates": tail_diagnostics.index.astype(str).tolist(),
    "laggedTailRisk": {
        "title": "Monthly Losses Versus Prior-Window 99% Tail Estimates",
        "yAxisName": "Monthly loss (% of starting value)",
        "series": {name: (tail_diagnostics[name] * 100).tolist() for name in tail_diagnostics},
    },
}, allow_nan=False))

At 99% confidence, 36 observations supply only 0.36 observations' worth of tail probability. Under the toolkit's equal-probability empirical ES convention, that tail lies within the single worst observed return, so ES equals that observed loss. Historical VaR uses linear interpolation and is a different finite-sample estimator. Neither reliably identifies an unobserved extreme tail from such a short sample.

The exceedance rate summarizes this realized evaluation history; it is not a known future probability. Overlapping windows and changing risk also make observations dependent. Covariance estimates, volatility models, and tail diagnostics answer related but distinct questions.

Chapter Summary

• Covariance estimates need a common horizon, consistent asset ordering, and a justified sample. PSD and numerical conditioning matter for portfolio construction.
• Sample covariance retains cross-sectional detail but can be noisy or singular. Constant correlation, factor structure, and shrinkage impose different restrictions.
• Window length determines the history used; return frequency determines annualization. A 60-month window still uses $\sqrt{12}$ for monthly volatility.
• Finite-history exponential dispersion, recursive EWMA, and GARCH conditional variance have different centering, initialization, and dynamic assumptions.
• Compare portfolio estimators on the same assets and evaluation dates, using only information available before each return. Report concentration and implementation assumptions alongside realized performance.
• Covariance does not describe the entire loss distribution. Tail diagnostics require their own definitions, sample-size assessment, and timing checks.

Next, Black–Litterman combines reference-implied expected returns with investor views and their uncertainty, keeping the covariance and allocation choices explicit.