# Analysis of bank leverage via dynamical systems and deep neural networks

Fabrizio Lillo Dipartimento di Matematica, Universitá di Bologna and Scuola
Normale Superiore, Pisa, Italy. Email address<EMAIL_ADDRESS>Giulia
Livieri Scuola Normale Superiore, Pisa, Italy. Email address:
<EMAIL_ADDRESS>Stefano Marmi Scuola Normale Superiore, Pisa, Italy.
Email address<EMAIL_ADDRESS>Anton Solomko Scuola Normale Superiore,
Pisa, Italy. Email address<EMAIL_ADDRESS>Sandro Vaienti Aix
Marseille Université, Université de Toulon, CNRS, CPT, 13009 Marseille,
France. Email address<EMAIL_ADDRESS>
###### Abstract

We consider a model of a simple financial system consisting of a leveraged
investor that invests in a risky asset and manages risk by using Value-at-Risk
(VaR). The VaR is estimated by using past data via an adaptive expectation
scheme. We show that the leverage dynamics can be described by a dynamical
system of slow-fast type associated with a unimodal map on $[0,1]$ with an
additive heteroscedastic noise whose variance is related to the portfolio
rebalancing frequency to target leverage. In absence of noise the model is
purely deterministic and the parameter space splits in two regions: (i) a
region with a globally attracting fixed point or a 2-cycle; (ii) a dynamical
core region, where the map could exhibit chaotic behavior. Whenever the model
is randomly perturbed, we prove the existence of a unique stationary density
with bounded variation, the stochastic stability of the process and the almost
certain existence and continuity of the Lyapunov exponent for the stationary
measure. We then use deep neural networks to estimate map parameters from a
short time series. Using this method, we estimate the model in a large dataset
of US commercial banks over the period 2001–2014. We find that the parameters
of a substantial fraction of banks lie in the dynamical core, and their
leverage time series are consistent with a chaotic behavior. We also present
evidence that the time series of the leverage of large banks tend to exhibit
chaoticity more frequently than those of small banks.

Keywords: leverage cycles, risk management, systemic risk, random dynamical
systems, unimodal maps, Lyapunov exponents, neural networks.

2020 Mathematics Subject Classiﬁcation: primary 91G80, 34F05; secondary 37H15,
62M45.

###### Contents

  1. 1 Introduction
  2. 2 From the structural model to the dynamical system
  3. 3 The mathematical model
    1. 3.1 Unimodal maps
    2. 3.2 Markov chains
    3. 3.3 Coupling with a stochastic process
    4. 3.4 The leverage model
    5. 3.5 Random transformations
  4. 4 Stationary measure
    1. 4.1 Existence
    2. 4.2 Uniqueness
  5. 5 Stochastic stability
  6. 6 Lyapunov exponent
    1. 6.1 Average Lyapunov exponent
    2. 6.2 Continuity of the Lyapunov exponent
    3. 6.3 Random Lyapunov exponent and random entropy
  7. 7 Numerical results
    1. 7.1 Bifurcation diagram
    2. 7.2 Lyapunov exponent
  8. 8 Estimating the map parameters via deep neural networks
  9. 9 Chaos in real bank leverage time series
    1. 9.1 Data set
    2. 9.2 Estimation via neural networks
    3. 9.3 Classification via the Chaos Decision Tree Algorithm
  10. 10 Conclusions
  11. 11 Acknowledgements
  12. A The Chaos Decision Tree Algorithm
    1. A.1 Simulations

##  1 Introduction

Leverage is one of the most important and controversial concepts in finance.
On one side, borrowing is essential in many economic activities, while, on the
other, it is intrinsically connected with risk. Concerning this last point,
the pro-cyclical nature of leverage has been highlighted and studied in the
recent literature (see, e.g., [30, 2, 35, 3, 1, 50], among many others).
Specifically, [2] and [3] argued that when assets are evaluated at fair value,
an increase in market prices of assets decreases the so-called “quasi-market
leverage ratio” – roughly the ratio of total assets to equity capital – and
this leaves room to build up debt for banks that operate through leverage or
Value at Risk111The Value at Risk measures the maximum loss that an asset
portfolio may suffer over a specific horizon and with a given level of
confidence. (VaR) targeting. The empirical findings in [2] and [3] show that
cycles of expansion (contraction) in the banks’ balance sheet size go hand in
hand with increases (decreases) in leverage (i.e., leverage is pro-cyclical);
a behavior that has been witnessed since the 1960s but exacerbated during
2007–2009 financial crisis. The creation of negative externalities in
financial markets because VaR model’s widespread use has been put forward in
e.g., [51, 26, 41]. It is shown that it can create market instability and
result in what has been called by [27] as endogenous risk, that is, the
systemic risk caused and amplified by the system itself rather than being the
result of an exogenous shock. Indeed, because of the imposed VaR capital
requirements, banks are forced to reduce their positions when the risk exceeds
these limits. Since the VaR of the trading portfolio increases when the
volatility goes up, banks are forced to reduce their positions rapidly, and
because of these fire sales, the price can drop abruptly. This leads to a new
drop in prices and likely increases volatility, which triggers (further)
portfolio reductions. This mechanism creates an exceptionally threatening
environment if many banks hold similar positions and use the same VaR model to
manage their risk since they are forced to sell the same assets
contemporaneously, leading to a destabilizing spiral. Finally, to implement
the VaR constraint (as well as any risk management mechanism), banks must
estimate both the riskiness of the investments and their dependencies. Since
the estimation of risk is typically done using historical data, additional
feedback is created between past and future risks, creating new threats for
the systemic stability of financial markets.

Understanding and modeling the leverage dynamics is therefore of paramount
importance. Some recent papers [8, 25, 44] have proposed stylized agent-based
models of financial institutions that lead to a dynamical system evolution for
the leverage. By using numerical methods, they also show that in some
parameters regime the leverage dynamics becomes chaotic via a cascade of
period doubling bifurcations222Other works that analyze systemic risk problems
through tools from dynamical system theory are, e.g. [23, 24, 52, 20, 9], to
cite only a few.. The map describing some of these models also contains a
noise component, describing the uncertainty naturally present in financial
institutions’ decision process. For example, optimal leverage depends on
portfolio risk, which is typically estimated with statistical methods on past
observations and these estimations are naturally modeled as random variables.
In this paper, we address three important questions related to this approach
(to be discussed below): (1) What are the mathematical properties (existence
and uniqueness of stationary measure, stochastic stability, Lyapunov exponent,
etc.) of the noisy deterministic maps emerging from these models? (2) Is it
possible to reliably estimate the map parameters from short time series as
those in publicly available datasets of banks balance sheets? (3) Is there
evidence of chaotic behavior in the leverage dynamics of real banks?

We consider a simple agent-based financial system where the mechanisms
described above are present. Our starting point is a simplified version of the
model proposed in [44], which in turn builds on [25], to the case of one bank
and one asset. As [44] shows numerically, the dynamics describing the
financial leverage displays a period-doubling bifurcation cascade resulting in
chaotic behaviour (as measured with the computation of the Lyapunov exponent).
This happens either when the parameter of the VaR constraint or the memory
used by the banks to estimate volatility from past data vary. These findings,
while suggestive, are however not rigorously proved in [44]. More importantly,
the mathematical properties of this family of models, i.e. where a
deterministic map is perturbed by a heteroscedastic additive noise (arising
from the coupling with a faster random dynamics), are not known in general. In
the simplified model, the system is composed of a representative leveraged
investor (a bank) that invests in a risky asset; the bank’s risk management
consists of two components. First, the bank estimates the future volatility
(the risk) of its investment in the risky asset using past market data.
Second, the bank uses the estimated volatility to set its desired leverage.
However, the bank faces a Value-at-Risk (VaR) capital requirement policy which
implies a constraint for the financial leverage
$\lambda_{t}:={A_{t}}/{E_{t}}$, where $A_{t}$ indicates the assets of the bank
at time $t$, whereas $E_{t}$ the liabilities. The bank is allowed a maximum
leverage $\bar{\lambda}_{t}$ which is a function of its own perceived risk.
The representative bank updates its expectation of risk at time intervals of
unitary length, say $(t,t+1]$ with $t\in\mathbbm{Z}$, and, accordingly, it
makes new decisions about the leverage. This process defines the slow
component of the model. Moreover, the model assumes that over the unitary time
interval $(t,t+1]$ representative bank re-balances its portfolio to target the
leverage without changing the risk expectations. The re-balancing takes place
in $\mathbbm{n}$ time sub-intervals within $(t,t+1]$. The time scale
$1/\mathbbm{n}$, with $\mathbbm{n}\in\mathbbm{N}$, characterizes the fast
component of the model. In particular, the slow variables evolve in time as a
function of averages over the fast variables. In summary, the considered model
is a discrete-time slow-fast dynamical system ([29] and [14]). Starting from
this model, in this paper we make three main contributions.

First, we show that the dynamics of leverage in our model follows, under
suitable approximation, a deterministic unimodal map on $[0,1]$ perturbed with
additive and heteroscedastic noise. The variance of the noise is related to
the frequency of portfolio rebalancing to target leverage. In particular, the
parameter space of this deterministic map has two regions: (i) a region where
the map has a globally attracting fixed point or a 2-cycle; (ii) the so-called
dynamical core region, where the map can exhibit chaotic behavior. In order to
study the mathematical properties of the map rigorously, we consider a more
general class of maps and describe the leverage dynamic utilizing a Markov
chain parametrized by the rebalance time $\mathbbm{n}$; we will study the
regime of finite $\mathbbm{n}$, as well as the limit for
$\mathbbm{n}\to\infty$. Although the stability of Markov chains is relatively
well studied (see e.g., [17] or [46] and references therein), some specific
properties of the stochastic kernel that defines our model do not allow us to
apply general results available; for instance, we do not know if our chain is
Harris recurrent. We instead exploit the unimodal dynamics of the
deterministic map. Perturbations of unimodal maps with uniform additive noise
were studied in [12], [13]. As far as we know, the Markov chains with the kind
of heteroscedastic noise we introduce are apparently new; see [32] for another
type of heteroscedastic nonlinear autoregressive process applied to financial
time series. To handle with them, we look at the Markov operator’s spectral
properties on suitable Banach spaces and prove the quasi-compactness of such
an operator. This allows us to obtain many rigorous results. In particular, we
get finitely many stationary measures with bounded variation densities. The
stationary measure’s uniqueness is achieved when the chain perturbs the
unimodal map which is either topologically transitive or admits an attracting
periodic orbit; such maps correspond to a major part of the parameter space.
From a financial point of view, should the stationary measure not be unique,
it would imply that, depending on the initial conditions, different banks
could experience completely different dynamics, corresponding to different
stationary measures. When this occurs in physical systems one speaks of phase
transitions and of coexistence of different mutually singular states. In our
case, this would imply, for instance, that policy measures could not be
universal. We also show the weak convergence of the unique stationary measure
to the invariant measure of the unimodal map. We point out that this step is
particularly delicate since the stochastic kernel becomes singular in the
limit of large $\mathbbm{n}$. It is well known that given a continuous Markov
chain which perturbs a given map $T$ as in our setting, one could construct a
sequence of random transformations close to $T$ and therefore replace the
deterministic orbit of $T$ with a random orbit given by the concatenation of
the maps randomly chosen in the sequence. This construction is formally
possible under general assumptions, but it is challenging to get
“representations by special classes of transformations” as Y. Kifer pointed
out in [37]. We provide an explicit construction of those transformations and
show their closeness with the unimodal map; we believe this inference from the
Markov chain to the random transformations is interesting and illustrates very
well the way the Markov chain moves randomly the states of the system. In
particular, our random maps are obtained by adding to the unimodal map $T$ an
additive term which also depends on the state variable, and this motivates the
attribute of heteroscedastic we gave to our noise. Once we dispose of the
stationary measure, we can define an average Lyapunov exponent by integrating
the logarithm of the derivative of the unimodal map $T$ with respect to such a
measure; this definition is suitable for the Markov chain approach. By
switching to random transformations, one can define the Lyapunov exponent of
the cocycle. We show numerically that the two approaches are asymptotically
equivalent and prove that the average exponent converges to the Lyapunov
exponent of $T$ in the limit of large $\mathbbm{n}$, as a consequence of the
stochastic stability. We finally show that the average Lyapunov exponent
depends continuously on the Markov chain parameters, and relate it to the
different chaotic behavior of the unperturbed unimodal map.

One of the purposes of the present work was to rigorously establish the
possibility of chaotic behaviour in leverage time series of banks, as well as
to detect it in financial datasets. For this reason we started from the
parametric slow-fast model of [44] and looked at the corresponding Lyapunov
exponent. We remind that for deterministic systems, the Lyapunov exponents
characterize the divergence of nearby orbits, allowing to distinguish between
regular and chaotic dynamics. In the presence of the noise, our model becomes
intrinsically stochastic, in particular the Markov chain will mix
exponentially fast, see section 6.3. Nevertheless the average Lyapunov
exponent still allows us to discriminate periodic and chaotic behaviours: it
is negative when we perturb a contracting map (and then the realizations of
the process fluctuate around the fixed point), and it becomes positive by
perturbing the dynamical core region. In both cases the distribution of the
realizations in the state space is governed by the unique absolutely
continuous stationary measure. We prove that both behaviours, chaotic and
periodic, are present in the parameter space and we actually observe both
regimes in the real data. This justifies our claims about chaotic dynamics of
the leverage.

Indeed, in our third contribution (to be detailed below) we will detect
possible chaos in banks’ leverage dynamics. As said, for the stochastic
stability and the Lyapunov exponent we mainly discuss the relation between
finite $\mathbbm{n}$ (thus noisy system) and infinite $\mathbbm{n}$ (purely
deterministic system). Thus our results indicate in which sense what we learn
for a noisy system is informative about the deterministic backbone. In our
empirical analysis we will do not study or use directly neither the stationary
measure nor the Lyapunov exponent (mainly because we have very short time
series), however the ”continuity” we observe from finite to infinite
$\mathbbm{n}$ suggests that the properties we observe empirically for finite
$\mathbbm{n}$ are informative of the underlying deterministic dynamics.

The paper’s second contribution concerns a methodology to infer the parameters
of the noisy map from short empirical time series. It is indeed interesting to
ask if there is evidence of chaotic behavior in the leverage dynamics of real
banks. We claim that applying the maximum likelihood estimation is not
feasible for two reasons. First, the likelihood function is highly non-convex
so that standard optimization methods may perform poorly. Second, although the
likelihood function for the process itself can be written explicitly, it may
happen that in many cases we observe only a certain iterate of the process,
e.g., we observe only one slow time scale portfolio decision event out of two.
Therefore we propose to use a powerful deep learning technique known as
Convolutional Neural Network (CNN) ([40]) to estimate the parameters of the
map. More precisely, our CNN takes as input the one-dimensional time-series
and gives as output the map’s corresponding parameters. We train the CNN via
extensive simulations of the model, considering different regions of the
parameter space. The robustness of the trained model is validated by its
prediction of parameters performance in a huge testing set. Results show the
merit of using our proposed CNN architecture to estimate the parameters.
Importantly, being based only on simulations, the NN-based approach can also
work for partial observations. Without attempting to review the literature of
parameter estimation of the dynamical system via NN, we only provide here a
few key references to position our contribution. [47] employed multi-layer NN
and recurrent networks to identify and control nonlinear deterministic
dynamical systems. Artificial NN have been used in [42] and [53] in a similar
framework. Batch and recursive prediction error estimation algorithm have been
derived for a NN model with a single hidden layer in [22] and [21] for the
identification of noisy discrete-time nonlinear dynamical systems.

The third contribution is the empirical analysis of real banks leverage time
series. Assuming the proposed unimodal map with heteroscedastic noise as data
generating process for the banks’ leverage, we estimate the parameters on
quarterly data of about $5,000$ US Commercial Banks provided by the Federal
Financial Institutions Examination Council (FFIEC) via the proposed CNN
architecture. We have at our disposal a time period going from March 2001 to
December 2014, for a total of 59 quarters. Remarkably, we find that the
parameters of a sizable fraction of banks lie in the map’s dynamical core and
that the large banks’ leverage tends to be more chaotic than one of the small
ones. As a robustness check, identifying chaotic/periodic behavior is tested
by following a non-parametric approach. In this latter case, the map is not
specified and estimators of indicators (such as the Lyapunov exponent [59]),
which assume different values in the two regimes, are used to discriminate
them from a finite length time series. We use a very recent algorithm dubbed
Chaos Decision Tree Algorithm [57] which combines several tools into an
automated processing pipeline that can detect the presence (or absence) of
chaos in noisy recordings, even for difficult edge cases. We apply the Chaos
Decision Tree Algorithm to our data set. Remarkably, the results corroborate
the CNN approach’s findings concerning the chaotic behavior for a significant
subset of typically large banks.

Outline of the paper. In Section 2 we present the financial model of a
representative bank managing its leverage. We show that the model leads to a
slow-fast deterministic-random dynamical system which can be recasted into a
unimodal deterministic map with heteroscedastic noise. In order to analyze it
rigorously, in Section 3 we recall some facts about unimodal maps and Markov
chains and then define the class of chains that we study. We also represent
our model in terms of random transformations. In Section 4 and 5 we show the
existence and uniqueness of an absolutely continuous stationary measure and
establish its convergence to the invariant measure of the deterministic map.
This allows us to define the Lyapunov exponent and prove its continuity with
respect to the model parameters in Section 6. We also discuss chaotic
indicators naturally arising from the random maps representation of the
process. The last part of the paper concerns numerical and empirical analyses.
Specifically, Section 7 presents some numerical investigations of the
bifurcation diagram and Lyapunov exponent of the map. Section 8 proposes an
estimation method of the map based on the use of deep neural networks and
Section 9 presents an empirical application to a large set of leverage time
series of US banks, showing evidence of chaotic behavior. Finally, in Section
10 we draw some conclusions and outline some potential extensions of our work.

##  2 From the structural model to the dynamical system

