Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Solving and Simulating

Once you have defined a Model and prepared your parameters, pylcm solves via backward induction and simulates forward.

This page covers the common workflow. Solver-specific return artifacts and collective dissolution routing are specified exactly in Runtime, results, and persistence and Collective regimes.

Solving

solution = model.solve(params=params, log_level="debug")
value_functions = solution.values

Backward induction returns one immutable SolutionResult. Its values mapping is indexed by period -> regime_name -> value_function_array; replay policies, diagnostics, metadata, and omission reasons remain in their labelled fields and addressed artifact stores.

To inspect values with their state coordinates, use solution.value_array(period=0, regime="working"). It returns an xarray DataArray with one named dimension per state, so .sel(education="college") picks a curve regardless of axis order. Discrete states carry categorical labels by default, and a state declared with values carries those numbers; pass use_categorical_codes=True for integer codes. Collective values also have a stakeholder dimension. solution.value_frame(...) returns the same values as a long DataFrame with state columns and V; it explicitly materializes the whole grid. Use model.lookup_policy for values between grid nodes, with the same interpolation as simulation.

Retention and replay

The default VALUES_AND_REPLAY retention is the safe choice for later simulation:

from lcm.solver_api import ResultRetention

solution = model.solve(
    params=params,
    log_level="debug",
    retention=ResultRetention.VALUES_AND_REPLAY,
)

result = model.simulate(
    params=params,
    initial_conditions=initial_conditions,
    solution=solution,
    log_level="debug",
)

VALUES drops replay artifacts and works only when every simulated decision can be recovered from values. ALL_PERSISTABLE_ARTIFACTS keeps only what the result carries on its own: continuation or solver-defined artifacts are kept only where their model-built authority marks them as independently verifiable, and the NNBEGM replay policy of an AdaptiveOuterMesh search, which is replayed against a solve-generated mesh the model instance holds privately, is omitted as NOT_PERSISTED. Diagnostics follow log_level, not retention.

Simulation validates the durable model fingerprint, exact solution-relevant canonical parameters, solver/plugin and replay-route versions, value coverage and schemas, and every required replay artifact before forward execution. An in-memory result also has to come from the originating model instance. A restored result may come from another process: compatibility rests on the durable model fingerprint and exact declared versions instead. The checks remain active at log_level="off". Omit solution to solve automatically; there are no separate value, policy, or dissolution inputs. See Runtime, results, and persistence for the artifact stores and compatibility rules.

The time coordinate kind participates in compatibility: an age model and a period model cannot exchange solutions just because their numeric coordinates happen to coincide. Temporal parameter compatibility uses the aligned values; changes confined to silently discarded surplus rows do not change the consumed parameters.

Saving and restoring a solution

Save the complete result, including its metadata, omissions, and every retained artifact that has model-verifiable persistence authority:

from pathlib import Path

from lcm import load_solution, save_solution

path = Path("solution.lcm")
save_solution(solution=solution, path=path)

restored = load_solution(path=path)
result = model.simulate(
    params=params,
    initial_conditions=initial_conditions,
    solution=restored,
    log_level="debug",
)

solution.save(path=path) is the equivalent convenience method. Saving uses an atomic sibling-file replacement, so a failed write does not publish a partial archive. Values and artifacts are independently lazy after loading:

from lcm.solver_api import LoadState

assert restored.values.load_state(period=0, regime="working") is LoadState.UNLOADED
V_working = restored.value(period=0, regime="working")
assert restored.values.load_state(period=0, regime="working") is LoadState.LOADED

Loading one value leaves every other value and replay entry unloaded. Pass verify_checksums=True to load_solution to verify the entire archive without materializing any entry. Loading requires the exact solution-format, labelled-result schema, and solver-interface versions; replay also requires exact route, plugin, and artifact-schema identities. pylcm rejects a mismatch rather than migrating it silently. load_legacy_solution(path=...) is the explicit migration reader for the old value-only HDF5 format, which cannot be passed to simulate as a complete solution.

