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, collective dissolution flags, diagnostics, metadata, and omission reasons remain in their labelled fields and addressed artifact stores.

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 and no applicable collective gate needs a dissolution flag. 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.

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 a required argument: it controls both console verbosity and the runtime-validation policy — how solve() / simulate() react to an invalid transition-probability ensemble or a NaN value function. Start every project at "debug" (validation runs and raises); ease to "warning" / "off" once the model is trusted.

Run every model at log_level="debug" at least once, with the parameters you intend to use, before trusting any output.

"off" skips runtime validation entirely, and some of what it skips cannot be reconstructed from the result. If a regime transition puts probability on a target that is not active in the next period, that mass is dropped and the surviving targets are renormalized — which would otherwise make the value function come back finite, plausible, and independent of the mass that went missing.

A gross departure from unit mass is caught in the arithmetic at every log level and turns the continuation into NaN. But that is a backstop, not a diagnosis: it names no regime, no target, and no age. Smaller mass errors pass it and still reverse the optimal action. Only the validation that "off" skips tells you what is wrong. See Debugging.

# Debug — validation runs and raises on the first failure
solution = model.solve(params=params, log_level="debug")

# Silent — no logging, no validation
solution = model.solve(params=params, log_level="off")

# Validation runs but only warns; the run continues
solution = model.solve(params=params, log_level="warning")

# Diagnostics + disk snapshots
solution = model.solve(params=params, log_level="debug", log_path="./debug/")

The full behaviour of every log_level × log_path combination:

`log_level`	`log_path`	Runtime validation	Console output	Snapshots to disk
`"off"`	(ignored)	not run	silent	none
`"warning"`	`None`	runs → failures warn	warnings	none
`"warning"`	set	runs → failures warn	warnings	one per warned failure, capped at `log_keep_n_latest`
`"progress"`	`None`	runs → failures warn	warnings + timing	none
`"progress"`	set	runs → failures warn	warnings + timing	one per warned failure, capped at `log_keep_n_latest`
`"debug"`	`None`	runs → failures raise	warnings + timing + V_arr stats	none
`"debug"`	set	runs → failures raise	warnings + timing + V_arr stats	one per solve and on raise, capped at `log_keep_n_latest`

log_path is optional at every level — snapshots are written only when it is set. In "warning" / "progress" mode, an invalid model produces warnings and a numerically meaningless result rather than an exception; use this to keep an estimation loop running, but read the warnings.

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. 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,
        ]
    ),
}

Every non-shock state must have an entry.
"regime_id" must be included, with integer codes from the regime_id_class.
All arrays must have the same length (= number of agents).
Shock states are drawn automatically.

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"],
    }
)

Codes come from model.stakeholder_names_to_ids — one vocabulary for the whole model, so a role means the same thing in every regime. It is empty for a model with no collective regime.
The entry is required for subjects who start in a collective regime that can reach one declaring a value-dependent transition with more than one route — its own, or one further along, since a row keeps its role across an ordinary regime transition. Omitting it there raises InvalidInitialConditionsError rather than defaulting to whichever partner was declared first, and a code that is valid model-wide but names no stakeholder of that regime is rejected too. A collective start that can never arrive at such a transition needs no entry.
A subject starting in a singleton regime occupies no role and needs no entry. It is given one on entering a collective regime, from the target_stakeholder of the StakeholderRoute it takes, and loses it again on landing in a singleton regime.

See Collective regimes for the full rules.

Further arguments¶

log_level: Required. Console verbosity and runtime-validation policy (same options and table as solve()); start at "debug". Initial-condition validation (states on-grid, regimes valid) follows this policy too — "off" skips it.
seed=None: Random seed for stochastic simulations (int). Collective dissolution gates consume their addressed replay artifacts from solution; the automatic-solve path retains and threads them itself.
log_path=None: Directory for diagnostic snapshots; optional at every level.
log_keep_n_latest=3: Maximum snapshot directories to retain.

Heterogeneous initial ages¶

"age" must always be provided in initial_conditions. Each value must be a valid point on the model’s AgeGrid, and each subject’s initial regime must be active at their starting age. 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, age, regime_name, value, plus all states and actions. 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.

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¶

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

# 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_labels=False)

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:

arrays/ — orbax checkpoint of the per-subject raw_results tree and flat_params.
V_arr/ — orbax checkpoint of the solution value-function arrays; orbax streams sharded leaves shard by shard rather than gathering them onto one device.
metadata.pkl — cloudpickle of regimes, ages, and the parameter scaffold.
simulated_data.arrow — a feather dump of to_dataframe, ready for downstream consumers that want the flat per-subject view without re-instantiating a SimulationResult.

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

# 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¶

Note

At float32 precision, GPU simulation results are not reproducible across process invocations. XLA compiles different fused kernels each time a process starts, changing float32 accumulation order and producing ~1e-3 value function differences for large models. This is a property of XLA’s GPU compiler, not a pylcm bug.

Key facts:

Float64 results are reproducible across processes.
Float32 results are deterministic within a single process (repeated calls give identical output).
--xla_gpu_deterministic_ops=true does not help — it only guarantees determinism within a single compiled program, not across separate compilations.
Smaller models may be reproducible at float32 because XLA chooses the same kernel fusion strategy; larger models are not.

If you need bitwise-reproducible results for testing or validation, use float64 precision (jax.config.update("jax_enable_x64", True)).