The stylized model of the leverage dynamics we are going to present is a
special case of the model of [44] (which in turn builds on [25]) restricted to
the case of a single (representative) financial institution and of a single
investment asset. We present below this model and show how, under suitable
approximation, the resulting dynamics of leverage follows a deterministic map
with additive and heteroscedastic noise. The mathematical properties of such
map are studied in the next sections.

In the model, a representative financial institution (hereafter a bank) takes
investment decisions at discrete times $t\in{\mathbb{Z}}$ (the slow time
scale). At each time the bank’s balance sheet is characterized by the asset
$A_{t}$ and equity $E_{t}$, which together define the leverage
$\lambda_{t}:=A_{t}/E_{t}$. The bank wants to maximize leverage (by taking
more debt) to increase profits, but regulation constrains the bank’s Value-at-
Risk (VaR) in such a way that

$\lambda_{t}=\frac{1}{\alpha\sigma_{e,t}},$

where $\alpha$ depends on the return distribution and VaR constraint333For
example, if returns are Gaussian and the probability of VaR is 5%, it is
$\alpha=1.64$.. $\sigma_{e,t}$ is the expected volatility at time $t$ of the
asset, which in this simple model is composed by a representative risky
investment. Thus at each time $t$ the bank recomputes $\sigma_{e,t}$ and
chooses $\lambda_{t}$. Then in the interval $[t,t+1]$ the bank trades the
risky investment to keep the leverage close to the target $\lambda_{t}$. The
trading process occurs on the points of a grid obtained by subdividing
$[t,t+1]$ in $\mathbbm{n}$ subintervals of length $1/\mathbbm{n}$ (the fast
time scale). The dynamics of the investment return can be written as

$r_{t+k/\mathbbm{n}}=\varepsilon_{t+k/\mathbbm{n}}+e_{t+(k-1)/\mathbbm{n}},\quad
k=1,2,\ldots,\mathbbm{n},$ (1)

where $\varepsilon_{t+k/\mathbbm{n}}$ and $e_{t+(k-1)/\mathbbm{n}}$ are,
respectively, the exogenous and endogenous component of the return. The former
is a white noise term with variance $\sigma^{2}_{\epsilon}$, while the latter
depends on the banks’ demand for the risky investment in the previous step.
For each bank, the demand for the risky investment at time $t+k/\mathbbm{n}$
is the difference between the target value of $A_{t}$ to reach $\lambda_{t}$
and its actual value. Since the bank’s asset is composed by the risky
investment, an investment return $r_{t+k/\mathbbm{n}}$ modifies $A_{t}$ and
the bank trades at each grid point to reach the target leverage. It is
possible to show (see [25, 44]) that to achieve this, at each time
$t+k/\mathbbm{n}$ the bank’s demand for the risky investment is

$D_{t+k/\mathbbm{n}}=(\lambda_{t}-1)A^{*}_{t+(k-1)/\mathbbm{n}}r_{t+k/\mathbbm{n}},$

where $A^{*}_{t+(k-1)}$ is the target asset size in the previous step. If
there are $M$ identical banks, the aggregated demand is
$MD_{t+k/\mathbbm{n}}$. The endogenous component of returns
$e_{t+k/\mathbbm{n}}$ is determined by the aggregated demand by the equation

$e_{t+k/\mathbbm{n}}=\frac{1}{\gamma}\frac{MD_{t+k/\mathbbm{n}}}{C_{t+k/\mathbbm{n}}},$

where $C_{t+k/\mathbbm{n}}=MA^{*}_{t+(k-1)/\mathbbm{n}}$ is a proxy of the
market capitalization of the risky asset and $\gamma$ is a parameter measuring
the investment liquidity. Using the above expression, it is

$e_{t+k/\mathbbm{n}}=\frac{\lambda_{t}-1}{\gamma}e_{t+(k-1)/\mathbbm{n}}=\phi_{t}e_{t+(k-1)/\mathbbm{n}}$

and thus in the period $[t,t+1]$ the return $r_{t+k/\mathbbm{n}}$ follows an
AR(1) process with autoregression parameter $\phi_{t}=(\lambda_{t}-1)/\gamma$
and idiosyncratic variance $\sigma^{2}_{\epsilon}$.

To close the model, we specify how the bank forms expectations $\sigma_{e,t}$
on future volatility at time $t$. We assume that bank uses adaptive
expectations, which implies that

$\sigma^{2}_{e,t}=\omega\sigma^{2}_{e,t-1}+(1-\omega)\hat{\sigma}^{2}_{e,t},$

where $\omega\in[0,1]$ is a parameter weighting between the expectation at
$t-1$ and the estimation $\hat{\sigma}^{2}_{e,t}$ of volatility obtained by
the return data in $[t-1,t]$. As done in practice, this is obtained by
estimating the sample variance of the returns in $[t-1,t]$, i.e.

$\hat{\sigma}^{2}_{e,t}=\widehat{\text{Var}}\left[\sum_{k=1}^{\mathbbm{n}}r_{t-1+k/\mathbbm{n}}\right]\\\
=\left(1+2\frac{\hat{\phi}_{t-1}(1-\hat{\phi}_{t-1}^{\mathbbm{n}})}{1-\hat{\phi}_{t-1}}-2\frac{(\mathbbm{n}\hat{\phi}_{t-1}-\mathbbm{n}-1)\hat{\phi}_{t-1}^{\mathbbm{n}+1}+\hat{\phi}_{t-1}}{\mathbbm{n}(1-\hat{\phi}_{t-1})^{2}}\right)\frac{\mathbbm{n}\hat{\sigma}_{\epsilon}^{2}}{1-\hat{\phi}^{2}_{t-1}},$
(2)

where the last expression gives the aggregated variance of an AR(1) process as
a function of the AR estimated parameters $\hat{\phi}_{t-1}$ and
$\hat{\sigma}^{2}_{\epsilon}$. In the following we will assume that these are
the MLE estimators. We remind that when $\mathbbm{n}$ is large,
$\hat{\phi}_{t-1}$ is a Gaussian distributed variable with mean $\phi_{t-1}$
and variance $(1-\phi^{2}_{t-1})/\mathbbm{n}$.

In conclusion, the leverage dynamics is described by the following equations:

$\begin{cases}&\lambda_{t}=\left(\omega\frac{1}{\lambda^{2}_{t-1}}+(1-\omega)\alpha^{2}\widehat{\text{Var}}[\sum_{k=1}^{\mathbbm{n}}r_{t-1+k/\mathbbm{n}}]\right)^{-1/2},\\\
&r_{s}=\phi_{t-1}r_{s-1/\mathbbm{n}}+\epsilon_{s},\qquad
s=t-1+k/\mathbbm{n},\quad k=1,2,\ldots,\mathbbm{n},\end{cases}$ (3)

Since slow variables evolve depending on averages of the fast variables, the
model is a slow-fast deterministic-random dynamical system. By using the
expression above for the variance, we can rewrite the equation for the slow
component only as

$\lambda_{t}=\left(\omega\frac{1}{\lambda^{2}_{t-1}}+(1-\omega)\alpha^{2}\hat{\sigma}^{2}_{e,t}\right)^{-1/2},$

where the estimator $\hat{\sigma}^{2}_{e,t}$ can be seen as a stochastic term
depending on $\lambda_{t-1}$ and whose variance goes to zero when
$\mathbbm{n}\to\infty$.

If $\mathbbm{n}$ is large, the above map reduces to

$\lambda_{t}=\left(\omega\frac{1}{\lambda^{2}_{t-1}}+\frac{(1-\omega)\alpha^{2}\mathbbm{n}\hat{\sigma}_{\epsilon}^{2}}{(1-\hat{\phi}_{t-1})^{2}}\right)^{-1/2},$

and using the relation $\phi_{t}=\frac{\lambda_{t}-1}{\gamma}$, the map
becomes

$\phi_{t}=-\frac{1}{\gamma}+\frac{1}{\gamma}\left(\frac{\omega}{(1+\gamma\phi_{t-1})^{2}}+\frac{(1-\omega)\alpha^{2}\mathbbm{n}\hat{\sigma}_{\epsilon}^{2}}{(1-\hat{\phi}_{t-1})^{2}}\right)^{-1/2}.$
(4)

When changing $\mathbbm{n}$ also $\sigma^{2}_{\epsilon}$ changes, since the
AR(1) can be seen as the discretization of a continuous time stochastic
process (namely an Ornstein-Uhlenbeck process). A simple scaling argument
shows that the quantity $\Sigma_{\epsilon}=\sigma^{2}_{\epsilon}\mathbbm{n}$
is instead constant and independent from the discretization step
$1/\mathbbm{n}$. In the limit $\mathbbm{n}\to\infty$, it is
$\hat{\phi}_{t}\to\phi_{t}$, thus the above map becomes purely
deterministic444This is the deterministic skeleton, whose properties are
discussed in detail in [44].. The map in this case has a fixed point
$\phi^{*}=\frac{1-\alpha\sqrt{\Sigma_{\epsilon}}}{1+\alpha\gamma\sqrt{\Sigma_{\epsilon}}}$.
By replacing this condition in (4) and assuming that the risky asset is very
liquid ($\gamma\gg 1$), the map becomes

$\phi_{t}\simeq\left(\frac{\omega}{\phi^{2}_{t-1}}+\left(\frac{1-\phi^{*}}{\phi^{*}}\right)^{2}\frac{1-\omega}{(1-\hat{\phi}_{t-1})^{2}}\right)^{-1/2}.$
(5)

Since in the large $\mathbbm{n}$ limit the ML estimator $\hat{\phi}_{t-1}$ is
a Gaussian variable with mean $\phi_{t-1}$ and variance
$(1-\phi^{2}_{t-1})/\mathbbm{n}$, we can write

$\hat{\phi}_{t-1}=\phi_{t-1}+\eta_{t-1},$

where $\eta_{t-1}\sim{\mathcal{N}}(0,(1-\phi^{2}_{t-1})/\mathbbm{n})$. If the
noise $\eta_{t-1}$ is small (i.e. $\mathbbm{n}$ is large), we can perform a
series expansion, obtaining

$\phi_{t}\simeq\frac{|\phi_{t-1}(1-\phi_{t-1})|}{\sqrt{b\phi^{2}_{t-1}+\omega(1-\phi_{t-1})^{2}}}(1+\zeta_{t-1}),$
(6)

where $b=b(\omega,\phi^{*})$ is given by

$b=(1-\omega)\left(\frac{1-\phi^{*}}{\phi^{*}}\right)^{2}$ (7)

and the noise term

$\zeta_{t-1}=\frac{-b\phi^{2}_{t-1}}{(1-\phi_{t-1})(b\phi^{2}_{t-1}+\omega(1-\phi_{t-1})^{2})}\eta_{t-1}.$

In this approximation, the map can be seen as deterministic map with additive,
but heteroscedastic, noise

$\phi_{t+1}=T(\phi_{t};\theta)+\sigma(\phi_{t};\theta)\epsilon_{t},$ (8)

where $\epsilon_{t}\sim{\cal N}(0,1)$ and $\theta$ is a vector of parameters.
In our setting $\theta=(b,\omega,\mathbbm{n})$ and the deterministic map $T$
does not depend on $\mathbbm{n}$. Specifically,

$T(\phi_{t};\theta)=\frac{|\phi_{t}(1-\phi_{t})|}{\sqrt{b\phi^{2}_{t}+\omega(1-\phi_{t})^{2}}}$
(9)

and

$\sigma(\phi_{t};\theta)=\frac{b\phi_{t}^{3}\sqrt{1-\phi_{t}^{2}}}{\sqrt{\mathbbm{n}}\bigl{(}b\phi^{2}_{t}+\omega(1-\phi_{t})^{2}\bigr{)}^{3/2}}.$
(10)

The map $T$ is represented in Fig. 1. The term under square root is always
positive, since so are $\omega$ and $b$.

In the next sections we develop a rigorous mathematical theory of additive
maps with heteroscedastic noise as in (8). In particular, we will study the
existence of a stationary measure, the stochastic stability, and the Lyapunov
exponent for this class of models, having in mind our main example of (9) and
(10) coming from the financial application. It is worth noting, however, that
our results remain valid for any noise $\epsilon_{t}$ in (8), not only
Gaussian; see Section 3.3. Then in Section 8 we will present an estimation
method for the map and in Section 9 we will present the empirical analysis on
real leverage time series.

Figure 1: Plot of the deterministic component $T(\phi)$, $\phi^{*}=0.73$,
$\omega=0.4$ ($b=0.082$).

##  3 The mathematical model

As we explained in the previous section, we model the dynamics of a financial
system where two time scales are present: a slow time scale where risk
expectations (thus portfolios) are updated at unit times $t$, and a fast time
scale, between $t-1$ and $t$, during which banks rebalance portfolios
$\mathbbm{n}$ times. Since slow variables evolve depending on averages of the
fast variables, the model is a slow-fast deterministic-random dynamical
system. The evolution of the slow component is described by equations (5) and
(8). Notice that since the distribution of $\phi_{t}$ only depends on
$\phi_{t-1}$, both processes are (continuous state) Markov chains. In this
section, we define some tools that we later use to analyse the model.

###  3.1 Unimodal maps

Our process is constructed by perturbing with a heteroscedastic additive noise
a deterministic unimodal map $T$ of the unit interval $I=[0,1]$. We now
describe the class of unimodal maps $T$ we will consider, with (9) being its
representative.

We will refer to the class of unimodal maps studied in [13]; they enjoy a
series of ergodic properties which allow us to establish rigorous results for
the problem we deal with.

We therefore consider unimodal maps $T\colon I\to I$ of class $C^{4}$ with
$T(0)=T(1)=0$ and with a non-degenerate critical point555The critical point
for (9) is $c=\bigl{(}1+\sqrt[3]{b/\omega}\bigr{)}^{-1}$. at $c$:
$T^{\prime}(c)=0$. The map $T$ is strictly increasing on $[0,c)$ and strictly
decreasing on $(c,1]$. Moreover, we suppose that $T$ satisfies the following
assumptions:

  1. (A1)

$T$ has negative Schwarzian derivative:
$S(T)=\dfrac{T^{\prime\prime\prime}}{T^{\prime}}-\dfrac{3}{2}\left(\dfrac{T^{\prime\prime}}{T^{\prime}}\right)^{2}<0$,

  2. (A2)

$\Delta:=T(c)<1$,

  3. (A3)

the critical point is _quadratic_ : $T^{\prime\prime}(c)\neq 0$,

  4. (A4)

$|T^{k}(c)-c|\geq e^{-\alpha k}$ for all $k\geq k_{0}$,

  5. (A5)

$|(T^{k})^{\prime}(\Delta)|\geq\lambda_{c}^{k}$ for all $k\geq k_{0}$,

where $k_{0}\geq 1,1<\lambda_{c}<2$, and $\alpha>0$ with
$e^{2\alpha}<\sqrt{\lambda_{c}}$ are fixed constants.

Figure 2: Partition of the parameter space for the unimodal map (9) according
to the classification (C1)–(C3).

According to the theory of unimodal maps (see for instance the review [56]),
we have the following classification:

  1. (C1)

If $\Delta\leq c$, then there exists a globally attracting fixed point.

  2. (C2)

If $c<T(\Delta)<\Delta$, then there is a globally attracting fixed point or a
2-cycle in $(c,\Delta)$.

  3. (C3)

If $T(\Delta)<c<\Delta$, we can reduce the study to the so-called _dynamical
core_ $[T(\Delta),\Delta]$, which is mapped onto itself and absorbs all
initial conditions (except $0$ which is a fixed point).

Following this classification, we will say that $T$ is periodic if there is a
globally attracting fixed point or a globally attracting cycle. Recall that a
map $T$ on a topological space $X$ is called _topologically transitive_ if for
all nonempty open sets $U,V\subset X$, there exists $n$ such that $T^{-n}U\cap
V\neq\emptyset$. We will say that $T$ is chaotic if, in addition to above
assumptions (A1)–(A5),

  1. (At)

$T$ is topologically transitive on the interval $[T(\Delta),\Delta]$.

In particular, any such $T$ satisfies the dynamical core condition (C3).
Notice also that $T(x)>x$ for any $x<T(\Delta)$, in particular,
$T^{\prime}(0)>1$.

The invariant sets (attractors) of a unimodal map have a variety of
structures, as it is stated by the following theorem of Blokh and Lyubich (we
quote the statement given in [56, Theorem 6]):

######  Theorem A ([16]).

Let $T:I\to I$ be an S-unimodal map with nonflat critical point (S means
$S(T)<0$). Then $T$ has a unique metric attractor $A$, such that the
$\omega$-limit set $\omega(x)=A$ for Lebesgue almost all $x\in I$. The
attractor $A$ is of one of the following types:

  1. 1.

an attracting periodic orbit;

  2. 2.

a Cantor set of measure zero;

  3. 3.

a finite union of intervals with a dense orbit.

In the first two cases, $A=\omega(c)$.

Associated to $T$ there is the _transfer operator_ (also called the _Perron-
Frobenius operator_) $L\colon L^{1}\to L^{1}$ which is the positive linear
operator defined by the duality relation666Without mention of the contrary,
all the $L^{p}$ spaces in the paper will be intended with respect to the
Lebesgue measure. The latter will be denoted as $dx$ or Leb.

$\int_{I}Lf\,g=\int_{I}f\,g\circ T,\quad f\in L^{1},g\in L^{\infty}.$

In order to get useful information from this operator, we need to restrict the
functional space where it acts; we choose here the Banach space $BV$ of
bounded variation functions on the unit interval equipped with the complete
norm

$\|f\|_{BV}=|f|_{TV}+\|f\|_{1},$

where $|f|_{TV}$ is the total variation of the function $f\in L^{1}$. For a
chaotic map $T$ it follows that it admits a unique absolutely continuous
invariant measure $\nu=\nu\circ T^{-1}$ with a density $h\in BV$ and supported
on the interval $[T(\Delta),\Delta]$ ([13, Section 5, Corollary 1]). Moreover,
$\nu$ is mixing with exponential decay of correlations on BV observable,
namely there exists $0<v<1$ and a constant $C>0$ such that