Log levels and runtime validation

log_level is required and controls console output. runtime_checks=True is the independent default on both solve() and simulate(): computed transition probabilities, NaN values and existing feasibility/replay checks are validated, and the first validation failure raises even with silent logging. Construction, parameter schema and solution-identity validation always run.

# Silent output with numerical checks enabled.
solution = model.solve(params=params, log_level="off", runtime_checks=True)

# Disable numerical checks explicitly for a trusted fast path.
solution = model.solve(params=params, log_level="warning", runtime_checks=False)

# Detailed output and optional disk snapshots.
solution = model.solve(params=params, log_level="debug", log_path="./debug/")
log_levelConsole output
"off"silent
"warning"warnings
"progress"warnings and timing
"debug"warnings, timing and value statistics

log_path is optional. A checked NaN failure saves a partial solve snapshot at any log level when a path is supplied. Debug logging additionally saves completed solves. runtime_checks=False skips numerical validation; it does not make an invalid result valid. Infinity-only value entries retain the existing diagnostic warning semantics: applications requiring every stored value finite must retain that stronger invariant.

Input schema and temporal alignment checks run at every log level. "off" does not permit missing required time labels, selected duplicate keys or a wrong coordinate kind. The manual-array warning in age mode is also independent of this logging policy.

simulate() adds one further record to the console output: the resolved execution plan it dispatched --- the forward route, the subject devices and their backend, the resolved planner axis widths by regime, the chunk count and admitted widths, and the budget mode with its effective device-memory bytes. A one-line summary appears at "progress" and "debug", the complete record at "debug" only, and neither at "warning" or "off". It is diagnostic only, and the same record is available as SimulationResult.plan_summary. See Inspect the resolved execution plan.

See Debugging for details on snapshots.

Simulating

result = model.simulate(
    params=params,
    initial_conditions=initial_conditions,
    solution=solution,
    log_level="debug",
)

Forward simulation using solved value functions. Each agent starts from the given initial conditions and makes optimal decisions at each period. Every subject’s starting time coordinate and regime must be an admissible root in model.graph.initial_nodes; a start elsewhere is rejected even where that problem’s value is solved. Returns a SimulationResult object. The complete SolutionResult is supplied through solution=....

Simulate without pre-solving

When solution is omitted, simulate() solves the model automatically before simulating. Use this when you don’t need the raw value function arrays:

result = model.simulate(
    params=params,
    initial_conditions=initial_conditions,
    log_level="debug",
)

Initial Conditions

From a DataFrame

The standard way to supply initial conditions is as a pandas DataFrame with one row per agent. Pass it directly to simulate():

import pandas as pd

df = pd.DataFrame(
    {
        "regime_name": ["working_life", "working_life", "retirement", "working_life"],
        "age": [25.0, 25.0, 25.0, 25.0],
        "wealth": [1.0, 5.0, 10.0, 20.0],
        "health": ["good", "bad", "bad", "good"],  # string labels, auto-converted
    }
)

result = model.simulate(
    params=params,
    initial_conditions=df,
    log_level="debug",
)

Discrete states (those backed by a DiscreteGrid) are mapped from string labels to integer codes automatically. See Working with DataFrames and Series for details.

As JAX arrays

You can also pass initial conditions directly as JAX arrays — useful for programmatic setups like grid searches or tests:

initial_conditions = {
    "age": jnp.array([25.0, 25.0, 25.0, 25.0]),
    "wealth": jnp.array([1.0, 5.0, 10.0, 20.0]),
    "health": jnp.array([0, 1, 1, 0]),  # integer codes for discrete states
    "regime_id": jnp.array(
        [
            RegimeId.working_life,
            RegimeId.working_life,
            RegimeId.retirement,
            RegimeId.working_life,
        ]
    ),
}

Household roles

