skillmodels three-estimator tutorial

skillmodels ships three estimators that share the same ModelSpec and the same parameter index:

CHS (Cunha-Heckman-Schennach 2010): square-root unscented Kalman MLE.
AF (Antweiler-Freyberger 2025): sequential MLE with Halton quadrature over the latent posterior, period by period.
AMN (Attanasio-Meghir-Nix 2020): three-stage estimator (EM on a mixture-of-normals → minimum distance → simulate-and-regress).

This tutorial walks CHS and AF 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. AMN appears in its role as the CHS seed: it cannot consistently fit the restricted CES standalone, so CHS uses AMN’s three stages to build start values (see the AMN section).

The notebook is committed with pre-rendered outputs (execute: false); rerun it to refresh the numbers if the API drifts. Single-machine wallclock (CPU, JAX CPU backend): roughly 1–2 hours, dominated by the CHS MLE; AF and the AMN seeding step are faster.

Setup¶

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).

import warnings

import optimagic as om
import pandas as pd
import plotly.graph_objects as go

from skillmodels import FactorSpec, ModelSpec, Normalizations
from skillmodels.af import AFEstimationOptions, estimate_af
from skillmodels.af.posterior_states import get_af_posterior_states
from skillmodels.chs import (
    CHSEstimationOptions,
    get_maximization_inputs,
)
from skillmodels.common.config import CNLSY_DATA_PATH
from skillmodels.common.individual_states import get_individual_states_from_params
from skillmodels.common.variance_decomposition import (
    decompose_measurement_variance,
    summarize_measurement_reliability,
)

warnings.filterwarnings("ignore", category=DeprecationWarning)
pd.options.display.float_format = "{:.3f}".format

An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.

data = pd.read_csv(CNLSY_DATA_PATH).set_index(["caseid", "period"])
print(f"individuals: {data.index.get_level_values(0).nunique()}")
print(f"periods: {sorted(data.index.get_level_values(1).unique())}")
print(f"columns: {list(data.columns)}")

individuals: 1403
periods: [0, 1, 2]
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_PERIODS = 3
MC_MEASURES = tuple(f"mc_{i + 1}" for i in range(6))
MN_MEASURES = ("mn_neg", "mn_pos", "mn_rotter")
SKILL_MEASURES = ("skill_math", "skill_recog", "skill_comp")
INV_MEASURES = ("inv_reads", "inv_museum", "inv_praised")
INCOME_MEASURE = "log_income_observed"


def _measurements(meas, active_periods=None):
    active = tuple(range(_N_PERIODS)) if active_periods is None else active_periods
    return tuple(meas if t in active else () for t in range(_N_PERIODS))


def _normalizations(
    meas, *, normalize_periods=None, active_periods=None, pin_first_intercept=True
):
    active = tuple(range(_N_PERIODS)) if active_periods is None else active_periods
    norm = tuple(range(_N_PERIODS)) if normalize_periods is None else normalize_periods
    loadings, intercepts = [], []
    for t in range(_N_PERIODS):
        if t in active and t in norm:
            loadings.append({meas[0]: 1.0})
            intercepts.append({meas[0]: 0.0} if pin_first_intercept else {})
        else:
            loadings.append({})
            intercepts.append({})
    return Normalizations(loadings=tuple(loadings), intercepts=tuple(intercepts))


factors = {
    "skills": FactorSpec(
        measurements=_measurements(SKILL_MEASURES),
        normalizations=_normalizations(SKILL_MEASURES, normalize_periods=(0, 1, 2)),
        transition_function="log_ces",
    ),
    # Cognitive (MC) and non-cognitive (MN) skills are time-invariant inputs to the
    # skills CES. `af_state_role="static_persistent"` tells the AF calendar adapter to
    # re-apply their period-0 measurement density as a static importance block at every
    # step (a no-op for CHS/AMN, which read the same spec).
    "MC": FactorSpec(
        measurements=_measurements(MC_MEASURES, active_periods=(0,)),
        normalizations=_normalizations(
            MC_MEASURES, active_periods=(0,), normalize_periods=(0,)
        ),
        transition_function="linear",
        has_production_shock=False,
        af_state_role="static_persistent",
    ),
    "MN": FactorSpec(
        measurements=_measurements(MN_MEASURES, active_periods=(0,)),
        normalizations=_normalizations(
            MN_MEASURES, active_periods=(0,), normalize_periods=(0,)
        ),
        transition_function="linear",
        has_production_shock=False,
        af_state_role="static_persistent",
    ),
    # Investment is endogenous and reconstructed from its own measurement equation
    # rather than drawn from the period-0 latent mixture, so it sets
    # `is_endogenous=True, has_initial_distribution=False`. This activates the AF
    # source/destination calendar adapter (period-t indicators measure investment in
    # period t).
    "investment": FactorSpec(
        measurements=_measurements(INV_MEASURES, active_periods=(0, 1)),
        normalizations=_normalizations(
            INV_MEASURES, active_periods=(0, 1), normalize_periods=(0, 1)
        ),
        transition_function="linear",
        is_endogenous=True,
        has_initial_distribution=False,
    ),
}