$\Bigl{|}\int L^{n}f\,gdx-\int fdx\int gdx\Bigr{|}\leq
Cv^{n}\|f\|_{BV}\|g\|_{\infty},$ (11)

see [13, Section 5, Corollary 3] and [58, Proposition 5.15].

Before continuing, it is useful to quote a sort of analog of the theorem of
Blokh and Lyubich given above, for what concerns invariant measures for the
map $T$. We give here the statement of Theorem 9 in [56].

######  Theorem B (see [45, Chapter V.1]).

Let $T$ be an S-unimodal map with nonflat critical point. If $T$ has a
periodic attractor, or a Cantor attractor, then $T$ admits a unique SRB
measure777We remind that an invariant measure $\mu$ is called a Sinai-Ruelle-
Bowen (SRB) measure if
$\mu=\lim_{n\to\infty}\frac{1}{n}\sum_{k=0}^{n-1}\delta_{T^{k}(x)}$ for
$\operatorname{Leb}$-a.e. $x\in[0,1]$, where $\delta_{x}$ is the Dirac mass at
$x$. supported on the attractor.
If $T$ admits an absolutely continuous invariant probability measure $\mu$,
then:

  1. 1.

$\mu$ is a SRB measure;

  2. 2.

the attractor $A$ of $T$ is an interval attractor;

  3. 3.

$\operatorname{supp}\mu=A$, in particular, $\mu$ is equivalent to the Lebesgue
measure on $A$.

######  Remark 3.1.

As pointed out in [13], (A4) ensures that $T$ has no periodic attractors and
is ergodic with respect to the measure $\mu$, which is the unique absolutely
continuous invariant probability measure for $T$. And (At) allows to prove
that $T$ is Bernoulli and therefore mixing. Therefore chaotic map fulfills
conditions 1–3 of Theorem B above. As it is pointed out in [56, Section 4]:
“The theorem [above] does not guarantee the existence of a natural measure in
the case of an interval attractor. Indeed there are uncountably many
parameters in the logistic family, for which the corresponding maps have
interval attractors and lack natural measures altogether, or have natural
measures with weird properties… But, at least in the logistic family, both
such singular phenomena and Cantor attractors are rare in the sense of
Lebesgue measure”.

In Section 6 we will address the question to compute the Lyapunov exponent
$\Lambda$ for the map $T$. The following theorem by G. Keller can be applied
to our situation:

######  Theorem C ([36]).

Let $T\colon I\to I$ be an S-unimodal map with nonflat critical point. Then
$T$ admits an absolutely continuous invariant probability measure if and only
if $\lim_{n\to\infty}\frac{1}{n}\log|DT^{n}(x)|=\Lambda>0$, for almost all
$x\in I$.

Of course if $T$ has a periodic attractor, the Lyapunov exponent will be
negative. The situation is different whenever $T$ has a Cantor attractor. For
an S-unimodal map with nondegenerate critical point that also has a Cantor
attractor, the Lyapunov exponent will be $0$, while there are families of
unimodal maps with critical point of sufficiently high order, which have
Cantor attractors with sensitive dependence on initial conditions, see [56,
Section 5].

###  3.2 Markov chains

Recall that a Markov chain $\\{X_{t}:t\in\mathbb{N}\\}$ on the interval $I$ is
given by _transition probabilities_

$P_{x}(A)=\mathbb{P}\\{X_{t+1}\in A\mid X_{t}=x\\}$

(the probability that a chain at $x\in I$ will be in a set888 All sets
considered will be assumed to be measurable. For brevity’s sake, we omit the
word ‘measurable’ everywhere in this text.  $A\subset I$ after one step) and
an initial distribution $\rho_{0}(A)=\mathbb{P}\\{X_{0}\in A\\}$. In the
particular case where all $P_{x}$, $x\in I$, and $\rho_{0}$ are absolutely
continuous (with respect to $\operatorname{Leb}$) and are given by densities
$p(x,\cdot)$ and $h_{0}$, respectively, we have

$P_{x}(A)=\int_{A}p(x,y)dy,\quad\rho_{0}(A)=\int_{A}h_{0}(y)dy.$

The map $p\colon I\times I\to\mathbb{R}_{+}$ (known as the _stochastic
kernel_) plays the role that the transition matrix does in the theory of
Markov chains with a finite state space. For $P_{x}$ to be a probability, it
should satisfy $\int p(x,y)dy=1$ for every $x\in I$.

Denote with $\mathcal{M}$ the space of (real-valued) Radon measures on $I$.
There is an associated operator $\mathcal{L}\colon\mathcal{M}\to\mathcal{M}$
(called the _Markov operator_ corresponding to $P$) acting by

$\mathcal{L}\rho=\int P_{x}d\rho(x),\quad\rho\in\mathcal{M},$ (12)

that is, $\mathcal{L}\rho(A)=\int P_{x}(A)d\rho(x)$ for every $A\subset I$,
or, equivalently,

$\int\varphi d\mathcal{L}\rho=\iint\varphi(y)dP_{x}(y)d\rho(x)$

for all $\varphi\in C^{0}$, where $C^{0}$ denotes the Banach space of
continuous functions on $I$ with the $\sup$ norm. We note that
$\mathcal{L}\colon L^{1}\to L^{1}$ is an isometry, where $L^{1}$ is intended,
from now on, with respect to the Lebesgue measure. If the chain is given by
the kernel $p$, formula (12) restricted to $L^{1}$ becomes

$(\mathcal{L}h)(y)=\int p(x,y)h(x)dx,\quad h\in L^{1}.$ (13)

If $\rho_{t}$ denotes the distribution of the random variable $X_{t}$, then
the distribution of $X_{t+1}$ is $\rho_{t+1}=\mathcal{L}\rho_{t}$. In other
words, fixing the distribution $\rho_{0}$ for $X_{0}$, the entire sequence of
future distributions can be obtained by iterating with $\mathcal{L}$.

A measure $\mu\in\mathcal{M}_{+}$ is said to be _stationary_ if

$\mathcal{L}\mu=\mu.$

Every stationary measure $\mu$ gives rise to a shift-invariant measure
$\mathbb{P}_{\mu}$ on the sequences space
$\Omega=\\{(x_{t})_{t\in\mathbb{N}}:x_{t}\in I$} of realizations of the
process, such that $\mathbb{P}_{\mu}(x_{t}\in A)=\mu(A)$ for all
$t\in\mathbb{N}$ (see, e.g., [39]). We say that $\mu$ is _ergodic_ if
$\mathbb{P}_{\mu}$ is ergodic. In the next section we will show that, under
some mild conditions, our model admits a unique (and thus ergodic) absolutely
continuous stationary probability $\mu$ with a density $h\in BV$. Then, by the
Ergodic Theorem (see, e.g., [15, Remark C4.1]), for every $f\in L^{1}(\mu)$,

$\frac{1}{n}\sum_{t=1}^{n}f(X_{t})\xrightarrow[n\to\infty]{}\int
fd\mu,\quad\text{$\mathbb{P}_{\mu}$-almost surely}.$

In particular, realizations of the process are distributed in the state space
according to the measure $\mu$.

###  3.3 Coupling with a stochastic process

We are now ready to define a Markov chain that describes our model. It is
obtained as a deterministic unimodal map $T\colon I\to I$ satisfying
assumptions of Section 3.1, coupled with a stochastic process, namely, by
perturbing $T$ with additive noise. Starting from (8) as our main motivation,
on the one hand, we consider a more general class and, on the other, impose
some mild technical restrictions that are necessary for rigorous analysis.

Since the noise varies in a neighborhood of 0, we will need to extend the
state space on the negative axes. We will see in a moment, however, that such
an extension is irrelevant for the asymptotic behavior of the perturbed
system, whose random trajectories spend all the time, but a relatively short
transient, on the positive unit interval.

We fix $T$ and parametrize the chain by the rebalance time $\mathbbm{n}$
(which is roughly inverse to the variance of the noise), consequently indexing
with $\mathbbm{n}$ the chain $(X_{t}^{(\mathbbm{n})})$, the transition
probabilities $P_{x}^{(\mathbbm{n})}$ and the stochastic kernel
$p_{\mathbbm{n}}(x,y)$ where it is necessary. We will be interested in the
limit for $\mathbbm{n}\to\infty$.

We also need to assume that the noise is compactly supported, in order for
trajectories of the process to stay bounded. Compared to the Gaussian noise in
(8), this is done by truncating the distribution tails that are exponentially
small for large $\mathbbm{n}$, see Section 3.4 for the main example.

Denote by

$\Gamma:=1-\Delta$

the gap between $T(c)$ and $1$. We now extend the domain of definition of $T$
to the larger interval $[-\Gamma,1]$ (which, by abuse of notation, will be
still denoted by $I$) so that $T$ is continuous at $0$ and on $[-\Gamma,0)$ is
$C^{4}$ smooth, positive and decreasing, with $T(-\Gamma)<\Delta$.999A similar
extension was considered in [13] to allow perturbations with additive noise;
in particular, it was supposed that $T(I)\subset\operatorname{int}(I)$ and
that $T$ admits an extension to some compact interval $J\supset I$, preserving
all the previous properties and satisfying $T(\partial J)\subset\partial J$.
Notice that with this extension the map $T$ is not anymore of class $C^{4}$ as
prescribed at the beginning of Section 3.1, but this regularity still persist
on the interval $(0,1)$ and this will be enough for the next considerations.

To construct the chain, we need to define transition probabilities. Let
$g_{x,\mathbbm{n}}(y)$ be a probability density supported on a compact
interval $[-s(x),s(x)]$. We assume that

  * •

$0<s(x)<\Gamma/2$ for $x\in(0,1)$,

  * •

$T(x)-s(x)>0$ for $x\in(0,1-\Gamma/2]$,

  * •

$T(x)-s(x)>x$ for $x$ small (in particular, $T^{\prime}(0)>1$).

We set for simplicity $s(x):=0$ for $x\leq 0$ (meaning
$P_{x}^{(\mathbbm{n})}=\delta_{T(x)}$); Lemma 4.4 below shows that this choice
does not affect the dynamics. We will also assume that both the mean and the
variance of $g_{x,\mathbbm{n}}$ decrease to $0$ as $\mathbbm{n}\to\infty$ and,
for every $\varepsilon>0$,

$\sup_{x\in[\varepsilon,1-\varepsilon]}|g_{x,\mathbbm{n}}|_{TV}<\infty.$ (14)

Fix any initial distribution $\rho_{0}\in\text{BV}$ and define transition
probabilities

$P_{x}^{(\mathbbm{n})}(A):=\int_{-s(x)}^{s(x)}{\mathbbm{1}}_{A}(T(x)+y)g_{x,\mathbbm{n}}(y)dy,$
(15)

which correspond to the stochastic kernel

$p_{\mathbbm{n}}(x,y)=g_{x,\mathbbm{n}}(y-T(x)).$ (16)

Informally speaking, the probability that the chain steps from $x$ to $A$ will
be high whenever $T(x)$ falls in $A$. Equivalently, we can write

$X_{t+1}^{(\mathbbm{n})}=T(X_{t}^{(\mathbbm{n})})+Y_{t+1},\;\text{ where
}\;Y_{t+1}\sim g_{x,\mathbbm{n}}.$ (17)

The values of $X_{t+1}^{(\mathbbm{n})}$ are spread in a neighborhood of $T(x)$
due to the addition of the random variable $Y_{t+1}$.

###  3.4 The leverage model

We will slightly modify (8) to satisfy the technical assumptions listed above.
The unimodal map is

$T(x)=\frac{|x(1-x)|}{\sqrt{bx^{2}+\omega(1-x)^{2}}},$ (18)

where the parameters $\phi^{*}$ and $\omega$ are such that $T$ satisfies
assumptions (A1)–(A5). (see Section 3.1). It always has negative Schwarzian
derivative, as we verified numerically. The condition $T(\Delta)<c<\Delta<1$
defines a nonempty subset of parameters (see Fig. 2). Notice that
$T^{\prime}(0)=1/\sqrt{\omega}>1$.

We want $g_{x,\mathbbm{n}}$ to be (truncated) normal with a variance close to
(10). For this, let us denote by $\mathcal{N}_{a}(0,\sigma)$ the smoothed
truncated normal distribution with the density
$g(y)=c_{a,\sigma}\chi_{a}(y)e^{-\frac{y^{2}}{2\sigma^{2}}},$ where
$c_{a,\sigma}$ is so that $\int g(y)dy=1$ and $\chi_{a}$ is a smooth bump
function supported on $[-a,a]$.101010The smoothness of the truncation function
is only used in the proof of Theorem 6.5.  For instance, we may set

$\chi_{a}(y)=\begin{cases}1,&\mbox{if }|y|\leq(1-\varepsilon)a,\\\
\Psi(\frac{y\pm(1-\varepsilon)a}{\varepsilon a}),&\mbox{if
}(1-\varepsilon)a<|y|\leq a,\\\ 0,&\mbox{if }|y|>a,\end{cases}$

where $\Psi(t)=e^{1-\frac{1}{1-t^{2}}}$ is the standard $C^{\infty}$ bump
function on $[-1,1]$. We set

$\sigma_{\mathbbm{n}}(x):=\frac{bx^{3}\sqrt{1-x^{2}}}{\sqrt{\mathbbm{n}}\bigl{(}bx^{2}+\omega(1-x)^{2}\bigr{)}^{3/2}}.$
(19)

Denote $\sigma(x):=\sigma_{1}(x)$ and
$\sigma_{\max}:=\max_{x\in[0,1]}\sigma(x)$. Set

$g_{x,\mathbbm{n}}(y):=c_{x,\mathbbm{n}}\chi_{s(x)}(y)e^{-\frac{y^{2}}{2\sigma^{2}_{\mathbbm{n}}(x)}},$
(20)

where

$s(x):=\frac{\sigma(x)}{\sigma_{\max}}\min\Bigl{\\{}\frac{\Gamma}{2},T(1-\frac{\Gamma}{2})\Bigr{\\}}\quad\text{and}\quad
c_{x,\mathbbm{n}}:=\biggl{(}\int_{-s(x)}^{s(x)}\chi_{s(x)}(y)e^{-\frac{y^{2}}{2\sigma^{2}_{\mathbbm{n}}(x)}}dy\biggr{)}^{-1}.$

We can then rewrite (17) as

$X_{t+1}^{(\mathbbm{n})}=T(X_{t}^{(\mathbbm{n})})+\sigma_{\mathbbm{n}}(X_{t})Z_{t+1},\quad
Z_{t+1}\sim\mathcal{N}_{b_{\mathbbm{n}}}(0,1),$ (21)

with
$b_{\mathbbm{n}}:=\frac{\sqrt{\mathbbm{n}}}{\sigma_{\max}}\min\\{\frac{\Gamma}{2},T(1-\frac{\Gamma}{2})\\}\to\infty$
as $\mathbbm{n}\to\infty$. The fact that, for $\mathbbm{n}$ fixed,
$g_{x,\mathbbm{n}}$ are rescaled copies of the same distribution will be used
in Section 3.5 to explicitly describe random maps associated to the process.
We get the following stochastic kernel:

$p_{\mathbbm{n}}(x,y)=c_{x,\mathbbm{n}}\chi_{s(x)}(y-T(x))e^{-\frac{(y-T(x))^{2}}{2\sigma_{\mathbbm{n}}^{2}(x)}}.$
(22)

Notice also that the support of $p_{\mathbbm{n}}(x,y)$ does not depend on
$\mathbbm{n}$ (see Fig. 3).

Figure 3: Support of the kernel (22).

Finally, (14) holds, because $|g_{x,\mathbbm{n}}|_{TV}=2c_{x,\mathbbm{n}}$ and
the latter is proportional to $1/\sigma_{\mathbbm{n}}(x)$, which is bounded on
$[\varepsilon,1-\varepsilon]$.

###  3.5 Random transformations

Our model was defined as a Markov chain. We now present a slightly different,
yet equivalent, point of view. Namely, we will pick up a family of maps
$T_{\eta}\colon I\to I$, $\eta\in[0,1]$, in such a way that

$\operatorname{Leb}\\{\eta:T_{\eta}(x)\in A\\}=P_{x}(A)$ (23)

for all $A\subset I$. We can then define a stochastic process

$\bar{x}_{t+1}=T_{\eta_{t+1}}(\bar{x}_{t}),$ (24)

where $\eta_{t}$ are independent and uniformly distributed in $[0,1]$. We can
write $\bar{x}_{t}=T_{\eta_{t}}\circ\cdots\circ T_{\eta_{1}}\bar{x}_{0}$,
where $(\eta_{t})_{t\in\mathbb{N}}$ is an i.i.d. stochastic process, i.e. the
process (24) follows the orbits under the concatenation of randomly chosen
maps from the family. One can show that the two processes are equivalent, see,
for instance, [37]. Conversely, starting with a family of maps $T_{\eta}$ one
can use (23) to define transition probabilities $P_{x}$ and thus a Markov
chain.

Rewriting (23) as $P_{x}(A)=\int_{0}^{1}\mathbbm{1}_{A}(T_{\eta}(x))d\eta$ and
plugging into (12) we get the _disintegration formula_ for the Markov
operator:

$\mathcal{L}=\int L_{\eta}d\eta,$ (25)

where $L_{\eta}$ are the transfer operators associated with $T_{\eta}$. In
particular, a measure $\mu$ is stationary for the Markov chain if and only if
it satisfies $\mu=\int L_{\eta}\mu d\eta$, i.e., for all $A\subset I$,

$\mu(A)=\int\mu(T_{\eta}^{-1}A)d\eta.$ (26)

Equation (26) is usually taken as the definition of a stationary measure for
the family of random maps. Every such measure corresponds to a product measure
that is invariant for the skew-product with the Bernoulli shift in the base
and the maps $T_{\eta}$ in the fibers; we refer to [4, Section 2] for details.

