Skip to content

Traders and exchange rates coupled through arbitrage as lattice electrodynamics

Exchange-rate potentials on links as lattice Maxwell fields, sourced by traders who move a good on sites.

  • Name: gauge_arbitrage_market
  • Field: economics/finance
  • Formalism: ordinary differential equations
  • References: Schwichtenberg 2020

Select this system

from skcomplex.systems import create_simulation, create_system, describe

system = create_system('gauge_arbitrage_market')
entry = describe('gauge_arbitrage_market')
trajectory = create_simulation('gauge_arbitrage_market').prepare().run().trajectory_

Factory

@register_system(
    field='economics/finance',
    tags=('gauge', 'arbitrage', 'exchange rates', 'lattice', 'maxwell', 'cell complex'),
    references=('Schwichtenberg 2020',),
    defaults={'until': 20.0, 'dt': 0.5, 'solver': RK4Solver(dt=0.02)},
    observables=_market_observables,
)
def gauge_arbitrage_market(
    size: int = 16,
    *,
    D: float = 0.5,
    q: float = 0.2,
    bump: float = 1.0,
    seed: int = 0,
) -> System:
    """Return exchange-rate potentials on links as lattice Maxwell fields, sourced by traders who move a good on sites.

    The potentials follow `A_t = E` and `E_t = -codifferential(dA) - J`, and the good follows `rho_t =
    -codifferential(J)` with the current `J = -D (d rho - q A rho)` along its covariant gradient, after Schwichtenberg
    (2020 §2.1 and §2.2.1). Here `E` is `F_0i` of the book, so `codifferential(E)` equals the excess of the good over
    its mean, as in its Gauss law, Eq. 2.37, with a uniform background of the mean density that keeps the periodic
    lattice neutral. The system departs from the book in two places. The book writes Eq. 2.39 with no metric, which
    gives `A_tt = +codifferential(dA) - J`, where every mode with a curl grows, so the system takes the sign of a wave.
    And Eq. 2.56 applies the covariant derivative twice, which does not conserve the good, so the system applies it once
    inside the divergence of the current. The book gives no values of `D` and `q` and sets `mu_0` by no number, so `D`,
    `q`, and `mu_0 = 1` are choices of this system.

    Args:
        size:
            The number of sites along both axes.

        D:
            The diffusivity of the good.

        q:
            The charge of the good.

        bump:
            The height of the random part of the initial density above a density of 1.

        seed:
            The seed of the initial density.

    Returns:
        system:
            The system.

    Examples:
        >>> from skcomplex.systems.economics.finance import gauge_arbitrage_market
        >>> [element.label for element in gauge_arbitrage_market(4).elements]
        ['sites', 'links']
    """
    space = CellComplexSpace(LatticeSpace((size, size))).prepare()
    density = 1.0 + bump * np.random.default_rng(seed).random(size * size)
    field = _initial_field(space.incidences_[0], density - density.mean())
    sites = GroupElement('sites', element=GoodSiteElement, space=space, on='sites', values={'rho': density})
    links = GroupElement(
        'links',
        element=ExchangeLinkElement,
        space=space,
        on='edges',
        values={'E': field},
        visual=VisualSpec('cell', {'state': 'A'}),
    )
    interactions = [
        RateInteraction('potential', target='links', state='A', rate=_potential_rate),
        RateInteraction('field', target='links', state='E', rate=_field_rate),
        TraderCurrentInteraction('current', target='sites', links='links', D=D, q=q),
    ]
    return System(elements=[sites, links], interactions=interactions)

Components

Elements:

  • GoodSiteElement, from Element
  • ExchangeLinkElement, from Element

Interactions:

  • RateInteraction
  • TraderCurrentInteraction, from NeighborInteraction