In a model with a collective regime, every simulated row carries a role: which stakeholder of the household that row is. The role belongs to the row, not to the run — one cohort holds both partners at once — and it decides which regime the row moves to when the household dissolves. Seed it with an "own_stakeholder" entry alongside the states:

initial_conditions = {
    "age": jnp.full(4, model.ages.values[0]),
    "wealth": jnp.array([1.0, 5.0, 10.0, 20.0]),
    "regime_id": jnp.full(4, model.regime_names_to_ids["couple"], dtype=jnp.int32),
    "own_stakeholder": jnp.array(
        [
            model.stakeholder_names_to_ids["f"],
            model.stakeholder_names_to_ids["m"],
            model.stakeholder_names_to_ids["f"],
            model.stakeholder_names_to_ids["m"],
        ],
        dtype=jnp.int32,
    ),
}

As a DataFrame the same column carries stakeholder labels, converted to codes like any other discrete column:

df = pd.DataFrame(
    {
        "regime_name": ["couple", "couple"],
        "age": [25.0, 25.0],
        "wealth": [1.0, 5.0],
        "own_stakeholder": ["f", "m"],
    }
)

See Collective regimes for the full rules.

Checking initial conditions without simulating

A proposed population can be checked before any solve or simulation. Both methods take the same initial_conditions (mapping or DataFrame) and params as simulate:

mask = model.initial_conditions_feasibility(initial_conditions=df, params=params)
feasible_rows = df[np.asarray(mask)]

model.validate_initial_conditions(initial_conditions=feasible_rows, params=params)

initial_conditions_feasibility returns a boolean array with one entry per row, True where at least one combination of the regime’s action-grid points satisfies every constraint jointly (or, in a regime without actions, where every state-only constraint holds). validate_initial_conditions raises InvalidInitialConditionsError exactly when that mask contains a False, with the diagnostics simulate(log_level="debug") would print. A malformed structure — unknown regime, missing or extra state, unequal lengths, off-grid age, inactive regime, unknown categorical label, a collective start without its own_stakeholder — raises from both methods rather than being counted as infeasible, so a population they accept is one simulate accepts. Neither method takes a log_level: they always raise. A constraint that depends on an age-specialized function while subjects start at different ages cannot be evaluated and raises UnsupportedOperationError.

Further arguments

Typed roots are consumed directly; integer seeds keep their Threefry streams. Raw uint32 keys from jax.random.PRNGKey and batched typed keys are refused. Derive scalar role roots with fold_in:

import jax

root = jax.random.key(123)
result = model.simulate(
    params=params,
    initial_conditions=initial_conditions,
    seed=jax.random.fold_in(root, 0),
    taste_shock_seed=jax.random.fold_in(root, 1),
    log_level="off",
)

Common taste shocks across counterfactuals

Pass the same taste_shock_seed to compare policies using the same standardized EV1 shocks, while choosing the ordinary seed independently:

baseline = baseline_model.simulate(
    params=baseline_params,
    initial_conditions=initial_conditions,
    seed=100,
    taste_shock_seed=200,
    log_level="debug",
)
counterfactual = counterfactual_model.simulate(
    params=counterfactual_params,
    initial_conditions=initial_conditions,
    seed=101,
    taste_shock_seed=200,
    log_level="debug",
)

A taste key identifies a computational period, the destination regime’s name and a subject’s original row position in the initial conditions. Age labels, policy parameters, horizon endpoints, subject chunks and padding do not identify the key. Keep people in corresponding input rows; DataFrame index labels do not provide subject identity. Renaming the destination changes its taste stream. Changing the policy can change choices and state paths even with common shocks.

Action grids do not enter the key address. Comparing individual standardized taste shocks still requires the same ordered discrete-action domain and draw shape: category labels are not matched across incompatible domains. Changing a taste-shock scale multiplies the same standardized draw, so equal scaled shocks additionally require equal scales.