Figure 4: Random maps (27) for $\phi^{*}=0.5$, $\omega=0.2$,
$\mathbbm{n}=10,100,1000$, $\eta=k/16$.

As follows from (21), the random maps for the main example have the form

$T_{\eta}(x)=T(x)+q_{\mathbbm{n}}(\eta)\sigma_{\mathbbm{n}}(x),$ (27)

where $q_{\mathbbm{n}}$ is the quantile function of the truncated normal
distribution $\mathcal{N}_{b_{\mathbbm{n}}}(0,1)$. Indeed, since
$q_{\mathbbm{n}}$ maps the uniform measure on $[0,1]$ to the truncated
standard Gaussian measure on $[-b_{\mathbbm{n}},b_{\mathbbm{n}}]$, we have

$\displaystyle P_{x}(A)=$
$\displaystyle\mathbb{P}\\{T(x)+\sigma_{\mathbbm{n}}(x)Z_{t}\in A\\}$
$\displaystyle=$
$\displaystyle\operatorname{Leb}\\{\eta:T(x)+\sigma_{\mathbbm{n}}(x)q_{\mathbbm{n}}(\eta)\in
A\\}$ $\displaystyle=$ $\displaystyle\operatorname{Leb}\\{\eta:T_{\eta}(x)\in
A\\}.$

We can equivalently rewrite (27) as

$T_{\eta}(x)=T(x)+\tilde{q}_{\mathbbm{n}}(\eta)\sigma(x),$ (28)

where $\tilde{q}_{\mathbbm{n}}$ is the quantile function of
$\mathcal{N}_{b_{1}}(0,\frac{1}{\mathbbm{n}})$. Notice that for different
$\mathbbm{n}$, the set $\\{T_{\eta}\\}_{\eta\in[0,1]}$ consists of the same
maps, however the ones close to $T=T_{\frac{1}{2}}$ are given bigger weights
for large $\mathbbm{n}$. More precisely, for every $\delta>0$, we have

$\sup_{\eta\in[\delta,1-\delta]}\sup_{x\in
I}|T_{\eta}(x)-T(x)|\leq\sigma_{\max}\sup_{\eta\in[\delta,1-\delta]}|\tilde{q}_{\mathbbm{n}}(\eta)|\xrightarrow[\mathbbm{n}\to\infty]{}0.$
(29)

##  4 Stationary measure

We are now going to establish the existence of a unique stationary measure for
the chain. This will be accomplished in the following steps: we first prove
the Lasota-Yorke inequality (30); as a consequence we will get a finite number
of ergodic absolutely continuous stationary measures whose supports are
mutually disjoint up to sets of zero Lebesgue measure. The uniqueness will be
proved by showing that all the previous components share a measurable set of
positive Lebesgue measure. Throughout this section $\mathbbm{n}$ is fixed and
we will omit it from notations.

###  4.1 Existence

We will first show that there are _finitely many_ ergodic stationary densities
of bounded variation. The following lemma will be useful in the sequel.

######  Lemma 4.1.

For any $\rho\in BV$, if
$C:=\operatorname*{ess\,sup}_{x\in\operatorname{supp}\rho}|p(x,\cdot)|_{TV}<\infty$,
then

$|\mathcal{L}\rho|_{TV}\leq
C\|\rho\|_{1}\quad\text{and}\quad\|\mathcal{L}\rho\|_{BV}\leq(C+1)\|\rho\|_{1}.$

###### Proof.

For the first inequality we have

$\displaystyle|\mathcal{L}\rho|_{TV}$ $\displaystyle=\sup\sum_{i}\left|\int
p(x,y_{i+1})\rho(x)dx-\int p(x,y_{i})\rho(x)dx\right|$
$\displaystyle\leq\sup\sum_{i}\int\left|p(x,y_{i+1})\rho(x)-p(x,y_{i})\rho(x)\right|dx$
$\displaystyle\leq\int|p(x,\cdot)|_{TV}\rho(x)dx\leq C\|\rho\|_{1}.$

The second inequality follows from the first one, since the Markov operator is
an isometry, i.e. $\|\mathcal{L}\rho\|_{1}=\|\rho\|_{1}$ for all $\rho\in
L^{1}$. ∎

We say that a stochastic kernel $p(x,y)$ has _uniformly bounded variations_ if
$|p(x,\cdot)|_{TV}\in L^{\infty}$, i.e. there is $C>0$ such that
$|p(x,\cdot)|_{TV}\leq C$ for almost every $x\in I$.

######  Proposition 4.2.

If the kernel $p$ has uniformly bounded variations, then the operator
$\mathcal{L}$ is quasi-compact and there exist finitely many ergodic
stationary measures with densities in $BV$ and, moreover, their supports are
mutually disjoint up to sets of zero Lebesgue measure.

###### Proof.

By Lemma 4.1, for every $n$,

$\|\mathcal{L}^{n}\rho\|_{BV}=|\mathcal{L}^{n}\rho|_{TV}+\|\mathcal{L}^{n}\rho\|_{1}\leq
C\|\mathcal{L}^{n-1}\rho\|_{1}+\|\rho\|_{1}=(C+1)\|\rho\|_{1}.$

In particular,

$\|\mathcal{L}\rho\|_{BV}\leq(C+1)\|\rho\|_{1}\leq\eta\|\rho\|_{BV}+(C+1)\|\rho\|_{1}$
(30)

for any $\eta<1$. This is the Lasota-Yorke inequality. The latter, plus the
fact that $BV$ is compactly embedded in $L^{1}$, implies that the peripheral
spectrum of $\mathcal{L}$ is discrete and therefore the chain will admit
_finitely many_ (at least one) absolutely continuous ergodic stationary
measures, with supports that are mutually disjoint up to sets of zero Lebesgue
measure. Moreover, the essential spectral radius is strictly smaller than the
spectral radius (spectral gap). These properties, which are consequences of
the Ionescu-Tulcea-Marinescu theorem, are summarized by saying that the
operator $\mathcal{L}$ acting on $BV$ is quasi-compact, see, e.g., [11, 19,
18] for an exhaustive presentation of these results and [10, Section 2.3] for
a specific application to random systems. ∎

######  Remark 4.3.

Let us mention that whenever the operator $\mathcal{L}$ is quasi-compact and
the largest eigenvalue, which is $1$ in our case, is simple and therefore
there is only one stationary measure with density in $BV$, then the norm of
$\|\mathcal{L}^{k}f\|_{BV}$ goes exponentially fast to zero when $k\to\infty$,
for $f\in BV$ and $\int fdx=0$ (exponential decay of correlations). This fact
will be extensively used in Section 6.

Since the variance (19) vanishes at 0 and 1, the kernel (22) is in fact
unbounded. However, we can still apply Proposition 4.2 under a suitable
restriction of the domain of $\mathcal{L}.$ We first state a general result
which allows us to confine the stationary measures. For any $\varepsilon>0$,
define the interval

$I_{\varepsilon,\Gamma}:=[\varepsilon,1-\Gamma/2].$

######  Lemma 4.4.

Under the assumptions of Section 3.3, there is $\varepsilon>0$ such that any
stationary measure $\mu$ has $\operatorname{supp}\mu\subset\\{0\\}\cup
I_{\varepsilon,\Gamma}$. If $\mu$ is continuous,
$\operatorname{supp}\mu\subset I_{\varepsilon,\Gamma}$.

###### Proof.

First, notice that any stationary measure is supported on the interval
$K_{\Gamma}:=[-\Gamma/2,1-\Gamma/2]$. Indeed, by invariance,
$\mu(K_{\Gamma}^{c})=\int P_{x}(K_{\Gamma}^{c})d\mu(x)=0$, because
$P_{x}(K_{\Gamma}^{c})=0$ for all $x$.

Fix $\varepsilon>0$ such that $T(x)-s(x)>x$ for $x\in(0,\varepsilon)$. By
choosing a smaller $\varepsilon$ if needed, we may also assume that
$T(x)-s(x)>\varepsilon$ for $x\in[\varepsilon,1-\Gamma/2]$. Then
$\inf\operatorname{supp}P_{x}>\min\\{x,\varepsilon\\}$ for every $x\in
K_{\Gamma}\setminus\\{0\\}$. This means that for any realization $(x_{t})$ of
the process, either all $x_{t}=0$ (clearly, $0$ is a fixed point, since
$P_{0}=\delta_{0}$) or _all but finitely many_ $x_{t}>\varepsilon$.

On the other hand, if $\mu([-\Gamma,\varepsilon)\setminus\\{0\\})>0$, then by
the Poincaré recurrence theorem, applied to the shift on
$(\Omega,\mathbb{P}_{\mu})$, $\mathbb{P}_{\mu}$-almost surely there would
exist a realization $(x_{t})$ with _infinitely many_ $0\neq
x_{t}<\varepsilon$, which is not possible, as we showed above. This finishes
the proof. ∎

######  Theorem 4.5.

The chain defined in Section 3.3 admits finitely many ergodic stationary
measures with densities in $BV$. Moreover, there is $\varepsilon>0$ such that
$\operatorname{supp}\mu\subset I_{\varepsilon,\Gamma}$ for any such measure
$\mu$.

###### Proof.

By Lemma 4.4, the density of any absolutely continuous stationary measure
belongs to the subspace $Y:=\\{h\in L^{1}\mid\operatorname{supp}h\subset
I_{\varepsilon,\Gamma}\\}$. From the first part of the proof of Lemma 4.4 it
also follows that $Y$ is $\mathcal{L}$-invariant. Moreover, the kernel (16)
has uniformly bounded variations when restricted to
$I_{\varepsilon,\Gamma}\times I_{\varepsilon,\Gamma}$. Indeed,
$|p_{\mathbbm{n}}(x,\cdot)|_{TV}=|g_{x,\mathbbm{n}}|_{TV}$ and the latter is
bounded on $I_{\varepsilon,\Gamma}$ by (14). We can therefore apply
Proposition 4.2. ∎

######  Remark 4.6.

It is worth noticing that the preceding result is completely independent of
the structure of the unimodal map $T$. In this respect we could consider maps
admitting attracting periodic points or Cantor sets of measure zero, but still
producing smooth stationary measures when perturbed with our additive noise.

###  4.2 Uniqueness

We begin with the following simple lemma that links the topological dynamics
of $T$ with the structure of any stationary measure.

######  Lemma 4.7.

For any stationary measure $\mu$ and any open set $U$, if $\mu(U)=0$, then
also $\mu(T^{-k}U)=0$ for all $k>0$.

###### Proof.

It is enough to show that $\mu(T^{-1}U)=0$, the result then follows by
induction. By invariance we have $0=\mu(U)=\int
P_{x}(U)d\mu(x)\geq\int_{T^{-1}U}P_{x}(U)d\mu(x)$. But for every $x\in
T^{-1}U$, $T(x)\in U\cap\operatorname{supp}P_{x}$, and hence $P_{x}(U)>0$.
Therefore the latter integral can only be zero if $\mu(T^{-1}U)=0$. ∎

With the help of the following lemma we will show that the support of any
stationary measure contains the support of the $T$-invariant measure (atomic
in the periodic case). Recall that $x\in\operatorname{supp}\mu$ iff $\mu(U)>0$
for any open $U\ni x$.

######  Lemma 4.8.

Let $A\subset I$ be such that (1) $T$ is topologically transitive on $A$ and
(2) $\bigcup_{k=0}^{\infty}T^{-k}U=I$ for any open set $U\supset A$. Then
$A\subset\operatorname{supp}\mu$ for any stationary measure $\mu$.

###### Proof.

Given an open set $U$ with $U\cap A\neq\emptyset$, by transitivity
$A\subset\bigcup_{k=0}^{\infty}T^{-k}U$, and therefore
$\bigcup_{k=0}^{\infty}T^{-k}U=I$. By Lemma 4.7, $\mu(U)>0$. ∎

Now we are ready to state the main result of this section.

######  Theorem 4.9.

If $T$ is either periodic or chaotic, then the chain defined in Section 3.3
admits a unique stationary measure $\mu$ with $BV$ density. Moreover,
$\operatorname{supp}\mu$ contains a neighbourhood of the periodic cycle if $T$
is periodic, or the interval $[T(\Delta),\Delta]$ if $T$ is chaotic.

###### Proof.

If $T$ is periodic, a globally attracting cycle $\mathcal{O}$ satisfies the
assumptions of Lemma 4.8, therefore $\mathcal{O}\subset\operatorname{supp}\mu$
for any stationary measure $\mu$. Let us show that $\operatorname{supp}\mu$
contains an open neighbourhood of $\mathcal{O}$. Recall that $p(x,y)>0$ if and
only if $y\in(T(x)-s(x),T(x)+s(x))$, and $T$ and $s$ are continuous. Given
$x_{0}\in\mathcal{O}$, let $x_{1}\in\mathcal{O}$ be such that
$x_{0}=T(x_{1})$. Since $(x_{1},x_{0})\in\\{(x,y)\mid p(x,y)>0\\}$ and the
latter set is open, we can find open sets $U\ni x_{1}$ and $V\ni x_{0}$ such
that $p(x,y)>0$ for all $x\in U$, $y\in V$. Also $\mu(U)>0$, because
$x_{1}\in\operatorname{supp}\mu$. Denoting $h$ the density of $\mu$, by
invariance we get $h(y)\geq\int_{U}p(x,y)d\mu(x)>0$ for all $y\in V$, i.e.
$V\subset\operatorname{supp}\mu$.

If $T$ is chaotic, then from (At) we know that the set
$I_{\Delta}=[T(\Delta),\Delta]$ is invariant for $T$ and that $T$ is
topologically transitive when restricted to $I_{\Delta}$. It also follows from
(C3) that $\bigcup_{k=0}^{\infty}T^{-k}I_{\Delta}=I$, so again we can apply
Lemma 4.8.

In both cases, we conclude that
$\operatorname{Leb}(\operatorname{supp}\mu_{1}\cap\operatorname{supp}\mu_{2})>0$
for any stationary measures $\mu_{1}$, $\mu_{2}$, and therefore by Proposition
4.2 they must coincide. Hence the stationary measure is unique. ∎

##  5 Stochastic stability

Once we consider random perturbations of a deterministic dynamics, an
important question is to investigate the stochastic stability of the system,
which means to determine if a sequence of stationary measure will converge, in
a sense to precise, to the invariant measure of the unperturbed map. In our
case the sequence of probability measures is given by
$\mu_{\mathbbm{n}}:=h_{\mathbbm{n}}dx$. These measures belong to the set of
Borel probability measures on the unit interval, which is a compact metric
space with the weak-* topology111111The _weak-* topology_ is given by the
family of seminorms $\|\rho\|_{\varphi}=\int\varphi d\rho$,
$\rho\in\mathcal{M}$, $\varphi\in C^{0}$.. There will be therefore at least
one subsequence $(\mu_{\mathbbm{n}_{k}})_{k\geq 1}$ converging to a
probability measure $\mu_{\infty}$ on $I$. Our objective is to prove that: (i)
$\mu_{\infty}$ is invariant, (ii) it is the same for any convergent
subsequence, if more than one, and (iii) it coincides with $\mu$. Whenever
that happens we will say that our random system is weakly stochastic stable.
This result could be strengthened by showing that
$\|h_{\mathbbm{n}}-h\|_{1}\to 0$, which is called the strong stochastic
stability; we are not able at the moment to get this result. Instead we now
give a sufficient condition to get the weak stochastic stability:

  1. (Aq)

There exist $q>1$ and $C_{q}>0$ such that for all $\mathbbm{n}\geq 1$ we have
$\|h_{\mathbbm{n}}\|_{q}\leq C_{q}$.

We will see in the next section that with the preceding assumption we can
prove the convergence of the Lyapunov exponent (Theorem 6.3) and then verify
it numerically, which is an indirect indication of the validity of (Aq).

######  Lemma 5.1.

For every $x\in I$, $P_{x}^{(\mathbbm{n})}$ converges to $\delta_{Tx}$ in the
weak-* topology as $\mathbbm{n}\to\infty$, i.e. $\int\varphi
dP_{x}^{(\mathbbm{n})}\to\varphi(Tx)$ for all $\varphi\in C^{0}(I)$.

###### Proof.

For arbitrary $\varepsilon>0$ we can split

$\int\varphi dP_{x}^{(\mathbbm{n})}=\int_{B_{\varepsilon}(Tx)}\varphi
dP_{x}^{(\mathbbm{n})}+\int_{B_{\varepsilon}(Tx)^{c}}\varphi
dP_{x}^{(\mathbbm{n})}.$

By Chebyshev’s inequality,
$P_{x}^{(\mathbbm{n})}(B_{\varepsilon}(Tx)^{c})\leq{\mathrm{Var}^{2}P_{x}^{(\mathbbm{n})}}/{\varepsilon^{2}}\to
0$ as $\mathbbm{n}\to\infty$, while $\varphi$ is bounded, so the second
integral can be made arbitrarily small for $\mathbbm{n}$ large. Consequently,
$P_{x}^{(\mathbbm{n})}(B_{\varepsilon}(Tx))\to 1$, and since $\varphi$ is
continuous, the first integral can be made arbitrarily close to $\varphi(Tx)$.
∎

The proof of (i) below follows a suggestion in [5, Theorem D].

######  Proposition 5.2.

Let $\mu$ be a weak-* limit measure of a sequence
$\mu_{\mathbbm{n}_{k}}=h_{\mathbbm{n}_{k}}dx$. If $h_{\mathbbm{n}_{k}}$
satisfy (Aq), then

  1. (i)

$\mu$ is absolutely continuous with density in $L^{q}$;

  2. (ii)

$\mu$ is invariant under $T$.

###### Proof.

(i) Let $\varphi\in C^{0}(I)$. By Hölder’s inequality, with $p=\frac{q}{q-1}$,

$\Bigl{|}\int\varphi d\mu\Bigr{|}=\Bigl{|}\lim_{k\to\infty}\int\varphi
h_{\mathbbm{n}_{k}}dx\Bigr{|}\leq\lim_{k\to\infty}\|h_{\mathbbm{n}_{k}}\|_{q}\|\varphi\|_{p}\leq
C_{q}\|\varphi\|_{p}.$

