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¶

period_to_regime_to_V_arr = model.solve(params=params, log_level="debug")

Performs backward induction using dynamic programming. Returns an immutable mapping of period -> regime_name -> value_function_array.

Optional solve artifacts¶

Most users need only value functions. Request additional artifacts explicitly:

value_functions, policies = model.solve(
    params=params,
    log_level="debug",
    return_simulation_policy=True,
)

value_functions, dissolution_flags = collective_model.solve(
    params=params,
    log_level="debug",
    return_dissolution_flags=True,
)

Simulation policies are for inspection; a fresh simulation obtains supported policy artifacts internally. Dissolution flags are required when a collective regime’s ValueDependentTransition has a gate that reads D_target — the dissolution flag of a collective target — during simulation.

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
period_to_regime_to_V_arr = model.solve(params=params, log_level="debug")

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

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

# Diagnostics + disk snapshots
period_to_regime_to_V_arr = 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,
    period_to_regime_to_V_arr=period_to_regime_to_V_arr,
    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.

Simulate without pre-solving¶

When period_to_regime_to_V_arr=None, 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,
    period_to_regime_to_V_arr=None,
    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,
    period_to_regime_to_V_arr=None,
    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).
period_to_regime_to_dissolution_flags=None: The flags returned by solve(return_dissolution_flags=True). Required only for a model whose gate reads D_target; a no-op for every other model.
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 period_to_regime_to_V_arr=None)
result = model.simulate(
    params=params,
    initial_conditions=initial_df,
    period_to_regime_to_V_arr=None,
    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)).