class GoodSiteElement(Element):
    """Site of a lattice with a density `rho` of a good.

    Args:
        label:
            The label of the site or the group on the sites.

        visual:
            How the site is drawn, or `None` for `default_visual_`.

        rho:
            The initial density of the good.

    Attributes:
        default_visual_:
            The visual of the site when it sets none.

    Examples:
        >>> from skcomplex.systems.economics.finance import GoodSiteElement
        >>> GoodSiteElement('site', rho=2.0).rho
        2.0
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        rho: Annotated[float, StateSpec()] = 1.0,
    ) -> None:
        super().__init__(label, visual=visual)
        self.rho = rho

class ExchangeLinkElement(Element):
    """Link of a lattice with an exchange-rate potential `A` and its rate of change `E`.

    Args:
        label:
            The label of the link or the group on the links.

        visual:
            How the link is drawn, or `None` for `default_visual_`.

        A:
            The initial exchange-rate potential.

        E:
            The initial value of `F_0i` of Schwichtenberg (2020 Eq. 2.22), the rate of change of the potential.

    Attributes:
        default_visual_:
            The visual of the link when it sets none.

    Examples:
        >>> from skcomplex.systems.economics.finance import ExchangeLinkElement
        >>> ExchangeLinkElement('link', A=0.5).A
        0.5
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        A: Annotated[float, StateSpec()] = 0.0,
        E: Annotated[float, StateSpec()] = 0.0,
    ) -> None:
        super().__init__(label, visual=visual)
        self.A = A
        self.E = E