Therefore the map $L^{p}\ni\varphi\mapsto\int\varphi d\mu\in\mathbb{R}$ is
continuous, since $C^{0}$ is dense in $L^{p}$, and therefore such a functional
will be in $L^{q}$, namely $\mu=hdx$, $h\in L^{q}$, $\|h\|_{q}\leq C_{q}$.

(ii) It suffices to prove that any test function $\varphi\in C^{0}(I)$
satisfies

$\int\varphi h_{\mathbbm{n}_{k}}dy-\int\varphi\circ Th_{\mathbbm{n}_{k}}dy\to
0.$

Since $h_{\mathbbm{n}_{k}}=\mathcal{L}h_{\mathbbm{n}_{k}}=\int
p_{\mathbbm{n}_{k}}(x,\cdot)h_{\mathbbm{n}_{k}}(x)dx$, by changing the order
of integration in the first integral and subtracting the second, we get

$\int h_{\mathbbm{n}_{k}}(x)\left[\int
p_{\mathbbm{n}_{k}}(x,y)\varphi(y)dy-\varphi(T(x))\right]dx.$ (31)

Since the function $\psi_{\mathbbm{n}_{k}}(x):=\int
p_{\mathbbm{n}_{k}}(x,y)\varphi(y)dy-\varphi(T(x))$ is uniformly bounded,

$\eqref{gh}\leq\|h_{\mathbbm{n}_{k}}\|_{q}\|\psi_{\mathbbm{n}_{k}}\|_{p}\leq
C_{q}\|\psi_{\mathbbm{n}_{k}}\|_{p}.$

By Lemma 5.1, $\psi_{\mathbbm{n}_{k}}(x)\to 0$ for every $x\in I$, and
therefore by dominated convergence $\|\psi_{\mathbbm{n}_{k}}\|_{p}\to 0$. ∎

######  Theorem 5.3.

Under Assumption (Aq), the chain defined in Section 3.3 is weakly stochastic
stable, i.e. the stationary probabilities converge to the unique $T$-invariant
probability in the weak-* topology as $\mathbbm{n}\to\infty$.

###### Proof.

Since $T$ admits a unique invariant measure, $\mu$ must be the same for all
convergent subsequences in Proposition 5.2, and therefore the entire sequence
$\mu_{\mathbbm{n}}$ converges to $\mu$. ∎

It follows from Proposition 5.2 that Assumption (Aq) cannot be satisfied in
the periodic case, since the limiting $T$-invariant measure is singular and
supported on the periodic orbit, so Theorem 5.3 only covers the chaotic case.
We will now give a proof in the periodic case under the following assumption:

  1. (As)

For all $\mathbbm{n}$ sufficiently large and all
$x\in\operatorname{supp}\mu_{\mathbbm{n}}$ we have
$|T^{\prime}(x)|\leq\tau<1$.

######  Proposition 5.4.

If $T$ is periodic and satisfies (As), then the chain defined in Section 3.4
is weakly stochastic stable.

###### Proof.

Let us first consider the case when $T$ has a globally attracting fixed point
$x_{0}$. We need to show that for any test function $\varphi\in C^{0}(I)$ we
have $\int\varphi(x)h_{\mathbbm{n}}(x)dx\to\varphi(x_{0})$ as
$\mathbbm{n}\to\infty$. Since $h_{\mathbbm{n}}$ is a fixed point of the random
transfer operator (25) and this operator is the dual of the random Koopman
operator $\varphi\mapsto\int\varphi\circ T_{\eta}d\eta$ (see, e.g., [4,
Section 2] for details), the previous weak limit leads to prove that the
following quantity

$\int_{I}\int_{[0,1]^{k}}\bigl{(}\varphi(T_{\eta_{k}}\circ\cdots\circ
T_{\eta_{1}}(x))-\varphi(x_{0})\bigr{)}h_{\mathbbm{n}}(x)\,d\bar{\eta}\,dx$
(32)

goes to 0 as $\mathbbm{n}\to\infty$, where $k$ is an arbitrary fixed number
and $\bar{\eta}=(\eta_{1},\ldots,\eta_{k})$.

Given $\varepsilon>0$, let $\zeta>0$ be such that
$|\varphi(x)-\varphi(x_{0})|<\varepsilon$ when
$|x-x_{0}|<\frac{2\zeta}{1-\tau}$. Fix $k$ such that, for all $\mathbbm{n}$,

$\sup_{x\in\operatorname{supp}\mu_{\mathbbm{n}}}|T^{k}(x)-x_{0}|<\zeta.$ (33)

Next, fix $\delta>0$ such that

$2\|\varphi\|_{\infty}(1-(1-2\delta)^{k})<\varepsilon.$ (34)

Finally, by (29), for all $\mathbbm{n}$ sufficiently large, we have

$\sup_{\eta\in[\delta,1-\delta]}\sup_{x\in I}|T_{\eta}(x)-T(x)|<\zeta.$ (35)

We now split the integral (32) in the $\bar{\eta}$ variable over the region
$E:=[\delta,1-\delta]^{k}$ and its complement. On $E^{c}$ the absolute value
of (32) is bounded by (34). Notice that the integral over $x$ takes place on
the support of $\mu_{\mathbbm{n}}$, where (As) holds. Moreover, since the map
$[0,1]\ni\eta\mapsto T_{\eta}\in C^{2}(I)$ is continuous, each $T_{\eta}$ maps
$\operatorname{supp}\mu_{\mathbbm{n}}$ to itself; see [6]. Therefore, for
$\bar{\eta}\in E$ and $x\in\operatorname{supp}\mu_{\mathbbm{n}}$, by (35) and
(As), we have

$|T_{\eta_{2}}\circ
T_{\eta_{1}}(x)-T^{2}(x)|\leq|T_{\eta_{2}}(T_{\eta_{1}}(x))-T(T_{\eta_{1}}(x))|+|T(T_{\eta_{1}}(x))-T^{2}(x)|<\zeta+\tau\zeta.$

By induction we easily get $|T_{\eta_{k}}\circ\cdots\circ
T_{\eta_{1}}(x)-T^{k}(x)|<\frac{\zeta}{1-\tau}$ and therefore, in view of
(33), $|T_{\eta_{k}}\circ\cdots\circ
T_{\eta_{1}}(x)-x_{0}|<\frac{2\zeta}{1-\tau}$ for all $\bar{\eta}\in E$ and
$x\in\operatorname{supp}\mu_{\mathbbm{n}}$. Then, by the choice of $\zeta$,
$|\varphi(T_{\eta_{k}}\circ\cdots\circ
T_{\eta_{1}}(x))-\varphi(x_{0})|<\varepsilon$ and the absolute value of (32)
over $E$ is therefore bounded by $\varepsilon$.

It is straightforward to modify the above proof for the case when $T$ has a
globally attracting periodic orbit of length $m>1$. One needs to replace $T$
with $T^{m}$, the latter will have $m$ attracting fixed points. The
corresponding random maps of the form
$T_{\bar{\eta}}=T_{\eta_{1}}\circ\cdots\circ T_{\eta_{m}}$ will be
parametrized by $\bar{\eta}=(\eta_{1},\cdots,\eta_{m})\in[0,1]^{m}$ endowed
with the Lebesgue measure. We leave the details to the reader. ∎

######  Remark 5.5.

It follows from the proof that Proposition 5.4 remains valid for the general
class of Markov chains defined in Section 3.3 whenever (29) holds, which is in
turn the case when $T$ and $g_{x,\mathbbm{n}}$ are sufficiently smooth.

######  Remark 5.6.

We conjecture that, if $T$ is periodic with the attracting periodic orbit
$\mathcal{O}$, then $h_{\mathbbm{n}}\to 0$ uniformly on compact sets $K\subset
I\setminus\mathcal{O}$ as $\mathbbm{n}\to\infty$. This property, that we
checked numerically, straightens the previous result. In particular,

  1. (Ac)

If $T$ is periodic and the critical point $c$ does not belong to the
attracting periodic orbit, then $h_{\mathbbm{n}}\to 0$ uniformly in a
neighbourhood of $c$ as $\mathbbm{n}\to\infty$.

##  6 Lyapunov exponent

###  6.1 Average Lyapunov exponent

We are interested in the existence of the _Lyapunov exponent_ for the slow
component, which in our case is defined $\mathbb{P}_{\mu}$-almost surely as
the limit

$\Lambda=\lim_{n\to\infty}\frac{1}{n}\sum_{t=0}^{n-1}\log|T^{\prime}(X_{t})|$
(36)

along the chain $(X_{t})_{t\geq 0}$. We now motivate such a choice. It is
twofold: first of all we want to reproduce the Lyapunov exponent of the
unperturbed map $T$ in the limit of zero noise, which we will get in Theorem
6.3; successively we want an indicator which kept memory of the underlying
slow dynamics played by the map $T$. We will, in particular, show that such an
exponent is negative for periodic $T$, even in presence of mixing stationary
measure.

We now return to (36); if the chain admits a unique stationary probability
$\mu$, then, by the ergodic theorem for Markov chains, the above limit equals

$\int\log|T^{\prime}|d\mu,$ (37)

assuming $\log|T^{\prime}|\in L^{1}(\mu)$.

######  Remark 6.1.

The Lyapunov exponent (37) was called the average Lyapunov exponent in [31,
48], and it was associated to the phenomenon of noise induced order, which
happens when the perturbed systems admit a unique stationary measure depending
on some parameter, say $\theta$, and the Lyapunov exponent depends
continuously on $\theta$ and exhibits a transition from positive to negative
values, see also [43] for an experimental evidence of this fact. We will
partially prove this phenomenon below by combining Corollary 6.4 and Theorem
6.5, and show it numerically in Section 7.2.

A unimodal map $T$ is said to have a _critical point of order $l$_ if there is
a constant $D$ such that $D^{-1}|x-c|^{l-1}\leq|T^{\prime}(x)|\leq
D|x-c|^{l-1}$. In this case it was proved in [49] that the invariant density
for $T$ is in $L^{q}$, with $q<\frac{l}{l-1}$. We assume in (A3) that $T$ has
a critical point of order 2. It is easy to check that (18) satisfies this
assumption. Consequently, $\log|T^{\prime}|$ is in $L^{p}$ for any $p\geq 1$.

######  Theorem 6.2.

If $T$ is periodic or chaotic, the limit (36) exists almost surely.

###### Proof.

The integral (37) is finite, because $\log|T^{\prime}|$ is in $L^{1}$ and the
unique stationary measure $\mu$ has bounded density, as we proved in Section
4. ∎

Once we know that the Lyapunov exponent exists almost surely, it is natural to
ask how it depends on the model parameters, for instance, the length
$\mathbbm{n}$ of the fast component series. We have the following

######  Theorem 6.3.

Suppose one of the following is satisfied: (a) $T$ verifies (Aq); (b) $T$ is
periodic and verifies (As) and (Ac). Then the Lyapunov exponent (36) converges
to the Lyapunov exponent of the deterministic map $T$ as
$\mathbbm{n}\to\infty$.

###### Proof.

(a) Denote with $\mu_{\mathbbm{n}}=h_{\mathbbm{n}}dx$ the unique stationary
measure associated to $\mathbbm{n}$ and with $\mu=hdx$ the unique invariant
measure for $T$. We need to show that

$\int\log|T^{\prime}|d\mu_{\mathbbm{n}}\xrightarrow[\mathbbm{n}\to\infty]{}\int\log|T^{\prime}|d\mu.$
(38)

Let $q>1$ be such that $h_{\mathbbm{n}},h\in L^{q}$ and set
$p:=\frac{q}{q-1}$. Since $\log|T^{\prime}|\in L^{p}$, for any $\varepsilon>0$
there is $\varphi_{\varepsilon}\in C^{0}$ such that
$\|\log|T^{\prime}|-\varphi_{\varepsilon}\|_{p}<\varepsilon$. Write

$\int\log|T^{\prime}|d\mu_{\mathbbm{n}}=\int(\log|T^{\prime}|-\varphi_{\varepsilon})h_{\mathbbm{n}}dx+\int\varphi_{\varepsilon}d\mu_{\mathbbm{n}}$

and

$\int\log|T^{\prime}|d\mu=\int(\log|T^{\prime}|-\varphi_{\varepsilon})hdx+\int\varphi_{\varepsilon}d\mu.$

Since $\log|T^{\prime}|-\varphi_{\varepsilon}\in L^{p}$ and
$h_{\mathbbm{n}},h\in L^{q}$, we have

$\int|\log|T^{\prime}|-\varphi_{\varepsilon}|h_{\mathbbm{n}}dx\leq\|\log|T^{\prime}|-\varphi_{\varepsilon}\|_{p}\|h_{\mathbbm{n}}\|_{q}\leq\varepsilon
C_{q},$

and the same inequality holds for the integral with respect to $\mu$. Finally,
from Theorem 5.3 we know that $\mu_{\mathbbm{n}}\xrightarrow[]{w^{*}}\mu$,
hence
$\int\varphi_{\varepsilon}d\mu_{\mathbbm{n}}\to\int\varphi_{\varepsilon}d\mu$
as $\mathbbm{n}\to\infty$.

(b) We know from Proposition 5.4 that $\int\varphi
d\mu_{\mathbbm{n}}\to\frac{1}{|\mathcal{O}|}\sum_{x\in\mathcal{O}}\varphi(x)$
for all $\varphi\in C^{0}(I)$ as $\mathbbm{n}\to\infty$, since the
$T$-invariant measure $\mu$ is atomic and supported on the attracting periodic
orbit $\mathcal{O}$. Let us first consider the case when the critical point
$c$ belongs to $\mathcal{O}$; the right-hand side of (38) is then $-\infty$.
Denoting $f_{m}(x):=\max\\{\log|T^{\prime}(x)|,-m\\}\in C^{0}(I)$, for every
$m$ we have

$\int\log|T^{\prime}|d\mu_{\mathbbm{n}}\leq\int
f_{m}d\mu_{\mathbbm{n}}\xrightarrow[\mathbbm{n}\to\infty]{}\frac{1}{|\mathcal{O}|}\sum_{x\in\mathcal{O}}f_{m}(x)\leq-\frac{m}{|\mathcal{O}|}+C,$

because $f_{m}\leq C:=\sup\log|T^{\prime}|$ and $f_{m}(c)=-m$. Therefore
$\int\log|T^{\prime}|d\mu_{\mathbbm{n}}\to-\infty$ as $\mathbbm{n}\to\infty$.

If $c\notin\mathcal{O}$, we can fix a neighbourhood $U\ni c$ given by (Ac) and
split

$\int\log|T^{\prime}|d\mu_{\mathbbm{n}}=\int_{U}\log|T^{\prime}|d\mu_{\mathbbm{n}}+\int_{U^{c}}\log|T^{\prime}|d\mu_{\mathbbm{n}}.$

The first term is bounded by $\|\log|T^{\prime}|\|_{1}\sup_{U}h_{\mathbbm{n}}$
and vanishes as $\mathbbm{n}\to\infty$ by (Ac), while the second one converges
to $\int\log|T^{\prime}|d\mu$ by Proposition 5.4 (approximate
$\log|T^{\prime}|\mathbbm{1}_{U^{c}}$ with a suitable continuous function). ∎

######  Corollary 6.4.

Under the assumptions of Theorem 6.3 and for $\mathbbm{n}$ large enough,
$\Lambda$ is positive if $T$ is chaotic, and negative if $T$ is periodic.

In a few cases the negativity of the Lyapunov exponent can be shown relatively
easily. For instance, denote $\\{x\in I\mid|T^{\prime}(x)|\leq 1\\}=[m,M]$ and
$\bar{\Delta}:=\sup_{x\in I}T(x)+s(x)$. If $T(x)-s(x)>\min\\{x,m\\}$ for all
$x\in(0,\bar{\Delta}]$, then, arguing as in the proof of Lemma 4.4, one can
show that any continuous stationary measure $\mu$ has
$\operatorname{supp}\mu\subset[m,\bar{\Delta}]$. So if, moreover,
$\bar{\Delta}\leq M$, then $\Lambda<0$. Following the classification given in
Section 3.1, let us consider the case $T(c)<c$, where the map $T$ exhibits a
globally attracting fixed point. In this case, the conditions above will be
satisfied if $s(x)$ is small enough, in other words, the stationary measure
will be supported in a neighbourhood of the fixed point, where
$|T^{\prime}(x)|\leq 1$. For other cases, we provide some numerical examples
in Section 7.2.

###  6.2 Continuity of the Lyapunov exponent

Denote by
$\Theta:=\\{\theta=(\phi^{*},\omega,\mathbbm{n})\in(0,1)^{2}\times(0,\infty)\mid\max
T_{\theta}<1\\}$ the (extended) parameter space. In order to prove the
continuity of the Lyapunov exponent, we will assume that $T_{\theta}(x)\in
C^{3}(\Theta\times[0,1])$ and $p_{\theta}(x,y)\in
C^{2}(\Theta\times(0,1)^{2})$. It is straightforward that our main example
defined in Section 3.4 satisfies this assumption. Let
$\tilde{\Theta}\subset\Theta$ be the set of parameters $\theta$ for which
there is a unique stationary measure $\mu_{\theta}$ with a density
$h_{\theta}\in BV$; we proved in Section 4 that this is the case if
$T_{\theta}$ is periodic or chaotic, but our numerical investigations confirm
that in fact $\operatorname{Leb}(\Theta\setminus\tilde{\Theta})=0$.

######  Theorem 6.5.

The mapping $\tilde{\Theta}\ni\theta\mapsto\Lambda_{\theta}\in\mathbb{R}$ is
continuous.

###### Proof.