The independent stream uses Threefry for integer seeds. A supplied typed root retains its PRNG implementation; with an RBG or unsafe-RBG root, each subject’s draw comes from its own addressed key under every subject width and padding. Exact comparisons require matching numerical precision, backend and JAX random configuration. Ordinary and taste innovations use distinct versioned address domains, and their roots remain unchanged. taste_shock_seed=None uses the ordinary root with the taste address domain. Simulation results expose the address version, PRNG implementations and precision through result.randomness; saved results retain that metadata. The optional seed changes neither solutions nor their compatibility fingerprints.

Heterogeneous initial ages

This section describes age models. A period model instead requires "period", with integer indices in [0, model.n_periods), and uses InitialNodes(by_period={0: "working"}) to declare admissible starts. Booleans, fractional values and out-of-range starts are invalid. The two initial-condition checking methods apply the same coordinate rules.

"age" must always be provided in initial_conditions. Each value must be a valid point on the model’s AgeGrid, and each subject’s (age, regime) pair must be a declared start in model.graph.initial_nodes. The most common case is that all subjects start at the initial age — just pass a constant array.

Subjects can start at different ages:

initial_conditions = {
    "age": jnp.array([40.0, 60.0]),
    "wealth": jnp.array([50.0, 50.0]),
    "regime_id": jnp.array(
        [
            model.regime_names_to_ids["working_life"],
            model.regime_names_to_ids["working_life"],
        ]
    ),
}

In the resulting DataFrame, each subject appears only from their starting age onward — earlier periods are omitted, not filled with placeholders.

Working with SimulationResult

Converting to DataFrame

df = result.to_dataframe()

Returns a pandas DataFrame with columns: subject_id, period, regime_name, value, plus all states and actions. An age model also supplies age; a period model has no age column, including after saving and loading. An NNBEGM regime adds nested_policy_fallback: True on a row means the off-grid nested policy read was refused, so the row carries the best admissible baseline instead. That baseline is chosen by the canonical Q, not by the action grid alone — the grid-argmax pair and every published replay branch are scored, the higher score is emitted, and the grid pair takes an exact tie. Inference must refuse whenever any entry is True. Discrete variables are pandas Categorical with string labels. In an EGM-family regime, value is the solved value representation read at the subject’s realized resources, not the utility plus continuation along the replayed path; on a coarse savings grid the two can differ slightly in either direction.

A model with a collective regime publishes two further things. An own_stakeholder column names the role each row occupies — in every regime, not only the collective ones, because a row that has left a household still has to say that it now occupies none. It is a Categorical over the declared stakeholder names, and a row in a singleton regime carries a missing entry. A collective regime also publishes one value_<stakeholder> column per stakeholder, since the household stores every partner’s own value at the shared maximizing action. Those columns sit alongside value, which is dropped only when no regime in the model publishes a scalar value — that is, when every regime is collective.

Additional targets

Functions used solely for reporting must be registered with output_only_functions= when building the model or regime. Functions already read by the model are available without that declaration. See Output-only functions.

Compute functions and constraints alongside the standard output:

# Specific targets
df = result.to_dataframe(additional_targets=["utility", "consumption"])

# All available targets
df = result.to_dataframe(additional_targets="all")

# The published target tuple can be passed directly
df = result.to_dataframe(additional_targets=result.available_targets)

# See what's available
result.available_targets  # ('consumption', 'earnings', 'utility', ...)

Each target is computed for regimes where it exists; rows from other regimes get NaN.

Integer codes instead of labels

df = result.to_dataframe(use_categorical_codes=True)

Returns discrete variables as raw integer codes instead of categorical labels.

Metadata

result.regime_names  # ('retirement', 'working_life')
result.state_names  # ('health', 'wealth')
result.action_names  # ('consumption', 'work')
result.n_periods  # 50
result.n_subjects  # 1000

Persistence

SimulationResult.save(directory=...) writes four sibling artifacts:

Two keywords select what that Feather table holds; neither affects the other three artifacts. save(directory=..., df_additional_targets=...) is passed through to to_dataframe: None (the default) writes only the base columns (states, actions, regime, period, V_arr, and age for age models), a sequence of target names bakes those DAG outputs in, and "all" includes every available target --- which can grow the file by an order of magnitude on a model with many DAG leaves. df_use_categorical_codes=False (the default) stores discrete variables as pandas Categorical columns; True stores integer codes.

save() consumes the in-memory result by clearing its value-function arrays and compiled regimes. Reload the saved directory before further access that needs either.

from pathlib import Path

from lcm import SimulationResult

# Save
result.save(directory=Path("my_results"))

# Load (reads arrays + V_arr + metadata; the arrow file is for downstream consumers)
loaded = SimulationResult.load(directory=Path("my_results"))

Raw data (advanced)

result.raw_results  # regime -> period -> PeriodRegimeSimulationData
result.flat_params  # processed parameter object
result.period_to_regime_to_V_arr  # value function arrays from solve()

Typical Workflow

import numpy as np
import pandas as pd
from lcm import Model

# 1. Define model (see previous pages)
model = Model(
    regimes={...}, ages=..., regime_id_class=..., edges=..., initial_nodes=...
)

# 2. Set parameters
params = {
    "discount_factor": 0.95,
    "interest_rate": 0.03,
    # Add the model-specific parameter branches from the template.
}

# 3. Prepare initial conditions as a DataFrame
initial_df = pd.DataFrame(
    {
        "regime_name": "working_life",
        "age": model.ages.values[0],
        "wealth": np.linspace(1, 50, 100),
    }
)

# 4. Simulate (solves automatically when solution is omitted)
result = model.simulate(
    params=params,
    initial_conditions=initial_df,
    log_level="debug",
)

# 5. Analyze
df = result.to_dataframe(additional_targets="all")
df.groupby("period")["wealth"].mean()

Float32 GPU Reproducibility

See Also

Given shock paths

Supply realised shocks to the forward simulation without changing the distribution agents use when solving:

result = model.simulate(
    params=params,
    solution=solution,
    initial_conditions=initial_conditions,
    shock_paths={"income_shock": common_income_path},
    seed=42,
    log_level="off",
)

Each path is a real array shaped (periods,) for a common history or (subjects, periods) for individual histories in the original initial-condition row order. Columns are computational periods, beginning at the model’s first age or period; a subject entering later reads its entry column. The given value overrides that subject’s initial shock and each subsequent realised draw. Continuous processes accept finite values in their distribution’s support, including off-node values; lognormal levels must be positive and uniform draws must lie inside the declared interval. Discrete shocks require exact category codes, also for a state declared with values, whose other boundaries use the values: a code replays the path exactly and losslessly. grid.values[code] maps a code to its value. Persistent processes consume the supplied realised level directly, rather than treating it as an innovation and recursing again.

Each value is checked against the support of every regime that may read it: the subject’s starting regime in its entry period, and afterwards every regime the simulation can reach from there, since later regime membership is itself drawn. A subject’s entry value is replaced only when its starting regime draws the shock. Cells no regime can read, such as those before a subject enters, are ignored and may hold NaN. Draws taken inside a transition are supported, including shocks a regime reads only as next_<shock>: a draw affecting the transition into period p comes from column p. A deterministic law reading next_<shock> sees the supplied draw before computing its own next state. Shocks omitted from the mapping retain the ordinary seeded stream, and solve expectations remain stochastic.

A named current shock declared in a ShockGroup takes a tagged record instead of a bare array: BaseInnovations(values=...) for independent base coordinates, mapped through the current law, or PhysicalValues(values=...) for realised values. A correlated group takes physical values only as a complete frame; a partial frame, or physical values mixed with base coordinates, is refused. Tagged entry values and their precedence are described in the shocks reference.

pylcm
Working with DataFrames and Series
pylcm
Choosing a solver