Feature: extend simulate() to residual covariances (~~) and HSGP terms

[drbenvincent]

Problem

simulate() currently raises NotImplementedError for:

• Residual covariances (Y1 ~~ Y2) — correlated equation errors
• hsgp() smooth terms — useful for flexible trends in MMM

Both compile in model() for cross-sectional data (HSGP); ~~ works cross-sectionally in estimation.

Scope

Phase 1: Cross-sectional ~~ in simulate() — draw correlated residuals via Cholesky of fixed covariance params in params dict.

Phase 2: Cross-sectional hsgp() in simulate() — clamp HSGP basis coefficients via pm.do(), draw smooth contribution.

Phase 3 (separate / deferred): ~~ and HSGP in scan-compiled panel models (currently limited in estimation too).

Acceptance criteria

• Remove NotImplementedError guards in simulate() for supported cases
• Tests mirroring tests/test_simulate_data.py validation patterns
• Clear errors for unsupported panel+HSGP combinations (already enforced in compile)

Related

• Parent meta: Meta: MMM generative workflows as a stress-test for pathmc panel + simulate APIs #426
• HSGP feature: Feature: Hilbert Space Gaussian Process (HSGP) term for smooth nonparametric effects #360
• Fourier/seasonality alternative: Generalize the HSGP basis machinery into a reusable basis-term abstraction (with fourier() as the second member) #371

Decision (2026-08-24)

Resolved during meta-issue run #426 pre-flight:

• Phase 1 accepts chol_{block} as the packed lower-triangular Cholesky vector (PyMC LKJCholeskyCov packing), matching what introspect already advertises. Block-member noisy draws propagate as mu_{var} deterministics into descendant equations within one draw pass (documented behavior, mirrors estimation). Phase-scoped docstring edits only: hsgp guard stays until Phase 2.
  427:Cold start at zero for temporal state (matches compile_scan_panel init semantics), documented in the API. params must supply every free RV including alpha{var} vectors under pooling='partial'; scalar values are broadcast to the required shape. Implementation may reuse the scan step functions directly rather than run_do_panel_unified's idata machinery — whichever is simpler.
  429:Template returns a plain dict[str, shape/dtype descriptor] built from the same zero-filled placeholder compile simulate() already performs internally; a typed dataclass wrapper is optional sugar, not required. Placeholder data needs real column names/dtypes but not plausible values.
  432:New module pathmc/scaling.py (AGENTS.md module list gains it); Scaling exported from pathmc/init. Methods: max / mean / fixed / divide-by-external-grid, computed per panel dimension group. simulate() applies inverse scaling only when scaling= is passed; default None = no scaling.

[drbenvincent]

Implemented on run/426-mmm-simulate-stress-test (commit 9f106ac), meta-run #426.

• Phase 1: simulate() draws correlated residuals per ~~ block from chol_{block} (packed LKJCholeskyCov form); block members propagate as mu_{var} deterministics into descendant equations (mirrors estimation; documented in the docstring).
• Phase 2: cross-sectional hsgp() flows through the existing pm.do/draw path — ell_/eta_/beta_hsgp_{lhs}_{var} are ordinary free RVs. Panel+HSGP rejection in compile untouched.
• Wrong-length chol_{block} raises ValueError naming problem + fix.
• Tests: test_residual_cov_raises → round-trip test asserting X-residualized empirical corr ≈ target (tol 0.05); test_simulate_with_hsgp_raises → positive test incl. zero-noise smooth-reproduction check.

Verification: uv run pytest tests/test_simulate_data.py tests/test_hsgp_compile.py tests/test_residual_cov.py tests/test_hsgp_do.py -q → 45 passed. Smoke: n=20000 block sim, empirical corr 0.6942 vs target 0.70.

Left open per issue scope: Phase 3 (panel ~~ / HSGP) stays deferred.