2021

Part of this work was done during an internship at ESA.

[2]Emmanuel Blazquez

1]Department of Computing Science, University of Alberta, Edmonton, Canada

2]Advanced Concepts Team, European Space Agency (ESA), Noordwijk, Netherlands

# Deep Learning of Dynamical System Parameters from Return Maps as Images

Connor James Stephens<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [

###### Abstract

We present a novel approach to system identification (SI) using deep learning
techniques. Focusing on parametric system identification (PSI), we use a
supervised learning approach for estimating the parameters of discrete and
continuous-time dynamical systems, irrespective of chaos. To accomplish this,
we transform collections of state-space trajectory observations into image-
like data to retain the state-space topology of trajectories from dynamical
systems and train convolutional neural networks to estimate the parameters of
dynamical systems from these images. We demonstrate that our approach can
learn parameter estimation functions for various dynamical systems, and by
using training-time data augmentation, we are able to learn estimation
functions whose parameter estimates are robust to changes in the sample
fidelity of their inputs. Once trained, these estimation models return
parameter estimations for new systems with negligible time and computation
costs.

###### keywords:

dynamical systems, system identification, machine learning, deep learning,
chaos

##  1 Introduction

The natural sciences are experiencing a watershed moment with the recent
success of sophisticated, data-driven methods for approaching problems that
have resisted traditional analytical and optimization-based methods such as
AlphaFold for protein structure prediction Jumper2021 , graph-based neural
network models for weather forecasting Lam2022  and the control of tokamak
plasmas using deep reinforcement learning Degrave2022 . These recent successes
have been driven by innovations in large-scale data-driven modeling and rapid
advancements in specialized computer hardware that have propelled the rise of
deep learning methods alongside domain-specific expertise and techniques from
scientific sub-fields. In line with this trend, we present a novel method for
solving parametric system identification (PSI) by learning a regression
function from return maps to system parameters.

System identification (SI) is the process of modeling and analyzing the
behavior of dynamical systems based on observed input and output
dataSoderstrom1989 . Accurately characterizing and even predicting the
behavior of a dynamical system is an invaluable tool in various areas spanning
engineering, physics, biology, and economics. It is particularly valuable in
control applications where accurate modeling of a system can result in
dramatic performance and robustness improvements Alvin2003 ; Bruggeman2007 ;
Ljung2010 ; Villaverde2014 . This work focuses on parametric system
identification, which is concerned with identifying the underlying parameters
of the target system given some parametrized class of dynamical systems
believed to contain, or at least closely approximate, the target system
Quaranta2020 . To understand how and why we solve PSI as a regression problem,
it is helpful to place this approach in opposition to the more traditional
setting of PSI as an optimization problem to be solved by means of meta-
heuristic algorithms such as genetic and evolutionary algorithmsKonnur2005 ;
Tao2007 ; Wang2011 . Given a parametric class of dynamical systems and data
from a target system, in the optimization setting one has a way to compute
simulated trajectories given parameter estimates, and the goal is to find
parameter estimates which match those of the target system. To accomplish
this, one typically combines a loss function which acts as a ’dissimilarity’
score between trajectories generated using candidate system parameters and
observed trajectories. An optimization scheme is then used to select new
candidate parameters for the system. This process is iterated until some
computational budget has been exhausted or the loss converges. The observed
losses of parameter estimates are then used to select a final estimate, e.g.
the parameter estimate with the lowest loss observed in the optimization
process.

One of the challenges of this approach is choosing the loss function that is
used to evaluate the quality of a given set of parameter estimates for this
solution method. Given equal-length simulated and observed trajectories, one
of the most common approaches is to minimize the average error between the
pairs of simulated and observed trajectory points Konnur2005 ; Tao2007 ;
Ho2010 ; Modares2010 ; Wang2011 . This can cause issues when dealing with
chaotic systems, which have the characteristic property that trajectories that
are initially neighboring in state space can diverge exponentially over time.
This means that such a loss function can be highly sensitive to small
measurement errors in the initial conditions of the target system. To address
this issue, recent work has developed loss functions that compare trajectories
as signals at the state space, as opposed to the time-domain level. where the
trajectories of chaotic systems tend to take a more structured form Jafari2014
; Mousazadeh2020 . These works have proposed several heuristic loss functions
for comparing the similarity of trajectories in state space Jafari2014 ;
Shekofteh2015 ; Shekofteh2019 ; Shekofteh2019a ; Mousazadeh2020 , however,
there is currently no clear choice as to which state-space loss function is
best suited for parameter estimation.

Even working in state space, current methods for PSI solve an optimization
problem using an iterative procedure that involves repeatedly simulating
trajectories using system parameter estimates, which for continuous-time
systems can be a time and computation-consuming process. Worse still, for
simple discrete maps, it is often feasible to compute the gradients of the
loss function with respect to the system parameters Mousazadeh2022 , however,
a numerical integrator is typically required for continuous-time dynamical
systems. This makes evaluating or estimating gradients significantly more
challenging. Due in part to this, most previous work has viewed their choice
of loss function as a non-differentiable function of the candidate parameters
and used various zeroth-order optimization methods such as particle-based
methods and genetic algorithms. See Quaranta2020  for a review of methods
based on these computational intelligence approaches. These potentially
expensive iterations of simulation and optimization limit the applicability of
these methods for real-time settings or processing large collections of
systems.

This work explores an alternative approach to PSI for chaotic systems, framing
it as a supervised learning problem of mapping state-space data to system
parameters. We were initially motivated by the observation that after seeing
sufficiently many examples of return maps, a human expert can begin to
identify patterns that hint about the parameters of the underlying dynamical
systems. Following this intuition, we transform collections of state-space
trajectory observations into image-like data to retain the state-space
topology of trajectories from dynamical systems and use these collections to
train convolutional neural networks to estimate the parameters of dynamical
systems from these images. In our approach, the regression loss is the
parameter estimation error, avoiding the use of heuristic loss functions.
Additionally, with almost all of the time and computational cost occurring
during data collection and model training, our estimation models return
parameter estimations for new samples with negligible time and computation
costs at inference time, neatly complementing the drawbacks of optimization-
based approaches.

In the remainder of this paper, we formalize the parametric system
identification problem as it has conventionally been approached and follow
this by presenting our supervised learning problem in Section 2. We then
detail our solution method for this new problem formulation, including
generating trajectory datasets, extracting useful features from them, and
using them to train a parameter estimation model in Section 3. We present the
results of applying our method to both discrete- and continuous-time dynamical
systems in Section 4 and discuss our findings in Section 5. We conclude by
discussing the opportunities and challenges which are introduced with this new
approach to analyzing dynamical systems in Section 6.

##  2 Background

###  2.1 Dynamical Systems

We consider the problem of estimating the parameters of a
$p^{\text{th}}$-order autonomous dynamical system, defined on some state-space
$\mathcal{X}\subseteq\mathbb{R}^{p}$ with $p\geq 1$ by parametric dynamics
$F(\cdot~{}\mathchar
24635\relax\;\mathbf{\mathbf{\theta}}):\mathbb{R}^{p}\to\mathbb{R}^{p}$ with
$\mathbf{\mathbf{\theta}}\in\Theta\subseteq\mathbb{R}^{d}$ for some $d\geq 1$
by the state equation

$\begin{cases}\dot{\mathbf{x}}(t)=F(\mathbf{x}(t)\mathchar
24635\relax\;\mathbf{\mathbf{\theta}}),t\in\mathbb{R}&\text{for a continuous
flow},\\\ \mathbf{x}_{k+1}=F(\mathbf{x}_{k}\mathchar
24635\relax\;\mathbf{\mathbf{\theta}}),k\in\mathbb{Z}&\text{for a discrete
map}.\end{cases}$ (1)

Dynamical systems frequently appear in the study of classical mechanics, with
the Hamiltonian of some physical system defining the dynamics $F$ on some
configuration manifold. For our purposes, the dynamics $F$ can be any function
of the state $\mathbf{x}$ and the parameters $\mathbf{\mathbf{\theta}}$ which
characterizes the evolution of the state of the system in time according to
equation 1.

In this section, we restrict our discussion of parameter estimation to
discrete dynamical systems since our approach for handling continuous-time
systems is to use Poincaré sections and their respective Poincaré maps to
represent them with discrete-time systems for analysisParker1991 .

####  2.1.1 Poincaré Maps

A Poincaré map is defined by a dynamical flow along with an oriented
hypersurface in state space, referred to as a Poincaré section. Some
references remove the requirement that the section is oriented, resulting in a
so-called ‘two-sided’ Poincaré map. This choice does not have significant
consequences for our purposes. A Poincaré map is the unique map associated
with the original dynamical system, which takes a point $P$ in the Poincaré
section to the next point $P^{\prime}$ in the Poincaré section where the flow
of the original dynamical system, originating from $P$ crosses the oriented
section in the same direction as at $P$. For this reason, Poincaré maps are
sometimes called ‘first recurrence’, or ‘return maps’. In this work, when we
refer to return maps, we refer either to a collection of state-space
observations from a discrete map or from the discrete map induced by a
continuous time dynamical system along with a particular choice of Poincaré
section. Poincaré maps allow us to present a unified method for parametric
system identification on both discrete maps and continuous time dynamical
systems.

###  2.2 Parametric System Estimation

In the parametric system identification problem, we are given the parametric
form of $F$ and are tasked with constructing an estimate
$\hat{\mathbf{\mathbf{\theta}}}^{*}$ on the basis of a collection of $n$
observed trajectories
$\mathbf{T}=(\mathbf{\tau}_{1},\mathbf{\tau}_{2},\ldots,\mathbf{\tau}_{n})$
where each trajectory
$\tau_{i}=(\mathbf{x}^{i}_{1},\mathbf{x}^{i}_{2},\ldots)$ is an indexed
collection of elements from $\mathcal{X}$ which follows
$\mathbf{x}^{i}_{k+1}=F(\mathbf{x}^{i}_{k}\mathchar
24635\relax\;\mathbf{\mathbf{\theta}}^{*})$. Most previous work on PSI
considers the case when $n=1$, but is straightforward to generalize for the
case of multiple sample trajectories.

The optimization approach to estimating $\mathbf{\mathbf{\theta}}^{*}$ from a
collection of observed trajectories of the target system, $\mathbf{T}$ is by
solving an optimization problem of the form

$\hat{\mathbf{\mathbf{\theta}}}^{*}=\operatorname*{\text{arg
min}}_{\mathbf{\mathbf{\theta}}\in\Theta}~{}\ell(\mathbf{\mathbf{\theta}},\mathbf{T}),$
(2)

where $\ell:\Theta\times\mathcal{X}^{n\times m}\rightarrow\mathbb{R}^{+}$ is a
loss function that maps parameters estimates to positive real numbers, given
observed trajectories $\mathbf{T}$, and where we assume for simplicity that
each of the $n$ observed trajectories in $\mathbf{T}$ consist of $m$ points.
The value of $\mathbf{\mathbf{\theta}}^{*}$ influences $\mathbf{T}$ through
the dynamics $F(\cdot\mathchar 24635\relax\;\theta^{*})$ imposed on the
trajectories. For this reason, previous work has typically defined $\ell$ by
composing a loss function $f$ that compares sets of trajectories with a
numerical solver that is used to produce simulated trajectories corresponding
to a parameter estimate. Concretely, to evaluate the loss of a parameter
estimate $\hat{\mathbf{\theta}}^{*}$ against the true parameter
$\mathbf{\theta}^{*}$, $\ell(\hat{\mathbf{\theta}}^{*},\mathbf{\theta}^{*})$
practitioners use a numerical solver to create a collection of trajectories
$\mathbf{T}(\hat{\mathbf{\theta}}^{*})=(\hat{\mathbf{\tau}}_{1},\hat{\mathbf{\tau}}_{1},\ldots,\hat{\mathbf{\tau}}_{n})$
with
$\mathbf{\tau}_{i}=(\hat{\mathbf{x}}^{i}_{1},\hat{\mathbf{x}}^{i}_{2},\ldots,\hat{\mathbf{x}}^{i}_{m}),\hat{\mathbf{x}}^{i}_{k}\in\mathcal{X}$
and $\hat{\mathbf{x}}^{i}_{k+1}=F(\hat{\mathbf{x}}^{i}_{k}\mathchar
24635\relax\;\hat{\mathbf{\theta}}^{*})$. They then evaluate a loss that
assigns values to pairs of collections of trajectories
$f(\mathbf{T}^{\prime},\mathbf{T})$. This composition of a numerical solver
and loss function over sets of trajectories defines an implicit loss function
on $\Theta$ given a collection of observed trajectories $\mathbf{T}$,

$\ell(\hat{\mathbf{\theta}},\mathbf{T})\coloneqq
f(\mathbf{T}(\hat{\mathbf{\theta}}),\mathbf{T}).$ (3)

The method proposed in Mousazadeh2022  does not follow this description,
though as mentioned earlier, this method is not applicable when numerical
methods are required to approximate trajectories, e.g., when they correspond
to a Poincaré map.

####  2.2.1 State-space Representations of Trajectories

Parameter estimation is especially challenging when dealing with dynamical
systems which exhibit chaotic behavior. This is because of the characteristic
sensitivity of these systems to their initial conditions. Even if the exact
parameters of a chaotic system are known, any measurement error in the
system’s initial conditions can result in a simulated trajectory that diverges
exponentially from the target observation set over time. Due to these issues,
choices of $f(\hat{\mathbf{T}}(\hat{\mathbf{\theta}}),\mathbf{T})$ which are
constructed around pairwise differences of $\hat{\mathbf{x}}^{i}_{k}$ and
$\mathbf{x}^{i}_{k}$ in the time-domain are often poorly-conditioned with
respect to $\hat{\mathbf{\theta}}$ as well as any measurement error in the
initial conditions
$\mathbf{x}^{1}_{1},\mathbf{x}^{2}_{1},\ldots,\mathbf{x}^{m}_{1}$ of the
target system. A choice of $f$ which appears frequently in the literature is
the temporal mean squared error (MSE)

$f_{\text{MSE}}(\hat{\mathbf{T}}(\hat{\mathbf{\theta}}),\mathbf{T})\coloneqq\dfrac{1}{mn}\sum_{i=1}^{n}\sum_{k=1}^{m}\|\hat{\mathbf{x}}^{i}_{k}-\mathbf{x}^{i}_{k}\|^{2}.$
(4)

To address these issues, recent work such as Jafari2014 ; Mousazadeh2020  has
investigated the use of state-space representations of trajectories to compare
simulated trajectories to the target data. The basic idea of this method is
outlined by Jafari et al. Jafari2014 . While chaotic systems appear to have
highly disordered behavior when considering trajectories in the time domain,
they are known to have a more structured topology in state-space, which for
our purposes refers to considering trajectories as unordered sets of points in
state-space as opposed to time-indexed points. Chaotic attractors manifest
this phenomenon, wherein collections of dynamical system trajectories that
appear stochastic in the time domain are constrained to a low-dimensional
manifold in state space. This property of chaotic systems makes using state-
space representations of trajectories from dynamical systems appealing for
designing the function $f$ in the right-hand side of equation 3.

The main contribution of Jafari2014  was the construction of a loss function
between return maps. Their loss function considers
$\mathbf{T}(\hat{\mathbf{\theta}})$ and $\mathbf{T}$ as structureless sets of
points in $\mathcal{X}$, and measures the average of the Euclidean distances
between each point in the simulated data with those of its nearest neighbor in
the target data, and vice versa.

###  2.3 Supervised Machine Learning for Parametric System Identification

So far, we have only discussed the optimization approach to parametric system
identification (equation 2). Our novel solution for PSI instead frames
parameter estimation as a supervised machine-learning problem (see e.g.,
Murphy2012 ). In our approach, we aim to identify a single function that maps
from observations of trajectories from a discrete dynamical system to an
estimate of the system’s parameter $g:\mathcal{X}^{n\times
m}\rightarrow\mathbb{R}^{d}$, where $g\in\mathcal{F}$ is an element of some
function class $\mathcal{F}$. We select $g$ through a data-driven optimization
process that attempts to minimize the parameter estimation error of $g$,
assuming access to a large dataset of input-output pairs
$\mathcal{D}=(\mathbf{T}_{i},\mathbf{\theta}_{i})_{i=1}^{N}$ with a wide range
of values $\theta_{i}\in\Theta$. The $i^{\text{th}}$ pair in this set consists
of a collection of trajectories $\mathbf{T}_{i}$ along with the parameter,
$\mathbf{\theta}_{i}$, which generated them. Concretely, given a parameterized
class of dynamical systems, we use real or simulated observations from systems
with known parameters to construct a dataset, $\mathcal{D}$. We then use
$\mathcal{D}$ to estimate an optimal parameter estimation function, $g^{*}$,
defined to minimize the mean squared prediction error, or another suitable
choice of loss function, over some joint probability distribution $\mathbb{P}$
of possible values of $\mathbf{\theta}^{*}$ and the resulting observations
$\mathbf{T}(\mathbf{\theta}^{*})$:

$g^{*}=\operatorname*{\text{arg
min}}_{g\in\mathcal{F}}~{}\mathbb{E}_{(\mathbf{\theta},\mathbf{T})\sim\mathbb{P}}\left[\|\mathbf{\theta}-g(\mathbf{T}(\mathbf{\theta}))\|^{2}\right].$
(5)

For example, $\mathbb{P}$ could be the joint distribution corresponding to
$\theta$ drawn uniformly over a bounded subset of $\mathbb{R}^{d}$ with
trajectories $\mathbf{T}(\mathbf{\theta})$ generated from $\mathbf{\theta}$ in
some consistent but possibly random fashion, e.g., applying a deterministic
numerical integration software to randomly sampled initial states.

Given an estimate of the minimizer of equation 5, $\hat{g}^{*}$, and a new
collection of trajectories $\mathbf{T}(\mathbf{\theta}^{\prime})$
corresponding to parameter $\mathbf{\theta}^{\prime}$, we take
$\hat{g}^{*}(\mathbf{T}(\mathbf{\theta}^{\prime}))$ to be our estimate
$\hat{\mathbf{\theta}}^{\prime}$. We have purposefully left the form of the
function class $\mathcal{F}$ vague in this section. These details are
clarified in Section 2.3 when we outline our supervised machine learning
approach.

##  3 Parameter Estimation via Deep Learning on Return Maps

Estimating system parameters directly from observed trajectories is a
challenging problem for the same reasons outlined in Section 2.2.1. Drawing
from the lessons of earlier work on PSI for chaotic systems, rather than
working with trajectories as time-series data, we restrict our attention to
functions $g$ which only consider the state-space information of the
trajectory collections $\mathbf{T}_{i}$ that are passed as input. Further, we
transform this data into a form that allows us to make use of existing CNN
architectures to construct a parameter estimation function. Specifically, we
discretize regions of state space into pixels and use this discretization
scheme to transform return maps into coarse, single-channel images, and use
these as the input ‘features’ for our models.

Figure 1: Return maps generated by sampling states from the uniform
distribution over the region $x,y\in[-2,2]\times[-2,2]$ and plotting 250
iterations of the Hénon Map (equation 7) from each initial state. The
different colors correspond to trajectories for different initial states.

Figure 1 conveys the central intuition for our work. Different parameter
values for the Hénon map result in visually distinct return maps. This
observation suggests that it might be possible to learn a regression function
from these images to parameter estimates. Our experiments show that it is
possible to learn such a regression function from data.

The main contribution of this work is the demonstration that our method of
predicting system parameters from images of return maps is viable, introducing
the possibility of further applications of deep learning on return maps. For
simplicity, we consider dynamical systems with two-dimensional return maps in
this work. However, conceptually, our method generalizes to higher dimensions
at the cost of the requirement of additional computation and data requirements
commensurate with the curse of dimensionality Murphy2012 .

###  3.1 Dataset Generation

Developing a supervised machine learning model for predicting a system’s
parameters necessitates a dataset to train the model. Dataset generation is
straightforward after specifying a parameterized dynamical system and the
range of parameter values of interest. To ensure that the trained model can
perform well on the entire range of parameter values, we first selected a
large collection of evenly spaced parameter values. For the $i^{\text{th}}$
parameter value $\mathbf{\theta}_{i}$, we sampled $n$ initial states
conditions and iterated the dynamical map $m-1$ times in the discrete case, or
made use of Heyoka numerical integration software Biscani2021  to evaluate a
total of $m$ crossing points for a specified crossing section. The result of
this generation process was a collection of input-output pairs,
$\mathbf{T}_{i}$ where the ‘input’ is $\mathbf{\theta}_{i}$, and the ‘output’
is a collection of discrete trajectories corresponding to the
$N_{\text{traj}}$ trajectories of length $m$. In the case of dynamical flows,
it can be simpler to propagate the system forward until a fixed system time.
This can lead to the trajectories starting from the different initial states
having different numbers of crossings through the section. This detail does
not introduce any issues for our method, but it is a consideration in
practice.

###  3.2 Feature Processing

Figure 2: Transforming subsampled return maps into images. Each row shows the
process of processing a return map for input to the model. The left column
shows return maps sampled in the same fashion as in Figure 1. The middle
column shows the output of our random data augmentation scheme, i.e. taking a
random subset of the trajectories and truncating them to a random length. The
right-most column shows examples of the single-channel images fed input to the
model. The images are obtained by discretizing a region of state-space, here
the region $\mathbf{x}\in[-4,4]^{2}$, into pixels and darkening each pixel
following equation 6.

As discussed in the previous section, there are inherent benefits to
considering state-space representations of trajectories. For each input-output
pair $(\mathbf{T}_{i},\mathbf{\theta}_{i})_{i=1}^{N}$ in our dataset, we take
the trajectories $\mathbf{T}_{i}$ and ‘flatten’ them into a state-space
representation, essentially overlaying coarse scatter-plots of each trajectory
in state-space. We define this transformation by selecting a state-space
region and splitting it into ‘pixels’, or a higher-dimensional discretization
of the space in a more general setting, e.g. ‘voxels’ in three dimensions. Our
experiments were performed with an axis-aligned uniform grid of $128\times
128$ pixels. Given the $i^{\text{th}}$ collection of trajectories
$\mathbf{T}_{i}$, the pixels were ‘shaded’ by defining the value of the pixel
in the $(h,w)\in\\{0,1,...127\\}^{2}$ position in the $i^{\text{th}}$ sample
by

$\mathbf{P}_{i,h,w}\coloneqq\alpha^{n_{i,h,w}},$ (6)

where $n_{i,h,w}$ is the number of points in $\mathbf{T}_{i}$ which lie in the
region of space ascribed to the pixel at position $(h,w)$ and $\alpha\in(0,1]$
is a transformation parameter which determines the exponential base that the
pixels ‘darken’ with. We will use $\mathbf{P}_{i}$ to refer to the collection
of pixels $\mathbf{P}_{i,h,w}$, or ’pixelized return map’. The input-output
pairs after this transformation are then
$(\mathbf{\theta}_{i},\mathbf{P}_{i})_{i=1}^{N}$. The result of this process
is shown in the rightmost column of Figure 2. We performed all experiments
with $\alpha=7/10$.

####  3.2.1 Data Augmentation

We experimented with an additional step during model training before creating
each pixelized return map. Data augmentation Shorten2019  is a common practice
in training machine learning models in which carefully chosen random
modifications are made to training samples to improve the generalizability of
models to new data. In short, data augmentation can reduce the sensitivity of
the estimation performance of trained models to certain changes to their
input. An example of this method in image recognition tasks is making small
random changes to the cropping of input photos. The rationale is that, for
example, the subject of a photo with small differences in cropping is the
same, and we would like our models to respect this invariance. An additional
benefit of data augmentation is that it acts as a regularization mechanism
that reduces the tendency of trained models to overfit to patterns in the
training dataset at the expense of worse performance on new data. In our
setting, we performed data augmentation to make our model robust to changes in
the number and length of the trajectories we flattened to form the images.
Specifically, given the $i^{\text{th}}$ input-output pair, we first sampled
$N_{\text{traj}}\sim\text{Uniform}(10,n)$. We then selected $N_{\text{traj}}$
of the $n$ trajectories in $\mathbf{T}_{i}$ uniformly at random. Next, we
sampled $N_{\text{steps}}\sim\text{Uniform}(10,m)$, and for each of the
$N_{\text{traj}}$ trajectories we took only the first $N_{\text{steps}}$
points. The resulting collection of $N_{\text{traj}}$, length
$N_{\text{steps}}$ trajectories was then used to create a pixelized return map
as discussed in Section 3.2. Figure 2 shows the process of creating pixelized
return maps with this data augmentation method. The choice of 10 as the lower
limit on the number of trajectories and steps in the augmentation step is
another hyperparameter. Although we do not explore other choices for this
value in this work, one should consider that choosing a value too close to $m$
or $n$ limits the diversity of the training data; this may reduce the
regularization effect. On the other hand, training the model on samples with
too few trajectories or with too few steps may result in a trained model which
fails to generalize to more densely-sampled return maps.

###  3.3 Parameter Estimation Model

Figure 3: A network architecture diagram for the ResNet18 neural network
He2016 . The convolutional (conv), batch-norm (BN) Ioffe2015a , and linear
layers contain ‘learnable’ parameters, which we optimize during the model
training process using stochastic gradient descent. The network’s name comes
from the fact that each input image passes through 17 convolutional layers and
one linear layer. We used the PyTorch Paszke2017  TorchVision
TorchVisionMaintainersAndContributors2016  library’s implementation of
ResNet18, only modifying the first convolutional layer to take single channel
images as input.

Our experiments were all performed with a ResNet18 neural network He2016 , a
small model from a widely used family of deep-learning models which feature
convolutional layers and residual connections designed for computer-vision
tasks. The structure of the network is depicted in Figure 3. We selected this
architecture as the model’s residual connections and batch-normalization
layers Ioffe2015a  lead to relatively easy training of the model with current
methods. In the context of the formalism that we introduced in Section 2.3, we
consider estimation functions $g\in\mathcal{F}$ which are the composition of
our transformation from collections of trajectories to images of return maps,
with a ResNet18 neural network with learnable parameters. Selecting a specific
function $g$ amounts to setting the learnable parameters of the ResNet18, and
possibly the parameters which describe the transformation from trajectories to
images.

####  3.3.1 Model Training

Given a dataset $\mathcal{D}$ of $N$ input-output pairs, we partitioned the
pairs into training ($\mathcal{D}_{\text{train}}$), validation
($\mathcal{D}_{\text{validation}}$) and testing ($\mathcal{D}_{\text{test}}$)
datasets containing respectively 65%, 15% and 20% of the pairs in
$\mathcal{D}$. This split was performed in a manner that ensured that each
partition contained a representative mix of samples from the full range of
parameter values $\Theta$.

As is standard in training supervised machine learning models,
$\mathcal{D}_{\text{train}}$ and $\mathcal{D}_{\text{validation}}$ were used
to optimize the parameters of each model. Specifically,
$\mathcal{D}_{\text{train}}$ was used to optimize the learnable parameters of
the ResNet18 model using stochastic mini-batch gradient descent with the Adam
optimizer Kingma2017  as our optimization routine. We used default parameters
for the optimizer other than the learning rate and value of the weight decay
parameter. When using data augmentation, we performed fresh random
augmentation on each pass over the training set, i.e., the model essentially
never saw the same sample twice in the training routine. When applicable, the
validation set $\mathcal{D}_{\text{validation}}$ had random data augmentation
applied only once at the start of training. We fixed the validation set
augmentation to provide a consistent set to evaluate the performance of each
model, reducing the variance of the validation error. We used the model’s
error on $\mathcal{D}_{\text{validation}}$ to optimize the hyper-parameters of
our optimization routine. We focused on tuning the Adam optimizer’s learning
rate and weight decay parameters, as well as the number of training samples
used to evaluate each step of the optimizer (commonly referred to as the
training batch size). After model training, we opted to select the model
version with the smallest validation loss observed during the training
process. The resulting model’s performance was then evaluated on
$\mathcal{D}_{\text{test}}$, providing an unbiased estimation of model
performance as the test set consists of samples not used in the model training
or model selection process.

###  3.4 Loss Function

By framing parametric system identification as a supervised machine learning
problem, we can develop an estimation model that directly optimizes a loss
function on the estimation error. Still, the appropriate loss function to
optimize depends on the intended purpose of the parameter estimates. For
example, simply choosing the mean squared error may be sufficient for systems
with a single parameter. However, for systems with parameters
$\mathbf{\theta}\in\mathbb{R}^{d}$ for $d>1$, we may require greater precision
when estimating some coordinates of $\mathbf{\theta}$ than others. One way to
encode this information is in the loss function used during model training.
For example, in Section 4.1 we use a weighted loss function to ensure that the
relative errors of our model’s parameter estimates are roughly equal for both
Hénon map parameters.

Figure 4: Schematic diagram showing the process of generating training data,
training the parameter prediction model, and finally using the trained model
for parameter estimation.

Figure 4 outlines the entire process of creating a dataset from a dynamical
system, training a supervised machine learning model, and deploying the model
for parameter estimation.

##  4 Experiments

We trained parameter estimation models on continuous- and discrete-time
dynamical systems and performed experiments to understand the sensitivity of
our method to the amount of available training data, as well as the impact of
our data augmentation process. We show results for two systems in this work.
In order to build intuition for our proposed method, we begin with the simpler
case of parameter estimation for a discrete map before showing experiments for
a continuous-time dynamical system. In what follows we will refer to models
that are trained with or without augmented data as augmented and non-augmented
models respectively.

All of our experiments were carried out using a single NVIDIA® GeForce® RTX
2080 Ti GPU and two Intel® Xeon® E5-2650L v4 CPUs with a clock rate of 1.70GHz
and 251GB of RAM.

###  4.1 Hénon Map

In our first experiment, we estimate the parameters of the Hénon Map Henon1976
, a well-studied discrete dynamical system that maps pairs of points
$(x_{k},y_{k})\in\mathbb{R}^{2}$ to $(x_{k+1},y_{k+1})$ according to

$\begin{cases}x_{k+1}=1-ax_{k}^{2}+y_{k}\\\ y_{k+1}=bx_{k},\end{cases}6$ (7)

where $a,b\in\mathbb{R}$ are the parameters which characterize the system.

For our experiments, we considered the Hénon map parameters
$a\in[0.05,.45],b\in[-1.1,1.1]$. We constructed a database of input-output
pairs by taking a large number of $a,b$ points spaced evenly in
$[0.05,.45]\times[-1.1,1.1]$. For the $i^{\text{th}}$ $a,b$ pair, we generated
trajectories by applying the Hénon Map with parameters $a,b$ $250$ times over
a uniform grid of $225(=15^{2})$ initial $x,y$ points over the region
$[-4,4]\times[-4,4]$. With this collection of parameters and trajectories
$\mathcal{D}$, we processed features as outlined in Section 3.2 to obtain
image-like data, which we then used to train a ResNet18 neural network to make
parameter estimates in $\mathbb{R}^{2}$. Given the different scales of the $a$
and $b$ parameters in this task, we chose to measure our parameter prediction
performance according to a weighted mean squared error loss

$\ell\left((\hat{a},\hat{b}),(a,b)\right)=\dfrac{1}{2}\left(\sigma_{a}^{2}(\hat{a}-a)^{2}+\sigma_{b}^{2}(\hat{b}-b)^{2}\right),$
(8)

where $\sigma_{a},\sigma_{b}=(25,4.\overline{54})$ were chosen so that
$\sigma_{a}\cdot(a-a_{\text{mid}})$ and $\sigma_{b}\cdot(b-b_{\text{mid}})$
take values in the interval $[-5,5]$, where
$(a_{\text{mid}},b_{\text{mid}})=(-.25,0)$ is the center of the parameter
region $[0.05,.45]\times[-1.1,1.1]$

####  4.1.1 Results

Figure 5: Inputs and outputs from a trained parameter estimation model for the
Hénon map taking randomly augmented return map images as input. The left
column shows the input to the network, as well as the ground truth values of
$a,b$. The center column shows high-quality return map images, not input to
the network, corresponding to the system’s parameters, obtained by iterating
the map for a uniform grid of initial conditions using actual values of $a,b$.
The right column shows the model outputs and a return map obtained in the same
fashion as the middle column but using the model’s estimated values
$\hat{a},\hat{b}$.

Figure 5 shows input and output examples from a test set for a trained
parameter estimation model. As well as displaying the output parameter
estimates of the network, the figure shows return maps that correspond to
these parameter estimates, providing a qualitative assessment of our model’s
estimation accuracy. The inputs to the estimation model were noisily sampled
in the same fashion as in our data augmentation process (Section 3.2.1). We
also show more refined return maps, generated without this noisy sampling,
that correspond to the true parameter values for a more direct comparison with
the return maps corresponding to the estimation model’s outputs. For a perfect
estimation model, the return maps in the middle and right columns would be
essentially identical. The example in the third row of Figure 5 illustrates
the challenge of estimating parameters from the sparsely sampled return maps
that our data augmentation process can create. While the ‘ground truth’ return
map and the return map corresponding to the parameter estimates visually
differ near the center of the image but are otherwise similar. The input to
the network lacks samples in this central region, which makes the model’s
estimate plausible given its input.

Figure 6: Test mean squared error of augmented and non-augmented models
trained on datasets of increasing numbers of sample pairs from the Hénon map.
Each model was trained for 25,000 optimization steps with the same
configuration of the Adam optimizer, i.e. the same learning rate, weight
decay, and batch size parameters. The five model versions with the lowest
validation error observed during training were selected for each dataset size.
We then evaluated these models on a test set. We evaluated model errors on
inputs that matched their respective training distributions. I.e., we only
applied data augmentation to the test samples for augmented models. Error bars
represent bootstrapped 95% confidence intervals for models trained with three
random seeds.

Figure 6 summarizes the results of an experiment investigating the impact of
training dataset size and our data augmentation process on parameter
estimation error on the Hénon map. We emphasize two aspects of this figure,
the first of which is that the test error of our estimation models decreases
in an approximately power-law fashion with the size of the training set. A
second aspect of Figure 6 is that the augmented models seem to have lower test
error than non-augmented models on smaller numbers of training samples until
this trend reverses on larger training datasets.

Figure 7: Test error for a pair of augmented and non-augmented models, trained
on a Hénon map dataset of 8172 input-output pairs for 25,000 optimizer steps.
The non-augmented model was trained on trajectories with length $m=250$,
whereas the augmented train data contained trajectories varying in length from
$m=10$ to $m=250$. We then evaluated each model on test datasets which
consisted of fixed-length trajectories, with and without test-time
augmentation. For this experiment, the testing data augmentation consisted of
taking a random subset of $n\in\\{10,11,\ldots,225\\}$ trajectories for each
input-output pair. Error bands represent bootstrapped 95% confidence intervals
for models trained with three random seeds.

Figure 7 further explores the generalization gains observed from applying data
augmentation to samples during model training. Both models evaluated in the
figure appear to generalize across inputs with trajectories with lengths
ranging from 200 to 1000 samples, achieving a roughly constant error across
different trajectory lengths. The clearest difference between models in Figure
7 is that the estimation error of the non-augmented model (diamond markers)
grows up to 100 times larger on randomly augmented test inputs compared with
the performance on inputs without this randomization (the dashed and solid
lines respectively). In comparison, the augmented model displays a
significantly smaller gap in performance on augmented versus non-augmented
test data and, notably, demonstrated improved performance when tested on
inputs without the random augmentation that the model was trained with.

Finally, we examined the limits of the generalization abilities of augmented
models. Figure 8 shows the results of an experiment in which we tested the
estimation performance of such a model on test pairs with between 50 and 1000
length trajectories and between 25 and 400 trajectories per sample. As in
Figure 7, we observed that the model generalized well to short and long
trajectory lengths, including ones longer than the maximum length of 225
samples in the training set. The model’s performance was more significantly
affected by the number of trajectories per sample, with performance improving
as the number of trajectories increased from 25 to 225 before a significant
deterioration in performance on samples with 324 and 400 trajectories in each
sample.

Figure 8: A heatmap (darker is better) displaying the logarithm of the test
error of the Hénon map augmented model model for 25,000 optimizer steps Figure
7 (dashed orange curve), on test datasets in which the number of trajectories
is varied, as well as the length of trajectories.

###  4.2 The Swinging Atwood’s Machine

The swinging Atwood’s machine (SAM) Tufillaro1984  is a Hamiltonian system
with nonlinear dynamics, consisting of an Atwood’s machine in which one of the
mechanical bobs is allowed to swing from its pulley along the two-dimensional
plane containing the bobs and their pulleys. The configuration manifold of the
system is two-dimensional, consisting of the distance between the swinging
pendulum and its pulley $r\in(0,\infty)$, and the angle formed between the
swinging bob and vertical, $\phi\in(-\pi,\pi]$.

This system exhibits chaotic motion for nearly all values of the mass ratio
$\mu\in(0,\infty)$ between the stationary and the swinging pendulum, taken as
the parameter of the system Casasayas1990 . For our experiments, we set the
gravitational acceleration $g$, and the system’s mechanical energy $E$ to
unity, as is common when studying this system. In our experiments, we used
examples with $\mu\in[1.5,15]$. For ease of comparison with our Hénon map
results, we used a rescaled loss function

$\ell\left(\hat{\mu},\mu\right)=\dfrac{1}{2}\sigma_{\mu}^{2}(\hat{\mu}-\mu)^{2},$
(9)

where $\sigma_{\mu}=.\overline{740}$ was chosen so that
$\sigma_{\mu}\cdot(\mu-\mu_{\text{mid}})$ takes values in the interval
$[-5,5]$, where $\mu_{\text{mid}}=8.25$ is the center of the parameter
interval $[1.5,15]$.

####  4.2.1 Choice of Poincaré Section

For our experiments we used a common Poincaré section for the SAM: the section
where the state of the system $(r,\phi,\dot{r},\dot{\phi})$ passes through
$\phi=0$ with $\dot{\phi}>0$. Figure 9 shows return maps for this section
corresponding to four different values of $\mu$, including $\mu=3$ for which
the system is non-chaotic Tufillaro1984 .

Figure 9: Return maps for the SAM system obtained using the Poincaré section
described in Section 4.2.1. The system was integrated forwards for 1000 units
of time for 200 initial conditions, drawn randomly from energetically-allowed
initial states.

Given that our choice of Poincaré section constrains $\phi=0$ and the
Hamiltonian of the system is time-independent, the system’s state at each
section crossing is uniquely determined by the values of $r$ and its conjugate
momentum $p_{r}$. We used two-dimensional return maps that only consider these
values, e.g. the return maps in Figure 9.

####  4.2.2 Data Generation

To create the dataset $\mathcal{D}$ for these experiments, we used 4,000
evenly spaced values for $\mu\in[1.5,15]$. Because the SAM is a Hamiltonian
system and we considered crossing points for which $\phi=0$, one can determine
that all crossing points of the system lie within fixed values of $r,\phi$ for
each value of the value of mass ratio parameter $\mu$ at a fixed mechanical
energy of 1. The roughly triangular boundary is visible in the different
portraits in Figure 9. Using this fact, for each value of $\mu$ we used Heyoka
to integrate 256 initial states randomly selected from the energetically
allowed regions of the state space for 1,000 units of system time. For each
initial state, we recorded the state of the system each time it crossed
through the Poincaré section specified in Section 4.2.1. This data generation
method resulted in trajectories of varying lengths, depending on how many
times each trajectory passed through the section within 1,000 units of time,
and so in order to perform data augmentation we cropped each trajectory down
from its original length to between 1 and 250 points.

####  4.2.3 Results

Figure 10: Test mean squared error of augmented and non-augmented models on
datasets of increasing numbers of sample pairs from the swinging Atwood’s
Machine system with the. Each model was trained for 25,000 optimization steps
using the same configuration of the Adam optimizer (learning rate, weight
decay, batch size). For each dataset size, the five model versions with the
lowest validation error observed during training were selected. We then
evaluated these models on a test set. We evaluated model errors on inputs that
matched their respective training distributions. I.e., we only applied data
augmentation to the test samples for augmented models. Error bars indicate
bootstrapped $95\%$ confidence intervals from models trained with three
different random seeds. Figure 11:  Test error for a pair of augmented and
non-augmented models, trained for 25,000 optimizer steps on the swinging
Atwood’s Machine consisting of 7680 input-output pairs. Both the non-augmented
and augmented training trajectories varied in length from 1 to 250 points. We
evaluated each model on test datasets which consisted of fixed-length
trajectories, with and without test-time augmentation, which, for this
experiment, consisted of taking a random subset of the trajectories for each
input-output pair. the testing data augmentation consisted of taking a random
subset of $n\in\\{10,11,\ldots,225\\}$ trajectories for each input-output
pair. Error bands represent bootstrapped 95% confidence intervals from models
trained with three different random seeds.

The experiments in this section closely mirror those for the Hénon map system
in Section 4.1.1. Figure 10 summarizes the results of an experiment
investigating the impacts of training dataset size and our data augmentation
process for the swinging Atwood’s Machine system. Similar to in Figure 6 we
observe that the test error of our estimation models decreases in an
approximately power-law fashion with the size of the training set, although we
observe that the performance of the augmented model appears to plateau. Prior
to this plateau, we see that the augmented models often perform better than
non-augmented models.

Figure 11 shows the results of our experiment looking at the generalization
benefit of applying data augmentation to samples during model training. In
contrast with Figure 7 we see that the performance of both models is worse on
inputs with shorter trajectory lengths. Once again we observe that there is a
smaller generalization gap for augmented model as evidenced by the vertical
gap between the lines with circular markers versus those with diamond markers.

##  5 Discussion

By framing PSI as a supervised machine learning problem, we circumvent the use
of loss functions that compare simulated trajectories to observations.
Instead, we directly learn estimation functions that minimize parameter
estimation errors, which is a natural objective for parametric system
identification. In this setting, the ability to obtain effective estimation
functions from potentially limited datasets becomes a central concern. In our
experiments, we observed an approximately power-law reduction in estimation
errors by increasing the number of training samples used. This is a common
phenomenon in statistical estimation problems Wainwright2019  and implies that
the availability of training data is an important factor in the success of our
method.

With that said, we observed that for small sample sizes, there were
significant improvements in estimation accuracy when using random data
augmentation during model training. While the reason for the empirical success
of data augmentation is an open area of research (see e.g. Shen2022 ), one
explanation for our observations is that data augmentation effectively
increases the number and variety of training samples used to optimize the
estimation model’s parameters, resulting in a model which generalizes better
to new samples. On the other hand, we also observed faster improvement in
model performance with the size of the training set for non-augmented models,
with the test MSE of these models overtaking their augmented counterparts on
both dynamical systems in this paper. On this last point, while the test error
of non-augmented models on larger training datasets surpassed that of non-
augmented models, we would argue that the augmented models still have more
favourable characteristics for a practical model, as we discuss next.

Figures 7 and 11 provide a more complete characterization of the performance
of augmented and non-augmented models trained on the largest datasets shown in
Figures 6 and 10 – 8172 and 7680 samples, respectively. In these line plots,
we observe that the augmented models have improved performance across input
distributions of varying quality, especially on samples with shorter
trajectories. This can be seen in both the smaller gap between test errors on
augmented and non-augmented test samples and in the observation that the best-
performing models on test inputs with the shortest trajectories were those
trained with random data augmentation (left-hand sides of Figures 7 and 11).
These results indicate that the use of these augmentation methods is currently
a useful tool for building practical data-driven estimation models.

All of our experiments were performed with $\alpha=.7$ and 128 by 128-pixel
images. These ‘hyperparameters’, as they are referred to in the machine
learning community, can all in principle be optimized on validation data, as
we did with the batch size and parameters of the optimization routine, at the
cost of increased computation time.

A consideration for the method we propose is that we rely on changes to system
parameter values translating into changes to the resulting return maps. If
large changes to parameter values result in only small changes to return maps
then we would expect our method to have a large estimation error. This
behavior should also occur with all previous methods that rely on loss
functions that compare observed and simulated trajectories in state space. In
addition, in most scenarios, this effect is likely benign, given that in such
a scenario a large estimation error may be acceptable since trajectory
predictions and system analysis relying on estimated parameters would remain
reasonably accurate due to the low sensitivity of the dynamical system to
changes in its parameter values.

Continuous-time dynamical systems pose an additional challenge, in that if the
chosen Poincaré section fails to adequately capture key features of
trajectories, then by unfortunate coincidence two different parameter values
may result in different dynamics but share the same Poincaré map. This issue
is separate from a possible issue with the collected data in which all
parameter values for a system have distinct Poincaré maps for a given Poincaré
section, but there are not enough data points collected for each parameter
value to disambiguate between parameter values in some cases. Instead, we
refer to a scenario where for distinct parameter values $\theta$ and
$\theta^{\prime}$, the dynamics of the system result in the same Poincaré map
but different trajectories away from the Poincaré section, in which case no
amount of data collected on such a Poincaré section can distinguish between
systems parameterized by $\theta$ versus $\theta^{\prime}$.

This consideration means that care should be taken to ensure that the choice
of Poincaré section used to study a dynamical system captures the dynamics of
interest. Our method straightforwardly accommodates using two or more return
maps for each input pair, for example, return maps corresponding to different
Poincaré maps for the same parameter value, by stacking each return map as an
image channel before input to the ResNet parameter estimation model. This
method effectively provides the estimation model with different ‘cross-
sections’ of trajectory dynamics to use for parameter estimates and may lead
to better estimates if one of the Poincaré sections is not informative for
some parameter values.

##  6 Conclusion

In this work, we have introduced a novel solution method for parametric system
identification (PSI), framing the problem of mapping sample trajectories to
system parameters as a supervised machine learning problem. Combining this
idea with recent approaches to PSI, which use state-space representations of
trajectories, we show that our approach is effective for parameter
identification on chaotic dynamical systems. Since we use return maps as input
to a supervised machine learning model, our method generalizes to continuous
time dynamical systems using Poincaré maps. Although training the estimation
models is a compute-intensive process, the resulting models provide fast,
inexpensive, and accurate parameter estimates that can be used to process
large collections of data or as subcomponents of a more complicated program.
While we focused on systems with two-dimensional in this work, our method
extends to more dimensions simply by replacing the two-dimensional pixelized
transformation and subsequent convolution operations in the ResNet with their
higher dimensional analogs.

This work opens the door to future exploration of machine learning for PSI. In
our current method, we effectively reduce the problem of PSI to a computer
vision task, and so there are likely performance improvements, both in terms
of estimation accuracy and in terms of reducing the number of samples required
to reach high accuracy by more carefully optimizing the estimation
architecture used. One possible avenue is to use the wide availability of
large, pre-trained deep learning models trained on image tasks with more
layers than our ResNet18 model and only optimize the final layers of the
network for the parameter estimation task. This ‘transfer learning’ Weiss2016
method can potentially reduce the number of samples required to train a new
estimation model.

A second, more open direction is to investigate more sophisticated methods for
structuring return maps as input to a deep learning model. Our approach of
flattening these maps into single-channel images introduces several extraneous
parameters. In addition, pixelating return maps potentially removes some of
the finer spatial structure present in the return map when using coarser
pixelation schemes. An alternative approach could be to explicitly represent
return maps as unstructured collections of points and use a graph neural
network Scarselli2009  as the parameter estimation model.

Deep learning and other data-driven modeling and optimization methods such as
reinforcement learning are becoming increasingly useful tools in the physical
sciences. While these approaches can be powerful, there are often significant
challenges associated with developing useful applications in different problem
domains. In this work, we take a step in the right direction, and we invite
the research community to explore the possibilities of using data-driven
methods and return maps to analyze complex dynamical systems.

Acknowledgments

## Declarations

### Funding

The authors declare that no funds, grants, or other support were received
during the preparation of this manuscript.

### Conflict of interest/Competing interests

The authors have no relevant financial or non-financial interests to disclose.

### Ethics approval

Not applicable.

### Consent to participate

Not applicable.

### Consent for publication

Not applicable.

### Availability of data and materials

See the subsection ‘Code Availability’.

### Code availability

We have included the code to replicate our experiments here. The code
repository contains the configuration files used to generate the SAM dataset
files, as well as the configuration files we used to train the models on the
SAM and Hénon map systems shown in this paper.

### Authors’ contributions

Emmanuel Blazquez provided project supervision and editorial support in
writing the manuscript. Connor James Stephens developed the code used to run
the experiments in this paper and wrote the manuscript.

## References

  *   * (1) Jumper, J., Evans, R., Pritzel, A., Green, T., Figurnov, M., Ronneberger, O., Tunyasuvunakool, K., Bates, R., Žídek, A., Potapenko, A., Bridgland, A., Meyer, C., Kohl, S.A.A., Ballard, A.J., Cowie, A., Romera-Paredes, B., Nikolov, S., Jain, R., Adler, J., Back, T., Petersen, S., Reiman, D., Clancy, E., Zielinski, M., Steinegger, M., Pacholska, M., Berghammer, T., Bodenstein, S., Silver, D., Vinyals, O., Senior, A.W., Kavukcuoglu, K., Kohli, P., Hassabis, D.: Highly accurate protein structure prediction with AlphaFold 596(7873), 583–589
  * (2) Lam, R., Sanchez-Gonzalez, A., Willson, M., Wirnsberger, P., Fortunato, M., Pritzel, A., Ravuri, S., Ewalds, T., Alet, F., Eaton-Rosen, Z., Hu, W., Merose, A., Hoyer, S., Holland, G., Stott, J., Vinyals, O., Mohamed, S., Battaglia, P.: GraphCast: Learning Skillful Medium-Range Global Weather Forecasting
  * (3) Degrave, J., Felici, F., Buchli, J., Neunert, M., Tracey, B., Carpanese, F., Ewalds, T., Hafner, R., Abdolmaleki, A., Diego de las Casas, Donner, C., Fritz, L., Galperti, C., Huber, A., Keeling, J., Tsimpoukelli, M., Kay, J., Merle, A., Moret, J.-M., Noury, S., Pesamosca, F., Pfau, D., Sauter, O., Sommariva, C., Coda, S., Duval, B., Fasoli, A., Kohli, P., Kavukcuoglu, K., Hassabis, D., Riedmiller, M.: Magnetic control of tokamak plasmas through deep reinforcement learning 602(7897), 414–419
  * (4) Söderström, T., Stoica, P.: System Identification. Prentice Hall
  * (5) Alvin, K.F., Robertson, A.N., Reich, G.W., Park, K.C.: Structural system identification: From reality to models 81(12), 1149–1176
  * (6) Bruggeman, F.J., Westerhoff, H.V.: The nature of systems biology 15(1), 45–50 17113776
  * (7) Ljung, L.: Perspectives on system identification 34(1), 1–12
  * (8) Villaverde, A.F., Banga, J.R.: Reverse engineering and identification in systems biology: Strategies, perspectives and challenges 11(91), 20130505
  * (9) Quaranta, G., Lacarbonara, W., Masri, S.F.: A review on computational intelligence for identification of nonlinear dynamical systems 99(2), 1709–1761
  * (10) Konnur, R.: Estimation of all model parameters of chaotic systems from discrete scalar time series measurements 346(4), 275–280
  * (11) Tao, C., Zhang, Y., Jiang, J.J.: Estimating system parameters from chaotic time series with synchronization optimized by a genetic algorithm 76(1), 016209
  * (12) Wang, L., Xu, Y., Li, L.: Parameter identification of chaotic systems by hybrid Nelder–Mead simplex search and differential evolution algorithm 38(4), 3238–3245
  * (13) Ho, W.-H., Chou, J.-H., Guo, C.-Y.: Parameter identification of chaotic systems using improved differential evolution algorithm 61(1), 29–41
  * (14) Modares, H., Alfi, A., Fateh, M.-M.: Parameter identification of chaotic dynamic systems through an improved particle swarm optimization 37(5), 3714–3720
  * (15) Jafari, S., Sprott, J.C., Pham, V.-T., Golpayegani, S.M.R.H., Jafari, A.H.: A New Cost Function for Parameter Estimation of Chaotic Systems Using Return Maps as Fingerprints 24(10), 1450134
  * (16) Mousazadeh, A., Shekofteh, Y.: Cost function based on the self-organizing map for parameter estimation of chaotic discrete-time systems 94, 103817
  * (17) Shekofteh, Y., Jafari, S., Sprott, J.C., Hashemi Golpayegani, S.M.R., Almasganj, F.: A Gaussian mixture model based cost function for parameter estimation of chaotic biological systems 20(2), 469–481
  * (18) Shekofteh, Y., Jafari, S., Rajagopal, K., Pham, V.-T.: Parameter Identification of Chaotic Systems Using a Modified Cost Function Including Static and Dynamic Information of Attractors in the State Space 38(5), 2039–2054
  * (19) Shekofteh, Y., Jafari, S., Rajagopal, K.: Cost function based on hidden Markov models for parameter estimation of chaotic systems 23(13), 4765–4776
  * (20) Mousazadeh, A., Shekofteh, Y.: Parameter estimation of chaotic dynamical systems using LS-based cost functions on the state space 96(1), 3
  * (21) Parker, T.S., Chua, L.O.: Practical Numerical Algorithms for Chaotic Systems. vol. 56
  * (22) Murphy, K.P.: Machine Learning: A Probabilistic Perspective. MIT Press
  * (23) Biscani, F., Izzo, D.: Revisiting high-order Taylor methods for astrodynamics and celestial mechanics 504(2), 2614–2628
  * (24) Shorten, C., Khoshgoftaar, T.M.: A survey on Image Data Augmentation for Deep Learning 6(1), 60
  * (25) He, K., Zhang, X., Ren, S., Sun, J.: Deep Residual Learning for Image Recognition. In: 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 770–778
  * (26) Ioffe, S., Szegedy, C.: Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift. In: Proceedings of the 32nd International Conference on Machine Learning, pp. 448–456. PMLR
  * (27) Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., Lerer, A.: Automatic differentiation in PyTorch
  * (28) TorchVision maintainers and contributors: TorchVision: PyTorch’s Computer Vision Library
  * (29) Kingma, D.P., Ba, J.: Adam: A Method for Stochastic Optimization
  * (30) Hénon, M.: A two-dimensional mapping with a strange attractor 50(1), 69–77
  * (31) Tufillaro, N.B., Abbott, T.A., Griffiths, D.J.: Swinging Atwood’s Machine 52(10), 895–903
  * (32) Casasayas, J., Nunes, A., Tufillaro, N.: Swinging Atwood’s Machine : Integrability and dynamics 51(16), 1693–1702
  * (33) Wainwright, M.J.: High-Dimensional Statistics: A Non-Asymptotic Viewpoint, 1st edn. Cambridge University Press
  * (34) Shen, R., Bubeck, S., Gunasekar, S.: Data Augmentation as Feature Manipulation
  * (35) Weiss, K., Khoshgoftaar, T.M., Wang, D.: A survey of transfer learning 3(1), 9
  * (36) Scarselli, F., Gori, M., Tsoi, A.C., Hagenbuchner, M., Monfardini, G.: The Graph Neural Network Model 20(1), 61–80


TOPOLOGY AND BEHAVIOUR OF AGENTS: CAPITAL MARKETS

O. Hudak

Faculty of Finance, Matej Bel University,

Cesta na amfiteater 1, SK-974 15 Banska Bystrica, Slovak Republic,

e-mail<EMAIL_ADDRESS>
## Abstract

On a capital market the social group is formed from traders. Individual
behaviour of agents is influenced by the need to associate with other agents
and to obtain the approval of other agents in the group. Making decisions an
individual equates own needs with those of the other agents from the group.
Any two agents from the group may interact. The interaction consists of the
exchange of information and it costs some money. We assume that agents
interact in such a way that they give reference to the origin of the
information if asked by other agents. Thus the agent may verify obtained
private information. Hudak recently used methods described by Rivier to study
social behaviour of such agents. He characterized the quantity which
corresponds to verification of information. Here we study a capital market and
social behaviour of agents relations. Quantity which characterizes
verification of information contributes to an aversion of an agent with
respect to a risk. The mix of investments of an agent in a given cell with an
average measure A of risk aversion in the cell is found from minimum of the
average per cell aim function $<FM>$. Absolute minimum corresponds to such a
state in which there is an optimal mix of the exchange of information for a
given expectations about the capital market. The crowd and personal
/$\approx<f>$/ contributions to the risk aversion of an agent are present in
the aversion constant A. We have discussed a stable and metastable states of
the market for different values of E, an expected return for a given
investment period, of EV, an expected risk for a given investment period, and
of b, a constant which characterizes contribution of the quantity $<f>$ to the
risk aversion. Characteristics of the variance of n for the distribution of
nonreducibile subgroups are found. Our model describes intermediary process
effects.

##  1 Introduction

In a social group its members are called agents. The social group is
characterized in general in accordance with [1] \- [5]: we assume that
individual behaviour of agents is influenced to some degree by the need to
associate with other agents and to obtain the approval of other agents in the
group which is characterized by a very large nonrational and emotional element
to decisions of agents. Making decisions an individual equates own needs with
those of the other agents from the group. We assume that any two agents from
the group may interact. The interaction consists of the exchange of
information and it costs some ”energy” e.i. there are costs associated with
the exchange of information. We assume that the information is well defined,
and we assume that agents interact in such a way that they give reference to
the origin of the information if asked by other agents. The information is
private [6]. Public information is not studied in this paper. Thus there is
present information asymmetry. We assume that the agent may verify obtained
information. It means that the agent is giving references to two and only two
(here we assume two and only two foe simplicity) with which he/she exchanged
information. These two agents may be then called by the law court to give
evidence about the information which they exchanged with the agent in a case
of a law trial. Unverifiability is itself a special kind of information
problem [6]. Then there exits a subgroup of interacting agents the interaction
of which has the following property: it is nonreducible [5].

A cell is such a configuration of a number of nonreducible subgroups in which
every two interacting agents belong to two nonreducibile subgroups /subgroups
are connected in this sense/ and which is closed. Such a cell may disappear
and may be created, may change number of nonreducible subgroups in a
reversible way. Because the structure, configuration of interactions between
agents in the group, forms a macroscopic structure we say that it is a
microreversibile process any process within a nonreducibile subgroup and
within a cell. Statistical equilibrium of the whole group is characterized by
a set of different subgroups of the type mentioned above and by a probability
that such a subgroup occurs. Thus we have probability distribution which
characterizes the group. Moreover there exists an equation of state which
enables to compare different macroscopic states of the group. The statistical
equilibrium due to microreversibility is characterized by the maximum of
entropy and by the minimum of costs of information exchange minus a return.

We used [5] methods of statistical physics to study social behaviour of
agents, mainly the presence of topological structure of interactions between
agents and its changes, which is the most important property of the group of
agents. There are three empirically observed dependencies of personal radius
which enabled us to characterize the quantities of cells, faces /nonreducibile
subgroups/, vertices /agents/ and bonds /interactions/ [2]. There exist
constrains, such as a fixed number V of agents in the group, a number E of
interactions within the group, a number F of subgroups which are nonreducibile
/faces/, and a number C of cells. Thus we have a structure which is equivalent
to random cellular networks. Such networks and their evolution were described
by Rivier, [3] and [4]. Rivier applied methods of statistical mechanics to
study these networks.

Three empirically observed dependencies of personal radius dependence on some
factors enabled us [5] to characterize the quantity F which characterizes
verification of information. In this paper we study capital market and social
behaviour of agents. Quantity F, and in general $<f>$ the average number of
nonreducibile subgroups per cell, is the quantity which characterizes
verification of information and thus contributes to aversion of agent with
respect to a risk. The mix of investments of an agent in a given cell with an
average measure A of risk aversion in the cell is found from extremes of the
function FM, the aim function. We find the minimum of the average /per cell/
function $<FM>$. An absolute minimum of the function $<FM>$ corresponds to
such a state in the group of agents in which there is an optimal mix of
exchange of information for a given expectation of return and a given
expectation of risk on the capital market. Finding a minimum of the aim
function $<FM>$ with respect to the parameter $<n>$ of the group of agents we
find the optimal mix of exchange of information. We assume that the number of
nonreducible subgroups F per cell $<f>$ contributes linearly to aversion
constant A, $(A=A_{0}+(<f>-4).b)$. The crowd /$A_{1}=A_{0}-4.b$/ and personal
/$\approx<f>$ / contributions to the risk aversion of an agent are thus
present. We use $\frac{<f>.<n>}{3}$ as an average number of agents per cell.

We studied in this paper capital market and social behaviour of agents. We
have formulated model describing behaviour of agents on the capital market. We
have found that a pyramidal structure appears in the group, $<f>=4$, then
$A=A_{0}$. This type of structure corresponds to hierarchical economy systems.
We have found that when the structure contains topologically only one cell
then $<n>=6$ and A tends to infinity. Better verification of information leads
to expectation of higher returns and thus the acceptable risk is larger. This
type of structure corresponds to market-based economy systems.

In practice both market- and hierarchy- based systems co-exist in modern
economy. There exist open markets and there exist large organizations. We have
found that this co-existence coresponds to a state which is the minimum of the
aim function is in the range $3<<n><6$. In our model the market-based system
coresponds to a state with $<n>=6$ and the hierarchy-based szstem to a state
with $<n>=3$.

The mix of investments of an agent in a given cell with an average measure of
risk aversion in the cell is found from extremes of the function FM, the aim
function. In our model we included intermediary property of firms - cells.
Thus we have formulated a general model for topology of exchange of
information and behaviour of agents on capital markets which may take into
account intermediary property. We find the minimum of the average /per cell/
function $<FM>$ for the simpler case in which there is zero intermediary. Thus
assume that in this paper in our model the individual investors and agents
from financial intermediary firms do not differ. More general case will be
studied elsewier.

The extreme of the aim function $<FM>$ is minimum for some conditions on the
market and in the group of agents. This minimum corresponds to $<n>$ between 3
to 6 as we mentioned above. This is a stable state. The extreme may be also a
maximum. Then there may be two minima of the function FM or one minimum. One
of them is that which corresponds to $<n>=3$ /the pyramidal hierarchical
structure/ which is either a metastable state either a stable state, depending
on conditions on the market and in the group of agents. The other one
corresponds to $<n>=6$ which is a stable state. We call efficient capital
markets markets which are of the second type. We have characterized variance
of n for the distribution of nonreducibile subgroups. This variance diverges
to infinity with $<n>$ tending to 6. Thus the market-based economy system with
$<n>$ near or equal to has very large variance of nonreducibile subgroup
number of agents.

##  2 Capital Market and Social Behaviour of Agents

We know [5] that a personal diameter r increases when F decreases. We
generalize this relation to every cell: we assume that personal diameter r
increases when $<f>$, the average number of nonreducubile groups per cell,
decreases, $<f>$ is an average number of faces per cell in language of
topology. Quantity F, and generalizing quantity $<f>$, is the quantity which
characterizes verification of information. The larger F, and generalizing
$<f>$, the more verified information in general, e.i. not only a specific
information exchanged by members of a subgroup. Aversion of an agent to risk
on capital markets also contains characterization of verification of
information. We assume that the number of nonreducible subgroups F per cell,
$<f>$, contributes linearly to this aversion constant A:

$A=A_{0}-4.b+<f>.b.$ (1)

Here $A_{0}$ is a risk aversion constant for $<f>=4$ when $A=A_{0}$, b is a
constant which characterizes contribution of the quantity $<f>$ to the risk
aversion. The larger quantity $<f>$ the better verification of information,
the larger expected return contribution of the investment and thus the larger
acceptable risk.

Note that when the pyramidal structure appears in the group, $<f>=4$, then
$A=A_{0}$. This type of structure corresponds to hierarchical economy systems,
for these systems see in [6]. The aversion constant A has minimum value as
concerning dependence on the quantity $<f>$ in this case. Then lower risk is
acceptable and correspondingly lower profit is acceptable also, especially for
low $A_{0}$ constant. The constant A does not contain contribution from the
verification of information. Due to very low risk inthis case price signals do
not work in these systems. Information on resources and aims and objectives,
is flowing /exchanging/ through the hierarchy to the decision makers.

Note that when the structure contains topologically only one cell then $<n>=6$
and A tends to infinity. Better verification of information leads to
expectation of higher returns and thus the acceptable risk is larger. This
type of structure corresponds to market-based economy systems, for these
systems see in [6]. They work via price signals. We interpret $(<f>.b)$ as an
average measure of risk aversion of a personal /individual/ contribution to
risk aversion different from the crowd contribution of this person $A_{0}-4.b$
in the cell, because F and $<f>$ are characteristics of verification of
information. Trader - agent supplies or demands goods and services /on the
capital market shares or bonds/ if the market price exceeds or undervalues
his/her own valuation of shares or bonds. To make own valuation of shares or
bonds the agent individually verifies information on it.

In practice, as noted in [6], both types of system involve interaction of
people: in the market system they interact as traders, in the hierarchy as
agents within an organization. Both market- and hierarchy- based systems co-
exist in modern economy [6]: we have open markets and we also have large
organizations such as the joint-stock company.

##  3 A financial intermediary

In market economies the key role of providing financial intermediaries and
transaction services play banks and other financial firms. Intermediation
process is understood well in a hypothetical market economy by considering how
it might function without financial intermediaries [7]. For this economy there
is assumed for simplicity that government receipts equal outlays and that the
exports of goods and services equal import. It is also assumed that all
savings is by households and all investment is by businesses. Household claim
on business may be in the form of debt or equity capital. This is minor
additional difference from reality [7]. The level of interest rates and the
associated valuation of equity in this economy is determined by supply and
demand for savings, the influence of any monetary authority is not taken into
account [7]. In the absence of intermediaries households would have to hold
their savings in the form of equity or debt claims on specific firms. They are
limited in their ability to diversify their holdings by the high costs of
obtaining information about many different companies and of dealing with a
large volume of small denomination securities, [7]. Thus each household would
have to expect the large risk due to particular firms in which it invested.
This could result in a loss of most or all of its savings. Such a large risk
would lead to a high rate of interest and cost of capital which a household
would have to take. Reduction of their risk would lead to households willing
to accept a lower rate of interest. The financial intermediaries play a key
role, [7], to reduce the risks faced by households /individual savers/, by
pooling their savings and using these to assemble diversified portfolios of
assets. Diversification is such that risks of bankruptcy of different firms do
not depend on the same economic conditions specific to the firm, industry,
geographic region or the entire economy. This diversification requires
specialized knowledge and expertise, and a large portfolio. Low, or even
negative, correlation among the returns on different assets in portfolio and
larger number of these assets reduce risk. Thus intermediaries can pay a lower
rate of interest to households. Competition among financial intermediaries
leads together with lower rate of interest payed to households to lower costs
of financing, which leads to lower costs of financing firm enterprises, [7].
This hypothetical market economy we will use in this paper to describe
intermediary process quantitatively.

In our model of a capital market every point /vertex/ corresponds to an agent
acting on the market. In the simplified model economy described above these
agents would be households. Financial intermediaries which lead to lower costs
will be described in our model as cells with a given risk aversion constant of
the firm to which risk aversion constants of agents from this firm - cell
should not be too far. Some other cells are cells with zero risk tolerance
cell constant, agents in these cells corespond to households /individual
savers - agents/. Thus we have in our model not only households but also
financial intermediaries description included.

##  4 The mix of investment and interaction of agents

We will assume that every agent has its objective to maximize its expected
utility of wealth [8]. If returns are normally distributed and the investor
has constant absolute risk aversion, then expected utility can be written as
[9]:

$EU=-exp{(-c(E-\frac{c}{2}.EV))}$ (2)

where the risk tolerance constant is $c=\frac{1}{A}$, A is the risk aversion
constant, E is the expected value of end-of-period wealth, EV is the expected
variance of end-of-period wealth. To make this utility as large as possible,
one maximizes [10]:

$E-\frac{c}{2}.EV.$ (3)

This leads to maximum expected utility of an agent, assuming that the risk
aversion is constant. Agents are interacting, some of them correspond to
households /individual investors/, some of them are from the financial
intermediary firms. In our model we will assume that every agent has the same
E / the expected value of end-of-period wealth, in percents/ and the same EV
/the expected variance of end-of-period wealth, in percents/. Thus c is a
number, and A is also a number. In reality our assumption is not true, there
exist dispersion of E and EV quantities. Assuming the same E and EV quantities
we assume that there exist a mean value of E and Ev for the dispersion.

We would like to consider not only maximization of expected utility of an
agent, but also optimization of interactions of agents in such a way in which
interacting agents corresponding to households and to intermediary financial
firms maximize their expected utility. While the first case corresponds to
maximization of (3), the second case is more general. To maximize the expected
utility of an agent in fact it is necessary to maximize the function:

$fM=c(E-\frac{c}{2}.EV)$ (4)

In this function there is a constant c. This constant is dependent on the
constant A which depends on the quantity $<f>$. Optimization of interactions
of agents in such a way in which interacting agents corresponding to
households and to intermediary financial firms maximize their expected utility
leads to optimization of (4) for every agent taking into account the structure
of a cell. The cell may represent agents - households, but also a financial
intermediary. The first one cell will be described in our model as a cell with
different aversion constants of agents from this cell, the second one will be
described in our model as a cell with almost the same aversion constant of
agents from this cell. To model this fact we will optimize the following
function:

$fM=\sum_{i=1}^{i=V_{C}}c_{i}(E-\frac{c_{i}}{2}.EV)-\sum_{i=1}^{i=V_{C}}\gamma(c_{i}-c_{B})^{2}$
(5)

Here $V_{C}$ is a number of agents in the cell C. The first sum is a sum of
terms of the type (3), the second term is a sum which describes how far are
constants $c_{i}$ of agents from a constant $c_{B}$ which is characteristic
for a given cell. The constant $\gamma$ is a positive or zero constant. If it
is zero, then the function fM describes a sum of functions fM for single
agents. This correspond to individual investors in the sense that their
aversion to risk constant is different, individual. If it is nonzero, then the
function fM describes a sum of functions fM for agents which are not
individual investors in the sense that their aversion to risk constant is not
too much different. These agents are agents from an intermediary financial
firm characterized by a constant $c_{B}$. If the constant $\gamma$ is very
large, then all the constants $c_{i}$ from this cell tend to the same value
$c_{B}$ characteristic for the firm.

To interpret the second sum in terms of expected returns and expected
variances let us rewrite the function fM (5) in the form:

$fM=\sum_{i=1}^{i=V_{C}}(c_{i}((E+2.c_{B}.\gamma)-\frac{c_{i}}{2}.(EV+2.\gamma))-c_{B}^{2}.\gamma)$
(6)

We see that the second term in (5) corresponds to higher total expected value
$E+2.c_{B}.\gamma$ in which besides an expected value E of the individual
investor a new contribution $2.c_{B}.\gamma$ appears. This contribution is due
to the firm financial intermediary, due to diversification which is such that
risks of bankruptcy of different businesses do not depend on the same economic
conditions specific to the firm, industry, geographic region or the entire
economy. We see that the expected return of an individual investor should be
higher to have the same expected value as total expected return /we are taking
always a return for a given period/ for the financial intermediary firm. This
diversification requires, as written above, specialized knowledge and
expertise, and a large portfolio. The higher risk tolerance constant $c_{B}$
of the firm leads to the higher total expected return.

We see further that the second term in (5) corresponds to higher total
expected risk $EV+2.\gamma$, in which besides an expected variance EV of the
individual investor a new contribution $2.\gamma$ appears. This contribution
is again due to firm financial intermediary, due to diversification which is
such that risks of bankruptcy of different businesses do not depend on the
same economic conditions specific to the firm, industry, geographic region or
the entire economy. This diversification requires, as written above,
specialized knowledge and expertise, and a large portfolio thus larger
expected risk may be accepted. We see that the expected risk of an individual
investor should be higher to have the same expected risk as total expected
risk for the financial intermediary firm. The stronger tendency /the higher
value of the constant $\gamma$ / in the firm to have the same risk tolerance
constant $c_{B}$ for all agents in the firm the higher total expected variance
- risk. The constant contribution $-c_{B}^{2}.\gamma$ gives a term in fM which
is proportional to the number of agents in the cell - firm. This contribution
is larger for larger risk tolerance constant $c_{B}$ /in absolute value/ and
is negative.

Optimization of (6) for given E, EV and $c_{B}$ leads to optimization of the
number of agents in the cell and to optimization of verification of
information which is characterized by f. Correspondingly optimization for the
whole group /agents on capital market/ for given E, EV and $c_{B}$ leads to
optimization of the number of agents in the cells and to optimization of
verification of information on the market, which is characterized by $<f>$.
Then structure of exchange of information in the group of agents will be such
that we obtain the highest value of the function fM.

We see further that the second term in (5) for $\gamma$ very large leads to
all risk tolerance constants in a cell to be the same. We see also that the
second term in (5) for $\gamma=0$, e.i. for vanishing $\gamma$, leads to all
risk tolerance constants in a cell to be different.

To obtain an aim function FM for such a cell, we will consider function fM
from (5) with negative sign.

The mix of investment of a single cell in the mean field approximation is
found from the function:

$<FM>=(-c.E_{B}+\frac{c^{2}}{2}EV_{B})\frac{<f>.<n>}{3}+\gamma.c_{B}^{2}.\frac{<f>.<n>}{3},$
(7)

where $E_{B}=E+2\gamma c_{B}$ is an expected return for a given investment
period modified by a contribution of the return from the cell /firm/ financial
intermediary, $EV_{B}=EV+2\gamma$ is an expected risk for a given investment
period modified by contribution of the risk from the cell /firm/ financial
intermediary, E and EV is an expected return and an expected risk respectively
for a given investment period. Here $c_{B}$ is a risk tolerance constant given
for the cell.

The mix of investment of a given cell /financial intermediary firm/ with an
average measure c of risk tolerance of a person /which is inversely equal to a
crowd contribution of a person plus a personal /individual/ contribution of
this person in the cell/ is found from the average /per cell/ function $<FM>$
above (7). Let us now study the case $\gamma=0$, which does not distinguish
between individual investors /for which $\gamma=0$ / and a financial
intermediary firm.

##  5 The case in which individual investors and a financial intermediary firm
are not distinguished

The case in which individual investors and a financial intermediary firm are
not distinguished we will study now. Thus we study the case in which financial
intermediary in the capital market is not considered, all agents of such
financial intermediary firm behave as individual investors. The case
$\gamma>0$, which does distinguish between individual investors /for which
$\gamma=0$ / and a financial intermediary firm / for which $\gamma>0$ / will
be studied elsewhere.

Thus the equation for the aim function FM is found from (7):

$<FM>=(-c.E+\frac{c^{2}}{2}EV)\frac{<f>.<n>}{3},$ (8)

where $\frac{<f>.<n>}{3}$ is an average number of agents per cell. From (1)
and (8) it has the form, note that $c=\frac{1}{A}=\frac{1}{A_{0}-4.b+<f>.b}$:

$<FM>=(-\frac{E}{A_{0}-4.b+<f>.b}+\frac{EV}{2.(A_{0}-4.b+<f>.b)^{2}})\frac{<f>.<n>}{3}$
(9)

and where for an equilibrium structure with a given number of cell C, of faces
F, of interactions E, and of agents V, see [3] and [4]:

$<f>=\frac{12}{(6-<n>)},$ (10)

and where $A_{1}=A_{0}-4.b$. From the equation for the aim function FM (9) we
will now find which number of cell C, of faces F and of interactions E it
minimizes, here number of agents V is given. It means that there is a
minimization of FM with respect to $<n>$, number of interaction in an
ireducibile subgroup /face/. Thus we are looking for such a structure of
agents on the capital market which gives the lowest value of the aim function
FM taking into account their interaction exchanging information. The case
which does distinguish between individual investors and a financial
intermediary firm is more complicated and will be studied elsewhere.

##  6 Extremes of the function FM

The function FM has extremes. Absolute minimum corresponds to such a state in
which there is an optimal mix. There exists an extreme of FM $<n>^{\prime}$
which is given by:

$<n>^{\prime}=6.\frac{(EV(A_{0}-2b)-2E(A_{0}-2b)^{2})}{(EV.A_{0}-2E(A_{0}-4b)(A_{0}-2b))}.$
(11)

Note that when contribution of nonreducibile subgroups is absent $b=0$, then
$<n>^{\prime}=6$. Thus only one cell exists with many nonreducibile subgroups
in this case. Note that when the crowd and personal contributions to the risk
aversion of an agent are present then $b>0$, $<n>^{\prime}\neq 6$ and more
than one cell maz exist in the group. There are several subgroups -cells. Note
that when the crowd risk aversion is zero, $A_{0}=4.b$, then:

$<n>^{\prime}=-3(1+\frac{4bE}{EV})$ (12)

which is less than 3, for E and EV positive. The quantity $<n>^{\prime}$ tends
to 0 when $A_{0}$ tends to $2.b$. Many cells exist in the group - pyramidal
structure may exist for $<n>^{\prime}=3$. Thus it is necessary to discuss
minima of the aim function FM.

##  7 The extreme is minimum

The extreme may be an absolute minimum for $A_{0}<2b$. This is the case when
the tendency to verify information is stronger. The extreme (11) is now
minimum if:

$2E(2b-A_{0})>EV>E\frac{A_{0}}{b}(2b-A_{0})$ (13)

and here the state with $<n>=6$ has lower value of the aim function than the
state $<n>=3$.

The extreme is also minimum if:

$E\frac{A_{0}}{b}(2b-A_{0})>EV\geq 2E\frac{A_{0}(2b-A_{0})}{4b-A_{0}},$ (14)

and here the state with $<n>=3$ has now lower value of the aim function than
the state $<n>=6$.

In both cases the extreme is a state in which some type of hierarchical
structures exists locally in the group, the group is not hierarchical as
whole. The extreme - minimum gives the average number of agents per
nonreducibie subgroup which is less than 6 and larger than 3.

The extreme (11) is also minimum if:

$EV>2E(2b-A_{0}).$ (15)

Here however the state with $<n>=6$ is minimum due to condition that
$3\leq<n>\leq 6$ is not fulfilled for $<n>^{\prime}$, this extreme is larger
than 6. In this case the minimum of the aim function with $<n>=6$ is a state
in which there is no type of hierarchical structure locally present in the
group.

The extreme (11) is also minimum if

$0\leq EV<E\frac{A_{0}}{b}(2b-A_{0}).$ (16)

Here the state with $<n>=3$ is minimum due to the condition that $<n>$ is from
the interval $3\leq<n>\leq$, the extreme has lower value than 3. In this case
the aim function with $<n>=3$ has minimum value in this interval. The state
$<n>=3$ corresponds to an hierarchical state.

##  8 The extreme is maximum

The extreme may be an absolute maximum for $A_{0}>2b$. In this case the state
$<n>=6$ has always lower value of the aim function than the state with
$<n>=3$.

The extreme (11) is maximum if, for $A_{0}>4b$. :

$EV>2E\frac{A_{0}(A_{0}-2b)}{A_{0}-4b}$ (17)

and then there exist two states which minimize function FM /e.i. maximize
return and minimize risk/.

The first one is with $<n>=3$. It gives the average number 3 of agents per
nonreducibile subgroup, this number means that there is a pyramidal
/hierarchical/ structure.

The second one is with $<n>=6$. It gives the average number 6 of agents per
nonreducibile subgroup, the value of the aim function FM is such that this
state is much more stable than the state above /the state corresponding to
$<n>=3$/ which is a metastable state. The number $<n>=6$ means that there is a
one cell structure. Efficient capital markets are markets which are of the
second type /strong stability, conservative and aggressive agents are present
on the market, return is high/ in our paper.

The extreme (11) is maximum if, for $A_{0}>4b$. :

$2E\frac{A_{0}(A_{0}-2b)}{A_{0}-4b}>EV>2E\frac{(A_{0}-4b)(A_{0}-2b)}{A_{0}}$
(18)

however then there exists one state which minimize function FM /e.i. maximize
return and minimize risk/ for $3\leq<n>\leq 6$. The state with $<n>=3$, which
was a metastable state is now not the metastable state. The minimum od the aim
function FM is for $<n>=6$. Efficient capital markets are markets which are of
this type.

The extreme (11) is maximum if, for $4b>A_{0}>2b$. :

$EV<2E\frac{A_{0}(A_{0}-2b)}{A_{0}-4b}$ (19)

This inequality cannot be fulfilled because right hand side is negative for
positive expected returns E and for positive expected risk EV.

The extreme (11) is maximum if, for $4b>A_{0}>2b$. :

$EV>2E\frac{A_{0}(A_{0}-2b)}{A_{0}-4b}$ (20)

This inequality is fulfilled always. Thus the state with $<n>=6$ is now a
state with minimum aim function. There is no metastable state in this case.

##  9 Aboav relation for the case of capital markets

Aboav relation [4] for the case of capital markets tells us:

$n.m(f,n)=5.f-11-K.(f-1-n),$ (21)

where K is a parameter of the group independent of f. Aboav’s law describes
how many, $/m(f,n)/$, agents are present in average in a nonreducibile
subgroup in a cell neighbouring to a cell with f nonreducibile subgroups with
agent average number n.

Variance of n for the distribution of agents in nonreducibile subgroups of a
cell with f nonreducibie subgroups is due to Weaire [4]:

$<(n-<n>)^{2}>=<n^{2}>+<n>^{2}$ (22)

This variance may be calculated from m(f,n) which is equal to n [4] to Weaire
and we obtain:

$<n^{2}>=<n.m(f,n)>=5.<f>-11-K.(<f>-1-<n>)$ (23)

$=\frac{60}{6-<n>}-11-K.(\frac{12}{6-<n>}-1-<n>)$.

We see that for states in which $<n>$ tends to 6, the variance diverges as
$\frac{12.(5-K)}{6-<n>}$. Note that $-1\leq K\leq 2$, [3] and [4].

Note that structures /stable states/ in both examples on capital market are
such that their entropy /informational/ is maximum /risk is minimum/ and their
return is maximum minimizing the function FM.

##  10 Conclusions

We studied capital market and social behaviour of agents. Quantity F, and in
general $<f>$, is the quantity which characterizes verification of information
and thus contributes to aversion of an agent with respect to a risk. We
generalized this relation to every cell: we assumed that personal diameter r
increases when $<f>$, the average number of nonreducubile groups per cell,
decreases, $<f>$ is an average number of faces per cell. Thus we assumed that
aversion of an agent to risk on capital markets also contains characterization
of verification of information. We assumed that the number of nonreducible
subgroups F per cell, $<f>$, contributes linearly to this aversion constant A.

When the pyramidal structure appears in the group, $<f>=4$, then $A=A_{0}$.
This type of structure corresponds to hierarchical economy systems. When the
structure contains topologically only one cell then $<n>=6$ and A tends to
infinity. Better verification of information leads to expectation of higher
returns and thus the acceptable risk is larger. This type of structure
corresponds to market-based economy systems. They work via price signals. We
interpret $A_{0}-4.b$ as a kind of measure of risk aversion of crowd and we
interpret $(<f>.b)$ as an average measure of risk aversion of a personal
/individual/ contribution to risk aversion different from the crowd
contribution of this person in the cell.

In practice both types of system involve interaction of people: in the market
system they interact as traders, in the hierarchy as agents within an
organization. Both market- and hierarchy- based systems co-exist in modern
economy: there are open markets and there are large organizations. This
corresponds to our state for which the aim function has its minimum for $<n>$
from the interval $3<<n><6$. In our model the market-based system corresponds
to a state with $<n>=6$ and the hierarchy-based to $<n>=3$, which are boundary
cases of this interval for values of $<n>$.

The mix of investments of an agent in a given cell with an average measure A
of risk aversion in the cell is found from extremes of the function FM, the
aim function. We find the minimum of the average /per cell/ function $<FM>$.
We assume that in our model the individual investors and agents from financial
intermediary firms do not differ. More general case will be studied elsewhere.

When $<f>=4$ then the pyramidal /hierarchical/ structure exists in which there
are only nonreducibie subgroups with $n=3$. Thus we consider this structure as
a structure with a uniform risk aversion constant $A_{0}$. Then a personal
contributions to the risk aversion constant of an agent are those
contributions which are corresponding to structures with $<n>>3$. Thus there
is a linear contribution to the risk aversion constant which is proportional
to $<f>-4$. This corresponds to nonuniform risk aversion constants in the
group of agents. This contains a personal contribution to the risk aversion
constant. We use $\frac{<f>.<n>}{3}$ as an average number of agents per cell.

The extreme of $<FM>$ is minimum for some conditions on the market, and in the
group of agents, see above. This minimum corresponds to $<n>$ between 3 to 6.
This is a stable state. The extreme may be a maximum. Then there are two
minima of the function FM, one of them is that which corresponds to $<n>=3$
/the pyramidal hierarchical structure/ which is a metastable state, and the
other one corresponds to $<n>=6$ /the structure with hexagons in average and
with one cell/ which is a stable state. Efficient capital markets are markets
which are of the second type /strong stability, conservative and aggressive
agents are present on the market, returns are high/ in our paper. Note that
under some conditions of agents and the expectations about the capital market
the hierarchy economy system may be more stable than the market-based economy
system. Aboav’s law describes how many, $/m(f,n)/$, agents are present in
average in nonreducibie subgroups in a cell neighbouring to a cell with f
nonredducibile subgroups with agent average number n. This enables us to
characterize variance of n for the distribution of nonreducibile subgroups.
This variance diverges to infinity with $<n>$ tending to 6.

The financial intermediaries play a key role to reduce the risks faced by
individual savers, by pooling their savings and using these to assemble
diversified portfolios of assets. This diversification requires specialized
knowledge and expertise, and a large portfolio. This will be studied using our
model with nonzero $\gamma$ constant corresponding to a firm with a given
tolerance risk constant.

## Acknowledgment

The paper represents a part of results of the VEGA project 1/0495/03.

## References

  * [1] T. Plummer, The Psychology of Technical Analysis, Rev. Ed., Probus Pub.Comp., Chicago-Cambridge, 1993
  * [2] D. Lewis, The Secret Language of Success, Carroll and Graf Pub. Inc., USA, 1989
  * [3] N. Rivier, Journal de Physique, C9 N12 T46 (1985) 155
  * [4] N. Rivier, Physica, 23D (1986) 129
  * [5] O. Hudak, Topology and Social Behaviour of Agents, http://arXiv.org/abs/cond-mat/0312723 , 2003
  * [6] I. Molho, The Economics of Information, Blackwell Pub., Oxford, Malden, 1997
  * [7] A. Grenspan, Commercial Banks and the Central Bank in a Market Economy, in Readings on Financial Institutions and Markets, P.S. Rose editor, 5th ed., R.D. Irwin Inc., Homewood, Boston, 1993, p. 294
  * [8] J. von Neumann and O. Morgenstern, Theory of Games and Economic Behaviour, 3rd. ed., Princeton University Press, Princeton, 1953
  * [9] J. Lintner, The Market Price of Risk, Size of Market and Investor s Risk Aversion, Journal of Business, April (1968)
  * [10] W.F. Sharpe, Integrated Aset Allocation, Financial Analysts Journal, September -October (1987)