diff --git a/.ai-instructions b/.ai-instructions
index f2744aab..68dc3f48 160000
--- a/.ai-instructions
+++ b/.ai-instructions
@@ -1 +1 @@
-Subproject commit f2744aabcc974fd101240c2e4393de232f45b061
+Subproject commit 68dc3f4867cb075b0e0d8eadf383800ea33f5995
diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index ca658f18..9e344a72 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -28,7 +28,7 @@ jobs:
- uses: actions/checkout@v6
- uses: prefix-dev/setup-pixi@v0.9.4
with:
- pixi-version: v0.66.0
+ pixi-version: v0.70.2
cache: true
cache-write: ${{ github.event_name == 'push' && github.ref_name == 'main' }}
environments: tests-cpu
@@ -54,7 +54,7 @@ jobs:
- uses: actions/checkout@v6
- uses: prefix-dev/setup-pixi@v0.9.4
with:
- pixi-version: v0.66.0
+ pixi-version: v0.70.2
cache: true
cache-write: ${{ github.event_name == 'push' && github.ref_name == 'main' }}
frozen: true
diff --git a/.gitmodules b/.gitmodules
index 10e09e46..66ab1130 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,4 +1,4 @@
[submodule ".ai-instructions"]
path = .ai-instructions
- url = git@github.com:OpenSourceEconomics/ai-instructions.git
+ url = https://github.com/OpenSourceEconomics/ai-instructions.git
branch = make-submodule
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 2877891b..680b6cc0 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -5,7 +5,7 @@ repos:
- id: check-hooks-apply
- id: check-useless-excludes
- repo: https://github.com/tox-dev/pyproject-fmt
- rev: v2.19.0
+ rev: v2.25.0
hooks:
- id: pyproject-fmt
- repo: https://github.com/lyz-code/yamlfix
@@ -47,7 +47,7 @@ repos:
hooks:
- id: yamllint
- repo: https://github.com/astral-sh/ruff-pre-commit
- rev: v0.15.6
+ rev: v0.15.18
hooks:
- id: ruff-check
args:
@@ -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
diff --git a/CLAUDE.md b/CLAUDE.md
index 5b6d3709..79fb7af3 100644
--- a/CLAUDE.md
+++ b/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.
diff --git a/GEMINI.md b/GEMINI.md
new file mode 100644
index 00000000..7dc4841d
--- /dev/null
+++ b/GEMINI.md
@@ -0,0 +1,2 @@
+@.ai-instructions/profiles/tier-a.md @.ai-instructions/modules/jax.md
+@.ai-instructions/modules/optimagic.md
diff --git a/docs/explanations/architecture.md b/docs/explanations/architecture.md
new file mode 100644
index 00000000..7ef049be
--- /dev/null
+++ b/docs/explanations/architecture.md
@@ -0,0 +1,105 @@
+# 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 / translog_af /
+│ │ log_ces_af / ...
+│ ├── 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
+│ ├── validate.py validate_af_model, kappa-scope / leakage checks
+│ ├── params.py per-period optimagic params + constraints
+│ ├── 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 posterior means from the chained sample
+│ └── inference.py compute_af_standard_errors (propagated
+│ influence-function score bootstrap)
+└── 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
+ │ (optional endogenous-investment control function)
+ ├── 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,
+ per-replicate seeds + nonconvergence filtering)
+```
+
+## 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.
diff --git a/docs/explanations/names_and_concepts.md b/docs/explanations/names_and_concepts.md
index 69e9e2d5..6db976a3 100644
--- a/docs/explanations/names_and_concepts.md
+++ b/docs/explanations/names_and_concepts.md
@@ -76,12 +76,57 @@ 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_algorithm**: the optimagic algorithm name passed to
+ `optimagic.minimize(algorithm=...)` (default `"fides"`; use
+ `"scipy_lbfgsb"` for MC sweeps).
+- **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) is
+ the paper's unweighted identity-metric criterion over per-component means and
+ full covariance matrices, and is currently the only implemented option.
+ `"optimal"` is reserved for a future Avar-weighted criterion and raises
+ `NotImplementedError`.
+- **investment_endogeneity**: apply the AMN (2020) eq. 7-8 / AF Sec. 3.5
+ investment control-function correction in Stage 3. Defaults to `False`. When
+ `True` and the model has an endogenous (investment) factor, a first-stage
+ investment equation is OLS-fit per period and its residual is added as an
+ additive `cf` covariate (coefficient `kappa_t`) to each state factor's
+ production regression; observed factors are then excluded from the production
+ function and act as instruments (at least one observed instrument is
+ required). The default stays `False` because `estimate_af` calls
+ `estimate_amn` for start values and the AF likelihood implements only
+ `kappa=0`; opt into the correction at the application call site. A no-op for
+ models without endogenous factors.
+
+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.
diff --git a/docs/explanations/notes_on_factor_scales.md b/docs/explanations/notes_on_factor_scales.md
index e938faf4..8bec8dcb 100644
--- a/docs/explanations/notes_on_factor_scales.md
+++ b/docs/explanations/notes_on_factor_scales.md
@@ -8,10 +8,16 @@ Wiswall and Agostinelli define a class of transition functions with Known Locati
Scale (KLS) that require fewer normalizations. Their critique potentially invalidates
certain empirical estimates from CHS, but not the general estimation approach.
-To get estimates that avoid renormalization issues, you can either:
+To reduce the risk of renormalization issues, you can either:
1. Use fewer normalizations with KLS transition functions, or
2. Use non-KLS transition functions with one normalization per period and factor
+skillmodels supports both schemes, but supporting a scheme is not the same as
+guaranteeing identification: the model checker performs syntactic checks on the
+normalizations, not a transition-specific rank or scale-invariance analysis.
+Whether a given scheme identifies the model still depends on the transition
+functions and data; see the cautions below.
+
As there is no natural scale for skills, neither approach is inherently better.
However, we prefer using flexible non-KLS transition functions with explicit
normalizations because:
diff --git a/docs/getting_started/tutorial.ipynb b/docs/getting_started/tutorial.ipynb
index b123b437..269397c2 100644
--- a/docs/getting_started/tutorial.ipynb
+++ b/docs/getting_started/tutorial.ipynb
@@ -5,63 +5,183 @@
"id": "0",
"metadata": {},
"source": [
- "# Skillmodels Quickstart\n",
+ "# skillmodels three-estimator tutorial\n",
"\n",
- "This tutorial demonstrates the basic workflow for estimating a latent factor model\n",
- "using skillmodels. We'll use Example 2 from the Cunha, Heckman, and Schennach (2010)\n",
- "replication files."
+ "skillmodels ships three estimators that share the same `ModelSpec` and the same parameter index:\n",
+ "\n",
+ "- **CHS** (Cunha-Heckman-Schennach 2010): square-root unscented Kalman MLE.\n",
+ "- **AF** (Antweiler-Freyberger 2025): sequential MLE with Halton quadrature over the latent posterior, period by period.\n",
+ "- **AMN** (Attanasio-Meghir-Nix 2020): three-stage estimator (EM on a mixture-of-normals → minimum distance → simulate-and-regress).\n",
+ "\n",
+ "This tutorial walks all three on the **CNLSY** dataset from the AF 2025 application: three waves (ages 7 / 9 / 11), a CES production function for the latent skill, and an endogenous investment factor.\n",
+ "\n",
+ "The notebook is committed with pre-rendered outputs (`execute: false`); rerun it to refresh the numbers if the API drifts. Single-machine wallclock at the time of writing: ~10–15 minutes total."
]
},
{
- "cell_type": "code",
- "execution_count": null,
+ "cell_type": "markdown",
"id": "1",
"metadata": {},
- "outputs": [],
"source": [
- "import pandas as pd\n",
+ "## Setup\n",
"\n",
- "from skillmodels import get_maximization_inputs\n",
- "from skillmodels.config import TEST_DATA_DIR\n",
- "from skillmodels.test_data.model2 import MODEL2"
+ "Load the long-format CNLSY measurements bundled with the package and build the CES `ModelSpec`. The model has four latent factors (`skills`, `MC`, `MN`, `investment`) plus one observed factor (`log_income`). Cognitive (`MC`) and non-cognitive (`MN`) skills are pinned to identity transitions: they enter the skills CES as additional inputs but evolve as constants over time. All factor weights inside the skills CES are estimated freely, subject to the simplex constraint (gammas non-negative, sum to one)."
]
},
{
- "cell_type": "markdown",
+ "cell_type": "code",
+ "execution_count": null,
"id": "2",
- "metadata": {},
- "source": [
- ""
- ]
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-05-14T07:13:11.517683Z",
+ "iopub.status.busy": "2026-05-14T07:13:11.517454Z",
+ "iopub.status.idle": "2026-05-14T07:13:13.015271Z",
+ "shell.execute_reply": "2026-05-14T07:13:13.014780Z"
+ }
+ },
+ "outputs": [],
+ "source": "import warnings\n\nimport optimagic as om\nimport pandas as pd\nimport plotly.graph_objects as go\n\nfrom skillmodels import FactorSpec, ModelSpec, Normalizations\nfrom skillmodels.af import AFEstimationOptions, estimate_af\nfrom skillmodels.af.posterior_states import get_af_posterior_states\nfrom skillmodels.amn import AMNEstimationOptions, estimate_amn\nfrom skillmodels.amn.posterior_states import get_amn_posterior_states\nfrom skillmodels.chs import (\n CHSEstimationOptions,\n get_maximization_inputs,\n)\nfrom skillmodels.common.config import CNLSY_DATA_PATH\nfrom skillmodels.common.individual_states import get_individual_states_from_params\nfrom skillmodels.common.variance_decomposition import (\n decompose_measurement_variance,\n summarize_measurement_reliability,\n)\n\nwarnings.filterwarnings(\"ignore\", category=DeprecationWarning)\npd.options.display.float_format = \"{:.3f}\".format"
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 2,
"id": "3",
- "metadata": {},
- "outputs": [],
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-05-14T07:13:13.016446Z",
+ "iopub.status.busy": "2026-05-14T07:13:13.016229Z",
+ "iopub.status.idle": "2026-05-14T07:13:13.027837Z",
+ "shell.execute_reply": "2026-05-14T07:13:13.027469Z"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "individuals: 1403\n",
+ "periods: [0, 1, 2]\n",
+ "columns: ['skill_math', 'skill_recog', 'skill_comp', 'mc_1', 'mc_2', 'mc_3', 'mc_4', 'mc_5', 'mc_6', 'mn_neg', 'mn_pos', 'mn_rotter', 'inv_reads', 'inv_museum', 'inv_praised', 'log_income_observed']\n"
+ ]
+ }
+ ],
"source": [
- "model = MODEL2\n",
- "\n",
- "# Show the structure\n",
- "print(\"Factors:\", list(model.factors.keys()))"
+ "data = pd.read_csv(CNLSY_DATA_PATH).set_index([\"caseid\", \"period\"])\n",
+ "print(f\"individuals: {data.index.get_level_values(0).nunique()}\")\n",
+ "print(f\"periods: {sorted(data.index.get_level_values(1).unique())}\")\n",
+ "print(f\"columns: {list(data.columns)}\")"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 3,
"id": "4",
- "metadata": {},
- "outputs": [],
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-05-14T07:13:13.029013Z",
+ "iopub.status.busy": "2026-05-14T07:13:13.028919Z",
+ "iopub.status.idle": "2026-05-14T07:13:13.036687Z",
+ "shell.execute_reply": "2026-05-14T07:13:13.036255Z"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "latent factors: ('skills', 'MC', 'MN', 'investment')\n",
+ "observed factors: ('log_income_observed',)\n",
+ "n_mixtures: 2\n"
+ ]
+ }
+ ],
"source": [
- "data = pd.read_stata(TEST_DATA_DIR / \"model2_simulated_data.dta\")\n",
- "data = data.set_index([\"caseid\", \"period\"])\n",
- "data.head()"
+ "_N_PERIODS = 3\n",
+ "MC_MEASURES = tuple(f\"mc_{i + 1}\" for i in range(6))\n",
+ "MN_MEASURES = (\"mn_neg\", \"mn_pos\", \"mn_rotter\")\n",
+ "SKILL_MEASURES = (\"skill_math\", \"skill_recog\", \"skill_comp\")\n",
+ "INV_MEASURES = (\"inv_reads\", \"inv_museum\", \"inv_praised\")\n",
+ "INCOME_MEASURE = \"log_income_observed\"\n",
+ "\n",
+ "\n",
+ "def _measurements(meas, active_periods=None):\n",
+ " active = tuple(range(_N_PERIODS)) if active_periods is None else active_periods\n",
+ " return tuple(meas if t in active else () for t in range(_N_PERIODS))\n",
+ "\n",
+ "\n",
+ "def _normalizations(\n",
+ " meas, *, normalize_periods=None, active_periods=None, pin_first_intercept=True\n",
+ "):\n",
+ " active = tuple(range(_N_PERIODS)) if active_periods is None else active_periods\n",
+ " norm = tuple(range(_N_PERIODS)) if normalize_periods is None else normalize_periods\n",
+ " loadings, intercepts = [], []\n",
+ " for t in range(_N_PERIODS):\n",
+ " if t in active and t in norm:\n",
+ " loadings.append({meas[0]: 1.0})\n",
+ " intercepts.append({meas[0]: 0.0} if pin_first_intercept else {})\n",
+ " else:\n",
+ " loadings.append({})\n",
+ " intercepts.append({})\n",
+ " return Normalizations(loadings=tuple(loadings), intercepts=tuple(intercepts))\n",
+ "\n",
+ "\n",
+ "factors = {\n",
+ " \"skills\": FactorSpec(\n",
+ " measurements=_measurements(SKILL_MEASURES),\n",
+ " normalizations=_normalizations(SKILL_MEASURES, normalize_periods=(0,)),\n",
+ " transition_function=\"log_ces\",\n",
+ " ),\n",
+ " \"MC\": FactorSpec(\n",
+ " measurements=_measurements(MC_MEASURES, active_periods=(0,)),\n",
+ " normalizations=_normalizations(\n",
+ " MC_MEASURES, active_periods=(0,), normalize_periods=(0,)\n",
+ " ),\n",
+ " transition_function=\"linear\",\n",
+ " has_production_shock=False,\n",
+ " ),\n",
+ " \"MN\": FactorSpec(\n",
+ " measurements=_measurements(MN_MEASURES, active_periods=(0,)),\n",
+ " normalizations=_normalizations(\n",
+ " MN_MEASURES, active_periods=(0,), normalize_periods=(0,)\n",
+ " ),\n",
+ " transition_function=\"linear\",\n",
+ " has_production_shock=False,\n",
+ " ),\n",
+ " \"investment\": FactorSpec(\n",
+ " measurements=_measurements(INV_MEASURES, active_periods=(0, 1)),\n",
+ " normalizations=_normalizations(\n",
+ " INV_MEASURES, active_periods=(0, 1), normalize_periods=(0,)\n",
+ " ),\n",
+ " transition_function=\"linear\",\n",
+ " ),\n",
+ "}\n",
+ "\n",
+ "# Pin the time-invariant MC and MN transitions to identity: self-coefficient 1,\n",
+ "# cross-factor 0, constant 0. The skills CES gammas stay free so the optimizer\n",
+ "# discovers which factors actually feed skill growth.\n",
+ "fixed_rows = []\n",
+ "for t in range(_N_PERIODS - 1):\n",
+ " for factor in (\"MC\", \"MN\"):\n",
+ " fixed_rows.append(((\"transition\", t, factor, factor), 1.0))\n",
+ " for other in (\"skills\", \"MC\", \"MN\", \"investment\"):\n",
+ " if other != factor:\n",
+ " fixed_rows.append(((\"transition\", t, factor, other), 0.0))\n",
+ " fixed_rows.append(((\"transition\", t, factor, \"constant\"), 0.0))\n",
+ "\n",
+ "fixed_idx = pd.MultiIndex.from_tuples(\n",
+ " [r[0] for r in fixed_rows], names=[\"category\", \"period\", \"name1\", \"name2\"]\n",
+ ")\n",
+ "fixed_params = pd.DataFrame({\"value\": [r[1] for r in fixed_rows]}, index=fixed_idx)\n",
+ "\n",
+ "model = ModelSpec(\n",
+ " factors=factors,\n",
+ " observed_factors=(INCOME_MEASURE,),\n",
+ " n_mixtures=2,\n",
+ ")\n",
+ "print(f\"latent factors: {tuple(model.factors)}\")\n",
+ "print(f\"observed factors: {model.observed_factors}\")\n",
+ "print(f\"n_mixtures: {model.n_mixtures}\")"
]
},
{
@@ -69,202 +189,2140 @@
"id": "5",
"metadata": {},
"source": [
- "## Getting Maximization Inputs\n",
- "\n",
- "The main entry point is `get_maximization_inputs()`. It takes a model specification\n",
- "and dataset, and returns everything needed to maximize the likelihood using optimagic:\n",
+ "## CHS: square-root unscented Kalman MLE\n",
"\n",
- "- `loglike`: The compiled log-likelihood function\n",
- "- `gradient`: The gradient of the log-likelihood\n",
- "- `loglike_and_gradient`: Combined function (more efficient)\n",
- "- `debug_loglike`: Uncompiled version for debugging\n",
- "- `params_template`: Parameter DataFrame with bounds and starting values\n",
- "- `constraints`: Parameter constraints for optimization"
+ "CHS estimation runs in two steps: `get_maximization_inputs(...)` compiles the jitted likelihood, gradients, and constraints from the spec + data; then we hand the bundle to `optimagic.maximize`. Defaults are CHS's `start_params_strategy=\"amn\"` (runs AMN first under the hood to seed the param template), `robust_bounds=True`, and `bounds_distance=1e-3`."
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 4,
"id": "6",
- "metadata": {},
- "outputs": [],
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-05-14T07:13:13.038140Z",
+ "iopub.status.busy": "2026-05-14T07:13:13.037984Z",
+ "iopub.status.idle": "2026-05-14T07:13:39.788024Z",
+ "shell.execute_reply": "2026-05-14T07:13:39.787579Z"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "params_template shape: (158, 3)\n",
+ "n_constraints: 32\n"
+ ]
+ },
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:339: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_sd, \"value\"] = float(result.meas_sds[local_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:336: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_load, \"value\"] = float(result.loadings[local_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:339: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_sd, \"value\"] = float(result.meas_sds[local_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:336: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_load, \"value\"] = float(result.loadings[local_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:339: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_sd, \"value\"] = float(result.meas_sds[local_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:336: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_load, \"value\"] = float(result.loadings[local_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:339: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_sd, \"value\"] = float(result.meas_sds[local_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:336: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_load, \"value\"] = float(result.loadings[local_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:336: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_load, \"value\"] = float(result.loadings[local_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:339: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_sd, \"value\"] = float(result.meas_sds[local_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:336: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_load, \"value\"] = float(result.loadings[local_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:339: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_sd, \"value\"] = float(result.meas_sds[local_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:336: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_load, \"value\"] = float(result.loadings[local_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:339: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_sd, \"value\"] = float(result.meas_sds[local_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:363: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc, \"value\"] = sd_factor\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:559: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc, \"value\"] = float(beta[col_idx])\n",
+ "/home/hmg/econ/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:579: PerformanceWarning: indexing past lexsort depth may impact performance.\n",
+ " params.loc[loc_sd, \"value\"] = shock_sd\n"
+ ]
+ }
+ ],
"source": [
- "max_inputs = get_maximization_inputs(model, data)\n",
- "print(\"Available keys:\", list(max_inputs.keys()))"
+ "chs_options = CHSEstimationOptions(\n",
+ " robust_bounds=True,\n",
+ " bounds_distance=1e-3,\n",
+ ")\n",
+ "max_inputs = get_maximization_inputs(\n",
+ " model_spec=model,\n",
+ " data=data,\n",
+ " chs_options=chs_options,\n",
+ " fixed_params=fixed_params,\n",
+ ")\n",
+ "print(f\"params_template shape: {max_inputs['params_template'].shape}\")\n",
+ "print(f\"n_constraints: {len(max_inputs['constraints'])}\")"
]
},
{
- "cell_type": "markdown",
+ "cell_type": "code",
+ "execution_count": 5,
"id": "7",
- "metadata": {},
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-05-14T07:13:39.789235Z",
+ "iopub.status.busy": "2026-05-14T07:13:39.789034Z",
+ "iopub.status.idle": "2026-05-14T07:19:28.096303Z",
+ "shell.execute_reply": "2026-05-14T07:19:28.095804Z"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "CHS success: True, loglike: -39620.92\n"
+ ]
+ }
+ ],
"source": [
- "## Parameter Template\n",
- "\n",
- "The `params_template` is a pandas DataFrame with:\n",
- "- A MultiIndex identifying each parameter (category, period, name1, name2)\n",
- "- Columns for `value` (to be filled with starting values), `lower_bound`, `upper_bound`"
+ "chs_result = om.maximize(\n",
+ " fun=max_inputs[\"loglike\"],\n",
+ " params=max_inputs[\"params_template\"],\n",
+ " algorithm=\"scipy_lbfgsb\",\n",
+ " fun_and_jac=max_inputs[\"loglike_and_gradient\"],\n",
+ " constraints=max_inputs[\"constraints\"],\n",
+ ")\n",
+ "chs_params = chs_result.params\n",
+ "print(f\"CHS success: {chs_result.success}, loglike: {chs_result.fun:.2f}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8",
- "metadata": {},
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-05-14T07:19:28.097610Z",
+ "iopub.status.busy": "2026-05-14T07:19:28.097505Z",
+ "iopub.status.idle": "2026-05-14T07:19:57.433481Z",
+ "shell.execute_reply": "2026-05-14T07:19:57.432980Z"
+ }
+ },
"outputs": [],
- "source": [
- "params_template = max_inputs[\"params_template\"]\n",
- "params_template.head(10)"
- ]
+ "source": "chs_filtered = get_individual_states_from_params(\n model_spec=model, data=data, params=chs_params\n)\nchs_states = chs_filtered[\"unanchored_states\"][\"states\"]\nchs_decomp = decompose_measurement_variance(\n model_spec=model, params=chs_params, filtered_states=chs_states\n)\nchs_reliability = summarize_measurement_reliability(chs_decomp)\nchs_reliability"
},
{
"cell_type": "markdown",
"id": "9",
"metadata": {},
- "source": [
- "## Choosing Starting Values\n",
- "\n",
- "Good starting values are important for optimization. As a rule of thumb:\n",
- "\n",
- "- If measurements are standardized, use 1.0 for free loadings and 0.0 for free intercepts\n",
- "- Start measurement and shock standard deviations slightly larger than expected\n",
- "- Initial state means can often start at 0\n",
- "\n",
- "Here we set reasonable defaults:"
- ]
+ "source": "## AF: sequential Halton-quadrature MLE\n\nAF estimates each period in turn: period 0 fits the joint mixture-of-normals + period-0 measurement system; each subsequent period takes the previous period's posterior and runs a period-specific MLE via `optimagic.minimize` with the algorithm in `AFEstimationOptions.optimizer_algorithm` (default `\"fides\"`; pass `\"scipy_lbfgsb\"` for Monte Carlo sweeps where a deterministic stopping rule matters).\n\nDefault `start_params_strategy=\"amn\"` runs the full AMN three-stage estimator upfront and uses its parameters to seed each AF period. The number of mixture components is read from `ModelSpec.n_mixtures` (set to 2 above), shared across all three estimators."
},
{
"cell_type": "code",
"execution_count": null,
"id": "10",
- "metadata": {},
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-05-14T07:19:57.434856Z",
+ "iopub.status.busy": "2026-05-14T07:19:57.434760Z",
+ "iopub.status.idle": "2026-05-14T07:20:56.673081Z",
+ "shell.execute_reply": "2026-05-14T07:20:56.672613Z"
+ }
+ },
"outputs": [],
- "source": [
- "params = params_template.copy()\n",
- "\n",
- "# Set starting values by category\n",
- "for category in params.index.get_level_values(\"category\").unique():\n",
- " if category == \"loadings\":\n",
- " params.loc[category, \"value\"] = 1.0\n",
- " elif category == \"controls\":\n",
- " params.loc[category, \"value\"] = 0.0\n",
- " elif category in (\"meas_sds\", \"shock_sds\") or category == \"initial_cholcovs\":\n",
- " params.loc[category, \"value\"] = 0.5\n",
- " elif category == \"initial_states\":\n",
- " params.loc[category, \"value\"] = 0.0\n",
- " elif category == \"mixture_weights\":\n",
- " params.loc[category, \"value\"] = 1.0\n",
- " elif category == \"transition\":\n",
- " # Set transition parameters to reasonable defaults\n",
- " params.loc[category, \"value\"] = 0.5\n",
- "\n",
- "params.head(10)"
- ]
+ "source": "af_options = AFEstimationOptions(\n n_halton_points=100,\n n_halton_points_shock=50,\n optimizer_algorithm=\"scipy_lbfgsb\",\n)\naf_result = estimate_af(\n model_spec=model,\n data=data,\n options=af_options,\n fixed_params=fixed_params,\n)\naf_lls = [pr.loglikelihood for pr in af_result.period_results]\nprint(\n f\"AF per-period success: {[bool(pr.success) for pr in af_result.period_results]}, \"\n f\"per-period log-likelihoods: {[f'{ll:.2f}' for ll in af_lls]}\"\n)"
},
{
- "cell_type": "markdown",
+ "cell_type": "code",
+ "execution_count": null,
"id": "11",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-05-14T07:20:56.674410Z",
+ "iopub.status.busy": "2026-05-14T07:20:56.674311Z",
+ "iopub.status.idle": "2026-05-14T07:21:01.658274Z",
+ "shell.execute_reply": "2026-05-14T07:21:01.657922Z"
+ }
+ },
+ "outputs": [],
+ "source": "af_posterior = get_af_posterior_states(\n af_result=af_result,\n model_spec=model,\n data=data,\n n_halton_points=200,\n)\naf_states = af_posterior[\"unanchored_states\"][\"states\"]\naf_decomp = decompose_measurement_variance(\n model_spec=model, params=af_result.params, filtered_states=af_states\n)\naf_reliability = summarize_measurement_reliability(af_decomp)\naf_reliability"
+ },
+ {
+ "cell_type": "markdown",
+ "id": "12",
"metadata": {},
- "source": [
- "## JAX Compilation\n",
- "\n",
- "Skillmodels uses JAX for just-in-time compilation and automatic differentiation.\n",
- "The first call to `loglike` or `gradient` triggers compilation, which takes a few\n",
- "seconds. Subsequent calls are very fast."
- ]
+ "source": "## AMN: three-stage mixture-of-normals\n\nAMN's three stages are EM on the augmented measurement vector (mixture-of-normals), minimum-distance recovery of the structural parameters, and a simulate-and-regress step on the fitted mixture for the production function. With `ModelSpec.n_mixtures=2` the Stage-1 EM fits a 2-component Gaussian mixture on the joint $(M, X)$ vector — the kind of non-Gaussian latent structure AMN is designed for."
},
{
"cell_type": "code",
"execution_count": null,
- "id": "12",
- "metadata": {},
+ "id": "13",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-05-14T07:21:01.659567Z",
+ "iopub.status.busy": "2026-05-14T07:21:01.659466Z",
+ "iopub.status.idle": "2026-05-14T07:21:26.698296Z",
+ "shell.execute_reply": "2026-05-14T07:21:26.695286Z"
+ }
+ },
"outputs": [],
- "source": [
- "loglike = max_inputs[\"loglike\"]\n",
- "gradient = max_inputs[\"gradient\"]\n",
- "loglike_and_gradient = max_inputs[\"loglike_and_gradient\"]"
- ]
+ "source": "amn_options = AMNEstimationOptions(\n em_max_iter=500,\n n_simulation_draws=50_000,\n seed=0,\n)\namn_result = estimate_amn(\n model_spec=model,\n data=data,\n options=amn_options,\n fixed_params=fixed_params,\n)\nprint(f\"AMN success: {amn_result.success}\")"
},
{
"cell_type": "code",
"execution_count": null,
- "id": "13",
- "metadata": {},
+ "id": "14",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-05-14T07:21:26.699696Z",
+ "iopub.status.busy": "2026-05-14T07:21:26.699564Z",
+ "iopub.status.idle": "2026-05-14T07:21:26.794622Z",
+ "shell.execute_reply": "2026-05-14T07:21:26.793933Z"
+ }
+ },
"outputs": [],
- "source": [
- "# First call includes compilation time\n",
- "loglike_value = loglike(params)\n",
- "print(f\"Log-likelihood at starting values: {loglike_value:.2f}\")"
- ]
+ "source": "amn_posterior = get_amn_posterior_states(amn_result=amn_result, data=data)\namn_states = amn_posterior[\"unanchored_states\"][\"states\"]\namn_decomp = decompose_measurement_variance(\n model_spec=model, params=amn_result.params, filtered_states=amn_states\n)\namn_reliability = summarize_measurement_reliability(amn_decomp)\namn_reliability"
},
{
"cell_type": "markdown",
- "id": "14",
+ "id": "15",
"metadata": {},
"source": [
- "## Constraints\n",
+ "## Cross-estimator comparison\n",
"\n",
- "Skillmodels automatically generates constraints from the model specification:\n",
- "- Fixed parameters (normalized loadings and intercepts)\n",
- "- Stagemap equality constraints\n",
- "- Bound constraints\n",
+ "The three estimators target the same likelihood under different approximations / objectives. Below we line up:\n",
"\n",
- "You can add additional constraints for your specific model."
+ "1. Period-0 measurement loadings for the `skill_*` indicators.\n",
+ "2. CES production-function gammas + φ for the period-0 → period-1 skills transition.\n",
+ "3. Signal-share scatter: do they agree on which measurements are high-noise?"
]
},
{
"cell_type": "code",
"execution_count": null,
- "id": "15",
- "metadata": {},
+ "id": "16",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-05-14T07:21:26.795949Z",
+ "iopub.status.busy": "2026-05-14T07:21:26.795823Z",
+ "iopub.status.idle": "2026-05-14T07:21:26.809737Z",
+ "shell.execute_reply": "2026-05-14T07:21:26.809254Z"
+ }
+ },
"outputs": [],
+ "source": "estimators = {\n \"CHS\": chs_params,\n \"AF\": af_result.params,\n \"AMN\": amn_result.params,\n}\n\n\ndef _free_loading(params, period, meas):\n loc = (\"loadings\", period, meas, \"skills\")\n if loc not in params.index:\n return None\n return float(params.loc[loc, \"value\"])\n\n\nloading_rows = []\nfor est, params in estimators.items():\n for meas in SKILL_MEASURES:\n value = _free_loading(params, period=0, meas=meas)\n if value is None:\n continue\n loading_rows.append({\"estimator\": est, \"measurement\": meas, \"loading\": value})\nloadings_df = pd.DataFrame(loading_rows)\nloadings_pivot = loadings_df.pivot(\n index=\"measurement\", columns=\"estimator\", values=\"loading\"\n)\nloadings_pivot"
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "id": "17",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-05-14T07:21:26.810916Z",
+ "iopub.status.busy": "2026-05-14T07:21:26.810787Z",
+ "iopub.status.idle": "2026-05-14T07:21:29.003649Z",
+ "shell.execute_reply": "2026-05-14T07:21:29.003242Z"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.plotly.v1+json": {
+ "config": {
+ "plotlyServerURL": "https://plot.ly"
+ },
+ "data": [
+ {
+ "name": "CHS",
+ "type": "bar",
+ "x": [
+ "skill_math",
+ "skill_recog",
+ "skill_comp"
+ ],
+ "y": {
+ "bdata": "AAAAAAAA8D/DIYazpav2Pyn68owbpfU/",
+ "dtype": "f8"
+ }
+ },
+ {
+ "name": "AF",
+ "type": "bar",
+ "x": [
+ "skill_math",
+ "skill_recog",
+ "skill_comp"
+ ],
+ "y": {
+ "bdata": "AAAAAAAA8D9fbAPZXff1P7wQMLJ+EvY/",
+ "dtype": "f8"
+ }
+ },
+ {
+ "name": "AMN",
+ "type": "bar",
+ "x": [
+ "skill_math",
+ "skill_recog",
+ "skill_comp"
+ ],
+ "y": {
+ "bdata": "AAAAAAAA8D+wS/CPTL/yP0IjsiJDuvI/",
+ "dtype": "f8"
+ }
+ }
+ ],
+ "layout": {
+ "barmode": "group",
+ "template": {
+ "data": {
+ "bar": [
+ {
+ "error_x": {
+ "color": "#2a3f5f"
+ },
+ "error_y": {
+ "color": "#2a3f5f"
+ },
+ "marker": {
+ "line": {
+ "color": "white",
+ "width": 0.5
+ },
+ "pattern": {
+ "fillmode": "overlay",
+ "size": 10,
+ "solidity": 0.2
+ }
+ },
+ "type": "bar"
+ }
+ ],
+ "barpolar": [
+ {
+ "marker": {
+ "line": {
+ "color": "white",
+ "width": 0.5
+ },
+ "pattern": {
+ "fillmode": "overlay",
+ "size": 10,
+ "solidity": 0.2
+ }
+ },
+ "type": "barpolar"
+ }
+ ],
+ "carpet": [
+ {
+ "aaxis": {
+ "endlinecolor": "#2a3f5f",
+ "gridcolor": "#C8D4E3",
+ "linecolor": "#C8D4E3",
+ "minorgridcolor": "#C8D4E3",
+ "startlinecolor": "#2a3f5f"
+ },
+ "baxis": {
+ "endlinecolor": "#2a3f5f",
+ "gridcolor": "#C8D4E3",
+ "linecolor": "#C8D4E3",
+ "minorgridcolor": "#C8D4E3",
+ "startlinecolor": "#2a3f5f"
+ },
+ "type": "carpet"
+ }
+ ],
+ "choropleth": [
+ {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ },
+ "type": "choropleth"
+ }
+ ],
+ "contour": [
+ {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ },
+ "colorscale": [
+ [
+ 0,
+ "#0d0887"
+ ],
+ [
+ 0.1111111111111111,
+ "#46039f"
+ ],
+ [
+ 0.2222222222222222,
+ "#7201a8"
+ ],
+ [
+ 0.3333333333333333,
+ "#9c179e"
+ ],
+ [
+ 0.4444444444444444,
+ "#bd3786"
+ ],
+ [
+ 0.5555555555555556,
+ "#d8576b"
+ ],
+ [
+ 0.6666666666666666,
+ "#ed7953"
+ ],
+ [
+ 0.7777777777777778,
+ "#fb9f3a"
+ ],
+ [
+ 0.8888888888888888,
+ "#fdca26"
+ ],
+ [
+ 1,
+ "#f0f921"
+ ]
+ ],
+ "type": "contour"
+ }
+ ],
+ "contourcarpet": [
+ {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ },
+ "type": "contourcarpet"
+ }
+ ],
+ "heatmap": [
+ {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ },
+ "colorscale": [
+ [
+ 0,
+ "#0d0887"
+ ],
+ [
+ 0.1111111111111111,
+ "#46039f"
+ ],
+ [
+ 0.2222222222222222,
+ "#7201a8"
+ ],
+ [
+ 0.3333333333333333,
+ "#9c179e"
+ ],
+ [
+ 0.4444444444444444,
+ "#bd3786"
+ ],
+ [
+ 0.5555555555555556,
+ "#d8576b"
+ ],
+ [
+ 0.6666666666666666,
+ "#ed7953"
+ ],
+ [
+ 0.7777777777777778,
+ "#fb9f3a"
+ ],
+ [
+ 0.8888888888888888,
+ "#fdca26"
+ ],
+ [
+ 1,
+ "#f0f921"
+ ]
+ ],
+ "type": "heatmap"
+ }
+ ],
+ "histogram": [
+ {
+ "marker": {
+ "pattern": {
+ "fillmode": "overlay",
+ "size": 10,
+ "solidity": 0.2
+ }
+ },
+ "type": "histogram"
+ }
+ ],
+ "histogram2d": [
+ {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ },
+ "colorscale": [
+ [
+ 0,
+ "#0d0887"
+ ],
+ [
+ 0.1111111111111111,
+ "#46039f"
+ ],
+ [
+ 0.2222222222222222,
+ "#7201a8"
+ ],
+ [
+ 0.3333333333333333,
+ "#9c179e"
+ ],
+ [
+ 0.4444444444444444,
+ "#bd3786"
+ ],
+ [
+ 0.5555555555555556,
+ "#d8576b"
+ ],
+ [
+ 0.6666666666666666,
+ "#ed7953"
+ ],
+ [
+ 0.7777777777777778,
+ "#fb9f3a"
+ ],
+ [
+ 0.8888888888888888,
+ "#fdca26"
+ ],
+ [
+ 1,
+ "#f0f921"
+ ]
+ ],
+ "type": "histogram2d"
+ }
+ ],
+ "histogram2dcontour": [
+ {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ },
+ "colorscale": [
+ [
+ 0,
+ "#0d0887"
+ ],
+ [
+ 0.1111111111111111,
+ "#46039f"
+ ],
+ [
+ 0.2222222222222222,
+ "#7201a8"
+ ],
+ [
+ 0.3333333333333333,
+ "#9c179e"
+ ],
+ [
+ 0.4444444444444444,
+ "#bd3786"
+ ],
+ [
+ 0.5555555555555556,
+ "#d8576b"
+ ],
+ [
+ 0.6666666666666666,
+ "#ed7953"
+ ],
+ [
+ 0.7777777777777778,
+ "#fb9f3a"
+ ],
+ [
+ 0.8888888888888888,
+ "#fdca26"
+ ],
+ [
+ 1,
+ "#f0f921"
+ ]
+ ],
+ "type": "histogram2dcontour"
+ }
+ ],
+ "mesh3d": [
+ {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ },
+ "type": "mesh3d"
+ }
+ ],
+ "parcoords": [
+ {
+ "line": {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ }
+ },
+ "type": "parcoords"
+ }
+ ],
+ "pie": [
+ {
+ "automargin": true,
+ "type": "pie"
+ }
+ ],
+ "scatter": [
+ {
+ "fillpattern": {
+ "fillmode": "overlay",
+ "size": 10,
+ "solidity": 0.2
+ },
+ "type": "scatter"
+ }
+ ],
+ "scatter3d": [
+ {
+ "line": {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ }
+ },
+ "marker": {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ }
+ },
+ "type": "scatter3d"
+ }
+ ],
+ "scattercarpet": [
+ {
+ "marker": {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ }
+ },
+ "type": "scattercarpet"
+ }
+ ],
+ "scattergeo": [
+ {
+ "marker": {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ }
+ },
+ "type": "scattergeo"
+ }
+ ],
+ "scattergl": [
+ {
+ "marker": {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ }
+ },
+ "type": "scattergl"
+ }
+ ],
+ "scattermap": [
+ {
+ "marker": {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ }
+ },
+ "type": "scattermap"
+ }
+ ],
+ "scattermapbox": [
+ {
+ "marker": {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ }
+ },
+ "type": "scattermapbox"
+ }
+ ],
+ "scatterpolar": [
+ {
+ "marker": {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ }
+ },
+ "type": "scatterpolar"
+ }
+ ],
+ "scatterpolargl": [
+ {
+ "marker": {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ }
+ },
+ "type": "scatterpolargl"
+ }
+ ],
+ "scatterternary": [
+ {
+ "marker": {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ }
+ },
+ "type": "scatterternary"
+ }
+ ],
+ "surface": [
+ {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ },
+ "colorscale": [
+ [
+ 0,
+ "#0d0887"
+ ],
+ [
+ 0.1111111111111111,
+ "#46039f"
+ ],
+ [
+ 0.2222222222222222,
+ "#7201a8"
+ ],
+ [
+ 0.3333333333333333,
+ "#9c179e"
+ ],
+ [
+ 0.4444444444444444,
+ "#bd3786"
+ ],
+ [
+ 0.5555555555555556,
+ "#d8576b"
+ ],
+ [
+ 0.6666666666666666,
+ "#ed7953"
+ ],
+ [
+ 0.7777777777777778,
+ "#fb9f3a"
+ ],
+ [
+ 0.8888888888888888,
+ "#fdca26"
+ ],
+ [
+ 1,
+ "#f0f921"
+ ]
+ ],
+ "type": "surface"
+ }
+ ],
+ "table": [
+ {
+ "cells": {
+ "fill": {
+ "color": "#EBF0F8"
+ },
+ "line": {
+ "color": "white"
+ }
+ },
+ "header": {
+ "fill": {
+ "color": "#C8D4E3"
+ },
+ "line": {
+ "color": "white"
+ }
+ },
+ "type": "table"
+ }
+ ]
+ },
+ "layout": {
+ "annotationdefaults": {
+ "arrowcolor": "#2a3f5f",
+ "arrowhead": 0,
+ "arrowwidth": 1
+ },
+ "autotypenumbers": "strict",
+ "coloraxis": {
+ "colorbar": {
+ "outlinewidth": 0,
+ "ticks": ""
+ }
+ },
+ "colorscale": {
+ "diverging": [
+ [
+ 0,
+ "#8e0152"
+ ],
+ [
+ 0.1,
+ "#c51b7d"
+ ],
+ [
+ 0.2,
+ "#de77ae"
+ ],
+ [
+ 0.3,
+ "#f1b6da"
+ ],
+ [
+ 0.4,
+ "#fde0ef"
+ ],
+ [
+ 0.5,
+ "#f7f7f7"
+ ],
+ [
+ 0.6,
+ "#e6f5d0"
+ ],
+ [
+ 0.7,
+ "#b8e186"
+ ],
+ [
+ 0.8,
+ "#7fbc41"
+ ],
+ [
+ 0.9,
+ "#4d9221"
+ ],
+ [
+ 1,
+ "#276419"
+ ]
+ ],
+ "sequential": [
+ [
+ 0,
+ "#0d0887"
+ ],
+ [
+ 0.1111111111111111,
+ "#46039f"
+ ],
+ [
+ 0.2222222222222222,
+ "#7201a8"
+ ],
+ [
+ 0.3333333333333333,
+ "#9c179e"
+ ],
+ [
+ 0.4444444444444444,
+ "#bd3786"
+ ],
+ [
+ 0.5555555555555556,
+ "#d8576b"
+ ],
+ [
+ 0.6666666666666666,
+ "#ed7953"
+ ],
+ [
+ 0.7777777777777778,
+ "#fb9f3a"
+ ],
+ [
+ 0.8888888888888888,
+ "#fdca26"
+ ],
+ [
+ 1,
+ "#f0f921"
+ ]
+ ],
+ "sequentialminus": [
+ [
+ 0,
+ "#0d0887"
+ ],
+ [
+ 0.1111111111111111,
+ "#46039f"
+ ],
+ [
+ 0.2222222222222222,
+ "#7201a8"
+ ],
+ [
+ 0.3333333333333333,
+ "#9c179e"
+ ],
+ [
+ 0.4444444444444444,
+ "#bd3786"
+ ],
+ [
+ 0.5555555555555556,
+ "#d8576b"
+ ],
+ [
+ 0.6666666666666666,
+ "#ed7953"
+ ],
+ [
+ 0.7777777777777778,
+ "#fb9f3a"
+ ],
+ [
+ 0.8888888888888888,
+ "#fdca26"
+ ],
+ [
+ 1,
+ "#f0f921"
+ ]
+ ]
+ },
+ "colorway": [
+ "#636efa",
+ "#EF553B",
+ "#00cc96",
+ "#ab63fa",
+ "#FFA15A",
+ "#19d3f3",
+ "#FF6692",
+ "#B6E880",
+ "#FF97FF",
+ "#FECB52"
+ ],
+ "font": {
+ "color": "#2a3f5f"
+ },
+ "geo": {
+ "bgcolor": "white",
+ "lakecolor": "white",
+ "landcolor": "white",
+ "showlakes": true,
+ "showland": true,
+ "subunitcolor": "#C8D4E3"
+ },
+ "hoverlabel": {
+ "align": "left"
+ },
+ "hovermode": "closest",
+ "mapbox": {
+ "style": "light"
+ },
+ "paper_bgcolor": "white",
+ "plot_bgcolor": "white",
+ "polar": {
+ "angularaxis": {
+ "gridcolor": "#EBF0F8",
+ "linecolor": "#EBF0F8",
+ "ticks": ""
+ },
+ "bgcolor": "white",
+ "radialaxis": {
+ "gridcolor": "#EBF0F8",
+ "linecolor": "#EBF0F8",
+ "ticks": ""
+ }
+ },
+ "scene": {
+ "xaxis": {
+ "backgroundcolor": "white",
+ "gridcolor": "#DFE8F3",
+ "gridwidth": 2,
+ "linecolor": "#EBF0F8",
+ "showbackground": true,
+ "ticks": "",
+ "zerolinecolor": "#EBF0F8"
+ },
+ "yaxis": {
+ "backgroundcolor": "white",
+ "gridcolor": "#DFE8F3",
+ "gridwidth": 2,
+ "linecolor": "#EBF0F8",
+ "showbackground": true,
+ "ticks": "",
+ "zerolinecolor": "#EBF0F8"
+ },
+ "zaxis": {
+ "backgroundcolor": "white",
+ "gridcolor": "#DFE8F3",
+ "gridwidth": 2,
+ "linecolor": "#EBF0F8",
+ "showbackground": true,
+ "ticks": "",
+ "zerolinecolor": "#EBF0F8"
+ }
+ },
+ "shapedefaults": {
+ "line": {
+ "color": "#2a3f5f"
+ }
+ },
+ "ternary": {
+ "aaxis": {
+ "gridcolor": "#DFE8F3",
+ "linecolor": "#A2B1C6",
+ "ticks": ""
+ },
+ "baxis": {
+ "gridcolor": "#DFE8F3",
+ "linecolor": "#A2B1C6",
+ "ticks": ""
+ },
+ "bgcolor": "white",
+ "caxis": {
+ "gridcolor": "#DFE8F3",
+ "linecolor": "#A2B1C6",
+ "ticks": ""
+ }
+ },
+ "title": {
+ "x": 0.05
+ },
+ "xaxis": {
+ "automargin": true,
+ "gridcolor": "#EBF0F8",
+ "linecolor": "#EBF0F8",
+ "ticks": "",
+ "title": {
+ "standoff": 15
+ },
+ "zerolinecolor": "#EBF0F8",
+ "zerolinewidth": 2
+ },
+ "yaxis": {
+ "automargin": true,
+ "gridcolor": "#EBF0F8",
+ "linecolor": "#EBF0F8",
+ "ticks": "",
+ "title": {
+ "standoff": 15
+ },
+ "zerolinecolor": "#EBF0F8",
+ "zerolinewidth": 2
+ }
+ }
+ },
+ "title": {
+ "text": "Period-0 skill loadings across estimators"
+ },
+ "yaxis": {
+ "title": {
+ "text": "loading"
+ }
+ }
+ }
+ }
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
"source": [
- "constraints = max_inputs[\"constraints\"]\n",
- "print(f\"Number of auto-generated constraints: {len(constraints)}\")"
+ "fig = go.Figure()\n",
+ "for est in (\"CHS\", \"AF\", \"AMN\"):\n",
+ " sub = loadings_df[loadings_df[\"estimator\"] == est]\n",
+ " fig.add_trace(go.Bar(name=est, x=sub[\"measurement\"], y=sub[\"loading\"]))\n",
+ "fig.update_layout(\n",
+ " title=\"Period-0 skill loadings across estimators\",\n",
+ " barmode=\"group\",\n",
+ " template=\"plotly_white\",\n",
+ " yaxis_title=\"loading\",\n",
+ ")\n",
+ "fig"
]
},
{
- "cell_type": "markdown",
- "id": "16",
- "metadata": {},
+ "cell_type": "code",
+ "execution_count": 13,
+ "id": "18",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-05-14T07:21:29.004727Z",
+ "iopub.status.busy": "2026-05-14T07:21:29.004623Z",
+ "iopub.status.idle": "2026-05-14T07:21:29.011824Z",
+ "shell.execute_reply": "2026-05-14T07:21:29.011388Z"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " gamma_skills | \n",
+ " gamma_investment | \n",
+ " phi | \n",
+ "
\n",
+ " \n",
+ " | estimator | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | CHS | \n",
+ " 0.893 | \n",
+ " 0.000 | \n",
+ " 1.337 | \n",
+ "
\n",
+ " \n",
+ " | AF | \n",
+ " 0.683 | \n",
+ " 0.271 | \n",
+ " -1.178 | \n",
+ "
\n",
+ " \n",
+ " | AMN | \n",
+ " 0.894 | \n",
+ " 0.000 | \n",
+ " 0.457 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "