Beta-Delta (Quasi-Hyperbolic) Discounting

This notebook shows how to implement naive beta-delta discounting in PyLCM using a custom aggregation function $H$ and SolveSimulateFunctionPair. It covers:

The beta-delta discount function
Why naive agents need different $H$ for solve vs. simulate
A 3-period consumption-savings model with analytical solutions
Verifying numerical results against closed-form solutions

Background: Beta-Delta Discounting¶

Standard exponential discounting uses a single discount factor $\delta$ to weight future utility:

U_t = u_t + \delta\, u_{t+1} + \delta^2\, u_{t+2} + \cdots

(1)

Quasi-hyperbolic (beta-delta) discounting (Laibson, 1997) introduces a present-bias parameter $\beta \in (0, 1]$ that discounts all future periods relative to the present:

U_t = u_t + \beta\delta\, u_{t+1} + \beta\delta^2\, u_{t+2} + \cdots

(2)

The discount weights are $\{1,\; \beta\delta,\; \beta\delta^2,\; \ldots\}$ . When $\beta = 1$ this reduces to exponential discounting. When $\beta < 1$ the agent is present-biased — they overweight current utility relative to any future period.

Note the key structure: $\beta$ appears once (not compounded). The weight $n$ periods ahead is $\beta\delta^n$ , not $(\beta\delta)^n$ .

PyLCM’s $H$ Function and the Naive Agent¶

PyLCM defines the value function recursively via an aggregation function $H$ :

V_t(s) = \max_a \; H\bigl(u(s, a),\; \mathbb{E}_t[V_{t+1}(s')]\bigr)

(3)

The default $H$ is $H(u, \mathbb{E}[V']) = u + \delta \cdot \mathbb{E}[V']$ .

A naive beta-delta agent believes their future selves will be exponential discounters, but at each decision point applies present bias $\beta\delta$ to the continuation value. This requires two different $H$ implementations:

Solve phase (backward induction): use exponential discounting with $\delta$ to compute the perceived continuation values $V^E$
Simulate phase (action selection): use $\beta\delta$ to choose actions, applying present bias to the perceived $V^E$

SolveSimulateFunctionPair wraps these two implementations so you can call simulate once with a single params dict.

Why not just use $H(u, V') = u + \beta\delta V'$ for both phases? Because the recursion would compound $\beta$ at every step, giving effective weights $(\beta\delta)^n$ instead of the correct $\beta\delta^n$ . Splitting the phases ensures $\beta$ is applied exactly once — at the action selection step.

Setup: A 3-Period Model¶

We use a minimal consumption-savings model:

3 periods: $t = 0, 1$ (decisions), $t = 2$ (terminal)
1 state: wealth $w$
1 action: consumption $c$
Utility: $u(c) = \log(c)$
Budget: $w' = w - c$ (interest rate $R = 1$ )
Constraint: $c \le w$
Terminal: $V_2(w) = \log(w)$ (consume everything)

import jax.numpy as jnp

from lcm import (
    AgeGrid,
    LinSpacedGrid,
    Model,
    Regime,
    SolveSimulateFunctionPair,
    categorical,
)
from lcm.typing import BoolND, ContinuousAction, ContinuousState, FloatND, ScalarInt

@categorical(ordered=False)
class RegimeId:
    working: int
    dead: int


def utility(consumption: ContinuousAction) -> FloatND:
    return jnp.log(consumption)


def terminal_utility(wealth: ContinuousState) -> FloatND:
    return jnp.log(wealth)


def next_wealth(
    wealth: ContinuousState, consumption: ContinuousAction
) -> ContinuousState:
    return wealth - consumption


def next_regime(age: float) -> ScalarInt:
    return jnp.where(age >= 1, RegimeId.dead, RegimeId.working)


def borrowing_constraint(
    consumption: ContinuousAction, wealth: ContinuousState
) -> BoolND:
    return consumption <= wealth


def exponential_H(
    utility: float,
    E_next_V: float,
    discount_factor: float,
) -> float:
    return utility + discount_factor * E_next_V


def beta_delta_H(
    utility: float,
    E_next_V: float,
    beta: float,
    delta: float,
) -> float:
    return utility + beta * delta * E_next_V

We build two models:

One with exponential_H for the standard exponential agent
One with SolveSimulateFunctionPair for the naive beta-delta agent — solving with exponential $H$ and simulating with beta-delta $H$

WEALTH_GRID = LinSpacedGrid(start=0.5, stop=50, n_points=200)
CONSUMPTION_GRID = LinSpacedGrid(start=0.001, stop=50, n_points=500)

dead = Regime(
    transition=None,
    states={"wealth": WEALTH_GRID},
    functions={"utility": terminal_utility},
    active=lambda age: age > 1,
)

# Standard exponential model
working_exp = Regime(
    actions={"consumption": CONSUMPTION_GRID},
    states={"wealth": WEALTH_GRID},
    state_transitions={"wealth": next_wealth},
    constraints={"borrowing_constraint": borrowing_constraint},
    transition=next_regime,
    functions={"utility": utility, "H": exponential_H},
    active=lambda age: age <= 1,
)

model_exp = Model(
    regimes={"working": working_exp, "dead": dead},
    ages=AgeGrid(start=0, stop=2, step="Y"),
    regime_id_class=RegimeId,
)

# Naive beta-delta model (different H per phase)
working_naive = Regime(
    actions={"consumption": CONSUMPTION_GRID},
    states={"wealth": WEALTH_GRID},
    state_transitions={"wealth": next_wealth},
    constraints={"borrowing_constraint": borrowing_constraint},
    transition=next_regime,
    functions={
        "utility": utility,
        "H": SolveSimulateFunctionPair(
            solve=exponential_H,
            simulate=beta_delta_H,
        ),
    },
    active=lambda age: age <= 1,
)

model_naive = Model(
    regimes={"working": working_naive, "dead": dead},
    ages=AgeGrid(start=0, stop=2, step="Y"),
    regime_id_class=RegimeId,
)

Analytical Solution¶

With $\log$ utility, the optimal consumption rule is $c_t = w_t / D_t$ , where $D_t$ is a “marginal propensity to save” denominator.

At $t = 1$ (one period before terminal), the agent maximizes:

\max_c \;\bigl[\log(c) + \beta\delta \cdot \log(w - c)\bigr]

(4)

First-order condition: $1/c = \beta\delta / (w - c)$ , giving $c_1 = w / (1 + \beta\delta)$ .

At $t = 0$ , the naive agent uses the perceived exponential continuation value $V^E_1$ , which has coefficient $(1 + \delta)$ on $\log w$ . So the agent maximizes:

\max_c \;\bigl[\log(c) + \beta\delta (1 + \delta) \log(w - c) + \text{const}\bigr]

(5)

giving $c_0 = w / (1 + \beta\delta(1 + \delta))$ .

	$D_1$	$D_0$
Exponential ( $\beta = 1$ )	$1 + \delta$	$1 + \delta(1 + \delta)$
Naive ( $\beta < 1$ )	$1 + \beta\delta$	$1 + \beta\delta(1 + \delta)$

BETA = 0.7
DELTA = 0.95
W0 = 20.0


def analytical_consumption(beta, delta, w0):
    """Return (c0, c1) from the closed-form solution."""
    bd = beta * delta
    d1 = 1 + bd
    d0 = 1 + bd * (1 + delta)
    c0 = w0 / d0
    c1 = (w0 - c0) / d1
    return c0, c1


c0_exp, c1_exp = analytical_consumption(1.0, DELTA, W0)
c0_naive, c1_naive = analytical_consumption(BETA, DELTA, W0)

print(f"{'Type':<15} {'c_0':>8} {'c_1':>8}")
print("-" * 33)
print(f"{'Exponential':<15} {c0_exp:8.4f} {c1_exp:8.4f}")
print(f"{'Naive':<15} {c0_naive:8.4f} {c1_naive:8.4f}")

Type                 c_0      c_1
---------------------------------
Exponential       7.0114   6.6608
Naive             8.7080   6.7820

The naive agent consumes more at $t = 0$ than the exponential agent — present bias pulls consumption forward. At $t = 1$ , the naive agent also consumes more because $\beta\delta < \delta$ .

Numerical Results¶

Exponential Agent ( $\beta = 1$ )¶

initial_conditions = {
    "age": jnp.array([0.0]),
    "wealth": jnp.array([W0]),
    "regime": jnp.array([RegimeId.working]),
}

result_exp = model_exp.simulate(
    params={"working": {"H": {"discount_factor": DELTA}}},
    initial_conditions=initial_conditions,
    period_to_regime_to_V_arr=None,
)

df_exp = result_exp.to_dataframe().query('regime == "working"')
print("Exponential agent:")
print(
    f"  c_0 = {df_exp.loc[df_exp['age'] == 0, 'consumption'].iloc[0]:.4f}  "
    f"(analytical: {c0_exp:.4f})"
)
print(
    f"  c_1 = {df_exp.loc[df_exp['age'] == 1, 'consumption'].iloc[0]:.4f}  "
    f"(analytical: {c1_exp:.4f})"
)

Starting solution

Age: 2  regimes=1  (116.9ms)

Age: 1  regimes=1  (123.8ms)

Age: 0  regimes=1  (57.4ms)

Solution complete  (300.4ms)

Starting simulation

Age: 0  regimes=1  (745.7ms)

Age: 1  regimes=1  (98.7ms)

Age: 2  regimes=1  (35.7ms)

Simulation complete  (1.0s)

Exponential agent:
  c_0 = 7.0149  (analytical: 7.0114)
  c_1 = 6.7143  (analytical: 6.6608)

Naive Agent¶

The naive agent uses model_naive, built with a SolveSimulateFunctionPair wrapping $H$ . During backward induction, exponential_H computes the perceived value function $V^E$ ; during simulation, beta_delta_H applies present bias for action selection.

The params dict is the union of both variants’ parameters — discount_factor for exponential_H, plus beta and delta for beta_delta_H. Each variant receives only the kwargs it expects.

result_naive = model_naive.simulate(
    params={
        "working": {
            "H": {"discount_factor": DELTA, "beta": BETA, "delta": DELTA},
        },
    },
    initial_conditions=initial_conditions,
    period_to_regime_to_V_arr=None,
)

df_naive = result_naive.to_dataframe().query('regime == "working"')
print("Naive agent:")
print(
    f"  c_0 = {df_naive.loc[df_naive['age'] == 0, 'consumption'].iloc[0]:.4f}  "
    f"(analytical: {c0_naive:.4f})"
)
print(
    f"  c_1 = {df_naive.loc[df_naive['age'] == 1, 'consumption'].iloc[0]:.4f}  "
    f"(analytical: {c1_naive:.4f})"
)

Starting solution

Age: 2  regimes=1  (33.3ms)

Age: 1  regimes=1  (60.3ms)

Age: 0  regimes=1  (58.7ms)

Solution complete  (154.7ms)

Starting simulation

Age: 0  regimes=1  (140.1ms)

Age: 1  regimes=1  (89.6ms)

Age: 2  regimes=1  (27.1ms)

Simulation complete  (261.3ms)

Naive agent:
  c_0 = 8.7183  (analytical: 8.7080)
  c_1 = 6.8145  (analytical: 6.7820)

Comparison¶

The small differences between numerical and analytical solutions are due to grid discretization.

import pandas as pd

comparison = pd.DataFrame(
    {
        "Agent": ["Exponential", "Naive"],
        "c_0 (numerical)": [
            df_exp.loc[df_exp["age"] == 0, "consumption"].iloc[0],
            df_naive.loc[df_naive["age"] == 0, "consumption"].iloc[0],
        ],
        "c_0 (analytical)": [c0_exp, c0_naive],
        "c_1 (numerical)": [
            df_exp.loc[df_exp["age"] == 1, "consumption"].iloc[0],
            df_naive.loc[df_naive["age"] == 1, "consumption"].iloc[0],
        ],
        "c_1 (analytical)": [c1_exp, c1_naive],
    }
)
comparison = comparison.set_index("Agent")
comparison.style.format("{:.4f}")

Summary¶

Beta-delta discounting in PyLCM requires no core modifications:

Define two $H$ functions — one exponential, one with present bias:

def exponential_H(utility, E_next_V, discount_factor):
    return utility + discount_factor * E_next_V

def beta_delta_H(utility, E_next_V, beta, delta):
    return utility + beta * delta * E_next_V

Wrap them in SolveSimulateFunctionPair so backward induction uses exponential $H$ while action selection uses beta-delta $H$ :

from lcm import SolveSimulateFunctionPair

functions={
    "utility": utility,
    "H": SolveSimulateFunctionPair(
        solve=exponential_H,
        simulate=beta_delta_H,
    ),
}

Provide the union of both variants’ parameters:

params = {"working": {"H": {"discount_factor": delta, "beta": beta, "delta": delta}}}

This split is essential: using $H(u, V') = u + \beta\delta V'$ for both phases would compound $\beta$ at every recursive step, giving effective weights $(\beta\delta)^n$ instead of the correct $\beta\delta^n$ .