Step 1 - Specify the Prior

Notation

Unlike the first track's Step 1, where $\mathbf{W}$ is fixed at its classical value and only the score distribution's $\mu_i,\sigma_i^2$ get priors, here $\mathbf{W}$ itself is a free parameter:

$$\mathbf{W} \sim \text{Normal}(0, 1), \qquad \mathbf{W} \in \mathbb{R}^{N \times 3}$$

$$\sigma \sim \text{HalfNormal}(1)$$

$\boldsymbol{\mu}$ is fixed at $\mathbf{0}$ rather than given a prior — the data this model will eventually be fit to ($\widetilde{\mathbf{Y}}_c$, Classical PCA's Step 5 output) is already column-demeaned, so there's no free mean left to estimate.

These priors on $\mathbf{W}$ and $\sigma$ are proper, not flat/improper. The eventual likelihood only depends on $\mathbf{W}\mathbf{W}^\top$, so $\mathbf{W}$ is non-identified up to rotation — combined with an improper prior, that non-identifiability would leave the posterior non-normalizable. No identifiability constraint is imposed on $\mathbf{W}$ at this stage; that decision is deferred to whichever step first attaches a likelihood and samples the model.

No likelihood is attached at this stage, so there is nothing to sample or plot yet — same starting point as the first track's Step 1.

Python

filename: step1_prior.py

import pymc as pm

'''
prior for the full-latent model: W (loadings, N x k) and sigma (noise scale) are free
parameters here, in contrast to bayesian/step1_prior.py where W is kept fixed at its
classical value. mu is fixed at 0 rather than given a prior, since the input data
(dataframes/step5_demeaning.py::demeaning()) is already column-demeaned.

W's prior must be proper (not flat), unlike the mu_i priors in bayesian/step1_prior.py.
the likelihood here will only depend on W @ W.T, so W is non-identified up to rotation -
combined with an improper prior that would give a non-normalizable posterior. no
identifiability constraint is imposed yet; that's deferred to the sampling step.
'''

def prior_model(n_maturities=71, n_components=3):
    with pm.Model() as model:
        W = pm.Normal("W", mu=0, sigma=1, shape=(n_maturities, n_components))
        sigma = pm.HalfNormal("sigma", sigma=1)
    return model