class RateInteraction(Interaction):
    """Rate `rate(ns)` and noise `noise(ns) dW` added to the highest time derivative of a state of a target.

    On a state of order 1 the rate is its first derivative, and on a state of order `k` it adds to the `k`-th
    derivative, so rates on one state superpose. With a divisor `d` the rate and the noise are divided by `d`, the name
    of a state or parameter of the target, such as an inertia, a friction, or a capacitance, which backends can read.

    A rate, a noise, or a divisor function takes the namespace of the interaction, with the namespaces `target` and
    `sources` of the bound elements, the time `t`, and the symbols of the parameters a subclass declares. A number
    stands for a constant, and a number on a vector state applies to every component. The noise is a Wiener process per
    component and member by default, a process of a `WienerProcessSpec` to share or correlate it across states and
    members, or a fractional Brownian motion of a Hurst index.

    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 rate changes.

        state:
            The name of the state of the target that the rate changes.

        rate:
            A function of the namespace that returns the rate, or a number.

        sources:
            The paths of the elements or groups the rate and the noise read but do not change.

        noise:
            A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.

        divisor:
            The name of the state or parameter of the target that divides the rate and the noise, a function of the
            namespace of the target that returns the divisor, such as a moment of inertia `m l^2`, or `None`.

        process:
            The Wiener processes that drive the noise, or `None` for one independent process per component and member.

        index:
            The process of `process` that drives the noise, or `None` for every process, one per component of a vector
            state.

        convention:
            The stochastic calculus of the noise, `'ito'`, `'stratonovich'`, or `'anti-ito'`.

        hurst:
            The Hurst index in `(0, 1)` of a fractional Brownian motion that drives the noise, a number or a function
            of the namespace that returns a parameter, or `None` for a Wiener process.

        conserved:
            Whether the noise is the divergence of a random flux on the links of the lattice of a field, so it keeps the
            total of the state.

    Raises:
        ValidationError:
            The state is not a name, a function is neither callable nor a number, the convention is unknown, the Hurst
            index is outside `(0, 1)`, or the noise options cannot combine.

    Examples:
        >>> from skcomplex.core.components.elements import CompartmentElement
        >>> from skcomplex.core.components.interactions import RateInteraction
        >>> from skcomplex.core.components.system import System
        >>> from skcomplex.core.engine.runtime import Simulation
        >>> decay = RateInteraction('decay', target='stock', state='n', rate=lambda s: -0.5 * s.target.n)
        >>> system = System(elements=[CompartmentElement('stock', n=8.0)], interactions=[decay])
        >>> trajectory = Simulation(system, dt=0.01).prepare().run(until=2.0).trajectory_
        >>> round(float(trajectory.values('stock.n')[-1]), 3)
        2.943
    """

    law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (RateLaw, NoiseLaw, FractionalNoiseLaw)

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('state', 'divisor'))],
        state: str = 'x',
        rate: Callable[[SimpleNamespace], object] | float = 0.0,
        sources: Annotated[Sequence[str], RoleSpec(variadic=True, elements=(ElementKind.ANY,))] = (),
        noise: Callable[[SimpleNamespace], object] | float | None = None,
        divisor: str | Callable[[SimpleNamespace], object] | None = None,
        process: Annotated[WienerProcessSpec | None, RoleSpec(noise=True)] = None,
        index: int | None = None,
        convention: str = 'ito',
        hurst: Callable[[SimpleNamespace], object] | float | None = None,
        conserved: bool = False,
    ) -> None:
        super().__init__(label)
        self.target = target
        self.state = state
        self.rate = rate
        self.sources = sources
        self.noise = noise
        self.divisor = divisor
        self.process = process
        self.index = index
        self.convention = convention
        self.hurst = hurst
        self.conserved = conserved

    def _check_functions(self) -> None:
        """Raise an error for a state that is not a name, or a function that is neither callable nor a number."""
        if not isinstance(self.state, str) or not self.state.isidentifier():
            error_msg = (
                f'{type(self).__name__} {self.label!r} changes the state {self.state!r}, which is not a name of a '
                'state.'
            )
            raise ValidationError(error_msg)
        for name in ('rate', 'noise'):
            value = getattr(self, name)
            if not (value is None or callable(value) or isinstance(value, int | float)):
                error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace or a number.'
                raise ValidationError(error_msg)
        if not (
            self.divisor is None
            or callable(self.divisor)
            or (isinstance(self.divisor, str) and self.divisor.isidentifier())
        ):
            error_msg = (
                f'The divisor of {self.label!r} is {self.divisor!r}. Give the name of a state or parameter, or a '
                'function of the namespace of the target.'
            )
            raise ValidationError(error_msg)

    def _check_noise(self) -> None:
        """Raise an error for an unknown convention, a Hurst index outside `(0, 1)`, or noise options that conflict."""
        if self.convention not in CONVENTIONS:
            error_msg = f'Unknown noise convention {self.convention!r}. Use one of {", ".join(CONVENTIONS)}.'
            raise ValidationError(error_msg)
        if not (
            self.hurst is None or callable(self.hurst) or (isinstance(self.hurst, int | float) and 0 < self.hurst < 1)
        ):
            error_msg = f'The Hurst index of {self.label!r} is {self.hurst!r}. Give a number in (0, 1) or a function.'
            raise ValidationError(error_msg)
        if self.index is not None and self.process is None:
            error_msg = f'{self.label!r} picks the process {self.index!r} but binds no processes. Give `process`.'
            raise ValidationError(error_msg)
        options = [
            name
            for name, given in (('hurst', self.hurst), ('process', self.process), ('conserved', self.conserved))
            if given
        ]
        if len(options) > 1:
            error_msg = f'The noise of {self.label!r} takes {" and ".join(options)}. Give at most one of them.'
            raise ValidationError(error_msg)

    def _check(self) -> None:
        """Raise an error for a state, functions, or noise options the interaction cannot take.

        Raises:
            ValidationError:
                The type declares states, the state is not a name, a function is neither callable nor a number, the
                convention is unknown, the Hurst index is outside `(0, 1)`, or the noise options cannot combine.
        """
        super()._check()
        self._check_functions()
        self._check_noise()

    def _evaluate(self, function: Callable[[SimpleNamespace], object] | float, namespace: SimpleNamespace) -> object:
        """Return the value of a function of the namespace, or a number as it is."""
        return function(namespace) if callable(function) else function

    def _broadcast(self, value: object, state: object) -> object:
        """Return a value with one component per component of a vector state, repeating a number."""
        if isinstance(state, FreeModuleElement) and not isinstance(value, FreeModuleElement | list | tuple):
            return vector(SR, [value] * len(state))
        if isinstance(value, list | tuple):
            return vector(SR, value)
        return value

    def _divide(self, value: object, namespace: SimpleNamespace) -> object:
        """Return a value divided by the divisor of the target, or the value without a divisor."""
        if self.divisor is None:
            return value
        divisor = (
            self.divisor(self._get_member(namespace, 'target'))
            if callable(self.divisor)
            else getattr(self._get_member(namespace, 'target'), self.divisor)
        )
        return operator.truediv(value, divisor)

    def _get_rate(self, namespace: SimpleNamespace) -> object:
        """Return the rate before the divisor, the hook a subclass overrides to derive it from its parameters."""
        return self._evaluate(self.rate, namespace)

    def _get_noise(self, namespace: SimpleNamespace) -> object:
        """Return the noise amplitude before the divisor, the hook a subclass overrides, or `None` for no noise."""
        return None if self.noise is None else self._evaluate(self.noise, namespace)

    def _get_noise_law(self, state: object, amplitude: object, namespace: SimpleNamespace) -> BasePart:
        """Return the noise law of an amplitude, by the process, the convention, and the Hurst index."""
        if self.hurst is not None:
            hurst = self._evaluate(self.hurst, namespace)
            return FractionalNoiseLaw(state, amplitude, hurst=hurst, convention=self.convention)
        process = None
        if self.process is not None:
            handle = self._get_member(namespace, 'process').dW
            process = handle if self.index is None else handle[self.index]
        return NoiseLaw(state, amplitude, convention=self.convention, process=process, conserved=self.conserved)

    def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
        """Return the rate law of the state and its noise law, both divided by the divisor.

        Args:
            namespace:
                The namespace of the symbols the laws can use, with the namespaces of the target and the sources.

        Returns:
            laws:
                The rate law, and the noise law when the interaction has noise.
        """
        state = getattr(self._get_member(namespace, 'target'), self.state)
        rate = self._divide(self._broadcast(self._get_rate(namespace), state), namespace)
        laws: list[BasePart] = [RateLaw(state, rate)]
        noise = self._get_noise(namespace)
        if noise is not None:
            laws.append(self._get_noise_law(state, self._divide(noise, namespace), namespace))
        return laws