# Pin the time-invariant MC and MN transitions to identity: self-coefficient 1,
# cross-factor 0, constant 0. The skills CES gammas stay free so the optimizer
# discovers which factors actually feed skill growth.
fixed_rows = []
for t in range(_N_PERIODS - 1):
    for factor in ("MC", "MN"):
        fixed_rows.append((("transition", t, factor, factor), 1.0))
        for other in ("skills", "MC", "MN", "investment"):
            if other != factor:
                fixed_rows.append((("transition", t, factor, other), 0.0))
        fixed_rows.append((("transition", t, factor, "constant"), 0.0))

fixed_idx = pd.MultiIndex.from_tuples(
    [r[0] for r in fixed_rows], names=["category", "period", "name1", "name2"]
)
fixed_params = pd.DataFrame({"value": [r[1] for r in fixed_rows]}, index=fixed_idx)

model = ModelSpec(
    factors=factors,
    observed_factors=(INCOME_MEASURE,),
    n_mixtures=2,
)
print(f"latent factors: {tuple(model.factors)}")
print(f"observed factors: {model.observed_factors}")
print(f"n_mixtures: {model.n_mixtures}")

latent factors: ('skills', 'MC', 'MN', 'investment')
observed factors: ('log_income_observed',)
n_mixtures: 2

CHS: square-root unscented Kalman MLE¶

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.

chs_options = CHSEstimationOptions(
    robust_bounds=True,
    bounds_distance=1e-3,
)
max_inputs = get_maximization_inputs(
    model_spec=model,
    data=data,
    chs_options=chs_options,
    fixed_params=fixed_params,
)
print(f"params_template shape: {max_inputs['params_template'].shape}")
print(f"n_constraints: {len(max_inputs['constraints'])}")

params_template shape: (214, 3)
n_constraints: 84

