Step 2 - Build the Probabilistic Model

Notation

Step 1 defined proper priors on $\mathbf{W}$ and $\sigma$, with $\boldsymbol{\mu}$ fixed at $\mathbf{0}$. This step attaches the likelihood, completing the model:

$$\widetilde{\mathbf{Y}}_c \sim \text{Normal}(\mathbf{0},\ \mathbf{W}\mathbf{W}^\top + \sigma^2\mathbf{I}), \qquad \widetilde{\mathbf{Y}}_c \in \mathbb{R}^{(T-12) \times N}$$

where $\widetilde{\mathbf{Y}}_c$ is the demeaned yield matrix from Classical PCA's Step 5, and each row (date) is modelled as an independent draw from the same $N$-dimensional Normal, sharing one candidate $(\mathbf{W}, \sigma)$ across all dates. Unlike the first track, where the covariance implicit in each component's likelihood is diagonal by construction (independent $z_i$), here $\mathbf{W}\mathbf{W}^\top + \sigma^2\mathbf{I}$ is a full $N \times N$ matrix — see Where the covariance matrix fits in for why that takes over $\mathbf{A}$'s role from classical Step 6.

The model built here is still unsampled and $\mathbf{W}$ is still non-identified up to rotation — Step 3 runs MCMC over it to get the posterior; identifiability is addressed there, not here.

Python

filename: step2_model.py

import pymc as pm
from dataframes.step5_demeaning import demeaning
from bayesian_latent.step1_prior import prior_model

'''
extends the step 1 prior with a multivariate Normal likelihood on the demeaned yield matrix:
each row (date) is an independent draw from Normal(0, W @ W.T + sigma^2 * I), the covariance
implied by the current candidate (W, sigma) - see bayesian_latent_pca.md's "Where the
covariance matrix fits in" for why this plays A's role from classical Step 6.
'''

def build_model(n_components=3):
    df = demeaning()
    value_cols = df.columns[1:]
    Y_c = df[value_cols].values
    n_maturities = Y_c.shape[1]
    model = prior_model(n_maturities=n_maturities, n_components=n_components)
    with model:
        W = model["W"]
        sigma = model["sigma"]
        cov = pm.math.dot(W, W.T) + sigma**2 * pm.math.eye(n_maturities)
        pm.MvNormal(
            "Y",
            mu=0,
            cov=cov,
            observed=Y_c,
        )
    return model