\documentclass[oneside,10pt]{book}

\usepackage{amsmath} % for "\cfrac" macro

\usepackage[export]{adjustbox}

\usepackage{relsize}

\usepackage{array}

\usepackage{enumitem}

\setcounter{tocdepth}{3}

\setcounter{secnumdepth}{3}

\pagestyle{plain}

\newcommand\Chapter[2]{

%\chapter[#1: {\itshape#2}]{#1\\[2ex]\Large\itshape#2}

\chapter[#1]{#1\\[2ex]\Large\itshape#2}

}

%\usepackage{fontspec}

%\setmainfont{Times New Roman} %Times New Roman

%\setmonofont{Consolas}

%\usepackage{selinput}

%\SelectInputMappings{Euro={€}}

%\usepackage[utf8]{inputenc}

\usepackage[titles]{tocloft}

\setlength{\cftbeforechapskip}{7pt} %%%%%%%%%% 6pt

\usepackage{longtable}

\usepackage{wrapfig}

\usepackage{amsmath} % need for subequations

\usepackage{amsfonts}

\usepackage{amssymb} % needed for mathbb OK

\usepackage{bigints}

\usepackage{graphicx} % need for figures

\usepackage{subfig}

\usepackage{verbatim} % useful for program listings

\usepackage{color} % use if color is used in text

%\usepackage{subfigure} % use for side-by-side figures

\usepackage{parskip}

\usepackage{float}

\usepackage{courier}

%\usepackage{artemisia} %%%

\usepackage{exercise}

\usepackage{sistyle}

\usepackage{textcomp}

%

%%%\usepackage[utf8]{luainputenc}

%\usepackage{luatextra}

%

%%\usepackage[utf8]{inputenc}

%%\usepackage[T1]{fontenc}

%%\usepackage{textcomp,upgreek}

%\usepackage{fontspec} %,xltxtra}

%\usepackage{unicode}

\usepackage[euler]{textgreek}

%%\DeclareUnicodeCharacter{3B8}{\ensuremath{\uptheta}}

%

\SIthousandsep{,}

%\usepackage{numprint}

\setlength\parindent{0pt}

\newtheorem{prop}{Proposition}

\renewcommand{\DifficultyMarker}{}

\newcommand{\AtBeginExerciseHeader}{\hspace{-21pt}} %-0.2pt

\renewcommand{\ExerciseHeader}{\AtBeginExerciseHeader\textbf{\ExerciseName~\ExerciseHeaderNB} \ExerciseTitle}

\renewcommand{\AnswerHeader}{\large\textbf{\AnswerName~\ExerciseHeaderNB}\smallskip\newline}

\setlength\AnswerSkipBefore{1em}

\usepackage{xspace}

\usepackage{imakeidx}

\makeindex

\usepackage[nottoc]{tocbibind}

\usepackage[colorlinks = true,

linktocpage=true,

pagebackref=true, % add back references to bibliography

linkcolor = red,

urlcolor = blue,

citecolor = red,

% refcolor =red,

anchorcolor = blue]{hyperref}

\definecolor{dkgreen}{rgb}{0,0.6,0}

\definecolor{gray}{rgb}{0.5,0.5,0.5}

\definecolor{gray2}{rgb}{0.35,0.35,0.35}

\definecolor{mauve}{rgb}{0.58,0,0.82}

\definecolor{index}{rgb}{0.88,0.32,0}

%------- source code settings

\usepackage{listings}

\lstset{frame=tb,

language=Python,

aboveskip=3mm,

belowskip=3mm,

showstringspaces=false,

columns=flexible,

basicstyle={\small\ttfamily},

numbers=none,

numberstyle=\tiny\color{gray},

keywordstyle=\color{blue},

commentstyle=\color{dkgreen},

stringstyle=\color{mauve},

breaklines=true,

breakatwhitespace=true,

tabsize=3

}

%-----------------------------------------------------------------

\usepackage{blindtext}

\usepackage{geometry}

\geometry{

a4paper,

total={170mm,257mm},

left=20mm,

top=20mm,

}

\setlength{\baselineskip}{0.0pt}

\setlength{\parskip}{3pt plus 2pt}

\setlength{\parindent}{20pt}

\setlength{\marginparsep}{0.0cm}

\setlength{\marginparwidth}{0.0cm}

\setlength{\marginparpush}{0.0cm}

\setlength{\tabcolsep}{4pt}

\renewcommand{\arraystretch}{1.4} %%%

\newtheorem{theorem}{Theorem}[section]

\newtheorem{lemma}[theorem]{Lemma}

\newtheorem{proposition}[theorem]{Proposition}

\newtheorem{corollary}[theorem]{Corollary}

\newenvironment{proof}[1][Proof]{\begin{trivlist}

\item[\hskip \labelsep {\bfseries #1}]}{\end{trivlist}}

\newenvironment{definition}[1][Definition]{\begin{trivlist}

\item[\hskip \labelsep {\bfseries #1}]}{\end{trivlist}}

\newenvironment{example}[1][Example]{\begin{trivlist}

\item[\hskip \labelsep {\bfseries #1}]}{\end{trivlist}}

\newenvironment{remark}[1][Remark]{\begin{trivlist}

\item[\hskip \labelsep {\bfseries #1}]}{\end{trivlist}}

\newcommand{\qed}{\nobreak \ifvmode \relax \else

\ifdim\lastskip<1.5em \hskip-\lastskip

\hskip1.5em plus0em minus0.5em \fi \nobreak

\vrule height0.75em width0.5em depth0.25em\fi}

\usepackage[symbols,nogroupskip,acronym]{glossaries-extra}

%\usepackage[xindy,symbols,nogroupskip,sort=def,acronym]{glossaries}

\makenoidxglossaries %%%%%%%%%%%%%%%

%\setlength{\glsdescwidth}{1.3\hsize}

\begin{document}

\hypersetup{linkcolor=blue}

%inserting a glossary entry in gloss: \gls{gls:keyword1} \\

\baselineskip=2\baselineskip

\thispagestyle{empty}

\hspace{0pt}

\vfill

%\hrulefill

\begin{center}

\rule{0.90\textwidth}{.4pt}

\end{center}

\begin{center}

{\Huge \bf{Gentle Introduction To Chaotic Dynamical Systems} }

\end{center}

\baselineskip=0.5\baselineskip

\addvspace{2cm}

\begin{center}

%\includegraphics[width=0.7\textwidth]{linear.png} \\

%\addvspace{1cm}

%\includegraphics[width=0.6\textwidth]{imgpyRiemannFinalOrbits-v2-small2.jpg}

\includegraphics[width=0.75\textwidth]{pillow.png}

\end{center}

\addvspace{2cm}

\begin{center}

\rule{0.90\textwidth}{.4pt}

\end{center}

\begin{center}

Vincent Granville, Ph.D. $|$ \href{https://mltechniques.com/}{www.MLTechniques.com} $|$ Version 2.0, May 2023

\end{center}

%\hrulefill

\hypersetup{linkcolor=red} % red %

\vfill

\hspace{0pt}

\pagebreak

\chapter*{Preface} %\clearpage

In less than 100 pages, the book covers all important topics about discrete chaotic dynamical systems and

related time series and stochastic processes, ranging from introductory to advanced, in one and two dimensions. State-of-the art methods and new results are

presented in simple English. Yet, some mathematical proofs appear for the first time in this book: for instance,

about the full autocorrelation function of the logistic map, the absence of cross-correlation between digit sequences

in a family of irrational numbers, and a very fast algorithm to compute the digits of quadratic irrationals. These are not just new important if not seminal theoretical developments: it leads to

better algorithms in random number generation (PRNG), benefiting applications such as

data synthetization, security, or heavy simulations. In particular, you will find an implementation of a very fast, simple PRNG based on

millions of digits of millions of quadratic irrationals, producing strongly random sequences superior in many respects to those available

on the market.

Without using measure theory, the invariant distributions of many systems are discussed in details, with numerous closed-form expressions for classic and new

maps, including the logistic, square root logistic, nested radicals, generalized continued fractions (the Gauss map), the ten-fold and dyadic maps, and more. The concept of bad seed, rarely discussed in the literature, is explored in details. It leads to singular fractal distributions with

no probability density function, and sets similar to the Cantor set. Rather than avoiding these monsters, you will be able to leverage them

as competitive tools for modeling purposes, since many evolutionary processes in economy, fintech, physics, population growth and so on, do not always behave nicely.

A summary table of numeration systems serves as a useful, quick reference on the subject. Equivalence between different maps is also discussed. In a nutshell, this book is dedicated to the study of two numbers: zero and one, with

a wealth of applications and results attached to them, as well as some of the toughest mathematical conjectures. It will appeal in particular to busy practitioners

in fintech, security, defense, operations research, engineering, computer science, machine learning, AI, as well as consultants

and professional mathematicians. For students complaining about how hard this topic is, and deterred by the amount of advanced mathematics, this book will help them get jump-started. While the mathematical level remains high in some sections, it is explained as

simply as possible, focusing on what is needed for the applications.

Numerous illustrations including beautiful representations of these systems (generative art), a lot of well documented Python code, and nearly 20 off-the-beaten-path exercises complementing the theory, will help you navigate through this fascinating field.

You will see how even the most basic systems offer such an incredible variety of configurations depending on a few parameters, allowing you

to model a very large array of phenomena. A surprising application -- a synthetic stock exchange and lottery -- is described in detail in chapter~\ref{pouti}, including full business model and legal aspects.

Finally, chapter~\ref{ch1}

also covers time-continuous processes including unusual clustered, reflective, constrained, and integrated Brownian-like processes, random walks and time series, with little math and jargon-free. In the end, my goal is to

get you to you use these systems fluently, and see them as gentle, controllable chaos. In short, what real life should be! Quantifying the amount

of chaos is also one of the topics discussed in the book.

\section*{About the author}

Vincent Granville is a pioneering data scientist and machine learning expert, co-founder of Data Science Central (acquired by TechTarget), founder of \href{https://mltechniques.com/}{MLTechniques.com}, former VC-funded executive, author and patent owner.

\begin{wrapfigure}{l}{0.14\textwidth}

\vspace{-1ex}

\includegraphics[width=0.98\linewidth]{vgr3.png}

%\caption{Caption1}

%\label{fig:wrapfig}

\end{wrapfigure}

\vspace{-2ex}\quad \\

\noindent Vincent’s past corporate experience includes Visa, Wells Fargo, eBay, NBC, Microsoft, and CNET.

Vincent is also a former post-doc at Cambridge University, and the National Institute of Statistical Sciences (NISS).

He published in {\em Journal of Number Theory}, {\em Journal of the Royal Statistical Society} (Series B), and {\em IEEE Transactions on Pattern Analysis and Machine Intelligence}. He is also the author of multiple books, available \href{https://mltechniques.com/resources/}{here}. He lives in Washington state, and enjoys doing research on stochastic processes, dynamical systems, experimental math and probabilistic number theory.

%\renewcommand{\baselinestretch}{0.97}\normalsize

\listoffigures

%\renewcommand{\baselinestretch}{1.00}\normalsize

%\listoftables

\hypersetup{linkcolor=red}

\renewcommand{\baselinestretch}{1.00}\normalsize

\tableofcontents

\renewcommand{\baselinestretch}{1.00}\normalsize

\chapter{Random Walks, Brownian Motions, and Related Stochastic Processes}\label{ch1}

I introduce these stochastic processes, routinely used by Wall Street quants, with a simple approach consisting of rescaling random walks to make them time-continuous, with a finite variance, based on the central limit theorem.

Stochastic processes have many applications, including in finance and physics. It is an interesting model to represent many phenomena.

One of the most simple examples is a random walk, and indeed easy to understand with no mathematical background. However, time-continuous stochastic processes are usually defined and studied using advanced and abstract mathematical tools such as measure theory, martingales, and filtration.

In this chapter, it is explained in simple terms, accessible to first-year college students with exposure to calculus and basic probability principles. At least initially, I discuss the one-dimensional case. After introducing standard Brownian motions, I dicuss other related processes

with interesting properties. It leads to a functional integral equation with exact solution, for the equilibrium distribution of the system: in this case, non-Gaussian unlike Brownian motions.

\section{From random walks to Brownian motions}\label{bmrwp}

The \textcolor{index}{Brownian motion}\index{Brownian motion} [\href{https://en.wikipedia.org/wiki/Brownian_motion}{Wiki}] is among the most basic time-continuous processes: it is simply an integrated \textcolor{index}{white noise}\index{white noise} [\href{https://en.wikipedia.org/wiki/White_noise}{Wiki}]. It is sometimes referred to as the \textcolor{index}{Wiener process}\index{Wiener process} [\href{https://en.wikipedia.org/wiki/Wiener_process}{Wiki}].

To build these processes, I use one of the most basic time-discrete stochastic processess: a \textcolor{index}{random walk}\index{random walk} [\href{https://en.wikipedia.org/wiki/Random_walk}{Wiki}]. The random walk is defined by $X_{k+1}$ equal to either $X_k + 1$ or $X_k -1$ with equal probabilities. Typically, one starts with $X_0=0$, and $k$ denotes the time.

Thus, a random walk is sequence of auto-correlated random variables indexed by time. Let $(U_k)$ with $k=1,2$ and so on be a sequence of independent binary random variables taking the values $+1$ or $-1$

with equal probability $\frac{1}{2}$. Then $X_n$ can be written as $X_n = U_1 + \dots + U_n$, and $-n\leq X_n\leq n$. Let's use $U_1$ as the representative for the family $(U_k)$ of random variables, since they all have the same distribution.

%xxx showing random walk, brownian motion

What happens if we change the time scale (horizontal axis) say from daily to hourly, or to every millisecond? We then also need to rescale the values (vertical axis) appropriately; otherwise the process exhibits massive oscillations in very short time periods. At the limit, if we consider infinitesimal time increments, the process becomes a continuous one. Much of the complex mathematics needed to define these continuous processes amount to performing the correct rescaling of the vertical axis, to make the limiting process meaningful.

You can define these time-continuous processes as the limit of their time-discrete version: using the correct rescaling is straightforward. Let us define $\{Y_t(n)\}$ as the same process as $\{X_k\}$, but with small time increments of $1/n$ instead of $1$. In other words, $Y_{k/n}(n) = X_k$. I just rescaled the horizontal time axis. I also introduced the notation $t=k/n$ for the more granular time. Note that $Y_t(n)$ can take on very large values, between $-n$ and $+n$ when $t = 1$. Thus we also need to rescale the vertical axis. We have:

$$

\text{Var}[Y_1(n)] =\text{Var}[X_n] = \sum_{t=1}^n \text{Var}[U_t] = n\cdot \text{Var}[U_1].

$$

The only way to make the right-hand side of the equation not depending on $n$ is to rescale the vertical axis as follows. Define

$$

Z_t(n) = \frac{Y_t(n)}{\sqrt{n}}, \quad Z_t = \lim_{n\rightarrow\infty} Z_t(n).

$$

Then

$$

\text{Var}[Z_1(n)] = \text{Var}[U_1], \quad \text{Var}[Z_t(n)] = t\cdot \text{Var}[U_1],

$$

which leads to

\begin{equation}

\text{Var}[Z_t] = t \cdot \text{Var}[U_1].\label{eq1}

\end{equation}

Also, because of the central limit theorem, by construction for any real value $t$, $Z_t$ has a Gaussian distribution, regardless of the distribution of $U_1$. The final process $Z_t$ is both time-continuous and continuous on the vertical axis, though nowhere differentiable. It looks like a fractal and it is known as a Brownian motion -- the standard time-continuous stochastic process -- from which many other processes are derived.

Note that if instead of using a binary random variable for $U_1$, you use a Gaussian one, then the limiting process

$Z_t$ is identical, but we are dealing with Gaussian variables throughout the construction, making it easier to study the covariance structure and other properties. It then becomes a simple exercise to derive the covariance between $Z_t$ and $Z_s$ for any $t,s$. The covariance can also be estimated using simulations. Finally, note that $Z_0 = 0$ and $\text{E}[U_1] = 0$.

%%

\section{General Properties}

The Brownian motion can also be viewed as a Gaussian stationary time series, characterized by its covariance or auto-correlation structure. It is also related to deterministic dynamical systems that exhibit a fractal behavior. Under appropriate transformations, many of these processes can be made equivalent.

One question is whether the above construction (the limit of a time-discrete random walk) covers all types of Brownian motions, or only a few particular cases. One way to investigate this is to check whether this construction can generate any kind of covariance structure that characterizes these processes. The answer is positive, making advanced mathematical theory unnecessary to build and study Brownian motions, as well as the numerous complex stochastic processes derived from this base process. However, if you allow the random variables $U_k$ used in our construction to {\em not} be independent, then you can build more sophisticated time-continuous stochastic processes, that are not Brownian motions.

All the stochastic processes introduced so far, whether time-discrete or time-continuous, share the following properties. In most cases, it is easy to turn a stochastic process into one that satisfies these properties, using simple transformations, as illustrated later in this section.\vspace{1ex}

\begin{itemize}

\item \textcolor{index}{Stationarity}\index{stationarity} [\href{https://en.wikipedia.org/wiki/Stationary_process}{Wiki}]: It means that there is no trend, drift, or more precisely, the fact that the properties of the process in question do not explicitly depend on the time parameter $t$. Such processes are usually characterized by their auto-correlation structure alone.

\item \textcolor{index}{Ergodicity}\index{ergodicity} [\href{https://en.wikipedia.org/wiki/Ergodicity}{Wiki}]: This means that one instance of the process is enough to derive all its properties. You don't need to make hundreds of simulations to study the process' properties or compute estimates: one simulation (also called {\em instance}) over a very long time period will do.

\item \textcolor{index}{Scale-invariant}\index{scale-invariant} [\href{https://en.wikipedia.org/wiki/Scale_invariance}{Wiki}] and fractal behavior: If you zoom in or out on any single realization of these processes, you will get a new process with the exact same properties and behavior, indistinguishable from the parent process. Two different time windows provide two versions of the process that are identical with respect to their statistical properties.

\item \textcolor{index}{Memoryless}\index{memoryless property} [\href{https://en.wikipedia.org/wiki/Memorylessness}{Wiki}]: The future observations depend on the present value only, not on past observations. This is sometimes referred to as

the \textcolor{index}{Markov property}\index{Markov property} [\href{https://en.wikipedia.org/wiki/Markov_property}{Wiki}].

\end{itemize}\vspace{1ex}

\noindent It is sometimes possible to transform a process so that it satisfies some of the above properties. For instance, if $X_k$ is a time series with a linear trend and discrete time increments, the differences $X_k - X_{k-1}$ may represent a stationary time series. Likewise, if

$X_k$ depends on $X_{k-1}$ and $X_{k-2}$, then the vector $(X_k, X'_k)$ with $X'_k =

X_{k-1}$ represents a bivariate memoryless time series: the next bivariate value depends only on the current one.

Another interesting topic with numerous applications is the distribution of records (minimum or maximum), and their arrival times. A summary of the most important results is found in my book on synthetic data~\cite{vgsynthetic}, in the section

entitled ``Record Distances Between a Point and its Vertex". It deals with 2D processes, but the results also apply to 1D random walks. The analysis of records is the main subject of \textcolor{index}{extreme value theory}\index{extreme value theory} [\href{https://en.wikipedia.org/wiki/Extreme_value_theory}{Wiki}].

Finally, two important results are the \textcolor{index}{arcsine law}\index{arcsine law} [\href{https://en.wikipedia.org/wiki/Arcsine_laws_(Wiener_process)}{Wiki}] and chance to return to the origin. The former states that for 1D random walks and Brownian motions, the proportion of the time that the process is positive follows an arcsine law. The latter states the following: The probability of a random walk returning to its origin is 1 in one or two dimensions but only 34\% in three dimensions:

this is \textcolor{index}{Pólya's theorem}\index{Pólya's theorem} [\href{https://mathworld.wolfram.com/PolyasRandomWalkConstants.html}{Wiki}].

\begin{Exercise} {\em Statistical estimations} -- Simulate $10^6$ realizations of a Brownian motion with $0\leq t \leq 1$, using the random walk construction described previously. Study the distribution of the following quantities, using computations averaged across all your simulations. In particular, what is the mean and variance, for the following quantities:

\begin{itemize}

\item Extreme values: $\min Z_t$, $\max Z_t$,

\item Proportion of the time when $Z_t > 0$ (note that $Z_0 = 0$),

\item Number of times when the sign of $Z_t$ changes.

\end{itemize}

Keep in mind that the $Z_t$'s are auto-correlated. Given a particular realization of a stochastic process, these statistics can be used to check if it is a Brownian motion or not. Another interesting exercise is to study the process in question if the variable $U_1$ does not have a variance, for instance if $U_1$ has a Cauchy distribution.

\end{Exercise}

\section{Integration, differentiation, moving averages}\label{movbc}

Let's use the construction scheme in section~\ref{bmrwp} to build a Brownian motion $\{Z_t\}$.

The underlying time-discrete random walk $\{X_k\}$ is referred to as the base process.

I introduce two transformations:\vspace{1ex}

\begin{itemize}

\item The cumulative or \textcolor{index}{integrated process}\index{integrated process} $\{ S_t \}$ derived from $\{ Z_t \}$,

\item The theoretical \textcolor{index}{moving average process}\index{moving average process} $\{ M_t \}$ derived from $\{ Z_t \}$.

\end{itemize}\vspace{1ex}

The inverse of integration is differentiation: the \textcolor{index}{differentiated}\index{differentiated process} $\{S_t\}$ is $\{Z_t\}$. In practice, the smoother process (integrated or moving average) is easier to study and sometimes displays patterns that can't be identified in the original process.

To build the processes in question, proceeed as follows.

In the construction of the Brownian motion described in section~\ref{bmrwp}, replace

$X_n = U_1 + \cdots + U_n$ by

$X_n = V_1 + \cdots + V_n$, where $V_k$ is described below for each transformation. \vspace{1ex}

\begin{itemize}

\item Integration: $V_k = U_1 + \cdots + U_k$.

\item Differention: $V_k = U_{k+1} - U_k$. If $\{ Z_t \}$ is a Brownian motion then the resulting process is a white noise: nowhere continuous, nowhere differentiable.

\item Moving average: $V_k = U_k + U_{k + 1} + \cdots + U_{k + h(k)}$ where $h(k)$ is as small as possible to make the resulting process continuous and differentiable everywhere.

\end{itemize}\vspace{1ex}

For moving averages applied to a Brownian motion, $h(k) = \lfloor \sqrt{k}\rfloor$ works. Here $\lfloor \cdot \rfloor$ stands for the integer part function. Does $h(k) = \lfloor \log k \rfloor$ work? This would make the resulting process far more similar to the original one, but maybe barely (if at all) continuous -- in other words, more chaotic than with $h(k) = \lfloor\sqrt{k} \rfloor$.

\begin{figure}[H]

\centering

\includegraphics[width=0.78\textwidth]{ma2.PNG} %0.86

\caption{Brownian motion (green), integrated (orange) and moving average (red)}

\label{fig:trc}

\end{figure}

As in section~\ref{bmrwp}, you need to use the correct rescaling of the vertical axis to obtain meaningful variances that do not depend on $n$.

I illustrate this for the integrated process. In this case we have:

$$

X_n = \sum_{k=1}^n V_k =\sum_{k=1}^n k\cdot U_{n-k+1},

$$

thus

$$

\text{Var}[X_n] = \sum_{k=1}^n k^2 \text{Var}[U_1] =\frac{n(n+1)(2n+1)}{6}\text{Var}[U_1]\sim \frac{n^3}{3}\text{Var}[U_1].

$$

The rescaling of the horizontal time axis is the same. Thus, the proper rescaling factor for the vertical axis, as $n$ tends to infinity, is

$X_n / \sqrt{n^3/3}$. This leads to

\begin{equation}

\text{Var}[S_t] = t^3 \cdot \text{Var}[U_1].\label{eq2}

\end{equation}

The same logic applies to compute $\text{Var}[M(t)]$. The details are left as an exercise. A more complicated exercise consists of computing the covariance between $S_t$ and $S_{t + s}$ for $s > 0$, and proving that $\{S_t\}$ is {\em not} a Brownian motion itself (being differentiable everywhere unlike Brownian motions.)

Figure~\ref{fig:trc} was produced with the Python code in this section. It represents one realization of a Brownian motion $\{Z_t\}$, together with its integration $\{S_t\}$ and moving average $\{M_t\}$. The time period is $0\leq t \leq 5$. The Python program \texttt{brownian.py} is also on my

GitHub directory, \href{https://github.com/VincentGranville/Stochastic-Processes/blob/master/Brownian.py}{here}.

The Brownian motion, integrated Brownian and moving average are denoted respectively as \texttt{X}, \texttt{S} and \texttt{M} in the code,

while \texttt{T} is the time.

Different metrics exist to characterize the smoothness stochastic process is. The

\textcolor{index}{Hurst exponent}\index{Hurst exponent} [\href{https://en.wikipedia.org/wiki/Hurst_exponent}{Wiki}] is a well known one, measuring the amount of long-term memory in your process. It is equal to $\frac{1}{2}$ for a Brownian motion or

\textcolor{index}{Brown noise}\index{Brown noise} [\href{https://en.wikipedia.org/wiki/Brownian_noise}{Wiki}], and to 0

for a \textcolor{index}{pink noise}\index{pink noise} [\href{https://en.wikipedia.org/wiki/Pink_noise}{Wiki}] and for a

one-dimensional \textcolor{index}{white noise}\index{white noise} [\href{https://en.wikipedia.org/wiki/White_noise}{Wiki}]. The higher the value, the smoother the process. When the parameter \texttt{smooth} in the Python code is increased,

the resulting simulation is smoother. Another methodology to generate integrated (smooth) processes is discussed in my book on synthetic data~\cite{vgsynthetic}, in the chapter dealing with linear algebra and auto-regressive time series.

%-----------------------------vince/riemann2and3.mp4

\begin{figure}%[H]

\centering

\includegraphics[width=0.9\textwidth]{linear.png}

\caption{Integrated Brownian (top left), Brownian (top right) and nowhere continuous (bottom)}

\label{fig:lollog1xx}

\end{figure}

%imgpy9979_2and3.PNG

%-------------------------

There are many different ways to simulate Brownian motions and related processes. Figure~\ref{fig:lollog1xx} is based on simple

\textcolor{index}{auto-regressive time series}\index{autoregressive models} [\href{https://en.wikipedia.org/wiki/Autoregressive_model}{Wiki}], properly scaled to make them time-continuous. The technique is described in the chapter ``Gentle Introduction to Linear Algebra -- Synthetic Time Series" in my book on synthetic data~\cite{vgsynthetic}. The originality here consists in choosing auto-regressive models where some of the roots of the \textcolor{index}{characteristic polynomial}\index{characteristic polynomial}

[\href{https://en.wikipedia.org/wiki/Characteristic_polynomial}{Wiki}] are multiple, in particular roots with

modulus equal to 1 (the ``largest" roots). The processes in the top part of Figure~\ref{fig:lollog1xx} are standard: roots with multiplicity 1 yield Brownian motions, multiplicity 2 yields integrated Brownian, multiplicity 3 yields doubly integrated Brownian and so on. In the bottom part, the underlying time series used in the simulation has $X_n$ depends on $X_{n-2}$, $X_{n-3}$ and so on, but not

on $X_{n-1}$. The resulting processes are non-Brownian. Indeed they look like a ``derivative" of a Brownian: a non-standard function that densely fills some peculiar domain, resulting in a shape with a \textcolor{index}{fractal dimension}\index{fractal dimension}

[\href{https://en.wikipedia.org/wiki/Fractal_dimension}{Wiki}] between 1 and 2: something intermediate between a curve and a surface.

The following Python script \texttt{Brownian.py} is used to produce Figure~\ref{fig:trc}. It is on

GitHub, \href{https://github.com/VincentGranville/Stochastic-Processes/blob/master/Brownian.py}{here}. \vspace{1ex}

%https://www.youtube.com/watch?v=W9jktqV3_Mc brownian wolfram alpha look for api

\begin{lstlisting}

import numpy as np

import matplotlib.pyplot as plt

n = 10000

m = 5*n

T = []

X = []

T.append(0.0)

X.append(0.0)

np.random.seed(1979)

for k in range(1,m):

u = np.random.uniform(0,1)

if u < 0.5:

X.append(X[k-1]-1)

else:

X.append(X[k-1]+1)

T.append(T[k-1] + 1/n)

S = []

S.append(0.0)

for k in range(1,m):

S.append(X[k]+S[k-1])

M = []

smooth = 2.5 # the larger, the smoother the moving average

M.append(0.0)

hn = int(smooth*np.sqrt(n))

for k in range(1,m):

sum = 0.0

for h in np.arange(-hn, hn+1):

idx = k + h

if idx >= m: # fix for index outside the array range

idx = m - 1 - (idx % n)

elif idx < 0: # fix for index outside the array range

idx = -idx

sum += X[idx]

sum /= (2*hn + 1)

M.append(sum)

for k in range(1,m):

X[k] = X[k]/(n**0.5)

S[k] = S[k]/(n**1.5)

M[k] = M[k]/(n**0.5)

axes = plt.axes()

axes.tick_params(axis='both', which='major', labelsize=8)

axes.tick_params(axis='both', which='minor', labelsize=8)

for axis in ['top','bottom','left','right']:

axes.spines[axis].set_linewidth(0.5)

plt.plot(T, X, linewidth = 0.4, color = 'green', alpha = 0.2) # Brownian motion

plt.plot(T, S, linewidth = 0.8, color = 'orange', alpha = 0.8) # integrated Brownian motion

plt.plot(T, M, linewidth = 0.8, color = 'red', alpha = 1.0) # moving average process

plt.axhline(y = 0.0, color = 'grey', linestyle = '--', linewidth = 0.4)

plt.show()

\end{lstlisting}

\section{Reflected random walks}\label{rflectr}

The goal here is to introduce the reader to a bounded process, where moves up or down are allowed or not based on the current state. It mimics an environment where constraints prevent the system from going too far up too quickly, or the other way around. Thus such processes do not exhibit massive explosions or implosions.

\begin{figure}%[H]

\centering

\includegraphics[width=0.85\textwidth]{ma4.PNG} %0.86

\caption{Reflected random walk with $a=b=\frac{1}{2}$}

\label{fig:ivf}

\end{figure}

Our \textcolor{index}{reflected random walk}\index{reflected random walk} [\href{https://en.wikipedia.org/wiki/Reflected_Brownian_motion}{Wiki}] is defined as follow. Start with $X_0=0$ and then

\begin{equation}

X_k = \begin{cases}X_{k-1} + U_k\cdot k^{-a} & \text{if } X_{k-1} < 0, \\

X_{k-1} - U_k\cdot k^{-a} & \text{if } X_{k-1} \geq 0,%

\end{cases}\label{poputres}

\end{equation}

%\begin{align}

%X(k) & = X(k-1) + \frac{U(k)}{k^a} \quad \text{ if } X(k-1) < 0, \nonumber \\

% X(k) & = X(k-1) - \frac{U(k)}{k^a} \quad \text{ if } X(k-1) \geq 0, \nonumber

%\end{align}

where $U_k = V^b_k$, the $V_k$'s are uniform independent deviates on $[0, 1]$ and $a\geq 0, b>0$ are two parameters.

Also define $Z_k=k^a X_k$.

The distribution of $Z_k = k^a X_k$ is stable over time, in contrast to Brownian motions that lack stationarity unless normalized. For the reflected random walk, the limiting distribution of $Z_k$ -- called the \textcolor{index}{invariant measure}\index{invariant measure} [\href{https://en.wikipedia.org/wiki/Invariant_measure}{Wiki}]

or \textcolor{index}{attractor distribution}\index{attractor distribution} in dynamical systems -- is not Gaussian. This is another indication that the process is not Brownian.

\begin{figure}[H]

\centering

\includegraphics[width=0.80\textwidth]{ma3.PNG} %0.86

\caption{Invariant measure (density function) of reflected random walk with $a=b=\frac{1}{2}$}

\label{fig:iv}

\end{figure}

However it satisfies

the \textcolor{index}{stochastic integral equation}\index{stochastic integral equation}~(\ref{seqa}), which has an exact solution. For this type of

\textcolor{index}{functional equations}\index{functional equation} [\href{https://en.wikipedia.org/wiki/Functional_equation}{Wiki}], the unknown is a probability distribution: here, the attactor distribution. The solution is also referred to as a

\textcolor{index}{fixed point}\index{fixed point algorithm} [\href{https://en.wikipedia.org/wiki/Fixed_point_(mathematics)}{Wiki}], in this case of infinite dimension. It may not be unique.

Let $Z$ be the limit of $Z_k$ as $k\rightarrow\infty$. What is the

distribution of $Z$? Let $F_Z$ and $f_Z$ denote respectively the CDF (cumulative distribution function) and density function attached to $Z$.

Likewise, $F_U$ and $f_U$ is the common CDF and density attached to the $U_k$'s. To find $f_Z$ you need to solve the following integral equation where $f_Z$ is the unknown:

\begin{equation}

f_Z(z) = \int_0^1 \Big[f_U(x+z) + f_U(x-z)\Big]\cdot f_Z(x) dz. \label{seqa}

\end{equation}

Figure~\ref{fig:iv} shows an approximation to the solution computed on $2\times 10^6$ iterates of Formula~(\ref{poputres}) with

$a=b=\frac{1}{2}$. In this case, the exact density function is known and equal to

$$

f_Z(z) = \frac{b+1}{2}\cdot \Big(1-|z|^{1/b}\Big) = \frac{1-F_U(z)}{2 \text{E}[U]}, \quad -1\leq z \leq 1.

$$

This is the correct solution if $b=1,b=\frac{1}{2}$, or $b\rightarrow 0$, regardless of $0<a<1$. It can be verified by plugging this solution in Formula~(\ref{seqa}). I haven't checked if the formula is still valid for other values of $b$. It is easy to empirically obtain the following result,

based on observations over a very large time period:

$$

\text{Var}[Z] = \frac{b+1}{3(3b+1)}.

$$

Of course $\text{E}[Z]=0$. Figure~\ref{fig:ivf} shows one realization of a reflected random walk (the first $2000$ values of $Z_k$) when $a=b=\frac{1}{2}$. The Python code to produce Figures~\ref{fig:ivf} and ~\ref{fig:iv} is on my GitHub repository

\href{https://github.com/VincentGranville/Stochastic-Processes/blob/master/brownian_reflective.py}{here}, and also in

section~\ref{dupuis}.

\subsection{Exercises}

The purpose of this section is to explore particular cases of reflective random walks, corresponding to specific values of $a$ and $b$.

In Exercise~\ref{knorr}, the formula for the integral equation is established.

\begin{Exercise} {\em Reflective random walks: special cases} -- Prove that if $0<a < 1$, then $X_k\rightarrow 0$ as

$k\rightarrow\infty$. Under the same condition, prove that the limiting distribution of $Z$

\begin{itemize}

\item always exists and its support domain is $[-1, 1]$,

\item is symmetric, with mean and median equal to 0,

\item does not depend on $a$, but only on $b$.

\end{itemize}

For instance, if $b =1$, the distribution of $Z$ is triangular regardless of $a$. If $a < 1$ and $b = 0$, (the non-stochastic case) prove that

$Z$ can only take on 3 values: $-1,+1$ and $0$ respectively with probability $\frac{1}{4},\frac{1}{4}$ and $\frac{1}{2}$.

Also show that when $b\rightarrow 0^{+}$, the distribution of Z converges to a uniform distribution on $[-1, 1]$.

When $b=0$ (the non-stochastic case), we also have the following cases:

\begin{itemize}

\item if $a=1$, $X_k\rightarrow 0$,

\item if $a = 3$, $X_k\rightarrow \zeta(3) -\frac{5}{4}\approx -0.048$,

\item if $a = 4$, $X_k\rightarrow \zeta(4) -\frac{9}{8}\approx -0.043$.

\end{itemize}

Here $\zeta(\cdot)$ is the \textcolor{index}{Riemann zeta function}\index{Riemann zeta function} [\href{https://en.wikipedia.org/wiki/Riemann_zeta_function}{Wiki}].

\end{Exercise}

\begin{Exercise}\label{knorr} {\em Establishing the integral equation} -- Prove that the density $f_Z(z)$ attached to the limiting distribution must

satisfy~(\ref{seqa}). \vspace{1ex} \\

{\bf Solution} \\

Based on Formula~(\ref{poputres}), we have

$$

P(X_k<z) = \int_{-1}^1 P(X_k < z | X_{k-1} = x) f_{X_{k-1}}(x) dx.

$$

Let $U$ be a random variable with the same distribution as any $U_k$, and let $k\rightarrow\infty$. Separating the cases $x<0$ and $x\geq 0$, we have:

$$

P(Z<z) = \int_0^1 P(U > x-z )f_Z(x)dx + \int_{-1}^0 P(U < z-x )f_Z(x)dx.

$$

Taking advantage of the symmetries of the problem, this can be further simplified to

$$

F_Z(z) = \frac{1}{2} + \int_0^1 \Big[F_U(x+z)-F_U(x-z)\Big] f_Z(x)dx.

$$

Finally, taking the derivative with respect to $z$ on both sides of the previous equality, we obtain the desired integral equation for $f_z$. \qed

We will encounter more of these integral equations with an exact solution in the next chapters. Many are related to dynamical systems that represent new numeration systems, leading to interesting invariant measures. But in general, these equations do not have an exact solution. They are typically solved numerically with an iterative algorithm converging to the solution.

\end{Exercise}

\subsection{Python code}\label{dupuis}

The Python code in this section is used to produce Figures~\ref{fig:ivf} and~\ref{fig:iv}. It is also on my GitHub repository

\href{https://github.com/VincentGranville/Stochastic-Processes/blob/master/brownian_reflective.py}{here}. I included it in this chapter because it illustrates the Seaborn and Statsmodels libraries: the former to produce empirical density plots, and the latter to compute empirical

distributions. \vspace{1ex}

\begin{lstlisting}

import numpy as np

import matplotlib.pyplot as plt

import seaborn as sns

from statsmodels.distributions.empirical_distribution import ECDF

m = 2000

a = 0.5

b = 0.5

T = []

X = []

Z = []

T.append(0.0)

X.append(0.0)

Z.append(0.0)

np.random.seed(1979)

for k in range(1,m):

u = np.random.uniform(0,1)

u = u**b

if X[k-1] < 0:

X.append(X[k-1] + u/(k**a))

else:

X.append(X[k-1] - u/(k**a))

Z.append((k**a) * X[k])

T.append(T[k-1] + 1)

axes = plt.axes()

axes.tick_params(axis='both', which='major', labelsize=8)

axes.tick_params(axis='both', which='minor', labelsize=8)

for axis in ['top','bottom','left','right']:

axes.spines[axis].set_linewidth(0.5)

plt.plot(T, Z, linewidth = 0.4, color = 'green', alpha = 1)

plt.axhline(y = 0.0, color = 'black', linestyle = '--', linewidth = 0.4)

plt.show()

ecdf = ECDF(Z[10:len(Z)])

sns.set_context("paper",font_scale=0.8, rc={"lines.linewidth": 0.8})

sns.displot(Z[10:len(Z)], kind="kde",linewidth=0.5) ### , bins=150)

plt.show()

\end{lstlisting}

\section{Constrained random walks}

%xxxx

%properties of RW / records distribution of records

%xxx random walk never comes back to the origin in 3D ???

%x to add to previous section

%arc sine law

%xxxx link to [or add entire chapter] to "subtle departures to randomness" where I study

% return to zero and so on. or link to extreme value theory in same book

%xxx add chapter on long term autocorrel time series

%xxx pillow article / article with solution to dynamical system

%xxx 2D numeration system

%xxxx appendic=x from stats book

% add chapter: PRNG based on quadratic irrationals

The standard symmetric 1D \textcolor{index}{random walk}\index{random walk} [\href{https://en.wikipedia.org/wiki/Random_walk}{Wiki}] fundamental to this section is a sequence $\{S_n\}$ with $n\geq 0$, starting at $S_0=0$, and recursively defined by

$S_{n}=X_{n}+S_{n-1}$, for $n>0$. Here $X_1,X_2$ and so on are independent random variables with $P[X_n=1]=P[X_n=-1]=\frac{1}{2}$.

Thus $\{S_n\}$ is a time-discrete stochastic process, and indeed the most basic one. In sections~\ref{azxa} and~\ref{azxb}, I drop the assumption

of independence, leading to modified random walks such as those described in~\cite{nkrn2018,lanwu2012}.

More general references include~\cite{gtm2021,peresbrown}.

With proper rescaling, a random walk becomes a time-continuous stochastic process $S_t$ called \textcolor{index}{Brownian motion}\index{Brownian motion} [\href{https://en.wikipedia.org/wiki/Brownian_motion}{Wiki}], with $t\in\mathbb{R}^+$. See the time series in gray in Figure~\ref{fig:walk}, displaying a particular instance: it shows the first $\num{50000}$ values of

$S_n$ in a short window, giving the appearance of a Brownian motion.

By contrast, each of the orange, red and gray time series represents one instance of a specific type of non-Brownian motion. Sections~\ref{azxa} and~\ref{azxb} focuses on these three types of processes, which are quasi, but not fully random.

\subsection{Three fundamental properties of pure random walks}\label{poyt}

The standard random walk $\{S_n\}$ (illustrated in gray in Figure~\ref{fig:walk}) is the base or reference process, used to build more sophisticated models. It has too many properties to list in this short chapter. However, the following are the most relevant to our discussion.

\begin{itemize}

\item \textcolor{index}{Law of the iterated logarithm}\index{law of the iterated logarithm}\index{iterated logarithm} [\href{https://en.wikipedia.org/wiki/Law_of_the_iterated_logarithm}{Wiki}].

In our context, it is stated as follows:

\begin{equation}

\lim \sup \frac{|S_n|}{\sqrt{2nv\log \log n}} = 1 \quad \text{as } n\rightarrow \infty.\label{lil12}

\end{equation}

Here, as per the

\textcolor{index}{Hartman–Wintner theorem}\index{Hartman–Wintner theorem} [\href{https://encyclopediaofmath.org/wiki/Law_of_the_iterated_logarithm}{Wiki}], $v=\text{Var}[X_1]=1$. See~\cite{peresbrown} pages 118--123 for a version adapted to Brownian motions.

\item Expected number of \textcolor{index}{zero crossings}\index{random walk!zero crossing} in $S_1,\dots,S_n$, denoted as $N_n$. Here a zero-crossing is an index $0<k\leq n$

such that $S_k=0$. For $n>0$, we have (see \href{https://math.stackexchange.com/questions/1684576/expected-of-returns-in-a-symmetric-simple-random-walk}{here}):

$$

\text{E}[N_{2n}]=-1+\frac{2n+1}{4^n} \binom{2n}{n} \sim \frac{2}{\sqrt{\pi}}\cdot \sqrt{n} \quad \text{as } n\rightarrow \infty.

$$

\item Distribution of \textcolor{index}{first hitting time}\index{random walk!first hitting time} to zero [\href{https://en.wikipedia.org/wiki/First-hitting-time_model}{Wiki}], or

first zero crossing after $S_0=0$, also called time of first return. The random variable in question is denoted as $T$. It is defined as follows:

$T=n$ (with $n>0$) if and only if $S_{n}=0$ and $S_k\neq 0$ if $0<k<n$. We have $P[T=n]=0$ if $n$ is odd, and $\text{E}[T]=\text{Var}[T]=\infty$. Yet, our random walks cross the X-axis infinitely many times. We also have the following

\textcolor{index}{probability generating function}\index{probability generating function} [\href{https://en.wikipedia.org/wiki/Probability-generating_function}{Wiki}] (see \href{https://math.stackexchange.com/questions/64919/biased-random-walk-and-pdf-of-time-of-first-return}{here}):

$$

\sum_{n=1}^\infty (2x)^{2n} P[T=2n] =1-\sqrt{1-4x^2} \quad \text{if } x\leq \frac{1}{4}.

$$

From there, one can obtain

\begin{align}

P[T=2n] & =\frac{1}{(2n-1)4^n}\binom{2n}{n}\sim \frac{1}{\sqrt{4\pi}}\cdot n^{-3/2} \quad \text{as } n\rightarrow \infty,\nonumber \\

\text{E}[T^{-1}] & = \int_{0}^{1/2} \frac{1-\sqrt{1-4x^2}}{x}dx = 1-\log 2. \nonumber

\end{align}

\end{itemize}

\noindent Note that $\text{E}[T^{-1}]$ is finite, while $\text{E}[T]$ is infinite. The fact that

$\text{E}[T]=\infty$ explains why the sequence $S_n$ can stay above or below the X-axis for incredibly long time periods, as shown

in Figure~\ref{fig:walk} for the gray curve.

%-----------------------------vince/riemann2and3.mp4

\begin{figure}%[H]

\centering

\includegraphics[width=0.8\textwidth]{walk.png}

\caption{Typical path $S_n$ with $0\leq n\leq \num{50000}$ for four types of random walks}

\label{fig:walk}

\end{figure}

%imgpy9979_2and3.PNG

%-------------------------

The above three statistics $|S_n|/\sqrt{2n\log\log n}, N_{2n}$ and $T^{-1}$can be used to design tests of randomness for

\textcolor{index}{pseudorandom number generators}\index{pseudo-random numbers}. Indeed, the \textcolor{index}{prime test}\index{pseudo-random numbers!prime test}\index{prime test (of randomness)} relies on a number theoretic version of the law of the iterated logarithm (LIL). See chapter ``High Quality Random Numbers for Simulations and Data Synthetization" in~\cite{vgsynthetic}.

The purpose is to detect very weak departures from randomness, even in sequences that are random enough to pass the classic LIL test, yet not fully random.

I now describe special types of modified random walks that lack true independence in the sequence $\{X_n\}$. In particular, I discuss why

they are special and of great interest, with a focus on applications.

\subsection{Random walks with more entropy than pure random signal}\label{azxa}

One way to introduce dependencies in the sequences is to increase the frequency of oscillations (and thus the entropy) in the gray curve in

Figure~\ref{fig:walk}. The gray curve represents a realization of a pure random walk.

To achieve this goal, you may want the sequence to violate the law of the iterated logarithm: you want to build a sequence that would satisfy a modified law of the iterated logarithm with $\sqrt{2n\log\log n}$ in Formula~(\ref{lil12}) replaced by (say) $n^{2/5}$.

To accomplish this, you need to add constraints when simulating the sequence in question. Yet you want to preserve quasi randomness:

the absence of drifts and auto-correlations in

the sequence $\{X_n\}$, even though there is some modest lack of independence. So modest indeed that most statistical tests would fail to catch it, even though it can made highly visible to the naked eye: see the red, and especially the blue curve in Figure~\ref{fig:walk}.

\subsubsection{Applications}

Such sequences can be used to generate \textcolor{index}{synthetic data} or to model barely constrained stochastic processes, such as stock price fluctuations in an almost perfect market. See also~\cite{pac203}. Another application is

to introduce an undetectable backdoor in some encryption systems without third parties (government or hackers) being able to notice it, depending on the strength of the dependencies. This type of backdoor can help the encryption company decrypt a message when requested by a

legitimate user who lost his key, even though the encryption company has no standard mechanism to store or retrieve keys (precisely to avoid government interference).

This assumes that there is a mapping between the type of weak dependencies introduced in a specific sequence, and the

type of algorithm (or the key) used to decrypt the sequence in question. The mapping can be made too loose for full decryption even by the parent company, but helpful to retrieve partial data, such as where the sequence originates from: in this case, the type of dependencies is a proxy for a signature. All that is needed is to add some extra bits so that the sequence has the desired statistical behavior.

Ironically, you need a very good, industrial-grade \textcolor{index}{pseudo-random number generator}\index{pseudo-random numbers} (PRNG) to generate almost perfectly random sequences. PRNG's that are not good enough -- such as the \textcolor{index}{Mersenne twister}\index{Mersenne twister}\index{pseudo-random numbers!Mersenne twister} -- may introduce irregularities that can interfere with the ones you want to introduce. This is discussed

in detail in chapter 10 about PRNG's, in~\cite{vgsynthetic}.

\subsubsection{Algorithm to generate quasi-random sequences}\label{qrrnd}

One way to generate such sequences is as follows: \vspace{1ex}

\noindent\textcolor{white}{00} $S=0$ \\

\textcolor{white}{00} {\bf For} $n = 1, 2,\dots$\\

\textcolor{white}{0000} Generate random deviate $U$ on $[0,1]$\\

\textcolor{white}{0000} $M=g(n)$\\

\textcolor{white}{0000} {\bf If} ($S< -M$ and $U < \frac{1}{2}-\epsilon$)

or ($S> M$ and $U < \frac{1}{2}+\epsilon$) or ($|S|\leq M$ and $U<\frac{1}{2}$) \\

\textcolor{white}{0000} {\bf Then} \\

\textcolor{white}{000000} $X_n=-1$ \\

\textcolor{white}{0000} {\bf Else} \\

\textcolor{white}{000000} $X_n=1$\\

\textcolor{white}{0000} $S=S+X_n$\\

\textcolor{white}{0000} $S_n=S$ \vspace{1ex}

\noindent Here $0<\epsilon<\frac{1}{2}$ and $\alpha>0$. The function $g(n)$ is positive and

growing more slowly than $\sqrt{n}$. Typically, $g(n)=\alpha n^\beta$ with $0\leq \beta\leq\frac{1}{2}$, or $g(n)=\alpha (\log n)^\beta$

with $\beta\geq 0$.

The Python code in section~\ref{paths} performs this simulation: choose the option \texttt{deviations='Small'}. You can customize the function

$g(n)$, denoted as \texttt{G} in the code. The option \texttt{mode='Power'} corresponds to $g(n)=\alpha n^\beta$,

while \texttt{mode='Log'} corresponds to $g(n)=\alpha (\log n)^\beta$.

%-----------------------------vince/riemann2and3.mp4

\begin{figure}%[H]

\centering

\includegraphics[width=0.9\textwidth]{iteratedlog1b.png}

\caption{$\delta_n=1-\text{Var}[S_{n+1}]+\text{Var}[S_n]$ for four types of random walks, with $0\leq n\leq\num{5000}$}

\label{fig:lollog1b}

\end{figure}

%imgpy9979_2and3.PNG

%-------------------------

\noindent Results are displayed in

Figure~\ref{fig:walk}. The color scheme is as follows:

\begin{itemize}

\item Gray curve: $\epsilon=0$, corresponding to a pure random walk.

\item Blue curve: $g(n)=\log n$, $\epsilon=0.05$.

\item Red curve: $g(n)=n^\beta$ with $\beta=0.35$, $\epsilon=0.05$.

\end{itemize}

The yellow curve represents a very different type of process, discussed in section~\ref{azxb}.

\subsubsection{Variance of the modified random walk}

The symmetric nature of the modified random walk $\{S_n\}$ defined in section~\ref{qrrnd} results in several

identities. Let $p_n(m)=P(S_n=m)$, with $-n\leq m \leq n$. Also, let $S_0=0$ and $p_0(0)=1$. Then $p_n(m)$ can be recursively computed using some modified version of the Pascal triangle recursion:

\begin{equation}

p_{n+1}(m)=\Big[\frac{1}{2}+\epsilon \cdot A_n(m-1)\Big]p_n(m-1)+\Big[\frac{1}{2}-\epsilon\cdot A_n(m+1)\Big]p_n(m+1),\label{zxzdc}

\end{equation}

where $A_n(m)=\chi[m<-g(n)]-\chi[m>g(n)]$. Here $\chi$ is the indicator function: $\chi(\omega)=1$ if $\omega$ is true, otherwise $\chi(\omega)=0$.

Some of the identities in question include:

$$

\sum_{m=-n}^n m \cdot p_n(m)=\text{E}[S_n]=0,\quad

\sum_{m=-n}^n A_n(m)p_n(m)=0,\quad

\sum_{m=-n}^n m^2 A_n(m)p_n(m)=0,

$$

$$

\sum_{m=-n}^n m^2 \Big[A_{n-1}(m-1)p_{n-1}(m-1)+A_{n-1}(m+1)p_{n-1}(m+1)\Big]=0.

$$

From these identities, it is easy to establish a recursion for the variance:

\begin{equation}

\text{Var}[S_{n+1}]=\text{Var}[S_n]+1-\delta_n, \quad \text{with }\delta_n=8\epsilon\cdot \sum_{m>g(n)} m \cdot p_n(m).\label{varzes}

\end{equation}

The sum for $\delta_n$ is finite since $p_n(m)=0$ if $m>n$. Of course,

$\text{Var}[S_0]=0$. Also, if $\epsilon=0$, the sequence is perfectly random: $\delta_n=0$, $\text{Var}[S_n]=n$ and

$S_n/\sqrt{n}$ converges to a normal distribution. In turn, the law of the iterated logarithm is satisfied. Conversely,

this is violated if $\epsilon>0$. Formula~(\ref{varzes})

combined with \textcolor{index}{Hoeffding's inequality}\index{Hoeffding inequality} [\href{https://en.wikipedia.org/wiki/Hoeffding\%27s_inequality}{Wiki}], may

provide some bounds for $\text{Var}[S_n]$.

\noindent Figure~\ref{fig:lollog1} shows $\delta_n$ for four types of modified random walks, using the following color scheme:

\begin{itemize}

\item Yellow: $g(n)=10,\epsilon=0.05$

\item Red: $g(n)=n^\beta,\beta=0.50, \epsilon=0.05$

\item Blue: $g(n)=n^\beta,\beta=0.45, \epsilon=0.05$

\item Purple: $g(n)=n^\beta,\beta=0.55, \epsilon=0.05$

\end{itemize}

The curve that coincides with the X-axis ($\delta_n = 0$) corresponds to $\epsilon=0$, that is, to pure randomness regardless of $g(n)$. It is not

colored in Figure~\ref{fig:lollog1}. Finally, the Python code in section~\ref{pypy1bv} computes $\text{Var}[S_n]$ exactly (not via simulations) using two different methods, proving that

Formula~(\ref{varzes}) is correct.

%-----------------------------vince/riemann2and3.mp4

\begin{figure}%[H]

\centering

\includegraphics[width=0.9\textwidth]{iteratedlog1.png}

\caption{Same as Figure~\ref{fig:lollog1b}, using a more aesthetic but less meaningful chart type}

\label{fig:lollog1}

\end{figure}

%imgpy9979_2and3.PNG

%-------------------------

\subsection{Random walks with less entropy than pure random signal}\label{azxb}

In section~\ref{azxa}, I focused on creating sequences with higher oscillation rates than dictated by randomness, resulting in lower amplitudes.

Doing the opposite -- decreasing the oscillation rate -- is more difficult. For instance, using $g(n)=n^\beta$ with $\beta>\frac{1}{2}$ won't work.

You can't do better than $\sqrt{n}$ because of the law of the iterated logarithm: boosting $\beta$ beyond the threshold $\frac{1}{2}$ is

useless.

\noindent A workaround is to use the following algorithm: \vspace{1ex}

\noindent\textcolor{white}{00} $S=0$ \\

\textcolor{white}{00} {\bf For} $n = 1, 2,\dots$\\

\textcolor{white}{0000} Generate random deviate $U$ on $[0,1]$\\

\textcolor{white}{0000} $M=g(n)$\\

\textcolor{white}{0000} {\bf If} ($-M< S< 0$ and $U < \frac{1}{2}+\epsilon$)

or ($0<S< M$ and $U < \frac{1}{2}-\epsilon$) or ($S=0$ and $U<\frac{1}{2}$) \\

\textcolor{white}{0000} {\bf Then} \\

\textcolor{white}{000000} $X_n=-1$ \\

\textcolor{white}{0000} {\bf Else} \\

\textcolor{white}{000000} $X_n=1$\\

\textcolor{white}{0000} $S=S+X_n$\\

\textcolor{white}{0000} $S_n=S$ \vspace{1ex}

\noindent The Python code in section~\ref{paths}, with the option \texttt{deviations='Large'}, performs this simulation. The yellow time series in

Figure~\ref{fig:walk} is a realization of such a modified random walk, in this case with $g(n)=\alpha n^\beta$, with

$\alpha=0.30,\beta=0.54$ and $\epsilon=0.01$. It is unclear if the yellow curve will ever cross again the horizontal axis after

$\num{50000}$ iterations, but it is expected to do so. To the contrary, the other three curves (gray, red, blue) are guaranteed to cross the

horizontal axis infinitely many times, even though the random variable $T$ measuring the spacing between two crossings

(referred to as the \textcolor{index}{hitting time}\index{random walk!first hitting time} in section~\ref{poyt}) has infinite expectation.

For pure random walks (the gray curve in Figure~\ref{fig:walk}), the average number of times that $S_k=0$ when $0<k\leq 2n$

is asymptotically equal to $\sqrt{4n/\pi}$ , as discussed in section~\ref{poyt}. One would expect this value to be about $178$

when $2n=\num{50000}$. For the gray curve, the observed value is $243$. Keep in mind that huge variations are expected between different realizations of the same random walk, due to the fact that $\text{E}[T]=\infty$. Indeed, averaged over three realizations, the value $243$ was down to

$185$. Also, a faulty pseudo-random number generator could easily lead to results that are off, in this case.

One would expect much larger

values for the ``non-random" red and blue curves. The observed

values are respectively $747$ and $1783$, based on a single realization in each case. Likewise, the yellow curve is expected to have a much smaller value: in Figure~\ref{fig:walk},

that value is $105$.

\subsection{Python code: computing probabilities and variances attached to $S_n$}\label{pypy1bv}

This Python code is related to section~\ref{qrrnd}, where you can find more details.

It computes the variance of $S_n$ for $n=1, 2$ and so on, using two different methods:

one based on the standard definition of the variance (denoted as \texttt{var1} in the code), and one based

on Formula~\ref{varzes}. The latter is denoted as \texttt{var2} in the code. Also, the variable \texttt{delta}

represents $\delta_n$. The output $\delta_n$ is featured in Figure~\ref{fig:lollog1b}. Finally the function \texttt{G} represents $g(n)$.

The code below is also available on Github, \href{https://github.com/VincentGranville/Machine-Learning/blob/main/Source\%20Code/brownian_var.py}{here}.

Program name: \texttt{brownian\_var.py}. \\

%\pagebreak %

\begin{lstlisting}

import math

epsilon=0.05

beta=0.45

alpha=1.00

nMax=5001

Prob={}

Exp={}

Var={}

Prob[(0,0)] =1

Prob[(0,-1)]=0

Prob[(0,1)] =0

Prob[(0,-2)]=0

Prob[(0,2)] =0

def G(n):

return(alpha*(n**beta))

def psi(n,m):

p=0.0

if m>G(n):

p=-1

if m<-G(n):

p=1

return(p)

Exp[0]=0

Var[0]=0

OUT=open("rndproba.txt","w")

for n in range(1,nMax):

Exp[n]=0

Var[n]=0

delta=0

for m in range(-n-2,n+3,1):

Prob[(n,m)]=0

for m in range(-n,n+1,1):

Prob[(n,m)]=(0.5+epsilon*psi(n-1,m-1))*Prob[(n-1,m-1)]\

+(0.5-epsilon*psi(n-1,m+1))*Prob[(n-1,m+1)]

Exp[n]=Exp[n]+m*Prob[(n,m)]

Var[n]=Var[n]+m*m*Prob[(n,m)]

if m>G(n-1) and m<n:

delta=delta+8*epsilon*m*Prob[(n-1,m)]

var1=Var[n]

var2=Var[n-1]+1-delta

string1=("%5d %.6f %.6f %.6f" % (n,var1,var2,delta))

string2=("%5d\t%.6f\t%.6f\t%.6f\n" % (n,var1,var2,delta))

print(string1)

OUT.write(string2)

OUT.close()

\end{lstlisting}

\subsection{Python code: path simulations}\label{paths}

This Python code performs all the simulations discussed in sections~\ref{qrrnd} and~\ref{azxb} and shown in Figure~\ref{fig:walk}. The option

\texttt{deviations='Small'} is discussed in detail in section~\ref{qrrnd}, while \texttt{deviations='Large'} is

explained in section~\ref{azxb}. The function \texttt{G} in the code corresponds to $g(n)$. Also, if you want to simulate a perfectly random walk, set $\epsilon$ (the parameter \texttt{eps} in the code) to zero. Finally, the code generates multiple realizations for any type of random walk. The number of realizations is determined by the parameter \texttt{Nsample}.

The code below is also available on Github, \href{https://github.com/VincentGranville/Machine-Learning/blob/main/Source\%20Code/brownian_path.py}{here}.

Program name: \texttt{brownian\_path.py}. \\

\begin{lstlisting}

import random

import math

random.seed(1)

n=50000

Nsample=1

deviations='Large'

mode='Power'

if deviations=='Large':

eps=0.01

beta=0.54

alpha=0.3

elif deviations=='Small':

eps=0.05

beta=0.35 #beta = 1 for log

alpha=1

def G(n):

if mode=='Power':

return(alpha*(n**beta))

elif mode=='Log' and n>0:

return(alpha*(math.log(n)**beta))

else:

return(0)

OUT=open("rndtest.txt","w")

for sample in range(Nsample):

print("Sample: ",sample)

S=0

for k in range(1,n):

x=1

rnd=random.random()

M=G(k)

if deviations=='Large':

if ((S>=-M and S<0 and rnd<0.5+eps) or (S<=M and S>0 and rnd<0.5-eps) or

(abs(S)>=M and rnd<0.5) or (S==0 and rnd<0.5)):

x=-1

elif deviations=='Small':

if (S<-M and rnd<0.5-eps) or (S>M and rnd<0.5+eps) or (abs(S)<=M and rnd<0.5):

x=-1

print(k,M,S,x)

S=S+x

line=str(sample)+"\t"+str(k)+"\t"+str(S)+"\t"+str(x)+"\n"

OUT.write(line)

OUT.close()

\end{lstlisting}

\section{Two-dimensional Brownian motions}

There are countless types of random walks or quasi-Brownian motions that are -- on purpose and by design -- not perfectly random. A good reference is the book by Mörters and Peres~\cite{peresbrown}, published in 2010. Here I discuss a 2B Brownian motion generated using some specific probability distributions. Depending on the parameters, these distributions may or may not have an infinite expectation or variance. Things start to get interesting when the expectation becomes infinite (and the Brownian motion is no

longer Brownian), resulting in a system exhibiting a strong clustering structure.

Note that all the simulations performed here consist of discrete random walks rather than

space-continuous Brownian motions.

They approach Brownian motions very well: when increments become infinitesimal, with proper rescaling the result is a continuous process

as illustrated in Figure~\ref{fig:lolbrown}.

In one dimension, we start with $S_0=0$ and $S_n=S_{n-1}+R_n\theta_n$, for $n=1,2$ and so on. If the $R_n$'s are independently and identically distributed (iid) with an exponential distribution of expectation $1/\lambda$ and $\theta_n=1$, then the resulting process is a stationary