Fix an exhaustion of $\Theta$ by nested compact sets $\Theta_{\iota}$ and set
$\tilde{\Theta}_{\iota}:=\Theta_{\iota}\cap\tilde{\Theta}$. It is enough to
prove that the mapping
$\tilde{\Theta}_{\iota}\ni\theta\mapsto\Lambda_{\theta}\in\mathbb{R}$ is
continuous on each $\tilde{\Theta}_{\iota}$, and from now on we fix one of
them. As we showed in Lemma 4.4, for each $\theta\in\Theta$,
$\operatorname{supp}\mu_{\theta}\subset
I_{\varepsilon_{\theta}}=[\varepsilon_{\theta},1-\varepsilon_{\theta}]$, and
since $\varepsilon_{\theta}$ can be shown to depend continuously on $\theta$,
we can find a single $\varepsilon>0$ that works for all
$\theta\in\Theta_{\iota}$.

Given $\theta,\theta^{\prime}\in\tilde{\Theta}_{\iota}$ we can write

$|\Lambda_{\theta}-\Lambda_{\theta^{\prime}}|=\Bigl{|}\int\log|T_{\theta}^{\prime}|h_{\theta}dx-\int\log|T_{\theta^{\prime}}^{\prime}|h_{\theta^{\prime}}dx\Bigr{|}\leq$
$\int\bigl{|}\log|T_{\theta}^{\prime}|\bigr{|}|h_{\theta}-h_{\theta^{\prime}}|dx+\int\bigl{|}\log|T^{\prime}_{\theta}|-\log|T_{\theta^{\prime}}^{\prime}|\bigr{|}h_{\theta^{\prime}}dx.$
(39)

To bound the second term in (39), first notice that, by Lemma 4.1,

$\|h_{\theta^{\prime}}\|_{\infty}\leq\|h_{\theta^{\prime}}\|_{BV}=\|\mathcal{L}_{\theta^{\prime}}h_{\theta^{\prime}}\|_{BV}\leq
C\|h_{\theta^{\prime}}\|_{1}=C,$

where

$C=1+\sup_{\theta\in\Theta_{\iota}}\sup_{x\in
I_{\varepsilon}}|p_{\theta}(x,\cdot)|_{TV}\leq
1+\sup_{\theta\in\Theta_{\iota}}\sup_{x,y\in
I_{\varepsilon}}\bigl{|}\frac{\partial p_{\theta}}{\partial
y}(x,y)\bigr{|}<\infty$

is finite, because $\frac{\partial p_{\theta}}{\partial y}(x,y)$ is continuous
and $\Theta_{\iota}\times I_{\varepsilon}^{2}$ is compact. The second term is
thus bounded by
$C\|\log|T^{\prime}_{\theta}|-\log|T_{\theta^{\prime}}^{\prime}|\|_{1}$ and,
by Lemma 6.7 below, goes to 0 as $\theta^{\prime}\to\theta$.

We now estimate the first term in (39). Since $\log|T^{\prime}_{\theta}|\in
L^{1}$, it is enough to bound $\|h_{\theta}-h_{\theta^{\prime}}\|_{\infty}$
which is again dominated by $\|h_{\theta}-h_{\theta^{\prime}}\|_{BV}$. By
invariance,

$\|h_{\theta}-h_{\theta^{\prime}}\|_{BV}=\|\mathcal{L}^{k}_{\theta}h_{\theta}-\mathcal{L}^{k}_{\theta^{\prime}}h_{\theta^{\prime}}\|_{BV}\leq\|\mathcal{L}^{k}_{\theta}(h_{\theta}-h_{\theta^{\prime}})\|_{BV}+\|(\mathcal{L}^{k}_{\theta}-\mathcal{L}^{k}_{\theta^{\prime}})h_{\theta^{\prime}}\|_{BV}.$

As we said in Remark 4.3, the Markov operator $\mathcal{L}_{\theta}$ enjoys
the exponential bound

$\|\mathcal{L}_{\theta}^{k}f\|_{BV}\leq
C_{\theta}\zeta_{\theta}^{k}\|f\|_{BV}$

for all $k>0$ and $f\in BV$ supported on $I_{\varepsilon}$ with $\int fdx=0$,
where the constants $C_{\theta}>0$, $0<\zeta_{\theta}<1$ depend on the
parameter $\theta$. Since $h_{\theta}-h_{\theta^{\prime}}$ has zero mean, we
therefore have

$\|\mathcal{L}^{k}_{\theta}(h_{\theta}-h_{\theta^{\prime}})\|_{BV}\leq
C_{\theta}\zeta_{\theta}^{k}\|h_{\theta}-h_{\theta^{\prime}}\|_{BV}.$

Then

$(1-C_{\theta}\zeta_{\theta}^{k})\|h_{\theta}-h_{\theta^{\prime}}\|_{BV}\leq\|(\mathcal{L}^{k}_{\theta}-\mathcal{L}^{k}_{\theta^{\prime}})h_{\theta^{\prime}}\|_{BV},$

and for $k$ sufficiently large, $C_{\theta}\zeta_{\theta}^{k}<1$. By a
standard trick, expanding a telescopic sum
$\mathcal{L}^{k}_{\theta}-\mathcal{L}^{k}_{\theta^{\prime}}=\sum_{j=1}^{k}\mathcal{L}_{\theta}^{k-j}(\mathcal{L}_{\theta}-\mathcal{L}_{\theta^{\prime}})\mathcal{L}_{\theta^{\prime}}^{j-1}$,
we get

$\|(\mathcal{L}^{k}_{\theta}-\mathcal{L}^{k}_{\theta^{\prime}})h_{\theta^{\prime}}\|_{BV}\leq\sum_{j=1}^{k}\|\mathcal{L}_{\theta}^{k-j}(\mathcal{L}_{\theta}-\mathcal{L}_{\theta^{\prime}})h_{\theta^{\prime}}\|_{BV}\\\
\leq\sum_{j=1}^{k}C_{\theta}\zeta_{\theta}^{k-j}\|(\mathcal{L}_{\theta}-\mathcal{L}_{\theta^{\prime}})h_{\theta^{\prime}}\|_{BV}\leq
C_{\theta}\frac{1}{1-\zeta_{\theta}}\|(\mathcal{L}_{\theta}-\mathcal{L}_{\theta^{\prime}})h_{\theta^{\prime}}\|_{BV}.$

Combining the above inequalities we come to the following estimate:

$\int\bigl{|}\log|T_{\theta}^{\prime}|\bigr{|}|h_{\theta}-h_{\theta^{\prime}}|dx\leq
M_{\theta}\|(\mathcal{L}_{\theta}-\mathcal{L}_{\theta^{\prime}})h_{\theta^{\prime}}\|_{BV},$

where
$M_{\theta}=\frac{C_{\theta}\|\log|T^{\prime}_{\theta}|\|_{1}}{(1-\zeta_{\theta})(1-C_{\theta}\zeta_{\theta}^{k})}$.
It therefore remains to bound
$\|(\mathcal{L}_{\theta}-\mathcal{L}_{\theta^{\prime}})h_{\theta^{\prime}}\|_{BV}$.
Since both $h_{\theta^{\prime}}$ and
$(\mathcal{L}_{\theta}-\mathcal{L}_{\theta^{\prime}})h_{\theta^{\prime}}$ are
supported on $I_{\varepsilon}$ and $\|h_{\theta^{\prime}}\|_{\infty}\leq C$,
we have

$\|(\mathcal{L}_{\theta}-\mathcal{L}_{\theta^{\prime}})h_{\theta^{\prime}}\|_{1}=\int_{I_{\varepsilon}}\int_{I_{\varepsilon}}|p_{\theta}(x,y)-p_{\theta^{\prime}}(x,y)|h_{\theta^{\prime}}(x)dxdy\leq
M_{1}\|\theta-\theta^{\prime}\|,$

where $M_{1}=C\sup_{\theta\in\Theta_{\iota}}\sup_{x,y\in
I_{\varepsilon}}\|\nabla_{\theta}p_{\theta}(x,y)\|$ is finite because
$\nabla_{\theta}p_{\theta}(x,y)$ is continuous and $\Theta_{\iota}\times
I_{\varepsilon}^{2}$ is compact. Similarly, arguing as in the proof of Lemma
4.1, we get

$|(\mathcal{L}_{\theta}-\mathcal{L}_{\theta^{\prime}})h_{\theta^{\prime}}|_{TV}\leq\int_{I_{\varepsilon}}|p_{\theta}(x,\cdot)-p_{\theta^{\prime}}(x,\cdot)|_{TV}h_{\theta^{\prime}}(x)dx\leq
M_{2}\|\theta-\theta^{\prime}\|,$

with $M_{2}=C\sup_{\theta\in\Theta_{\iota}}\sup_{x,y\in
I_{\varepsilon}}\|\nabla_{\theta}\frac{\partial p_{\theta}}{\partial
y}(x,y)\|<\infty$. Therefore, the first term in (39) is bounded by
$M_{\theta}(M_{1}+M_{2})\|\theta-\theta^{\prime}\|$. This finishes the proof.
∎

######  Remark 6.6.

Clearly, the above proof works if we replace $\log|T^{\prime}|$ with any
continuous function. Therefore, the mapping
$\tilde{\Theta}\ni\theta\mapsto\mu_{\theta}\in\mathcal{M}$ is continuous with
respect to the weak-* topology on $\mathcal{M}$, i.e.
$\tilde{\Theta}\ni\theta\mapsto\int\varphi\mu_{\theta}\in\mathbb{R}$ is
continuous for any $\varphi\in C^{0}(I)$. Theorem 6.5 is more delicate,
however, because $\log|T^{\prime}_{\theta}|$ is neither continuous nor bounded
and also depends on $\theta$.

######  Lemma 6.7.

$\|\log|T^{\prime}_{\theta}|-\log|T_{\theta^{\prime}}^{\prime}|\|_{1}\to 0$ as
$\theta^{\prime}\to\theta$.

###### Proof.

First notice that the critical point $c_{\theta}$ of the map $T_{\theta}$
depends continuously on the parameter $\theta\in\Theta$. Indeed, since
$T^{\prime}_{\theta}$ is continuous on $\Theta\times I$, the set
$\\{(\theta,c_{\theta})\\}=(T^{\prime}_{\theta})^{-1}(\\{0\\})$ is closed, and
then the map $\Theta\ni\theta\mapsto c_{\theta}\in I$ is continuous by the
closed graph theorem.

The functions $\log|T^{\prime}_{\theta}|$ and
$\log|T^{\prime}_{\theta^{\prime}}|$ have logarithmic singularities at
$c_{\theta}$ and $c_{\theta^{\prime}}$ respectively. We will show that these
singularities cancel out as $c_{\theta^{\prime}}$ approaches $c_{\theta}$. As
in the proof of Theorem 6.5, we may assume that
$\theta,\theta^{\prime}\in\Theta_{\iota}$, where $\Theta_{\iota}$ is compact.

Let $\alpha:=\sup_{\theta\in\Theta_{\iota}}|T^{\prime\prime\prime}_{\theta}|$.
As a direct consequence of the mean value theorem,

$(|T^{\prime\prime}_{\theta}(c_{\theta})|-2\alpha\delta)|x-c_{\theta}|\leq|T^{\prime}_{\theta}(x)|\leq(|T^{\prime\prime}_{\theta}(c_{\theta})|+2\alpha\delta)|x-c_{\theta}|$

for all $|x-c_{\theta}|\leq 2\delta$, and the same inequality holds if we
replace $\theta$ with $\theta^{\prime}$. Set
$D_{\delta}^{\pm}:=|T^{\prime\prime}_{\theta}(c_{\theta})|\pm 3\alpha\delta$;
both $D_{\delta}^{+}$ and $D_{\delta}^{-}$ are positive, since $T_{\theta}$
has quadratic critical point ($T^{\prime\prime}_{\theta}(c_{\theta})<0$). If
$\theta^{\prime}$ is sufficiently close to $\theta$, then
$|c_{\theta}-c_{\theta^{\prime}}|<\delta/2$ and
$|T^{\prime\prime}_{\theta}(c_{\theta})-T^{\prime\prime}_{\theta}(c_{\theta^{\prime}})|<\delta$,
and for all $|x-c_{\theta}|\leq\delta$ we then simultaneously have

$D_{\delta}^{-}|x-c_{\theta}|\leq|T^{\prime}_{\theta}(x)|\leq
D_{\delta}^{+}|x-c_{\theta}|,$
$D_{\delta}^{-}|x-c_{\theta^{\prime}}|\leq|T^{\prime}_{\theta^{\prime}}(x)|\leq
D_{\delta}^{+}|x-c_{\theta^{\prime}}|,$

and hence

$\bigl{|}\log|T^{\prime}_{\theta}|-\log|T_{\theta^{\prime}}^{\prime}|\bigr{|}\leq\log\frac{D_{\delta}^{+}}{D_{\delta}^{-}}+\Bigl{|}\log\frac{|x-c_{\theta}|}{|x-c_{\theta^{\prime}}|}\Bigr{|}.$

Given $\varepsilon>0$ and $\theta\in\Theta_{\iota}$, we first fix $\delta>0$
such that $\log\frac{D_{\delta}^{+}}{D_{\delta}^{-}}<\frac{\varepsilon}{3}$
and then let $\theta^{\prime}\to\theta$. The integral of the second term is
elementary and vanishes as $c_{\theta^{\prime}}\to c_{\theta}$, so
$\int_{B_{\delta}(c_{\theta})}\bigl{|}\log\frac{|x-c_{\theta}|}{|x-c_{\theta^{\prime}}|}\bigr{|}dx<\frac{\varepsilon}{3}$,
provided $\theta^{\prime}$ and $\theta$ are sufficiently close, and

$\int_{B_{\delta}(c_{\theta})}\bigl{|}\log|T^{\prime}_{\theta}|-\log|T_{\theta^{\prime}}^{\prime}|\bigr{|}dx<\frac{2\varepsilon}{3}.$
(40)

Let us show that
$\log|T^{\prime}_{\theta}|-\log|T_{\theta^{\prime}}^{\prime}|\to 0$ uniformly
on $B_{\delta}(c_{\theta})^{c}$. Denote $\beta:=\inf_{x\in
B_{\delta}(c_{\theta})^{c}}|T^{\prime}_{\theta}(x)|>0$. We have
$\|T^{\prime}_{\theta}-T^{\prime}_{\theta^{\prime}}\|_{\infty}\leq
M\|\theta-\theta^{\prime}\|$, where $M=\sup_{\theta\in\Theta_{\iota},x\in
I}\|\nabla_{\theta}T^{\prime}_{\theta}(x)\|$. Therefore $\inf_{x\in
B_{\delta}(c_{\theta})^{c}}|T^{\prime}_{\theta^{\prime}}(x)|>\frac{\beta}{2}$
if $\|\theta-\theta^{\prime}\|<\frac{\beta}{2M}$. Consequently,
$\bigl{|}\log|T^{\prime}_{\theta}|-\log|T_{\theta^{\prime}}^{\prime}|\bigr{|}<\frac{2M}{\beta}\|\theta-\theta^{\prime}\|$,
and for $\|\theta-\theta^{\prime}\|<\frac{\varepsilon\beta}{6M}$ we have

$\int_{B_{\delta}(c_{\theta})^{c}}\bigl{|}\log|T^{\prime}_{\theta}|-\log|T_{\theta^{\prime}}^{\prime}|\bigr{|}dx<\frac{\varepsilon}{3}.$
(41)

Combining (41) with (40) we finally get
$\|\log|T^{\prime}_{\theta}|-\log|T_{\theta^{\prime}}^{\prime}|\|_{1}<\varepsilon$.
∎

###  6.3 Random Lyapunov exponent and random entropy

In (36) we used the derivative of the deterministic map only. Alternatively,
if we define the process using the random transformations (24), we are led to
compute the Lyapunov exponent of the cocycle given by the derivative computed
along the random orbit, namely we define the random Lyapunov exponent (RLE)
$\bar{\Lambda}$ as

$\bar{\Lambda}:=\lim_{n\to\infty}\frac{1}{n}\log|D(T_{\eta_{n}}\circ\cdots\circ
T_{\eta_{1}})(x)|,$ (42)

for almost every sequence $(\eta_{k})\in[0,1]^{\mathbb{N}}$ with respect to
the measure $\operatorname{Leb}^{\otimes\mathbb{N}}$ (see Section 3.5), and
almost every $x\in I$ with respect to the stationary measure
$\mu_{\mathbbm{n}}$. By using the notation introduced in Section 3.5:
$\bar{x}_{k}(x):=T_{\eta_{k}}\circ\cdots\circ T_{\eta_{1}}(x)$, with
$\bar{x}_{0}(x)=x$, formula (42) is equal to

$\bar{\Lambda}=\lim_{n\to\infty}\frac{1}{n}\sum_{k=1}^{n}\log|T_{\eta_{k}}^{\prime}(\bar{x}_{k-1})|,$
(43)

again for $\operatorname{Leb}^{\otimes\mathbb{N}}$-a.e.
$(\eta_{k})\in[0,1]^{\mathbb{N}}$ and $\mu_{\mathbbm{n}}$-a.e. $x\in I$. By
using the ergodic theorem for random transformations, see [37] or [7, Section
3.1], we have that

$\bar{\Lambda}=\int\log|T^{\prime}_{\eta}(x)|d\mu_{\mathbbm{n}}(x)d\eta.$

Notice that if we compare this random exponent with $\Lambda$, we see that the
difference between the two is bounded by

$|\Lambda-\bar{\Lambda}|\leq\int\bigl{|}\log|T^{\prime}_{\eta}(x)|-\log|T^{\prime}(x)|\bigr{|}d\mu_{\mathbbm{n}}(x)d\eta.$

By using expression (28) for $T_{\eta}$, we can bound the error term in a more
explicit way as

$|\Lambda-\bar{\Lambda}|\leq\int\left|\log\Bigl{|}1+\frac{\tilde{q}_{\mathbbm{n}}(\eta)\sigma^{\prime}(x)}{T^{\prime}(x)}\Bigr{|}\right|h_{\mathbbm{n}}(x)dxd\eta.$