class TraderCurrentInteraction(NeighborInteraction):
    """Trader current `J = -D (rho_t - rho_s - q A (rho_s + rho_t) / 2)` along the links, down the covariant gradient.

    The covariant derivative is `D_mu = d_mu - A_mu` of Schwichtenberg (2020 Eq. 2.55), with the charge `q` of its
    Eq. 4.24, and the current is minus the diffusion constant times this derivative of the good, as in Eq. 2.56. The
    current moves the good from the source of a link to its target (footnote 18 of §2.1), and it is the source of the
    field of the link with the sign of the Ampere-Maxwell law of Eq. 2.39, `E' = ... - J`, so the Gauss law holds as the
    good moves. These signs give the potential a mass, `A'' = ... - D q rho A`, for a positive density and charge, so
    the fields stay bounded.

    Args:
        label:
            The label of the interaction, unique among the parts of its system.

        target:
            The path of the group on the sites.

        links:
            The path of the group on the edges.

        D:
            The diffusivity of the good.

        q:
            The charge of the good.

    Examples:
        >>> from skcomplex.systems.economics.finance import TraderCurrentInteraction
        >>> TraderCurrentInteraction('current', target='sites', links='links', q=0.5).q
        0.5
    """

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        links: Annotated[str, RoleSpec()],
        D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
        q: Annotated[float, ParameterSpec()] = 0.2,
    ) -> None:
        super().__init__(label)
        self.target = target
        self.links = links
        self.D = D
        self.q = q

    def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
        """Return the rates of the good at both ends of a link and of the electric field of the link.

        Args:
            namespace:
                The namespace of an edge, with the ends `source` and `target` and the member `edge`.

        Returns:
            laws:
                The rates.
        """
        source, target, edge = namespace.source, namespace.target, namespace.edge
        middle = (source.rho + target.rho) / 2
        current = -namespace.D * (target.rho - source.rho - namespace.q * edge.A * middle)
        return [RateLaw(source.rho, -current), RateLaw(target.rho, current), RateLaw(edge.E, -current)]

Equations

links.A[i]: diff(A_links(t), t) == E_links