Add Antweiler-Freyberger (2025) iterative quadrature estimator

.pre-commit-config.yaml

@@ -65,6 +65,9 @@ repos:
     rev: 0.9.1
     hooks:
       - id: nbstripout
+        # The getting-started tutorial ships pre-rendered (`execute: false` in
+        # myst.yml); stripping its outputs would leave the docs site empty.
+        exclude: ^docs/getting_started/tutorial\.ipynb$
         args:
           - --extra-keys
           - metadata.kernelspec metadata.language_info.version metadata.vscode

CLAUDE.md

@@ -259,6 +259,10 @@ When writing new public-facing code, always accept and return `period`. Convert
 
 ## Testing
 
-- pytest with markers: `wip`, `unit`, `integration`, `end_to_end`
+- pytest with markers: `wip`, `unit`, `integration`, `end_to_end`, `long_running`
 - Test files mirror source structure in `tests/`
 - Memory profiling available via pytest-memray (Unix only)
+- MATLAB AF CES / translog reproduction tests live in the parent workspace at
+  `../matlab_ces_repro/` (alongside `sim_repro/`), not in this library. They depend on
+  reference data at `/home/hmg/sciebo/Skill estimation/` and the CNLSY xls bundled
+  beside them. Run from the workspace root.

docs/explanations/architecture.md

@@ -0,0 +1,100 @@
+# Package Architecture
+
+Skillmodels hosts three estimators under one model specification. The package
+layout reflects that:
+
+```
+src/skillmodels/
+├── common/                  Estimator-agnostic machinery
+│   ├── model_spec.py          ModelSpec, FactorSpec, AnchoringSpec, Normalizations
+│   ├── types.py               ProcessedModel, Dimensions, Labels,
+│   │                          EndogenousFactorsInfo, ParsingInfo, ...
+│   ├── process_model.py       process_model(spec) -> ProcessedModel
+│   ├── process_data.py        long-format data -> internal arrays
+│   ├── params_index.py        4-level MultiIndex used by all estimators
+│   ├── parse_params.py        flat vector <-> structured params
+│   ├── constraints.py         get_constraints, FixedConstraintWithValue,
+│   │                          collect_fixed_locs, project_to_probability_constraints
+│   ├── selector.py            select_by_loc, align_index_names
+│   ├── transition_functions.py  linear / translog / log_ces / ...
+│   ├── transitions.py         apply_anchored_transition (sigma-points-agnostic)
+│   ├── anchoring.py           anchor / unanchor states
+│   ├── state_ranges.py        create_state_ranges
+│   ├── simulate_data.py       simulate_dataset, simulate_policy_effect
+│   ├── variance_decomposition.py  signal/noise decomposition
+│   └── diagnostic_plots.py    plot_residual_boxplots, plot_likelihood_contributions
+├── chs/                     Cunha-Heckman-Schennach Kalman MLE
+│   ├── options.py             CHSEstimationOptions
+│   ├── kalman_filters.py      square-root unscented Kalman filter
+│   ├── likelihood.py          jitted log-likelihood
+│   ├── likelihood_debug.py    non-jitted variant with debug arrays
+│   ├── maximization_inputs.py get_maximization_inputs(...)
+│   ├── filtered_states.py     get_filtered_states(...)
+│   └── process_debug_data.py  post-process Kalman debug arrays
+├── af/                      Antweiler-Freyberger sequential Halton MLE
+│   ├── types.py               AFEstimationOptions, AFEstimationResult, ...
+│   ├── estimate.py            estimate_af(...) -- top-level orchestration
+│   ├── initial_period.py      period-0 mixture + measurement system MLE
+│   ├── transition_period.py   period-t transition + measurement-system MLE
+│   ├── likelihood.py          jitted period-specific log-likelihoods
+│   ├── halton.py              quadrature nodes / weights
+│   ├── batching.py            obs-batching for the autodiff chunking
+│   ├── posterior_states.py    conditional-distribution materialisation
+│   ├── inference.py           compute_af_standard_errors (cluster bootstrap)
+│   └── jaxopt_backend.py      on-device L-BFGS-B alternative
+└── amn/                     Attanasio-Meghir-Nix 2020 (three-stage)
+    ├── types.py               AMNEstimationOptions, ...
+    ├── estimate.py            estimate_amn(...) -- top-level orchestration
+    ├── mixture_em.py          Stage 1: EM on the augmented mixture
+    ├── minimum_distance.py    Stage 2: structural recovery
+    ├── simulate_and_regress.py Stage 3: synthetic-panel regression
+    ├── moments.py             Spearman + Bartlett start-values
+    ├── start_values.py        get_spearman_start_params, pool_equality_groups
+    ├── posterior_states.py    simulate factor paths from fitted mixture
+    └── inference.py           compute_amn_standard_errors (cluster bootstrap)
+```
+
+## How the layers interact
+
+Every estimator reads the same `ModelSpec` and produces the same canonical
+params DataFrame (4-level MultiIndex
+`(category, period, name1, name2)`). The differences live entirely below the
+spec:
+
+- **CHS** consumes `process_model(spec) -> ProcessedModel`, then plugs that
+  into the Kalman recursion. `CHSEstimationOptions` is passed at call time
+  to `get_maximization_inputs(spec, data, chs_options=...)`.
+- **AF** also calls `process_model`, but uses `ProcessedModel` only for the
+  parameter index, labels, and transition info. The Kalman filter is not
+  invoked; period-specific Halton designs replace the predict step.
+- **AMN** likewise calls `process_model` for the index/labels, then runs its
+  three-stage pipeline. The result re-uses the same params DataFrame format
+  so the AMN output can seed CHS or AF estimation when desired.
+
+`process_model` itself is structural: it takes only the spec and produces
+shapes, labels, transition info, and an `EndogenousFactorsInfo`. It does not
+carry any estimator-specific tuning. Each estimator's options class
+(`CHSEstimationOptions`, `AFEstimationOptions`, `AMNEstimationOptions`) is
+passed in at call time.
+
+## Why this split
+
+The package grew organically: CHS was the original codebase; AF and AMN were
+later additions. Earlier iterations stored CHS-only options on `ModelSpec`,
+which made the spec leak CHS assumptions into a notionally agnostic container.
+The split into `common/`, `chs/`, `af/`, `amn/` makes the scope of each piece
+explicit at the import site:
+
+- `from skillmodels import ModelSpec` — pure structural description.
+- `from skillmodels.chs import CHSEstimationOptions, get_maximization_inputs`
+  — CHS-specific.
+- `from skillmodels.af import estimate_af, AFEstimationOptions` — AF-specific.
+- `from skillmodels.common.variance_decomposition import decompose_measurement_variance`
+  — works for any estimator, given pre-computed filtered states.
+
+The architectural principle: a function lives in `common/` iff it does not
+import from `chs/`, `af/`, or `amn/`. Anything that does belongs in the
+relevant subpackage. There is one practical exception:
+`CHSEstimationOptions` is defined in `chs/options.py` but the
+`process_model` orchestration in `common/` doesn't read it (it reads the
+structural `ModelSpec.n_mixtures` field instead), so the layering is clean.