Since $\sigma^{\prime}$ is bounded over $I_{\varepsilon,\Gamma}$,
$\log|T^{\prime}|\in L^{p},p\geq 1$, and finally $\|h_{\mathbbm{n}}\|_{q}\leq
C_{q}$ for all $\mathbbm{n}$, we have

$|\Lambda-\bar{\Lambda}|\leq
C_{q}\int\left\|\log\Bigl{|}1+\frac{\tilde{q}_{\mathbbm{n}}(\eta)\sigma^{\prime}}{T^{\prime}}\Bigr{|}\right\|_{p}d\eta.$

We expect that this error converges to zero for large $\mathbbm{n}$, since the
quantile function converge to zero almost everywhere. Table 1 in Section 7.2
shows that, in fact, the error remains very small even for small values of
$\mathbbm{n}$.

It is known that whenever the random Lyapunov exponent (42) is positive, then
it equals the random entropy, which is the random generalization of the
Kolmogorov-Sinai entropy: this equality is called the entropy formula. Roughly
speaking, the random entropy computes the limit of the Shannon entropy of the
join partition generated by the successive application of the backward images
of the random maps on an initial generating partition, aka the entropy rate.
We defer to [37, Theorem 1.3] for the precise definition and [7, Theorem 3.2]
for the connection with the RLE. What is important for us is that the random
entropy coincides with the much easier object which is the RLE when the latter
is positive.

In Corollary 6.4 we proved that the Lyapunov exponent (36) may be negative,
and there is strong numerical evidence that the RLE (43) is also negative for
certain parameters; see Section 7.2. By assuming that the RLE is negative or
zero, we then get that the random entropy is also zero, by using another
important result connecting entropy and Lyapunov exponents, namely the
Margulis-Ruelle inequality (see for instance [54]), which states that the
metric entropy is bounded by the maximum between zero and the sum of the
positive Lyapunov exponents. The random version of this inequality, which we
use, has been proved by Kifer in [37, Theorem 1.4].

We should now stress the interesting fact that although the entropy is zero,
the Markov chain _mixes exponentially fast_ , as we pointed out in Remark 4.3.
This means that for any observables $f\in L^{1}$, $g\in BV$ the correlations
$\int(\mathcal{L}_{\mathbbm{n}}^{k}f)(x)g(x)dx$ converge to $\int
f(x)d\mu_{\mathbbm{n}}(x)\int g(x)dx$ exponentially fast when $k\to\infty$.
This result can be stated in a more suggestive way by relying on random
transformations. By using the notations introduced in Section 3.5, we can in
fact rewrite the previous correlation in terms of composition of randomly
chosen maps and say that there exist $0<v<1$ and $C>0$, depending only on the
system, such that, for all $k\geq 0$,

$\left|\iint f(x)g(T_{\eta_{k}}\circ\dots\circ
T_{\eta_{1}}(x))d\overline{\eta}dx-\int fd\mu_{\mathbbm{n}}\int
g(x)dx\right|\leq Cv^{k}\|f\|_{1}\|g\|_{BV}.$

This exponential decay of correlations is a consequence of the spectral gap
prescribed by the quasi-compactness of the Markov operator proved in
Proposition 4.2, and of the uniqueness of the absolutely continuous stationary
measure, see [10] for details; of course these properties hold even when the
Lyapunov exponent is negative or zero.

##  7 Numerical results

In this section we describe and discuss some numerical experiments in support
of our rigorous theoretical investigations. Specifically, we present the
bifurcation diagram associated with the unimodal map (18) and compute the
corresponding Lyapunov exponents for both the deterministic and the stochastic
version of the map, see Subsections 7.1 and 7.2, respectively.121212The code
to reproduce all the figures and tables of this section are available from the
corresponding author upon reasonable request.

###  7.1 Bifurcation diagram

In this subsection we analyse the dynamics of the unimodal map (18). The
bifurcation diagram of a dynamical system shows how the asymptotic
distribution of a typical orbit varies as a function of a parameter. For our
map either the memory parameter $\omega$ or the parameter $\phi^{*}$ can be
employed as bifurcation parameter. Fig. 7 shows the bifurcation diagram as a
function of $\phi^{*}$ (the bifurcation diagram as a function of $\omega$
looks similar). The choice of the parameters for this plot corresponds to a
vertical segment in the parameter space (see Fig. 2) with $\omega=0.5$ and
$\phi^{*}$ varying in a neighbourhood of the dynamical core area. As explained
in Section 3.1, Theorem A, the invariant set of the unimodal map (18) could be
an attracting periodic orbit, a Cantor set of measure zero or a finite union
of intervals with a dense orbit, depending on the parameters $\omega$ and
$\phi^{*}$. Specifically, there is an attracting fixed point or a 2-cycle
outside the dynamical core region, whereas in the dynamical core the situation
is more complex as small parameter variations can change the dynamics from
chaotic to periodic and back, as we see in Fig. 7. To identify more precisely
the signature of a chaotic behaviour, in the next subsection we compute the
Lyapunov exponent as a function of $\phi^{*}$.

Figure 5: Bifurcation diagram for $T$ in the dynamical core region
($\omega=0.5$).

Figure 6: Lyapunov exponents for deterministic and stochastic maps
($\omega=0.5$).

$\omega$$\phi^{*}$

(a)

$\phi^{*}$

(b)

Figure 7:  Contour plot of the average Lyapunov exponents for (a)
deterministic and (b) stochastic maps ($\mathbbm{n}=1$).

###  7.2 Lyapunov exponent

The Lyapunov exponent for the deterministic map (18) is positive if and only
if $T$ admits an absolutely continuous invariant measure, see Theorem C in
Section 3.1. Fig. 7 shows the estimated Lyapunov exponent in the same slice of
the parameter space as in Fig. 7. Notice that the exponent becomes a smooth
function of $\phi^{*}$ when we add even a very small noise, in agreement with
the results of Section 6.2. Fig. 7 shows the contour plot of the Lyapunov
exponent as a function of $\phi^{*}$ and $\omega$ both for the deterministic
map (18) and for the stochastic process described in Section 3.4. The right
panel shows a clear noise-induced regularization phenomenon. In fact, for the
stochastic version of the map the intricate fine structure in the parameter
dependence of the Lyapunov exponent disappears and is replaced by a smooth
dependence.

To provide a numerical exemplification of the stochastic stability and of
Theorem 6.3, we computed average (36) and random (43) Lyapunov exponents, as
well as the Lyapunov exponent for the deterministic map (18), for different
values of $\phi^{*}$, $\omega$ and $\mathbbm{n}$. The results are presented in
Table 1. Within each row, the two subrows are ALE and RLE, respectively: the
two agree very well, in most cases up to the precision of the numerical
computation. In both cases we sampled 128 independent realizations of the
process, each of length $10,000$.

$\phi^{*}$ | $\omega$ | DC | Per. | Lyapunov exponent
---|---|---|---|---
$\mathbbm{n}$ | Det.
$1$ | $10^{3}$ | $10^{6}$ | $10^{9}$
0.845 | 0.557 | yes | no | 0.287 | 0.287 | 0.349 | 0.341 | 0.340
0.286 | 0.286 | 0.349 | 0.341
0.795 | 0.390 | yes | no | 0.345 | 0.346 | 0.389 | 0.398 | 0.400
0.345 | 0.346 | 0.389 | 0.399
0.904 | 0.627 | yes | no | 0.557 | 0.557 | 0.560 | 0.550 | 0.552
0.558 | 0.557 | 0.560 | 0.550
0.821 | 0.439 | yes | yes | 0.378 | 0.375 | $-0.051$ | $-0.159$ | $-0.158$
0.378 | 0.375 | $-0.052$ | $-0.159$
0.944 | 0.826 | yes | yes | 0.296 | 0.297 | 0.049 | $-0.123$ | $-0.122$
0.297 | 0.296 | 0.052 | $-0.123$
0.766 | 0.323 | yes | yes | 0.320 | 0.324 | 0.286 | $-0.076$ | $-0.046$
0.320 | 0.325 | 0.286 | $-0.076$
0.258 | 0.837 | no | yes | $-0.243$ | $-0.248$ | $-0.248$ | $-0.248$ | $-0.248$
$-0.286$ | $-0.248$ | $-0.248$ | $-0.248$
0.908 | 0.804 | no | yes | $-0.284$ | $-0.285$ | $-0.365$ | $-0.362$ | $-0.362$
$-0.286$ | $-0.287$ | $-0.366$ | $-0.362$
0.541 | 0.227 | no | yes | $-0.619$ | $-0.441$ | $-0.380$ | $-0.380$ | $-0.380$
$-0.578$ | $-0.446$ | $-0.380$ | $-0.380$
Table 1: Average and random Lyapunov exponents for different values of
$\mathbbm{n}$ compared to the Lyapunov exponent for the deterministic map. DC
stands for dynamical core, Per. for periodic, and Det. for the Lyapunov
exponent of the deterministic map.

##  8 Estimating the map parameters via deep neural networks

We now consider the problem of estimating the parameters of the map from
(short) time series. This is motivated by the fact that in Section 9 we
empirically investigate a dataset of US commercial banks leverage. We will
consider the time series of bank’s leverage as realizations of the process
described in Section 3.4 and will estimate for each bank the model parameters.
Each time series is very short, being composed by only 59 points.

Given the random nature of the map, one could use maximum likelihood
estimation to estimate the parameters. However, this approach is not feasible
for two reasons. First, the likelihood function is highly non-convex, so that
standard optimization methods may perform poorly. Second, although the
likelihood function for the process itself can be written explicitly, in many
cases it may happen that the observed time series are systematically
undersampled and this prevents an explicit calculation of the likelihood
function. For example, we may observe only one slow time scale, corresponding
to portfolio rebalancing, out of two, or even out of three (i.e. the bank’s
risk assessment and portfolio composition may be updated more frequently than
our quarterly observations, for instance at a monthly frequency). If we
observe, for instance, only the second iterate of the process,

$\phi_{t+2}=T(T(\phi_{t};\theta);\theta)+\sigma(\phi_{t};\theta)\epsilon_{t})+\sigma(\phi_{t+1};\theta)\epsilon_{t+1},\quad
t\in\mathbb{Z},$

the transition probabilities $p(\phi_{t+2}|\phi_{t};\theta)$,
$t\in\mathbb{Z}$, are no longer Gaussian (as it would be the case if we
observe the first iterate). Hence, there is apparently no effective formula
for the likelihood function.

For this reason, to estimate the parameters of the map, we propose to use a
convolutional neural network (CNN) consisting of a sequence of convolutional
layers followed by a sequence of dense, or fully connected, layers; we refer
to [33, Chapter 9] for detailed exposition of CNNs. In order to deal with the
possibility that the observed time-series are realizations of certain iterates
of the process, we separately optimize two CNN architectures to be used
sequentially. First, we optimize a CNN (henceforth denoted by $\mathsf{CNN1}$)
for estimating the number of iterates between observations: it takes as input
time series of length 59 and outputs the corresponding value $k$ of the map’s
iterate that generated the time-series. Second, for each value of $k$, we
optimize a CNN (henceforth denoted by $\mathsf{CNN2}(k)$), having the exact
same inputs, to output the corresponding parameters $(\phi^{*},\omega)$ that
generated the time-series. Once the parameters $(\phi^{*},\omega)$ are
estimated, the variance of the noise (and therefore $\mathbbm{n}$) can be
estimated by standard methods. To train $\mathsf{CNN1}$ and $\mathsf{CNN2}(k)$
we used a training set of one million samples simulated from the model in
Section 3.4 with values of the parameters
$\theta=(\phi^{*},\omega,\mathbbm{n})$ which uniformly span the parameter
space. For both steps, when simulating the series, the initial state of the
system was taken randomly from a uniform distribution on $[0,1]$. This is
especially important because of the relatively short length of the series at
hand. Therefore, being based only on simulations, the NN approach, contrary to
the maximum likelihood one, can work also for partial observations.

The architectures of $\mathsf{CNN1}$ and $\mathsf{CNN2}(k)$ are schematized in
Fig. 8. In general, a convolutional layer is composed of $n_{f}$ filters and
each filter is associated with one kernel that is applied to a small moving
window of the time-series; for instance, in our first convolutional layer
$n_{f}=128$ and all the windows are of width $2$. The outputs of one
convolutional layer are connected to the next layer. The weights of these
connections constitute the NN parameters to be optimized. After seven
convolutions, the output is passed to a sequence of dense layers, which
concludes the NN. We use the rectified linear unit (ReLU) function as
activation function. Now, we describe the experimental setup.

Model CNN1: "convolutional_categorical_model"
_________________________________________________________________ Layer (type)
Output Shape Param #
================================================================= reshape
(Reshape) (None, 59, 1) 0
_________________________________________________________________ conv1d_1
(Conv1D) (None, 58, 128) 384
_________________________________________________________________ conv1d_2
(Conv1D) (None, 29, 64) 16448
_________________________________________________________________ conv1d_3
(Conv1D) (None, 15, 64) 8256
_________________________________________________________________ conv1d_4
(Conv1D) (None, 8, 64) 8256
_________________________________________________________________ conv1d_5
(Conv1D) (None, 4, 64) 8256
_________________________________________________________________ conv1d_6
(Conv1D) (None, 2, 64) 8256
_________________________________________________________________ conv1d_7
(Conv1D) (None, 1, 64) 8256
_________________________________________________________________ flatten
(Flatten) (None, 64) 0
_________________________________________________________________ dense_1
(Dense) (None, 128) 8320
_________________________________________________________________ dense_2
(Dense) (None, 64) 8256
_________________________________________________________________ dense_3
(Dense) (None, 3) 195
================================================================= Trainable
params: 74,883 Model CNN2: "convolutional_model" ...
_________________________________________________________________ dense_3
(Dense) (None, 2) 130
================================================================= Trainable
params: 74,818

Figure 8: Architectures of the $\mathsf{CNN1}$ model used to estimate the
iterate $k$ and the $\mathsf{CNN2}(k)$ model used to estimate the parameters
$(\phi^{*},\omega)$ for each $k$. The two models differ only in the output
layer.

The implementation is carried out in Python. To generate training and testing
data we simulate one million samples. $\mathsf{CNN1}$ is optimized with the
stochastic gradient descent method by using the Adam algorithm [38], the
categorical cross-entropy as loss function, and the accuracy as metric, with
$L^{2}$ regularization of weights. Instead, to optimize $\mathsf{CNN2}$ we use
the Mean Squared Error (MSE) both as loss function and as metric. The batch
size is $32$ in both cases. The seven convolutional and three dense layers
have a total of $74,818$ trainable parameters.

The CNN models show a good performance. We tested our methods on a testing set
of 100,000 out of sample time series. Fig. 9 shows the accuracy of
$\mathsf{CNN1}$ to estimate the iterates on test data. We choose $k=1,2,3$
because of our empirical application of Section 9. The MSE of
$\mathsf{CNN2}(k)$ on the test set is about 0.001 for each $k$. Since both
$\phi^{*}$ and $\omega$ are uniformly distributed in $[0,1]$, the MSE is quite
small and the NN effective.

Figure 9:  Accuracy of the $\mathsf{CNN1}$ model used to estimate the
iterates.

##  9 Chaos in real bank leverage time series

In this section we perform an empirical analysis of a large set of bank’s
leverage time series. We first describe the the data set. Next, through
$\mathsf{CNN1}$ and $\mathsf{CNN2}(k)$, we estimate the observed iterate $k$
as well as the parameters $(\phi^{*},\omega)$ of the model and discuss the
results, investigating the relation between the estimated parameters and the
bank’s size. Finally, to perform an independent analysis supporting our
conclusion, we apply the Chaos Decision Tree Algorithm [57] to these time-
series and compare the resulting classification with the one obtained with NN
estimates.

###  9.1 Data set

We use the data set of US Commercial Banks and Savings and Loans Associations
provided by the Federal Financial Institutions Examination Council (FFIEC).
For the sake of completeness, we provide here a description of it, referring
to [28] and references therein for further details. A Commercial Bank is
defined131313See http://www.ffiec.gov/nicSearch/FAQ/Glossary.html. officially
by the FFIEC as: “[…] a financial institution that is owned by stockholders,
operates for a profit, and engages in various lending activities”. Commercial
banks quarterly fill the Consolidated Report of Condition and Income
(generally referred to as Call Report) as required by the FFIEC. A Savings and
Loan Associations, instead, is a financial institution that accepts deposits
primarily from individuals and channels its funds primarily into residential
mortgage loans. Starting from the first quarter of 2012 they are required to
file the same reports of Commercial Banks, thus they are included in the data
set since then. The data provided by the Call reports are publicly available
since 1986 although the level of details required has increased over time. To
have a good compromise between the fine structure of data and a reasonably
populated statistics we follow [28] and consider the time period going from
March 2001 to December 2014, for a total of 59 quarters. Also we consider only
the financial institutions that are present in the data set in all the
quarters for a total of $5,031$ banks. The financial leverage $\lambda_{t}$ of
each institution at time $t$ is defined as the ratio between the sum of its
assets and its equity at time $t$. In particular, the latter is given by
$E_{t}=A_{t}-L_{t}$ where $L_{t}$ represents the liabilities and $A_{t}$ the
assets of the bank, thus $\lambda_{t}={A_{t}}/{E_{t}}$.

###  9.2 Estimation via neural networks

In order to estimate the parameters of the map on the just-described data set,
we need to fix the value of the liquidity parameter $\gamma$; remind that we
consider the linear transformation $\phi_{t}=(\lambda_{t}-1)/\gamma$. In this
work, we assume that the liquidity parameter $\gamma$ of the risk investment
is the same for all the banks in our data set. Admittedly, this is a
simplifying assumption, coherent with the so-called assumption of statistical
equivalence for risky investments (see also [44]), which allows for an
analytical tractability of the model; a complete exploration of a relaxation
of this hypothesis is beyond the scopes of the present paper and is,
therefore, left for future work. In order to fix its value, we exclude $662$
time series (out of $5,031$) that contain outliers, which we define to be
values that are two standard deviations away from the mean. We then set
$\gamma$ to the maximum over the remaining $4,369$ series, obtaining
$\gamma=15.969$.