/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:431: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc, "value"] = float(features[rows, f_idx].mean())
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:434: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[w_loc, "value"] = float(rows.mean())
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:552: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_sd, "value"] = float(result.meas_sds[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:549: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_load, "value"] = float(result.loadings[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:552: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_sd, "value"] = float(result.meas_sds[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:549: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_load, "value"] = float(result.loadings[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:552: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_sd, "value"] = float(result.meas_sds[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:549: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_load, "value"] = float(result.loadings[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:552: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_sd, "value"] = float(result.meas_sds[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:549: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_load, "value"] = float(result.loadings[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:552: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_sd, "value"] = float(result.meas_sds[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:549: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_load, "value"] = float(result.loadings[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:552: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_sd, "value"] = float(result.meas_sds[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:549: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_load, "value"] = float(result.loadings[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:552: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_sd, "value"] = float(result.meas_sds[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:549: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_load, "value"] = float(result.loadings[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:576: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc, "value"] = sd_factor

chs_result = om.maximize(
    fun=max_inputs["loglike"],
    params=max_inputs["params_template"],
    algorithm="scipy_lbfgsb",
    fun_and_jac=max_inputs["loglike_and_gradient"],
    constraints=max_inputs["constraints"],
)
chs_params = chs_result.params
print(f"CHS success: {chs_result.success}, loglike: {chs_result.fun:.2f}")

CHS success: True, loglike: -39542.42

chs_filtered = get_individual_states_from_params(
    model_spec=model, data=data, params=chs_params
)
chs_states = chs_filtered["unanchored_states"]["states"]
chs_decomp = decompose_measurement_variance(
    model_spec=model, params=chs_params, filtered_states=chs_states
)
chs_reliability = summarize_measurement_reliability(chs_decomp)
chs_reliability

/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:431: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc, "value"] = float(features[rows, f_idx].mean())
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:434: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[w_loc, "value"] = float(rows.mean())
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:552: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_sd, "value"] = float(result.meas_sds[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:549: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_load, "value"] = float(result.loadings[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:552: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_sd, "value"] = float(result.meas_sds[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:549: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_load, "value"] = float(result.loadings[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:552: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_sd, "value"] = float(result.meas_sds[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:549: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_load, "value"] = float(result.loadings[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:552: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_sd, "value"] = float(result.meas_sds[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:549: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_load, "value"] = float(result.loadings[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:552: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_sd, "value"] = float(result.meas_sds[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:549: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_load, "value"] = float(result.loadings[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:552: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_sd, "value"] = float(result.meas_sds[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:549: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_load, "value"] = float(result.loadings[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:552: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_sd, "value"] = float(result.meas_sds[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:549: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc_load, "value"] = float(result.loadings[local_idx])
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/amn/start_values.py:576: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc, "value"] = sd_factor

AF: sequential Halton-quadrature MLE¶

AF 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).

Default 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.

af_options = AFEstimationOptions(
    n_halton_points=100,
    n_halton_points_shock=50,
    optimizer_algorithm="scipy_lbfgsb",
)
af_result = estimate_af(
    model_spec=model,
    data=data,
    options=af_options,
    fixed_params=fixed_params,
)
af_lls = [pr.loglikelihood for pr in af_result.period_results]
print(
    f"AF per-period success: {[bool(pr.success) for pr in af_result.period_results]}, "
    f"per-period log-likelihoods: {[f'{ll:.2f}' for ll in af_lls]}"
)

/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/estimate.py:345: UserWarning: Restricted-CES scale system ('MC', 'MN', 'skills') period 0: 3 independent scale pins (loadings) across the system, but the CES restrictions identify the relative scales from a SINGLE primitive anchor. 2 of them are testable restrictions, not normalizations.
  validate_af_model(model_spec, fixed_params, constraints)
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/estimate.py:345: UserWarning: Factor 'skills': built-in transition 'log_ces' enumerates parameters over ALL factors including observed factors ('log_income_observed',), so they enter the production function with free coefficients. The AF model assumes observed factors affect skills only through the investment equation. Use a production-factors-only transition ('translog_af' or 'log_ces_af'), or pin every observed-factor transition coefficient to 0.0 via `fixed_params`, to avoid changing the production estimand.
  validate_af_model(model_spec, fixed_params, constraints)
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/estimate.py:345: UserWarning: Factor 'MC': built-in transition 'linear' enumerates parameters over ALL factors including observed factors ('log_income_observed',), so they enter the production function with free coefficients. The AF model assumes observed factors affect skills only through the investment equation. Use a production-factors-only transition ('translog_af' or 'log_ces_af'), or pin every observed-factor transition coefficient to 0.0 via `fixed_params`, to avoid changing the production estimand.
  validate_af_model(model_spec, fixed_params, constraints)
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/estimate.py:345: UserWarning: Factor 'MN': built-in transition 'linear' enumerates parameters over ALL factors including observed factors ('log_income_observed',), so they enter the production function with free coefficients. The AF model assumes observed factors affect skills only through the investment equation. Use a production-factors-only transition ('translog_af' or 'log_ces_af'), or pin every observed-factor transition coefficient to 0.0 via `fixed_params`, to avoid changing the production estimand.
  validate_af_model(model_spec, fixed_params, constraints)

/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/params.py:357: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[idx, "lower_bound"] = bounds_distance
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/params.py:364: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc, "value"] = val
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/params.py:365: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc, "lower_bound"] = val
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/params.py:366: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc, "upper_bound"] = val
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/params.py:372: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc, "value"] = val
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/params.py:373: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc, "lower_bound"] = val
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/params.py:374: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[loc, "upper_bound"] = val
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/initial_period.py:429: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[idx, "value"] = obs_sds.get(parts[0], meas_sd * 0.5)
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/initial_period.py:431: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[idx, "value"] = 0.0
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/initial_period.py:377: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[idx, "value"] = max(obs_sd * 0.5, 0.01)
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/initial_period.py:383: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[idx, "value"] = 1.0
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/initial_period.py:392: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[idx, "value"] = 0.0

/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/transition_period.py:1079: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[idx, "value"] = 0.5
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/transition_period.py:1090: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[idx, "value"] = max(obs_sd * 0.5, 0.01)
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/transition_period.py:1096: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[idx, "value"] = 1.0

/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/transition_period.py:1079: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[idx, "value"] = 0.5
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/transition_period.py:1090: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[idx, "value"] = max(obs_sd * 0.5, 0.01)
/home/hmga/skillmodels-applications/skillmodels/src/skillmodels/af/transition_period.py:1096: PerformanceWarning: indexing past lexsort depth may impact performance.
  params.loc[idx, "value"] = 1.0

AF per-period success: [True, True, True], per-period log-likelihoods: ['-14.93', '-21.16', '-22.86']

af_posterior = get_af_posterior_states(
    af_result=af_result,
    model_spec=model,
    data=data,
    n_halton_points=200,
)
af_states = af_posterior["unanchored_states"]["states"]
af_decomp = decompose_measurement_variance(
    model_spec=model, params=af_result.params, filtered_states=af_states
)
af_reliability = summarize_measurement_reliability(af_decomp)
af_reliability

AMN: three-stage mixture-of-normals¶

AMN’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.

On this model AMN does not run standalone: estimate_amn deliberately refuses the restricted CES (log_ces) because its Stage-3 regression omits the Freyberger (2025) primitive-scale recovery, so a standalone fit would be inconsistent. Its role here is instead to seed CHS — CHSEstimationOptions defaults to start_params_strategy="amn", so the CHS fit above already ran AMN’s three stages to build its start values. For a standalone AMN fit of a CES production function, use log_ces_general (AMN fits the transformed general CES and recovers the primitive scales via recover_primitive_ces_scales); see the AMN how-to for AMN on its native mixture-of-normals ground.

# AMN cannot consistently fit the restricted CES (`log_ces`) standalone, so here it
# seeds CHS instead: `CHSEstimationOptions` defaults to `start_params_strategy="amn"`,
# i.e. the CHS fit above already ran AMN's three stages to build its start values.
seed_strategy = chs_options.start_params_strategy
print(f"CHS start-value strategy: {seed_strategy!r} (AMN seeds CHS)")

CHS start-value strategy: 'amn' (AMN seeds CHS)

Cross-estimator comparison¶

The three estimators target the same likelihood under different approximations / objectives. Below we line up:

Period-0 measurement loadings for the skill_* indicators.
CES production-function gammas + φ for the period-0 → period-1 skills transition.
Signal-share scatter: do they agree on which measurements are high-noise?

estimators = {
    "CHS": chs_params,
    "AF": af_result.params,
}


def _free_loading(params, period, meas):
    loc = ("loadings", period, meas, "skills")
    if loc not in params.index:
        return None
    return float(params.loc[loc, "value"])


loading_rows = []
for est, params in estimators.items():
    for meas in SKILL_MEASURES:
        value = _free_loading(params, period=0, meas=meas)
        if value is None:
            continue
        loading_rows.append({"estimator": est, "measurement": meas, "loading": value})
loadings_df = pd.DataFrame(loading_rows)
loadings_pivot = loadings_df.pivot(
    index="measurement", columns="estimator", values="loading"
)
loadings_pivot

fig = go.Figure()
for est in ("CHS", "AF"):
    sub = loadings_df[loadings_df["estimator"] == est]
    fig.add_trace(go.Bar(name=est, x=sub["measurement"], y=sub["loading"]))
fig.update_layout(
    title="Period-0 skill loadings: CHS vs AF",
    barmode="group",
    template="plotly_white",
    yaxis_title="loading",
)
fig

def _safe_float(params, loc):
    if loc not in params.index:
        return float("nan")
    return float(params.loc[loc, "value"])


ces_rows: list[dict[str, float | str]] = []
for est, params in estimators.items():
    row: dict[str, float | str] = {"estimator": est}
    for col in ("skills", "investment"):
        row[f"gamma_{col}"] = _safe_float(params, ("transition", 0, "skills", col))
    row["phi"] = _safe_float(params, ("transition", 0, "skills", "phi"))
    ces_rows.append(row)
ces_df = pd.DataFrame(ces_rows).set_index("estimator")
ces_df

decomp_by_est = {"CHS": chs_decomp, "AF": af_decomp}
signal_share = pd.DataFrame(
    {est: d["fraction_signal"] for est, d in decomp_by_est.items()}
)
signal_share = signal_share.dropna()
fig = go.Figure()
fig.add_trace(
    go.Scatter(
        x=signal_share["CHS"],
        y=signal_share["AF"],
        mode="markers",
        name="AF vs CHS",
        marker={"size": 8},
    )
)
lo = float(signal_share.min().min())
hi = float(signal_share.max().max())
fig.add_trace(
    go.Scatter(
        x=[lo, hi],
        y=[lo, hi],
        mode="lines",
        line={"dash": "dash", "color": "gray"},
        showlegend=False,
    )
)
fig.update_layout(
    title="Signal share by measurement (CHS reference vs AF)",
    template="plotly_white",
    xaxis_title="signal share (CHS)",
    yaxis_title="signal share (AF)",
    height=520,
)
fig

Next steps¶

How to estimate AF covers the AF API in depth.
How to estimate AMN shows AMN on its native ground (a synthetic 2-mixture DGP where the mixture-of-normals advantage is visible).
How to compare estimators extends this tutorial with 95% confidence intervals (CHS analytic OPG / inverse-score; AF a propagated influence-function score bootstrap; AMN a cluster bootstrap) and overlaid posterior-factor trajectories.
Estimator prerequisites compares what data features and model constructs each estimator supports.
Architecture maps the common/chs/af/amn subpackage layout.