docs/explanations/names_and_concepts.md

@@ -76,12 +76,47 @@ of factors are arbitrary).
 
 ## Estimation Options
 
-The `EstimationOptions` dataclass controls numerical aspects:
+Each estimator has its own options dataclass, passed at call time rather than
+embedded in `ModelSpec`. The three classes share no fields — what counts as a
+tuning knob differs between estimators.
+
+`CHSEstimationOptions` (from `skillmodels.chs`) controls the Kalman MLE:
 
 - **robust_bounds**: Tightens parameter bounds to avoid numerical issues
 - **bounds_distance**: How much stricter to make bounds (zeroed if robust_bounds is
   false)
-- **n_mixtures**: Number of mixture components in the distribution
 - **sigma_points_scale**: Controls spread of sigma points in unscented Kalman filter
 - **clipping_\***: Parameters for soft-clipping the log-likelihood to prevent
   infinities
+- **start_params_strategy**: How to seed the `params_template`. `"amn"` (default)
+  runs the full AMN three-stage estimator and uses its parameters as the start;
+  `"spearman"` uses moment-based start values; `"none"` leaves entries as NaN
+  for the caller to fill in.
+
+`AFEstimationOptions` (from `skillmodels.af`) controls the sequential MLE:
+
+- **n_halton_points**, **n_halton_points_shock**: quadrature counts.
+- **n_mixture_components**: number of components in the latent-factor mixture.
+- **optimizer_backend**: `"auto"` (default), `"optimagic"`, or `"jaxopt"`. Auto
+  picks `"jaxopt"` if a JAX GPU is visible and the model has no probability or
+  equality constraints; otherwise `"optimagic"`.
+- **optimizer_algorithm**: the optimagic algorithm name used when the backend
+  is `"optimagic"`. Ignored under `"jaxopt"`.
+- **initialization_strategy**: `"amn"`, `"spearman"`, or `"constant"`. Same
+  meaning as in CHS.
+
+`AMNEstimationOptions` (from `skillmodels.amn`) controls the three-stage
+pipeline:
+
+- **n_mixture_components**: Stage-1 EM components.
+- **em_max_iter**, **em_tol**, **em_n_init**, **em_reg_covar**: Stage-1 EM
+  numerical knobs.
+- **n_simulation_draws**: Stage-3 synthetic-panel size.
+- **minimum_distance_weighting**: Stage-2 weighting; `"identity"` (default) or
+  `"optimal"`.
+- **investment_endogeneity**: include the control-function residual in Stage 3
+  for endogenous-investment models.
+
+The shared structural field — number of mixture components in the latent
+distribution — lives directly on `ModelSpec.n_mixtures`, since it changes the
+model itself rather than the optimizer.