Brock-Mirman stochastic growth
The Brock-Mirman economy, the optimal growth of a representative agent whose productivity is Markov.
- Name:
brock_mirman_growth - Field:
economics/macroeconomics - Formalism: discrete map
- References: Caiani et al. 2016 §1.1.4; Brock and Mirman 1972
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('brock_mirman_growth')
entry = describe('brock_mirman_growth')
trajectory = create_simulation('brock_mirman_growth').prepare().run().trajectory_
Factory
@register_system(
field='economics/macroeconomics',
tags=('optimal growth', 'representative agent', 'markov chain', 'productivity shocks', 'discrete map'),
references=('Caiani et al. 2016 §1.1.4', 'Brock and Mirman 1972'),
defaults={'until': 200},
observables=lambda system: [
ObservableSpec(name, function, element='economy')
for name, function in (
('productivity', StochasticGrowthElement.get_productivity),
('output', StochasticGrowthElement.get_output),
('consumption', StochasticGrowthElement.get_consumption),
)
],
)
def brock_mirman_growth(
K: float = 0.2,
*,
alpha: float = 1 / 3,
beta: float = 0.95,
z_low: float = 0.9,
z_high: float = 1.1,
persistence: float = 0.9,
regime: int = 0,
) -> System:
"""Return the Brock-Mirman economy, the optimal growth of a representative agent whose productivity is Markov.
A representative consumer maximises `E_0 sum_t beta^t ln C_t` subject to `K_{t+1} = Y_t - C_t`, with full
depreciation and the output `Y_t = z_t K_t^alpha` of a representative firm, whose productivity `z_t` follows a
Markov chain. The solution is the linear consumption rule `C_t = (1 - alpha beta) Y_t`, `C = 0.68 Y` for
`alpha = 1/3` and `beta = 0.95`, so the capital follows `K_{t+1} = alpha beta z_t K_t^alpha` and `ln K` is a
stationary autoregression around `(ln(alpha beta) + E ln z) / (1 - alpha)`. This is the standard model that Caiani
et al. (2016, Section 1.1.4, Eqs. 1.3 to 1.9) compare with an agent-based economy. The source leaves the chain open,
so the productivity takes two levels that stay with the probability `persistence`.
Args:
K:
The initial capital.
alpha:
The capital share of the technology.
beta:
The discount factor of the household.
z_low:
The low productivity.
z_high:
The high productivity.
persistence:
The probability that the productivity keeps its level from one period to the next.
regime:
The initial regime, 0 for the low and 1 for the high productivity.
Returns:
system:
The system, with the capital `economy.K` and the regime `economy.regime`. The observables of the system are
`economy.productivity`, `economy.output`, and `economy.consumption`.
Examples:
>>> from skcomplex.systems.economics.macroeconomics import brock_mirman_growth
>>> brock_mirman_growth(beta=0.9).elements[0].beta
0.9
"""
chain = ((persistence, 1 - persistence), (1 - persistence, persistence))
return System(
elements=[
StochasticGrowthElement(
'economy',
K=K,
regime=regime,
alpha=alpha,
beta=beta,
z_low=z_low,
z_high=z_high,
),
],
interactions=[
MapInteraction(
'saving',
target='economy',
state='K',
function=lambda s: s.target.alpha * s.target.beta * StochasticGrowthElement.get_output(s.target),
),
TransitionInteraction('productivity', target='economy', state='regime', matrix=lambda s: chain),
],
)
Components
Elements:
StochasticGrowthElement, fromElement
Interactions:
MapInteractionTransitionInteraction
class StochasticGrowthElement(Element):
"""Representative economy with capital, a productivity regime, and a Cobb-Douglas technology `Y = z K^alpha`.
The productivity `z` is `z_low` in regime 0 and `z_high` in regime 1, and a Markov chain of the system moves the
regime. With logarithmic utility and full depreciation, the optimal consumption is `C = (1 - alpha beta) Y` and
the capital of the next period is the output saved, `K' = alpha beta z K^alpha`.
Args:
label:
The label of the economy.
visual:
How the element is drawn, or `None` for `default_visual_`.
K:
The capital.
regime:
The productivity regime, 0 for low and 1 for high.
alpha:
The capital share of the technology.
beta:
The discount factor of the household.
z_low:
The productivity of regime 0.
z_high:
The productivity of regime 1.
Examples:
>>> from skcomplex.systems.economics.macroeconomics import StochasticGrowthElement
>>> StochasticGrowthElement('economy', beta=0.9).beta
0.9
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
K: Annotated[float, StateSpec(domain=RealSet.open(0, oo))] = 0.2,
regime: Annotated[float, StateSpec(domain=NN)] = 0,
alpha: Annotated[float, ParameterSpec(domain=RealSet.open(0, 1))] = 1 / 3,
beta: Annotated[float, ParameterSpec(domain=RealSet.open(0, 1))] = 0.95,
z_low: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.9,
z_high: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.1,
) -> None:
super().__init__(label, visual=visual)
self.K = K
self.regime = regime
self.alpha = alpha
self.beta = beta
self.z_low = z_low
self.z_high = z_high
@staticmethod
def get_productivity(s: SimpleNamespace) -> Expression:
"""Return the productivity of the current regime.
Args:
s:
The namespace of the economy.
Returns:
z:
The productivity.
"""
return s.z_low + (s.z_high - s.z_low) * s.regime
@staticmethod
def get_output(s: SimpleNamespace) -> Expression:
"""Return the output `z K^alpha`.
Args:
s:
The namespace of the economy.
Returns:
output:
The output.
"""
return StochasticGrowthElement.get_productivity(s) * s.K**s.alpha
@staticmethod
def get_consumption(s: SimpleNamespace) -> Expression:
"""Return the optimal consumption `(1 - alpha beta) Y`.
Args:
s:
The namespace of the economy.
Returns:
consumption:
The consumption.
"""
return (1 - s.alpha * s.beta) * StochasticGrowthElement.get_output(s)
class MapInteraction(Interaction):
"""Map of a state of discrete time, whose next value is a function, or the state plus it.
In the mode `'replace'` the next value of the state is `f`, and in the mode `'increment'` it is the state plus `f`.
Under the `'synchronous'` schedule every member of a bound group updates at once, from the states at the start of
the step. Under the `'sequential'` schedule the members update one after the other, in their order, and under the
`'random-sequential'` schedule in a new random order every step, each member reading the values the members before
it set. Under the `'member'` schedule of the agent form, every member updates in its own activation, together with
the rules of its element and the other member interactions, in the order the system lists them.
In the `'symbolic'` form, `function(namespace)` returns a symbolic expression, or a vector of them for a vector
state, where `namespace.target` holds the states and parameters of the target, `namespace.sources` the namespaces
of the sources, and `namespace.t` the time, and may hold draws such as `normal('noise')`. In the `'agent'` form,
the target is an agent group and `function(agent, system)` returns the value as a number, or as an array of the
shape of a vector state, so it can read the agent, its neighbours, and the history and storage of the system view.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state the map sets.
state:
The name of the state of the target.
function:
The function `f(namespace)` in the symbolic form, or `f(agent, system)` in the agent form.
sources:
The paths of the elements, or members of groups such as `cells[0]`, that the function reads.
mode:
`'replace'` or `'increment'`.
schedule:
`'synchronous'`, `'sequential'`, `'random-sequential'`, or `'member'` for the agent form.
form:
`'symbolic'` for update laws of the states of the target, or `'agent'` for a rule over an agent group.
phase:
The phase of the step the agent form runs in, or `None` for the default phase.
Raises:
ValidationError:
The mode, the schedule, or the form is unknown, the symbolic form has the schedule `'member'`, the
function is not callable, or the agent form returns a value that does not match the state, raised by
`run` of the simulation.
Examples:
>>> from skcomplex.core.components.elements import Element
>>> from skcomplex.core.components.interactions import MapInteraction
>>> from skcomplex.core.components.spaces import DiscreteTime
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.specs import StateSpec
>>> from typing import Annotated, SupportsFloat
>>> class ValueElement(Element):
... def __init__(self, label, x: Annotated[float, StateSpec()] = 0.5):
... super().__init__(label)
... self.x = x
>>> logistic = MapInteraction(
... 'map', target='value', state='x', function=lambda s: 3 * s.target.x * (1 - s.target.x)
... )
>>> system = System(elements=[ValueElement('value')], interactions=[logistic], time=DiscreteTime())
>>> Simulation(system).prepare().run(until=2).trajectory_.values('value.x').tolist()
[0.5, 0.75, 0.5625]
>>> agents = MapInteraction('map', target='walkers', state='x', function=lambda a, m: a.x + 1, form='agent')
>>> type(agents.prepare().run(SimpleNamespace()).laws_[0]).__name__
'RuleLaw'
>>> from skcomplex.core.components.elements import GroupElement
>>> class WeightsElement(Element):
... def __init__(self, label, visual=None, *, w: Annotated[object, StateSpec(vector=True)] = 0.0):
... super().__init__(label, visual=visual)
... self.w = w
>>> group = GroupElement('weights', element=WeightsElement, size=1, values={'w': [1.0, 2.0]})
>>> doubling = MapInteraction(
... 'double', target='weights', state='w', function=lambda a, m: [2 * x for x in a.w], form='agent'
... )
>>> trajectory = Simulation(System(elements=[group], interactions=[doubling])).prepare().run(until=2)
>>> trajectory.trajectory_.agents('weights', 2).values['w'].tolist()
[[4.0, 8.0]]
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (RuleLaw, UpdateLaw)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('state',))],
state: str,
function: Callable[..., object],
sources: Annotated[Sequence[str], RoleSpec(variadic=True, elements=(ElementKind.ANY,))] = (),
mode: str = 'replace',
schedule: str = 'synchronous',
form: str = 'symbolic',
phase: str | None = None,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.function = function
self.sources = sources
self.mode = mode
self.schedule = schedule
self.form = form
self.phase = phase
def _check(self) -> None:
"""Raise an error for an unknown mode, schedule, or form, or a function that cannot be called.
Raises:
ValidationError:
The mode, the schedule, or the form is unknown, the symbolic form has the schedule `'member'`, or the
function is not callable.
"""
super()._check()
self._check_variant('mode', self.mode, MODES, (BaseSpec,))
for name, value, variants in (
('schedule', self.schedule, SCHEDULES),
('form', self.form, FORMS),
):
if value not in variants:
error_msg = f'The {name} of map {self.label!r} is {value!r}. Use one of {", ".join(variants)}.'
raise ValidationError(error_msg)
if self.schedule == 'member' and self.form != 'agent':
error_msg = (
f'The map {self.label!r} updates every member in its activation, which only the agent form does.'
)
raise ValidationError(error_msg)
if not callable(self.function):
error_msg = f'The function of map {self.label!r} must be callable, but it is {self.function!r}.'
raise ValidationError(error_msg)
def _apply_mode(self, current: Value, value: Value) -> Value:
"""Return the next value of a state from its current value and the value of the function, by the mode."""
if self.mode == 'increment':
return current + value
if self.mode == 'replace':
return value
raise self._get_variant_error(self.mode, 'mode', '_apply_mode')
def _get_vector(self, agent: AgentViewProtocol, system: SystemViewProtocol) -> Value:
"""Return the next value of a vector state of an agent."""
current = np.asarray(getattr(agent, self.state), dtype=np.float64)
value = np.asarray(self.function(agent, system), dtype=np.float64)
if value.shape != current.shape:
error_msg = (
f'The function of map {self.label!r} returns the shape {value.shape} for the state {self.state!r} of '
f'the shape {current.shape}. Return an array of the shape of the state.'
)
raise ValidationError(error_msg)
return self._apply_mode(current, value)
def _get_value(self, agent: AgentViewProtocol, system: SystemViewProtocol) -> Value:
"""Return the next value of the state of an agent, a number or the array of a vector state."""
if np.ndim(getattr(agent, self.state)) > 0:
return self._get_vector(agent, system)
found = self.function(agent, system)
if not isinstance(found, SupportsFloat):
error_msg = f'The function of map {self.label!r} returns {found!r} for an agent. Return a number.'
raise ValidationError(error_msg)
value = float(found)
return self._apply_mode(float(getattr(agent, self.state)), value)
def _set_value(self, agent: AgentViewProtocol, value: Value) -> None:
"""Set the state of an agent to its next value."""
setattr(agent, self.state, value)
def _update_member(self, agent: AgentViewProtocol, system: SystemViewProtocol) -> None:
"""Set the state of one member to its next value."""
if agent.alive:
self._set_value(agent, self._get_value(agent, system))
def update(self, system: SystemViewProtocol) -> None:
"""Set the state of every agent of the target to its next value, by the schedule.
Args:
system:
The system view of the interaction.
"""
self._update_agents(system, self.schedule, self._get_value, self._set_value)
def _get_order(self) -> object:
"""Return the order of the members of the symbolic form: none, their own order, or a uniform draw."""
if self.schedule == 'synchronous':
return None
if self.schedule == 'sequential':
return SR(0)
return uniform(re.sub(r'\W', '_', f'order_{self.label}'))
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the update of the state to its next value, or the rule of the agent form.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target` and `namespace.sources`.
Returns:
laws:
The update law of the state, or the rule that runs the map of the agents.
"""
if self.form == 'agent' and self.schedule == 'member':
return [RuleLaw(self._update_member, phase=self.phase, members=self._get_role('target'))]
if self.form == 'agent':
return [RuleLaw(self.update, phase=self.phase)]
value = getattr(self._get_member(namespace, 'target'), self.state)
mapped = cast(Value, self.function(namespace))
next_value = self._apply_mode(value, mapped)
return [UpdateLaw(value, next_value, increment=self.mode == 'increment', order=self._get_order())]
class TransitionInteraction(Interaction):
"""Move of a discrete state of the target from one value to another, or between all its values by a matrix.
A member whose state holds `source` moves to `destination` at a `rate`, as a jump in continuous time, with a
`probability` in every step of discrete time, or in the step in which a `condition` holds. With `over='neighbors'`
the rate and the probability are those of one neighbour, written with `namespace.target.other` for the neighbour:
the rates of the neighbours add up, and their probabilities combine as `1 - prod(1 - p)`, the probability that at
least one of them moves the member. With a matrix instead, the state takes the values `0` to `n - 1` and steps by
the rows of a transition matrix in discrete time, the row of its current value. An `emission` matrix then draws the
`observed` state of the target from the row that the next value of the state selects, as a hidden Markov source.
With a residence time `after`, a member moves by a probability or a condition only once it has held the source
value for `after` steps, at the earliest at the end of the step that completes them.
In the `'symbolic'` form, every function reads `namespace.target`, the states and parameters of the target, and the
time `namespace.t`. In the `'agent'` form, the target is an agent group and every function takes the agent and
the system view, `f(agent, system)`, and returns a number, a truth value, or the rows of a matrix, and a rate `r`
moves an agent with probability `1 - exp(-r)` in a step. With `over='cell'`, every member acts on the other members
of its cell: it moves each one that holds the source value with the probability, the rate, or the condition
`f(agent, other, system)` of that member `agent` and the acting member `other`, each by its own draw.
The `schedule` of the agent form orders the members. Under `'synchronous'` every member decides from the states at
the start of the step, under `'sequential'` the members act one after the other in their order, and under
`'random-sequential'` in a new random order every step, each reading the values the members before it set. Under
`'member'` every member acts in its own activation, together with the rules of its element and the other member
interactions, in the order the system lists them, so one rule of a member runs its movement, transfers,
transitions, and life cycle.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state moves.
state:
The name of the discrete state of the target.
source:
The value the state moves from, a collection of such values, or a function of the value of the state that
returns whether the state moves from it, such as `lambda value: value >= 0`, a condition in the symbolic
form. A state that moves from several values moves to the destination from each.
destination:
The value the state moves to.
rate:
The function of the rate of the move, or `None`.
probability:
The function of the probability of the move in a step, or `None`.
condition:
The function of the condition under which the state moves in a step, or `None`.
matrix:
The function of the transition matrix, one row and one column per value of the state, or `None`.
emission:
With a matrix, the function of the emission matrix, one row per value of the state and one column per value
of the observed state, or `None`.
observed:
The name of the state of the target that the emission matrix draws, or `None`.
over:
`'member'` for a rate or a probability of the member, `'neighbors'` for that of one neighbour, or `'cell'`
for the agent form whose members act on the members of their cell.
neighborhood:
The neighbourhood of `over='neighbors'`, or `None` for the neighbours of the network or the lattice.
after:
The number of steps a member holds the source value before it can move, or `None`.
schedule:
The order of the members of the agent form, `'synchronous'`, `'sequential'`, `'random-sequential'`, or
`'member'`.
form:
`'symbolic'` for laws of the states of the target, or `'agent'` for a rule over an agent group.
phase:
The phase of the step the agent form runs in, or `None` for the default phase.
Raises:
ValidationError:
The transition has not exactly one kind, the source and the destination are equal, an emission has no
matrix or no observed state, `over` is unknown or goes with a form or a kind it cannot take, the
residence time is negative or goes with a rate, a matrix, or members acting on their cell, the form or the
schedule is unknown, the symbolic form has a schedule, or a function is not callable.
Examples:
>>> from skcomplex.core.components.interactions import TransitionInteraction
>>> recovery = TransitionInteraction(
... 'recovery', target='hosts', state='state', source=1, destination=0, rate=lambda s: 0.5
... )
>>> recovery.over
'member'
>>> agents = TransitionInteraction(
... 'recovery', target='hosts', state='state', source=1, destination=2, probability=lambda a, m: 0.1,
... after=3, form='agent'
... )
>>> type(agents.prepare().run(SimpleNamespace()).laws_[0]).__name__
'RuleLaw'
>>> separation = TransitionInteraction(
... 'separation', target='workers', state='employer', source=lambda value: value >= 0, destination=-1,
... probability=lambda a, m: 0.05, form='agent'
... )
>>> separation.destination
-1
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (JumpLaw, UpdateLaw, RuleLaw, TransitionLaw, EmissionLaw)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('state', 'observed'))],
state: str = 'state',
source: SourceValues = 0.0,
destination: float = 1.0,
rate: Callable[..., SupportsFloat] | None = None,
probability: Callable[..., SupportsFloat] | None = None,
condition: Callable[..., object] | None = None,
matrix: Callable[..., Sequence[Sequence[object]]] | None = None,
emission: Callable[..., Sequence[Sequence[object]]] | None = None,
observed: str | None = None,
over: str = 'member',
neighborhood: BaseNeighborhood | None = None,
after: float | None = None,
schedule: str = 'synchronous',
form: str = 'symbolic',
phase: str | None = None,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.source = source
self.destination = destination
self.rate = rate
self.probability = probability
self.condition = condition
self.matrix = matrix
self.emission = emission
self.observed = observed
self.over = over
self.neighborhood = neighborhood
self.after = after
self.schedule = schedule
self.form = form
self.phase = phase
def _check_kind(self) -> str:
"""Raise an error for a transition without exactly one callable kind, or with an emission it cannot take."""
kinds = [kind for kind in KINDS if getattr(self, kind) is not None]
if len(kinds) != 1 or not callable(getattr(self, kinds[0])):
error_msg = (
f'The transition {self.label!r} has the forms {kinds}. Give exactly one of {", ".join(KINDS)}, as a '
'function of the namespace.'
)
raise ValidationError(error_msg)
unpaired = (self.emission is None) != (self.observed is None)
if unpaired or (self.emission is not None and self.matrix is None):
error_msg = (
f'The transition {self.label!r} has the emission {self.emission!r} of the state {self.observed!r}. '
'Give both with a matrix, or neither.'
)
raise ValidationError(error_msg)
return kinds[0]
def _check_form(self) -> None:
"""Raise an error for a transition without exactly one callable kind, or with arguments its kind ignores."""
kind = self._check_kind()
if self.over not in OVERS:
error_msg = f'The transition {self.label!r} runs over {self.over!r}. Use one of {", ".join(OVERS)}.'
raise ValidationError(error_msg)
if self.form not in FORMS:
error_msg = f'The form of transition {self.label!r} is {self.form!r}. Use one of {", ".join(FORMS)}.'
raise ValidationError(error_msg)
if self.schedule not in SCHEDULES or (self.form == 'symbolic' and self.schedule != 'synchronous'):
error_msg = (
f'The schedule of transition {self.label!r} is {self.schedule!r}. Use one of {", ".join(SCHEDULES)}, '
'and the agent form for another than `synchronous`.'
)
raise ValidationError(error_msg)
if self.over == 'cell' and (self.form != 'agent' or kind == 'matrix' or self.after is not None):
error_msg = (
f'The members of transition {self.label!r} act on their cell, which only the agent form does, by a '
'rate, a probability, or a condition and without a residence time.'
)
raise ValidationError(error_msg)
if self.over == 'neighbors' and (kind not in {'rate', 'probability'} or self.form == 'agent'):
error_msg = (
f'The transition {self.label!r} runs over neighbours, which only a rate or a probability of the '
'symbolic form does. In the agent form, read the neighbours of the agent in the function.'
)
raise ValidationError(error_msg)
single = not (callable(self.source) or isinstance(self.source, Collection))
if self.matrix is None and single and self.source == self.destination:
error_msg = f'The transition {self.label!r} moves from {self.source} to itself. Give distinct values.'
raise ValidationError(error_msg)
if self.after is not None and (kind not in {'probability', 'condition'} or not self.after >= 0):
error_msg = (
f'The transition {self.label!r} has the residence time {self.after!r}. Give a number of steps at least '
'0 with a probability or a condition.'
)
raise ValidationError(error_msg)
def _check_transition(self) -> None:
"""Raise an error for arguments the transition cannot take."""
self._check_form()
def _check(self) -> None:
"""Raise an error for a transition without exactly one kind, or with arguments its kind ignores.
Raises:
ValidationError:
The transition has not exactly one kind, the source and the destination are equal, an emission has
no matrix or no observed state, `over` is unknown or goes with a condition, a matrix, or the agent
form, the residence time is negative or goes with a rate or a matrix, the form is unknown, or a
function is not callable.
"""
super()._check()
self._check_transition()
def _get_moving(self, value: Expression, source: SourceValues) -> Expression:
"""Return 1 for a member whose state holds a source value, and 0 otherwise."""
if callable(source):
return where(SR(source(value)), 1, 0)
if isinstance(source, Collection):
return sum((where(value == item, 1, 0) for item in set(source)), SR(0))
return where(value == source, 1, 0)
def _get_jump(self, value: Expression, source: SourceValues, destination: float, rate: Expression) -> JumpLaw:
"""Return the jump of the state from a source value to the destination at a rate."""
if callable(source) or isinstance(source, Collection):
return JumpLaw(rate * self._get_moving(value, source), {value: destination - value})
return JumpLaw(rate * self._get_moving(value, source), {value: destination - source})
def _get_step(self, value: Expression, source: SourceValues, destination: float, moves: Expression) -> UpdateLaw:
"""Return the update of the state that moves a member from a source value where a condition holds."""
if callable(source) or isinstance(source, Collection):
return UpdateLaw(value, where(self._get_moving(value, source) > 0, where(moves, destination, value), value))
return UpdateLaw(value, where(value == source, where(moves, destination, source), value), source=source)
def _holds_source(self, value: float) -> bool:
"""Return whether the value of the state of an agent is a source value."""
if callable(self.source):
return bool(self.source(value))
if isinstance(self.source, Collection):
return value in set(self.source)
return value == self.source
def _combine_rates(self, rate: Expression) -> Expression:
"""Return the rate of the member, or the sum of the rates of its neighbours."""
return neighbors(rate, over=self.neighborhood) if self.over == 'neighbors' else rate
def _combine_probabilities(self, probability: Expression) -> Expression:
"""Return the probability of the member, or the probability that at least one neighbour moves it."""
if self.over == 'member':
return probability
return 1 - exp(neighbors(log(1 - probability), over=self.neighborhood))
def _get_draw(self) -> Expression:
"""Return the uniform draw of the transition in a step."""
return uniform(re.sub(r'\W', '_', f'transition_{self.label}'))
def _get_chain(self, namespace: SimpleNamespace, value: Expression) -> list[BasePart]:
"""Return the transition law of the state by its matrix, with the emission law of the observed state."""
laws: list[BasePart] = []
if self.matrix is not None:
laws.append(TransitionLaw(value, self.matrix(namespace), label=self.label))
if self.emission is not None and self.observed is not None:
observed = getattr(self._get_member(namespace, 'target'), self.observed)
laws.append(EmissionLaw(observed, value, self.emission(namespace), label=self.label))
return laws
def _get_held(self, agent: AgentViewProtocol, system: SystemViewProtocol) -> int:
"""Return the number of whole steps an agent has held its value, counted in the storage of the system."""
value = float(getattr(agent, self.state))
record = system.storage.get(f'{self.label}.{agent.id}')
held = record[1] + 1 if isinstance(record, tuple) and record[0] == value else 0
system.storage[f'{self.label}.{agent.id}'] = (value, held)
return held
def _draw_row(self, row: Sequence[object], system: SystemViewProtocol) -> float:
"""Return a value drawn by the probabilities of a row of a matrix."""
weights = np.asarray(row, dtype=float)
return float(system.rng.choice(weights.size, p=weights / weights.sum()))
def _evaluate[T](
self,
function: Callable[..., T],
agent: AgentViewProtocol,
other: AgentViewProtocol | None,
system: SystemViewProtocol,
) -> T:
"""Return a function of an agent, or of an agent and the member of its cell that acts on it."""
return function(agent, system) if other is None else function(agent, other, system)
def _moves_agent(
self,
agent: AgentViewProtocol,
system: SystemViewProtocol,
other: AgentViewProtocol | None = None,
) -> bool:
"""Return whether an agent that holds the source value moves in this step, moved by another one or alone."""
if self.condition is not None:
return bool(self._evaluate(self.condition, agent, other, system))
if self.rate is not None:
probability = 1.0 - math.exp(-float(self._evaluate(self.rate, agent, other, system)))
elif self.probability is not None:
probability = float(self._evaluate(self.probability, agent, other, system))
else:
probability = 0.0
return bool(system.rng.random() < probability)
def _get_next(self, agent: AgentViewProtocol, system: SystemViewProtocol) -> float | None:
"""Return the next value of the state of an agent, or `None` when it keeps its value."""
value = float(getattr(agent, self.state))
if self.matrix is not None:
return self._draw_row(self.matrix(agent, system)[int(value)], system)
held = self._get_held(agent, system) if self.after is not None else 0
if not self._holds_source(value) or (self.after is not None and held + 1 < self.after):
return None
return float(self.destination) if self._moves_agent(agent, system) else None
def _get_moves(self, agent: AgentViewProtocol, system: SystemViewProtocol) -> list[tuple[AgentViewProtocol, float]]:
"""Return the members an agent moves in this step, itself or the members of its cell, with their values."""
if self.over != 'cell':
value = self._get_next(agent, system)
return [] if value is None else [(agent, value)]
return [
(other, float(self.destination))
for other in agent.agents_here(agent.group)
if other is not agent
and self._holds_source(float(getattr(other, self.state)))
and self._moves_agent(other, system, agent)
]
def _apply(self, moves: list[tuple[AgentViewProtocol, float]], system: SystemViewProtocol) -> None:
"""Set the moved states, and draw the observed state of an emission."""
for agent, value in moves:
setattr(agent, self.state, value)
if self.emission is not None and self.observed is not None:
setattr(agent, self.observed, self._draw_row(self.emission(agent, system)[int(value)], system))
def _move_member(self, agent: AgentViewProtocol, system: SystemViewProtocol) -> None:
"""Move the state of one member, or of the members of its cell, from the states the members before it set."""
if agent.alive:
self._apply(self._get_moves(agent, system), system)
def move(self, system: SystemViewProtocol) -> None:
"""Move the state of every agent of the target that moves in this step, by the schedule.
Args:
system:
The system view of the interaction.
"""
agents = list(system.agents(self._get_role('target')))
if self.schedule == 'synchronous':
moves = [move for agent in agents for move in self._get_moves(agent, system)]
self._apply(moves, system)
return
order = system.rng.permutation(len(agents)) if self.schedule == 'random-sequential' else range(len(agents))
for index in order:
self._move_member(agents[int(index)], system)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the jump, the update, the transition or emission law of the state, or the rule of the agent form.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target`.
Returns:
laws:
The law of the move, or the rule that moves the agents.
"""
if self.form == 'agent' and self.schedule == 'member':
return [RuleLaw(self._move_member, phase=self.phase, members=self._get_role('target'))]
if self.form == 'agent':
return [RuleLaw(self.move, phase=self.phase)]
value = getattr(self._get_member(namespace, 'target'), self.state)
if self.matrix is not None:
return self._get_chain(namespace, value)
if self.rate is not None:
rate = self._combine_rates(SR(self.rate(namespace)))
return [self._get_jump(value, self.source, self.destination, rate)]
if self.probability is not None:
probability = self._combine_probabilities(SR(self.probability(namespace)))
moves = self._get_draw() < probability
else:
moves = SR(self.condition(namespace)) if self.condition is not None else SR(0)
if self.after is not None:
moves = both(residence(value) + 1 >= self.after, moves)
return [self._get_step(value, self.source, self.destination, moves)]
Equations
economy.K: K_economy(t + 1) == ((z_high_economy - z_low_economy)*regime_economy(t) + z_low_economy)*alpha_economy*beta_economy*K_economy(t)^alpha_economy
economy.regime: regime_economy(t + 1) == where(categorical__transition_regime_economy_productivity*(where(regime_economy(t) == 0, 0.9, 0.09999999999999998) + where(regime_economy(t) == 0, 0.09999999999999998, 0.9)) >= where(regime_economy(t) == 0, 0.9, 0.09999999999999998), 1, 0)