Since the time series that we analyze contain quarterly data and portfolio
decisions may be made more frequently, it is natural to assume that the
observed time series are realizations of certain iterates of the process; we
assume $k\in\\{1,2,3\\}$. Fig. 10(a) displays the output of $\mathsf{CNN1}$.
It turns out that only a small percentage (about $1\%$) of the banks in our
data set rebalance their portfolios at a quarterly frequency. Most banks seem
to rebalance either every six weeks ($k=2$, about $55\%$) or every month
($k=3$, about $43\%$). One may ask if the portfolio re-balancing frequency is
related to the size of the bank (defined as the average across the 59 quarters
of the sum of the dollar-amount of all the type of assets detained by it), for
example because larger banks manage more actively their portfolio. Fig. 10(b)
shows the box plots of the logarithm of the size of the banks for $k=1,2,3$.
We observe that there is not a statistically significant difference among
them.

(a) Number of banks by iterate.

(b) Log bank size by iterate.

Figure 10: Results on the iterate of the process.

Once the number of iterates $k$ has been identified, we proceed to the study
of the chaotic behaviour of the time series. We divide the banks in the three
groups identified by $k$ and employ $\mathsf{CNN2}(k)$ in order to estimate
the parameters $(\phi^{*},\omega)$. In Fig. 11(a) we plot the estimates of
$\phi^{*}$ against those of $\omega$; pairs belonging to the dynamical core
region are displayed in red, whereas those falling outside the dynamical core
region are displayed in blue. Interestingly, the percentage of banks for which
the estimates $(\phi^{*},\omega)$ are in the dynamical core region is about
$12\%$. Moreover, Fig. 11(b) indicates that $k$ is very often equal to two for
these banks.

(a) Parameters distribution by DC.

(b) Parameters distribution by iterate.

Figure 11:  Estimated parameters for iterates in $\\{1,2,3\\}$.

We now ask if bank size is related with the fact that the estimated pair
$(\phi^{*},\omega)$ is or not in the dynamical core region. Fig. 12(a) shows
the probability density functions of the logarithm of the banks size, by
considering separately banks inside and outside the dynamical core, and Fig.
12(b) displays the corresponding probability-probability plot. To test that
the difference between the distribution of banks sizes in and outside the
dynamical core region is statistically significant we perform the Kolmogorov-
Smirnov test of the null hypothesis that the two samples have the same
distribution. The statistics of the test is $9.31\times 10^{-2}$ corresponding
to a $p$-value of $5.9\times 10^{-4}$. This latter value shows that the two
subsamples have different distributions.

Summarizing, we have found that the parameters of a sizable fraction of banks
lie in the dynamical core region and that the dynamics of the leverage of the
larger banks tends to be more frequently in the dynamical core than that of
the smaller banks.

(a) Probability density functions.

(b) Probability-probability plot.

Figure 12: Results on the parameters of the process and bank size.

###  9.3 Classification via the Chaos Decision Tree Algorithm

Finally, we perform an independent analysis on the bank’s leverage time series
by making use of the recently proposed Chaos Decision Tree Algorithm (CDTA)
[57], described in detail in Appendix A. This is a non-parametric method which
classifies an input time series as chaotic, periodic, or
stochastic141414Notice that the definitions of _chaotic_ and _periodic_ for
the CDTA differ from the ones given in Section 3.1, see Appendix A.. We
perform this analysis for two reasons. First, we know that chaotic behavior
can be present only for series generated by our map with parameters in the
dynamical core. Thus we test whether the series classified as chaotic by CDTA
have estimated parameters in the dynamical core. The second reason is to count
how many banks in the dynamical core are identified as chaotic or periodic by
CDTA. The Appendix also contains the results of the application of CDTA to
data simulated by our map for different time series length, level of noise
$\mathbbm{n}$, and number of iterates $k$.

Applying CDTA, we find that $64\%$ of the banks are classified as stochastic,
$\sim 12\%$ as periodic and $\sim 23\%$ as chaotic. The consistency between
the classification made by CDTA and the partition ‘dynamical core’ and ‘not
dynamical core’ of the parameters space $(\phi^{*},\omega)$ found by the NN
model can be assessed by looking at Table 2. We find that a large fraction of
series outside the dynamical core are classified as stochastic by CDTA, while
a third of banks in the dynamical core is classified as chaotic. This fraction
is significantly smaller for banks outside the dynamical core. Thus, despite
the agreement is not perfect, we find a reasonable consistency between the
conclusions of the two methods and, more importantly, find significant (and
independent) support to the conclusion that a sizable fraction of bank time
series are described by a chaotic dynamics.

| Periodic | Chaotic | Stochastic
---|---|---|---
Non dynamical core | 382 (9.98%) | 648 (16.93%) | 2798 (73.09%)
Dynamical core | 107 (20.34%) | 176 (33.46%) | 243 (46.20%)
Table 2: Number of banks by classes.

Finally, the findings reported in the previous subsection suggest a positive
relation between the size of a bank and the probability that the dynamics of
the (corresponding) leverage time series is chaotic. To verify this
observation by using the CDTA classification, we first rank the banks in
quintiles according to their size and within each quintiles we compute the
percentage of banks that are detected to be stochastic, periodic and chaotic.
Table 3 collects the results. In a nutshell, banks having a larger size have,
on average, a larger percentage of leverage time series detected as chaotic. A
$\chi^{2}$-test applied to contiguous quintiles rejects the hypothesis of
independence of the CDTA classification from the quintile, indicating that the
difference in frequencies across quintiles are statistically significant. Thus
also the CDTA analysis confirms that larger banks are more likely
characterized by chaotic time series of leverage.

Statistics | q1 | q2 | q3 | q4 | q5
---|---|---|---|---|---
Chaotic (%) | 18.2 | 20.2 | 21.7 | 23.6 | 29.1
Periodic (%) | 12.9 | 11.6 | 12.9 | 13 | 11.8
Stochastic (%) | 68.6 | 67.3 | 65 | 63 | 58.9
Table 3: Fraction of banks classified as chaotic, periodic, or stochastic by
CDTA conditionally to the decile of the bank size.

##  10 Conclusions

Most risk management practices (as, for example, Value-at-Risk) assume that
prices are not affected by actions of other financial institutions that are
managing the risk of their portfolio. In other words, these practices assume
that risk is exogenous. In reality, in the presence of limited liquidity,
coordinated and homogeneous risk management can create market instability and
result in what is known as endogenous risk. This has the potential to amplify
market instabilities and create crashes through the well-known feedback
between leverage, risk, and asset prices. An additional, and less considered,
feedback between past and future risks is present because financial
institutions use historical data to estimate both the riskiness of their
investments and their correlations. This creates new threats for the systemic
stability of financial markets. Studying how these two feedbacks affect the
leverage dynamics is of paramount importance for understanding systemic risk.

In this paper we consider a stylized model where both feedbacks are present.
We showed that the dynamics of the bank’s leverage is described by a unimodal
map on $[0,1]$ perturbed with additive and heteroscedastic noise. The
perturbed system can be described in two equivalent ways as a stationary
Markov chain or in terms of random transformations. In both cases a
fundamental object is the stationary measure of the process which allow us to
properly define and state all the statistical properties of the system. We are
able to construct such a measure and to prove its uniqueness. Moreover we
show, under a few assumptions, the stochastic stability of the perturbed
system, namely the weak convergence of the stationary measure to the invariant
measure for the unimodal map in the zero noise limit. We also define an
average Lyapunov exponent, still in terms of the stationary measure, as a
sensible indicator of the slow motion, and prove its continuity with respect
to the parameters defining the system. We show that, depending on the
parameters, the average Lyapunov exponent can be either negative or positive,
leading to two qualitatively different (periodic- and chaotic-like) dynamics
of the leverage.

We then estimate the parameters of the map via a method based on deep neural
networks, whose efficiency was tested in a large testing set. Assuming the
proposed unimodal map with heteroscedastic noise as data generating process
for the banks leverage, we estimated the parameters on quarterly data of about
5,000 US Commercial Banks via the proposes CNN architecture. By investigating
the period from March 2001 to December 2014, for a total of 59 quarters, we
found that the parameters of a sizable fraction of banks lie in the dynamical
core region of the parameter space and that the large banks’ leverage tends to
be more chaotic than the one of small ones. The latter finding was
corroborated also by using a non-parametric approach.

We believe that the proposed methodologies may offer revealing perspectives
for future works. For instance, it would be interesting to extend the employed
mathematical techniques to study a model in which more than one asset and one
bank are present in the system.

##  11 Acknowledgements

This research was supported by the research project ‘Dynamics and Information
Research Institute - Quantum Information, Quantum Technologies’ within the
agreement between UniCredit Bank and Scuola Normale Superiore. F.L. and S.M.
acknowledge partial support by the European Program scheme
‘INFRAIA-01-2018-2019: Research and Innovation action’, grant agreement
#871042 ’SoBigData++: European Integrated Infrastructure for Social Mining and
Big Data Analytics’. S.V. is grateful to the Centro di Ricerca Matematica
Ennio de Giorgi, Scuola Normale Superiore, Laboratoire International Associé
LIA LYSM, the INdAM (Italy) and the UMI-CNRS 3483, Laboratoire Fibonacci
(Pisa), where this work has been initiated and completed under a CNRS
delegation, for the support. S.V. thanks C. Gonzalez-Tokman for useful
discussion about Section 6.3.

##  Appendix A The Chaos Decision Tree Algorithm

The Chaos Decision Tree Algorithm [57] is a non-parametric chaos-detection
tool which has been developed with the goal of being especially robust to
measurement noise. It provides an automated processing pipeline which has been
showed to be able to detect the presence (or absence) of chaos in noisy
recordings, even for difficult edge cases. We use it in our work to identify
periodic/chaotic time series without any reference to our model. It is meant
to provide an independent check on the existence of chaotic behaviour in
leverage time series and to support the evidence that it may depend on the
bank’s size.

The algorithm classifies a time-series as either stochastic or periodic, or
chaotic. We now briefly explain how the algorithm works. The first step is to
test if data are stochastic. This is done via surrogate-based approach by
comparing the permutation entropy of the original time-series to the
permutation entropy of random surrogates of that time-series by using a
combination of Amplitude Adjusted Fourier Transform surrogates and Cyclic
Phase Permutation surrogates. If the permutation entropy of the original time-
series falls within either surrogate distribution, the time-series is
classified as stochastic. If the permutation entropy falls outside the
surrogate distribution, then the algorithm proceeds to de-noise the inputted
signal by using the Schreiber’s noise-reduction algorithm [55]. Notice that
the calculation of the permutation entropy relies on two parameters: the
permutation order and the time-lag. The time-lag has been set to 1 as
suggested in [57]. The choice of the order of the permutation is made in order
to maximize the detection of chaotic series in our model.

To test the CDTA algorithm on our model, we first consider the deterministic
map and generate $100$ chaotic and $100$ periodic time-series of length 59
from the dynamical core. We then apply CDTA to these series for different
values of the permutation order ($\in\\{3,\ldots,8\\}$). While the periodicity
accuracy is maximized for a value equal to 8 ($92\%$ of the periodic series
are correctly detected as such and the remaining $8\%$ are labeled chaotic),
the chaos detection accuracy is maximized for a choice of the permutation
order equal to 5 ($65\%$ of the series are correctly detected as chaotic and
the remaining periodic). Because of the purpose of this paper, we fix the
permutation order to 5, but we have checked that the conclusions of our data
analysis (in particular Table 4) would have been the same with a different
choice.

At this point, the algorithm checks for signal oversampling and, if the data
are over-sampled, the algorithm iteratively down-samples the data until they
are no longer over-sampled. Finally, CDTA performs the 0-1 chaos test [34] on
the input data. Ref. [57] points out that the 0-1 chaos test has been modified
from the original one to be less sensitive to noise. Then, it suppresses the
correlations arising from quasi-periodicity, and normalizes the standard
deviation of the test signal. The value for the parameter that suppresses
signal correlations is chosen based on ROC analyses. The modified 0-1 test
provides a single statistic, K, which approaches 1 for chaotic systems and
approaches 0 for periodic systems. The algorithm sets up a cutoff for K based
on the length of the time-series. If K is greater than the cutoff, the data
are classified as chaotic, and if they are smaller than or equal to the
cutoff, they are classified as periodic.

###  A.1 Simulations

Iterate | Series | $\mathbbm{n}$ | Dynamical core | Not dynamical core
---|---|---|---|---
length | S (%) | P (%) | C (%) | S (%) | P (%) | C (%)
1 | 59 | 5 | 33.5 | 5.14 | 61.4 | 57.4 | 4.81 | 37.7
20 | 29.3 | 4.28 | 66.5 | 69.9 | 3.28 | 26.8
100 | 24.3 | 6.99 | 69.6 | 88.1 | 2.94 | 8.97
295 | 5 | 2.2 | 1.7 | 96.1 | 22.5 | 6.24 | 71.2
20 | 0.1 | 1.9 | 98 | 43.8 | 10.4 | 45.8
100 | 0 | 2.3 | 97.7 | 73.8 | 8.35 | 17.9
590 | 5 | 0 | 0.7 | 99.3 | 13.1 | 6.08 | 80.9
20 | 0 | 0.4 | 99.6 | 33.6 | 8.5 | 57.9
100 | 0 | 0.4 | 99.6 | 66.4 | 8.04 | 25.6
1180 | 5 | 0 | 0.1 | 99.9 | 10.9 | 3.44 | 85.6
20 | 0 | 0 | 100 | 27.7 | 5.57 | 66.8
100 | 0 | 0 | 100 | 60.2 | 5.29 | 34.5
2 | 59 | 5 | 75.7 | 2.17 | 22.2 | 83.8 | 1.33 | 14.9
20 | 80.1 | 1.65 | 18.2 | 92.8 | 0.26 | 6.91
100 | 86.6 | 1.43 | 11.9 | 96.6 | 0.56 | 2.81
295 | 5 | 39.4 | 0 | 60.6 | 40.5 | 3.26 | 56.2
20 | 38.6 | 0.6 | 60.8 | 70 | 3.7 | 26.3
100 | 21 | 1.2 | 77.8 | 83.9 | 2.96 | 13.1
590 | 5 | 27.6 | 0 | 72.4 | 25.6 | 3.23 | 71.2
20 | 10.6 | 0 | 89.4 | 52.3 | 4.24 | 43.4
100 | 4.8 | 0.6 | 94.6 | 74.7 | 2.88 | 22.4
1180 | 5 | 11 | 0 | 89 | 13.7 | 2.22 | 84.1
20 | 0.2 | 0 | 99.8 | 39.4 | 2.82 | 57.7
100 | 0.2 | 0 | 99.8 | 64.1 | 2.24 | 33.7
3 | 59 | 5 | 85.3 | 0.47 | 14.2 | 83.9 | 1.33 | 14.9
20 | 84.5 | 1.17 | 14.3 | 92.8 | 0.26 | 6.91
100 | 89.5 | 1.43 | 8.79 | 96.6 | 0.56 | 2.81
295 | 5 | 36 | 0 | 64 | 40.5 | 3.26 | 56.2
20 | 32.6 | 0.4 | 67 | 70 | 3.7 | 26.3
100 | 23 | 1.6 | 75.4 | 83.9 | 2.96 | 13.1
590 | 5 | 20.6 | 0 | 79.4 | 25.6 | 3.23 | 71.2
20 | 9.6 | 0 | 90.4 | 52.3 | 4.24 | 43.4
100 | 9.4 | 0.8 | 89.8 | 74.7 | 2.88 | 22.4
1180 | 5 | 6 | 0 | 94 | 13.7 | 2.22 | 84.1
20 | 2 | 0 | 98 | 39.4 | 2.82 | 57.7
100 | 2.8 | 0.8 | 96.4 | 64.1 | 2.24 | 33.7
Table 4: Percentage of time-series detected by the Chaos Decision Tree
Algorithm as stochastic (S), periodic (P) and chaotic (C) in the dynamical
core and its complement as a function of the time-series length and
$\mathbbm{n}$ for different iterates of the map.

In this section we present some numerical investigations showing how CDTA
performs when simulating (noisy) time-series from the map described in Section
3.4. Specifically, we run the following two numerical experiments:

  1. (i)

First, we simulate time-series from the dynamical core area (i.e. time-series
for which the pairs $(\phi^{*},\omega)$ satisfy condition (C3):
$T(\Delta)<c<\Delta<1$). More precisely, we simulate $500$ samples of
different length and level of noise, which is captured by the variable
$\mathbbm{n}$.

  2. (ii)

Second, we simulate time-series from outside the dynamical core area (the map
$T$ thus satisfies (C1) or (C2)). The remaining simulation setting coincides
with that in (i).

The procedures explained in (i)–(ii) are also repeated when $T$ is replaced by
its $k$-th iterate $T^{k}$, $k=2,3$; see discussion in Section 8. Table 4
collects the results. We observe that when both $\mathbbm{n}$ and the time
series length are large, CDTA classifies almost all the time series in the
dynamical core as chaotic, while those outside it are never classified as
such. This is quite independent on $k$. By decreasing either $\mathbbm{n}$ or
the time series length, the classification is less precise and this effect is
stronger for larger values of $k$. In the regime of length comparable with our
empirical data ($N=59$), roughly a third (for $k=1$) or up to $85\%$ (for
$k=3$) of the time series in the dynamical core are classified as stochastic,
showing the limits of the CDTA when the time series are short and/or the noise
is large.

## References

  * [1] Viral V. Acharya and Stephen G. Ryan. “Banks’ Financial Reporting and Financial System Stability”. In: Journal of Accounting Research 54.2 (2016), pp. 277–340.
  * [2] Tobias Adrian and Hyun Song Shin. “Liquidity and leverage”. In: Journal of Financial Intermediation 19.3 (2010), pp. 418–437.