branching processes. In detail: the new method does not require processing
spike-time recordings by attaching time bins and analysing the resulting step
process. The new method directly works with spike time recordings and
calculates ratios of the moments of the inter-spike interval. Once the ratios
are calculated, the parameters of the process can be read off a map, the
moment-ratio map, Fig. 11. The user can generate this map to the desired
precision by using the provided python codePausch (2020). Some conceptual
drawbacks were pointed out in the setup of the model, Sec. II.1. The new
approach has some significant success, of which several are highlighted in the
following.

On the bare level of inter-spike intervals, the new method predicts moments,
the coefficient of variation and its skewness to high accuracy. Their
approximation errors were analysed in Fig. 21, which showed small percentage
errors. This shows that the method and process represent the observed spike
time series very well.

In comparison to other methods, the approach avoids creating a branching-like
processes and allows analysing inter-spike intervals directly. In particular,
the proposed methodology does not require time-binning and the defining
parameters (degree of criticality $\frac{r}{s}$ and relative spontaneous
creation $\frac{\gamma}{s}$) are dimensionless. Therefore, it is less
sensitive compared to a number of commonly used data processing
methods.Priesemann  _et al._ (2009); Wilting and Priesemann (2018) More
precisely, it allows identifying how close the system is to criticality with a
time-scale independent measure and can therefore be regarded advantageous over
criticality estimates based on either time-bin-dependent power-law
distributions of avalanche sizes or on time-bin-dependent auto-correlation
functions.Priesemann  _et al._ (2009); Wilting and Priesemann (2018) In
particular, the proposed new method applied to experimental data indicates
that neuronal circuits are further away from criticality than most articles
suggest.Beggs and Plenz (2003, 2004); Beggs (2008); Priesemann  _et al._
(2009, 2013); Wilting and Priesemann (2018); Wilting  _et al._ (2018)

The most common way of confirming criticality in the brain has been the
identification of a power law distribution of avalanche sizes. However, that
approach is a notoriously difficult one in statistical physics and is prone to
introduce biases.Goldstein  _et al._ (2004); Priesemann and Shriki (2018)
These difficulties are avoided in this article by looking at a moment-ratio
map of the inter-spike intervals. It avoids generating a branching-like time
series and thus also avoids fitting power laws. In addition, the effects of
binning on a pumped branching process strongly resemble its observed effects
when applied to experimental data, Figs. 15 and 16. Consequently and more
importantly, it indicates that neuronal avalanche dynamics in the used in vivo
data sets Mizuseki  _et al._ (2009) is not scale free and not at criticality.
The used in vitro data sets Wagenaar  _et al._ (2006) show that criticality
increases for a plated developing neural network over time and that external
input is significantly lower compared to the in vivo data.

The external signal immigration (‘pumping’) into the locally measured neuronal
circuit is also determined in this article. It suggests that for the in vivo
data, avalanche overlap is relatively common and cannot be ignored. In
particular, the common separation-of-timescales (STS) assumption is not
justified for the used in vivo data Mizuseki  _et al._ (2009). It is therefore
also expected that the avalanche size distribution does not follow a power law
distribution with exponent $-3/2$ even if it were at criticality.di Santo  _et
al._ (2017) In contrast for the in vitro data Wagenaar  _et al._ (2006),
external input was rare and the STS assumption is valid.

Another advantage of the presented model is that relaxation behavior can be
predicted analytically. The presented methodology would allow to investigate
whether different steady states, such as wakefulness or deep sleep Priesemann
_et al._ (2013), relate to changes in external input $\frac{\gamma}{s}$ or to
internal adjustment of branching efficiency $\frac{r}{s}$.

The presented model and data analysis also opens several routes for future
research. Including topological constraints in a pumped branching process is a
major challenge which will most likely have significant impact on more
detailed comparisons with experiments. In addition, when considering
subsampling, Sec. III.6, the data indicated beyond-mean-field subsampling
which would require detailed knowledge of the connectome for more accurate
models. Furthermore, better understanding finite size effects, which are
unintentionally introduced due to the size of used probes, might clarify some
of the anomalous statistics observed. In a different direction, future work
could also try to understand mechanism that cause inter-avalanche
correlations.

## Acknowledgements

The author wishes to thank Gunnar Pruessner, Ignacio Bordeu Weldt, Rosalba
Garcia Milan, Benjamin Walter, and Viola Priesemann for helpful discussions.
He received financial support through an EPSRC PhD scholarship and an EPSRC
Doctoral Prize Fellowship at Imperial College London.

##  Appendix A Derivation of spontaneous creation action

Let $P(N(t)=n)$ denote the probability that there are $n$ particles in the
system at time $t$. Then, a spontaneous creation of particles can be described
by a master equation van Kampen (1992) as follows

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}P(N(t)=n)=\gamma\Bigl{(}P\bigl{(}N(t)=n-1\bigr{)}-P\bigl{(}N(t)=n\bigr{)}\Bigr{)}.$
(71)

Such a master equation can be transformed into an action of a Doi-Peliti field
theory Täuber (2014); Garcia-Millan  _et al._ (2018); Pausch (2019), which
results in the action $\mathcal{A}_{c}$ shown in Eq. (2).

For comparison, the branching/extinction process with offspring distribution
$p_{k}$ and occurrence rate $s$ is described the following master equation

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}P(N(t)=n)=$ $\displaystyle
s\sum\limits_{k=0}^{\infty}p_{k}(n-k+1)P(N(t)=n-k+1)$ $\displaystyle-
snP(N(t)=n),$ (72)

derived in Garcia-Millan  _et al._ (2018).

Although written as a single equation, master equations (71) and (72)
represent one ordinary differential equations for each $n\in\mathbb{N}_{0}$.
Thus, they actually describe an infinite system of coupled ordinary
differential equations.

##  Appendix B $\langle\phi^{m}\rangle=\mathbb{E}[(N(t))_{m}]$

In Doi-Peliti field theories, the $m$th factorial moment naturally appears as
observable and takes the form $\phi^{m}$. This is explained at length in
Pausch (2019) and a brief outline is given here.

The actions of Doi-Peliti field theories are derived from Master equations
which take the form

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}P(N(t))=\mathcal{L}P(N(t)),$ (73)

where $P(N(t))=(P(N(t)=0),P(N(t)=1),\dots)^{T}$ is a vector and $\mathcal{L}$
is a linear transformation. Eq. (73) is a system of infinitely many coupled
ordinary differential equations. It can be transformed into a single partial
differential equation for the probability generating function, defined as

$\displaystyle\mathcal{M}(z,t)=\sum\limits_{n=0}^{\infty}P(N(t)=n)z^{n},$ (74)

by using derivatives w.r.t. $z$ or multiplying by $z$. For example, Eq. (71)
is transformed into

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{M}(z,t)=\gamma\left(z-1\right)\mathcal{M}(z,t),$
(75)

and Eq. (72) is transformed into

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{M}(z,t)=s\left(\sum\limits_{n=0}^{\infty}p_{k}z^{k}-z\right)\frac{\mathrm{d}}{\mathrm{d}z}\mathcal{M}(z,t).$
(76)

From this PDE perspective, path integrals are just solving PDEs perturbatively
around $z=0$. However in most references, $z$ and $\mathrm{d}/\mathrm{d}z$ are
represented as ladder operators $a^{\dagger}$ and $a$, respectively, and
interpreted as particle creators and annihilators.

A useful step is introducing a new variable $\widetilde{z}$ by setting
$z=\widetilde{z}+1$, and work with
$\widetilde{\mathcal{M}}(\widetilde{z},t)=\mathcal{M}(z,t)$. This is called
the Doi-shift. The path integral solves now the corresponding PDE for
$\widetilde{\mathcal{M}}$ perturbatively around $\widetilde{z}=0$, i.e.

$\displaystyle\langle\phi^{m}(t)\rangle=\frac{\mathrm{d}^{m}}{\mathrm{d}z^{m}}\widetilde{\mathcal{M}}(\widetilde{z},t)\bigr{|}_{\widetilde{z}=0}=\mathbb{E}[(N(t))_{m}],$
(77)

where $(N)_{m}=N(N-m)\dots(N-m+1)$ is the falling factorial.

##  Appendix C Steady-state probability distribution

Using the fact that the probability generating function is equal to the
factorial moment generating function, the probability distribution of the
particle number $N$ can be deduced from the factorial moments
$\langle\phi^{k}\rangle$:

$\displaystyle P(N=\ell)=$
$\displaystyle\frac{1}{\ell!}\frac{\mathrm{d}^{\ell}}{\mathrm{d}z^{\ell}}M_{N}(z)\Biggr{|}_{z=0}=\sum\limits_{k=\ell}^{\infty}\frac{(-1)^{k-\ell}}{(k-\ell)!}\left(\frac{q_{2}}{r}\right)^{k}\left(\frac{\gamma}{q_{2}}\right)^{(k)}$
$\displaystyle=$
$\displaystyle\frac{1}{\ell!}\left(\frac{q_{2}}{r}\right)^{\ell}\sum\limits_{m=0}^{\infty}\frac{(-1)^{m}}{m!}\left(\frac{q_{2}}{r}\right)^{m}\left(\frac{\gamma}{q_{2}}\right)^{(m+\ell)}$
$\displaystyle=$
$\displaystyle\frac{1}{\ell!}\left(\frac{q_{2}}{r}\right)^{\ell}\left(\frac{\gamma}{q_{2}}\right)^{(\ell)}\sum\limits_{m=0}^{\infty}\frac{(-1)^{m}}{m!}\left(\frac{q_{2}}{r}\right)^{m}\left(\frac{\gamma}{q_{2}}+\ell\right)^{(m)}$
$\displaystyle=$
$\displaystyle\frac{1}{\ell!}\left(\frac{q_{2}}{r}\right)^{\ell}\left(\frac{\gamma}{q_{2}}\right)^{(\ell)}\sum\limits_{m=0}^{\infty}\frac{1}{m!}\left(\frac{q_{2}}{r}\right)^{m}\left(-\frac{\gamma}{q_{2}}-\ell\right)_{m}$
$\displaystyle=$
$\displaystyle\frac{1}{\ell!}\left(\frac{q_{2}}{r}\right)^{\ell}\left(\frac{\gamma}{q_{2}}\right)^{(\ell)}\left(1+\frac{q_{2}}{r}\right)^{-\frac{\gamma}{q_{2}}-\ell}$
(78)

where in the second to last line, the umbral Taylor expansion of
$(1+q_{2}/r)^{-\gamma/q_{2}-\ell}$ around $q_{2}/r=0$ is identified.
$(x)^{(\ell)}=x(x+1)\dots(x+\ell-1)$ is the rising factorial.

Let’s check that this probability distribution is normalized:

$\displaystyle\sum\limits_{\ell=0}^{\infty}P(N=\ell)=$
$\displaystyle\frac{1}{\left(1+\frac{q_{2}}{r}\right)^{\frac{\gamma}{q_{2}}}}\sum\limits_{\ell=0}^{\infty}\frac{1}{\ell!}\left(\frac{q_{2}}{r+q_{2}}\right)^{\ell}\left(\frac{\gamma}{q_{2}}\right)^{(\ell)}$
$\displaystyle=$
$\displaystyle\frac{1}{\left(1+\frac{q_{2}}{r}\right)^{\frac{\gamma}{q_{2}}}}\sum\limits_{\ell=0}^{\infty}\frac{1}{\ell!}\left(\frac{-q_{2}}{r+q_{2}}\right)^{\ell}\left(\frac{-\gamma}{q_{2}}\right)_{\ell}$
$\displaystyle=$
$\displaystyle\frac{1}{\left(1+\frac{q_{2}}{r}\right)^{\frac{\gamma}{q_{2}}}}\left(1-\frac{q_{2}}{r+q_{2}}\right)^{-\frac{\gamma}{q_{2}}}=1$
(79)

##  Appendix D Initializing an empty system

In the operator picture, where $z=a^{\dagger}$ and
$\frac{\mathrm{d}}{\mathrm{d}z}=a$, states of the system are described as Fock
states $|n\rangle=z^{n}$. Their dual states over the $L^{2}$ inner product are
compactly written as $\langle n|$. Before time $t_{0}=0$, the system is
allowed to evolve freely resulting in state
$|\mathcal{M}(t_{0}^{-})\rangle=\sum_{n=0}^{\infty}P(N(t_{0}^{-})=n)|n\rangle$.
But at $t_{0}=0$, it is projected out using
$\langle\sun|=\sum_{n=0}^{\infty}\langle n|$. Because $\langle
m|n\rangle=\delta_{mn}$, the result is
$\langle\sun|\mathcal{M}(t_{0}^{-})\rangle=1$. Then, the system is
reinitialized with zero particles by attaching the empty state $|0\rangle$.
Thus, given that the system was empty at time $t_{0}=0$, the $n$th factorial
moment at time $t$ equals

$\displaystyle\mathbb{E}[(N)_{n}(t)|N(0)=0]=\langle\sun|a^{\dagger
n}a^{n}e^{\mathcal{A}t}|0\rangle\langle\sun|\mathcal{M}(t_{0})\rangle$ (80)

In the field theory, the projection with $\langle\sun|$ and reinitialization
with $|0\rangle$ is expressed as the observable
$e^{-\widetilde{\phi}(0)\phi(0)}$. See Pausch (2019) for more details on how
to derive the field-theoretic expressions from operator expressions.

##  Appendix E Inter-event time distribution

In steady state, the branching process’ particle number probability
distribution is a negative binomial (or Pólya) distribution with the general
form

$\displaystyle
P(N=\ell)=\frac{1}{\ell!}\frac{\Gamma(u+\ell)}{\Gamma(u)}p^{u}(1-p)^{\ell},$
(81)

with $u=\gamma/q_{2}$, $p=r/(r+q_{2})$. Its moment generating function equals

$\displaystyle\mathcal{M}_{N}(x)=\left(\frac{p}{1-(1-p)e^{x}}\right)^{u}\qquad\text{for
}x<-\ln(1-p)$ (82)

This distribution can be regarded as a prior distribution to a conditional
exponential distribution $f_{e}(t|\ell)$ with rate $\lambda=\ell s+\gamma$:

$\displaystyle f_{e}(t|\ell)=(s\ell+\gamma)e^{-(s\ell+\gamma)t}.$ (83)

The marginal event time distribution $f_{e}(t)$ can be recovered as

$\displaystyle f_{e}(t)=$
$\displaystyle\sum\limits_{\ell=0}^{\infty}f_{e}(t|\ell)P(N=\ell)$ (84)
$\displaystyle=$ $\displaystyle e^{-\gamma
t}\left(\gamma-\frac{\partial}{\partial t}\right)\mathcal{M}_{N}(-st)$
$\displaystyle=$ $\displaystyle\gamma e^{-\gamma
t}\left(\frac{(r+q_{2})(1+e^{-st})}{r+q_{2}(1-e^{-st})}\right)\left(\frac{r}{r+q_{2}(1-e^{-st})}\right)^{\frac{\gamma}{q_{2}}}$

##  Appendix F Verification of moments of the inter-event and inter-spike
interval

It is important to verify analytical derivations as much as possible and
reasonable and to make such verifications available. Here, the moments of the
inter-event time and the inter-spike time are verified using simulations in
the parameter space that is relevant for comparison with experimental data.
Fig. 18 shows the first three moments of the inter-event interval over a range
of $\frac{r}{s}$, the degree of criticality. As expected, intervals become
shorter as criticality $\frac{r}{s}=0$ is approached.

Figure 18: Rescaled first, second and third moment of the time between events
over $r/s\in[0,1]$. Events include both extinctions and creations. The
parameter values are $\gamma/s\in\\{0.1,0.6,1.0\\}$. Symbols: simulation
results. Lines: Analytical prediction. For the first moment, all lines
collapse after rescaling.

Fig. 19 shows the first three moments of the inter-spike interval over a range
of $\frac{r}{s}$ values. Analogously to the inter-event intervals, the inter-
spike intervals tend to zero as criticality is approached.

Although the moments of intervals between events and spikes appear to be very
similar, they show an important difference in the limit of a pure Poisson
process $\frac{r}{s}\rightarrow 1$. In this limit, spike times will be equal
to spontaneous creation times which occur with rate $\gamma$. Event intervals
however still depend on the extinction time scale $sp_{0}\rightarrow s$.

Figure 19: Rescaled first, second and third moment of the time between spikes,
i.e. between creations of particles, over $r/s\in[0,1]$. The parameter values
are $\gamma/s\in\\{0.1,0.6,1.0\\}$. Symbols: simulation results. Lines:
Analytical prediction. For the first moment, all lines collapse after
rescaling.

##  Appendix G Derivation of the moments of inter-spike intervals given
initial state

The result in Eq. (66) can be found by starting to calculate it explicitly for
$m=0,1,2,\dots$ from which the general rule can be derived.

Given the system is in state $n=0$, the next creation event time is
$\sim\text{Exp}(\gamma)$ distributed. If the system is in state $n=1$, the
next creation event time has the following moments

$\displaystyle\mathbb{E}[T_{c}^{n}|N=1]=$
$\displaystyle\int\limits_{0}^{\infty}\biggl{(}t^{n}(s+\gamma)e^{-(s+\gamma)t}\frac{sp_{2}+\gamma}{s+\gamma}+$
(85)
$\displaystyle+\int\limits_{0}^{\infty}(t+t^{\prime})^{n}(s+\gamma)e^{-(s+\gamma)t}\frac{sp_{0}}{s+\gamma}\gamma
e^{-\gamma t^{\prime}}\mathrm{d}t^{\prime}\biggr{)}\mathrm{d}t$
$\displaystyle=(sp_{2}+\gamma)\left(-\frac{\mathrm{d}}{\mathrm{d}\gamma}\right)^{n}\frac{1}{s+\gamma}+sp_{0}\gamma\left(-\frac{\mathrm{d}}{\mathrm{d}\gamma}\right)^{n}\frac{1}{\gamma(s+\gamma)}.$

If the system is in state $n=2$, the moments are

$\displaystyle\mathbb{E}[T^{n}|N=2]=$
$\displaystyle\int\limits_{0}^{\infty}t^{n}(2sp_{2}+\gamma)e^{-(2s+\gamma)t}\mathrm{d}t+$
(86)
$\displaystyle+\int\limits_{0}^{\infty}\int\limits_{0}^{\infty}(t+t^{\prime})^{n}2p_{0}e^{-(2s+\gamma)t}(sp_{2}+\gamma)e^{-(s+\gamma)t^{\prime}}\mathrm{d}t^{\prime}\mathrm{d}t+$
$\displaystyle+\iiint\limits_{0^{3}}^{\infty^{3}}(t+t^{\prime}+t^{\prime\prime})^{n}2(sp_{0})^{2}\gamma
e^{-(2s+\gamma)t-(s+\gamma)t^{\prime}-\gamma
t^{\prime\prime}}\mathrm{d}t^{\prime\prime}\mathrm{d}t^{\prime}\mathrm{d}t$
$\displaystyle=(2sp_{2}+\gamma)\left(-\frac{\mathrm{d}}{\mathrm{d}\gamma}\right)^{n}\frac{1}{2s+\gamma}+$
$\displaystyle+2sp_{0}(sp_{2}+\gamma)\left(-\frac{\mathrm{d}}{\mathrm{d}\gamma}\right)^{n}\frac{1}{(s+\gamma)(2s+\gamma)}+$
$\displaystyle+2sp_{0}sp_{0}\gamma\left(-\frac{\mathrm{d}}{\mathrm{d}\gamma}\right)^{n}\frac{1}{\gamma(s+\gamma)(2s+\gamma)}$

This can be generalized by induction for any condition $N=m$. If the system is
in state $N=m$, then after an exponentially distributed amount of time with
rate $sm+\gamma$, either a particle is created, i.e. a spike occurs, or a
particle goes extinct. After the the extinction event, the previous result
$\mathbb{E}[T_{s}^{n}|N=m-1]$ can be used and the result is shown in Eq. (66).

##  Appendix H Phase space boundary for $\gamma\rightarrow 0$

Considering Eq. (66), in the limit $\gamma\rightarrow 0$, the inter-spike
interval moments conditioned on initial state $N=m$ diverge as follows

$\displaystyle\lim\limits_{\gamma\rightarrow 0}\mathbb{E}[T_{s}^{n}|N=m]\sim$
$\displaystyle\frac{n!p_{0}^{m}}{\gamma^{n}}$ (87)

On the other hand, the limit of $f(N=m)$ for $\gamma\rightarrow 0$ is

$\displaystyle\lim\limits_{\gamma\rightarrow 0}f(N=m)=$
$\displaystyle\,\lim\limits_{\gamma\rightarrow
0}\frac{r^{u+1}\Gamma(u+m)q_{2}^{m}(sm+\gamma)}{m!\Gamma(u)\gamma(r+s)(r+q_{2})^{u+m}}$
$\displaystyle=$ $\displaystyle\,\frac{srq_{2}^{m-1}}{(r+s)(r+q_{2})^{m}},$
(88)

which implies for $f(m-1\uparrow m)$

$\displaystyle\lim\limits_{\gamma\rightarrow 0}f(m-1\uparrow
m)=\frac{srq_{2}^{m-2}}{(r+s)(r+q_{2})^{m-1}}2p_{2}.$ (89)

Combining Eq. (87) and Eq. (89) gives

$\displaystyle\lim\limits_{\gamma\rightarrow
0}\mathbb{E}[T_{s}^{n}]\sim\frac{n!}{\gamma^{n}}\frac{2r}{(r+s)}.$ (90)

This implies for the coefficient of variation $c_{V}$

$\displaystyle\lim\limits_{\gamma\rightarrow 0}c_{V}=\sqrt{\frac{s}{r}}.$ (91)

Furthermore, the following ratios can be calculated:

$\displaystyle x=$ $\displaystyle\lim\limits_{\gamma\rightarrow
0}\frac{\mathbb{E}[T_{s}^{3}]}{\mathbb{E}[T_{s}]^{3}}-6=6\left(\frac{(r+s)^{2}}{4r^{2}}-1\right)$
(92) $\displaystyle\Leftrightarrow$
$\displaystyle\,\frac{(r+s)}{2r}=\sqrt{(x+6)/6}$ (93) $\displaystyle y=$
$\displaystyle\lim\limits_{\gamma\rightarrow
0}\frac{\mathbb{E}[T_{s}^{4}]}{\mathbb{E}[T_{s}^{2}]^{2}}-6$ (94)
$\displaystyle=$
$\displaystyle\,6\left(\frac{(r+s)}{2r}-1\right)=6\left(\sqrt{(x+6)/6}-1\right),$
(95)

which describes the phase space boundary in the moment-ratio map.

##  Appendix I Verification of the moment-ratio map

Figure 20: Moment-ratio map. Blue lines: analytically calculated
$\frac{\gamma}{s}$-level sets. Green lines: analytically calculated
$\frac{r}{s}$-level sets. Symbols: Monte Carlo simulations results. The area
below the thick blue line for $\frac{\gamma}{s}=0$ is not in the phase space
of pumped branching processes.

In Fig. 11A, data points from experiments were plotted onto a map of spike
moments, which was then transformed into a map in
$\frac{r}{s}$-$\frac{\gamma}{s}$-space. In this section, the analytics in this
figure are verified using simulations which are shown in Fig. 20.

Although the formulas in Eqs. (66), (67) for the level sets are analytical,
calculating specific values requires special functions and high precision
numerical evaluation, which makes accessing this map difficult for areas in
the top left of Fig. 20, i.e. for small $X$ and large $Y$. A program for
calculating moments of inter-spike intervals induced by pumped branching
processes is available at Pausch (2020). The area below the thick blue line
for $\frac{\gamma}{s}=0$ is not in the phase space of the spike time series
induced by pumped branching processes.

##  Appendix J Approximation Errors

Figure 21: Approximation errors for in vivo data sets from Mizuseki  _et al._
(2009) and in vitro data sets from Wagenaar  _et al._ (2006). A: Signed
approximation errors of the coefficient of variation $c_{V}$ and skewness
$\tilde{\mu}_{3}$. B: Signed approximation error of the 2nd and 3rd moments of
the inter-spike interval $T_{s}$. All approximation errors appear to indicate
that the estimation is unbiased over all data sets.

The moment-ratio map in Fig. 11A is used to identify values for the degree of
criticality $\frac{r}{s}$ and relative spontaneous creation $\frac{\gamma}{s}$
of experimental data, which then are used to identify the time scale $s$ using
the first moment of inter-spike interval $T_{s}$ alone. There are many reasons
why these parameter values have errors. Firstly, the pumped branching model
and its induced spike process is an idealization of the true process in the
brain. Second, experimental measurements always contain errors due to used
instruments and choice of data processing. Furthermore, the identification of
the parameters is done algorithmically using a gradient descent method and has
inherent inaccuracies that are chosen to limit run time. In addition, it is
based on calculation using 128bit double precision. Although 128bit double
precision seems to be an unusual restriction, it is a serious one in this case
because close to criticality, the probabilities for the system to be in a
specific state $N=n$ tend to zero. This implies that the calculation of the
moments of $T_{s}$ involves summing over many near-zero probabilities and thus
the 128bit double precision becomes a limiting factor and a source of errors.

The effect of all of these sources of errors result in approximation errors of
the true inter-spike interval distribution. The magnitude of these errors
determines the accuracy and predictive power of the model. For an observable
$\mathcal{O}$, the signed approximation error is defined as

$\displaystyle\delta_{\mathcal{O}}=\frac{\mathcal{O}_{\text{data}}-\mathcal{O}_{\text{predicted}}}{\mathcal{O}_{\text{predicted}}}.$
(96)

Fig. 21 shows the signed approximation errors of the coefficient of variation
$c_{V}$ and the skewness $\tilde{\mu}_{3}$ (Panel A), as well as the 2nd and
3rd moments of the inter-spike interval (Panel B). Except for the skewness,
the approximation errors lie in a range of a few percent and seem to be
unbiased over all data sets. As the skewness and coefficient of variation are
both depending on the mean and standard deviation, a systematic relationship
between their errors is not surprising.

##  Appendix K Description of data collection procedure

The first batch of data sets Mizuseki  _et al._ (2009) was downloaded from
crcns.org, where it is called data set hc-2. The following is a quote from the
description document of the data set:

’Three male Long-Evans rats (rat ID; ec13, ec14, ec16, 250-400 g) were
implanted with a 4-shank or 8-shank silicon probe in layer CA1 of the right
dorsal hippocampus. The individual silicon probes were attached to respective
micromanipulators and moved independently. Each shank had 8 recording sites
(160 $\mu$m2 each site; 1-3 M$\Omega$ impedance). These recordings sites were
staggered to provide a two-dimensional arrangement (20 $\mu$m vertical
separation; Fujisawa et al., 2008). The shanks were aligned parallel to the
septo-temporal axis of the hippocampus (45 degrees parasagittal), positioned
centrally at anteroposterior=3.5 mm from bregma and mediolateral=2.5 mm from
midline. Two stainless steel screws inserted above the cerebellum were used as
indifferent and ground electrodes during recordings. All protocols were
approved by the Institutional Animal Care and Use Committee of Rutgers
University.’

The second batch of data sets Wagenaar  _et al._ (2006) was downloaded from
neurodatasharing.bme.gatech.edu/development-data/. All data sets labelled as
’dense’ are used in this article. The following is paraphrasing and a short
summary of the full description in Wagenaar  _et al._ (2006): For the data
collection, multi-electrode arrays were used with 59 electrodes of 30$\mu$m
diameter. The electrodes were arranged in a square grid with 200$\mu$m spacing
between them, measured center-to-center. Neurons and glia were obtained from
rat embryos. Quote: ”[…] timed-pregnant Wistar rats were sacrificed using
CO-32 inhalation, according to NIH approved protocols, at day 18 of gestation.
Embryos were removed and decapitated, and the anterior part of the cortices
(including somatosensory, motor, and association areas) were dissected out. At
the rostral edge, the boundary with the olfactory bulb was used as a landmark;
at the caudal edge, the third ventricle and the boundary with the lateral horn
of the hippocampus were used as landmarks. Striatum and hippocampus were not
included. Cortices from several embryos from the same litter were combined,
[…]”. For the maintenance of the cultures, quote: ”Half of the medium was
replaced approximately every five days in most experiments […]. To test
whether feeding schedule affected activity, all medium was replaced every
seven days in some experiments (N = 3). This did not result in significantly
different activity patterns compared with sister cultures. Feeding always took
place after the day’s recording session, to allow at least 12 hours for
transient effects to disappear before the next recording.”

## References

  * Hodgkin and Rushton (1946) A. Hodgkin and W. Rushton,  _The electrical constants of a crustacean nerve fiber,_Proc. R. Soc. Lond. B 133, 444–479 (1946).
  * Stein (1965) R. B. Stein,  _A theoretical analysis of neuronal variability,_Biophys. J. 5, 173 (1965).
  * Pfeiffer and Kiang (1965) R. Pfeiffer and N. Y.-S. Kiang,  _Spike discharge patterns of spontaneous and continuously stimulated activity in the cochlear nucleus of anesthetized cats,_Biophys. J. 5, 301 (1965).
  * Lamarre  _et al._ (1971) Y. Lamarre, M. Filion, and J. Crodeau,  _Neuronal discharge of the ventrolateral nucleus of the thalamus during sleep and wakefulness in the cat,_Exp. Brain Res. 12, 480 (1971).
  * Wilbur and Rinzel (1983) W. J. Wilbur and J. Rinzel,  _A theoretical basis for large coefficient of variation and bimodality in neuronal interspike interval distributions,_J. theor. Biol. 105, 345 (1983).
  * Softky and Koch (1993) W. R. Softky and C. Koch,  _The highly irregular firing of cortical cells is inconsistent with temporal integration of random epsps,_J. Neurosci. 13, 334 (1993).
  * Dayan and Abbott (2001) P. Dayan and L. Abbott,  _Theoretical Neuroscience_ (Massachusetts Institute of Technology, 2001).
  * Beggs and Plenz (2003) J. M. Beggs and D. Plenz,  _Neuronal avalanches in neocortical circuits,_J. Neurosci. 23, 11167 (2003).
  * Beggs and Plenz (2004) J. M. Beggs and D. Plenz,  _Neuronal avalanches are diverse and precise activity patterns that are stable for many hours in cortical slice cultures,_J. Neurosci. 24, 5216 (2004).
  * Beggs (2008) J. M. Beggs,  _The criticality hypothesis: how local cortical networks might optimize information processing,_Phil. Trans. R. Soc. A 366, 329 (2008).
  * Priesemann  _et al._ (2009) V. Priesemann, M. H. Munk, and M. Wibral,  _Subsampling effects in neuronal avalanche distributions recorded in vivo,_BMC Neurosci. 10, 1 (2009).
  * Priesemann  _et al._ (2013) V. Priesemann, M. Valderrama, M. Wibral, and M. L. V. Quyen,  _Neuronal avalanches differ from wakefulness to deep sleep – evidence from intracranial depth recordings,_PLoS Comput. Biol. 9 (2013), 10.1371/journal.pcbi.1002985.
  * Wilting  _et al._ (2018) J. Wilting, J. Dehning, J. P. Neto, L. Rudelt, M. Wibral, J. Zierenberg, and V. Priesemann,  _Operating in a reverberating regime enables rapid tuning of network states to task requirements,_Front. Syst. Neurosci. 12 (2018), 10.3389/fnsys.2018.00055.
  * Wilting and Priesemann (2018) J. Wilting and V. Priesemann,  _Inferring collective dynamical states from widely unobserved systems,_Nat. commun. 9, 1 (2018).
  * Wilting and Priesemann (2019a) J. Wilting and V. Priesemann,  _Between perfectly critical and fully irregular: A reverberating model captures and predicts cortical spike propagation,_Cereb. Cortex 29, 2759 (2019a).
  * Wilting and Priesemann (2019b) J. Wilting and V. Priesemann,  _25 years of criticality in neuroscience – established results, open controversies, novel concepts,_Curr. Opin. Neurobiol. 58, 105 (2019b).
  * Stephanie R. Miller and Plenz (2019) S. Y. Stephanie R. Miller and D. Plenz,  _The scale-invariant, temporal profile of neuronal avalanches in relation to cortical $\gamma$-oscillations, _Sci. Rep. 2019, 1 (2019).
  * Tuckwell (1979) H. C. Tuckwell,  _Synaptic transmission in a model of stochastic neural activity,_J.theor. Biol. 77, 65 (1979).
  * Christodoulou and Bugmann (2001) C. Christodoulou and G. Bugmann,  _Coefficient of variation vs. mean interspike interval curves: What do they tell us about the brain?_Neurocomputing 38–40, 1141 (2001).
  * Moreno-Bote (2014) R. Moreno-Bote,  _Poisson-like spiking in circuits with probabilisticsynapses,_PLOS Comput. Biol. 10, 1 (2014).
  * Lengler and Steger (2017) J. Lengler and A. Steger,  _Note on the coefficient of variations of neuronal spike trains,_Biol. Cybern. 111, 229–235 (2017).
  * Mijatović  _et al._ (2018) G. Mijatović, T. Lončar-Turukalo, E. Procyk, and D. Bajić,  _A novel approach to probabilistic characterisation of neuralfiring patterns,_J. Neurosci. Methods 305, 67 (2018).
  * Zierenberg  _et al._ (2020) J. Zierenberg, J. Wilting, V. Priesemann, and A. Levina,  _Description of spreading dynamics by microscopic network models and macroscopic branching processes can differ due to coalescence,_Phys. Rev. E 101 (2020), 10.1103/PhysRevE.101.022301.
  * Bellac (1991) M. L. Bellac,  _Quantum and Statistical Field Theory_ (Oxford Science Publications, 1991).
  * Garcia-Millan  _et al._ (2018) R. Garcia-Millan, J. Pausch, B. Walter, and G. Pruessner,  _Field-theoretic approach to the universality of branching processes,_ Phys. Rev. E 98, 1 (2018).
  * Paczuski  _et al._ (1996) M. Paczuski, S. Maslov, and P. Bak,  _Avalanche dynamics in evolution, growth, and depinning models,_Phys. Rev. E 53, 414 (1996).
  * di Santo  _et al._ (2017) S. di Santo, P. Villegas, R. Burioni, and M. A. Muñoz,  _Simple unified view of branching process statistics: Random walks in balanced logarithmic potentials,_Phys. Rev. E 95, 1 (2017).
  * Goldstein  _et al._ (2004) M. Goldstein, S. Morris, and G. Yen,  _Problems with fitting to the power-law distribution,_Eur. Phys. J. B 41, 255 (2004).
  * Priesemann and Shriki (2018) V. Priesemann and O. Shriki,  _Can a time varying external drive give rise to apparent criticality in neural systems?_PLoS Comput. Biol. 14 (2018), 10.1371/journal.pcbi.1006081.
  * Lombardi  _et al._ (2012) F. Lombardi, H. Herrmann, C. Perrone-Capano, D. Plenz, and L. DeArcangelis,  _Balance between excitation and inhibition controls the temporal organization of neuronal avalanches,_Phys. Rev. Lett. 108, 1 (2012).
  * Lombardi  _et al._ (2014) F. Lombardi, H. J. Herrmann, D. Plenz, and L. DeArcangelis,  _On the temporal organization of neuronal avalanches,_Front. Syst. Neurosci. 8, 1 (2014).
  * Petermann  _et al._ (2009) T. Petermann, T. C. Thiagarajan, M. A. Lebedev, M. A. L. Nicolelis, D. R. Chialvo, , and D. Plenz,  _Spontaneous cortical activity in awake monkeys composed of neuronal avalanches,_Proc. Natl. Acad. Sci. 106, 15921 (2009).
  * Martinello  _et al._ (2017) M. Martinello, J. Hidalgo, A. Maritan, S. di Santo, D. Plenz, and M. A. Muñoz,  _Neutral theory and scale-free neural dynamics,_Phys. Rev. X 7, 1 (2017).
  * Pruessner (2012) G. Pruessner,  _Self-Organised Criticality_ (Cambridge University Press, 2012).
  * Harris (1964) T. E. Harris,  _The Theory of Branching Processes_ (RAND, 1964).
  * Mizuseki  _et al._ (2009) K. Mizuseki, A. Sirota, E. Pastalkova, and G. Buzsáki,  _Theta oscillations provide temporal windows for local circuit computation in the entorhinal-hippocampal loop,_Neuron 64, 267 (2009), link to data set: http://dx.doi.org/10.6080/K0Z60KZ9, protocols were approved by the Institutional Animal Care and Use Committee of Rutgers University.
  * Zierenberg  _et al._ (2018) J. Zierenberg, J. Wilting, and V. Priesemann,  _Homeostatic plasticity and external input shape neural network dynamics,_Phys. Rev. X 8, 1 (2018).
  * Das and Levina (2019) A. Das and A. Levina,  _Critical neuronal models with relaxed timescale separation,_Phys. Rev. X 9, 1 (2019).
  * Yu  _et al._ (2017) S. Yu, T. l. Ribeiro, C. Meisel, S. Chou, A. Mitz, R. Saunders, and D. Plenz,  _Maintained avalanche dynamics during task-induced changes of neuronal activity in nonhuman primates,_eLife 6, 1 (2017).
  * Fontenele  _et al._ (2019) A. J. Fontenele, N. A. P. de Vasconcelos, T. Feliciano, L. A. A. Aguiar, C. Soares-Cunha, B. Coimbra, L. D. Porta, S. Ribeiro, A. J. Rodrigues, N. Sousa, P. V. Carelli, and M. Copelli,  _Criticality between cortical slices,_Phys. Rev. Lett. 122, 1 (2019).
  * Harris (1974) T. Harris,  _Contact interactions on a lattice,_Ann. Probab. 2, 969 (1974).
  * Nicoletti  _et al._ (2020) G. Nicoletti, S. Suweis, and A. Maritan,  _Scaling and criticality in a phenomenological renormalization group,_Phys. Rev. Res. 2, 1 (2020).
  * Doi (1976) M. Doi,  _Second quantization representation for classical many-particle system,_ J. Phys. A: Math. Gen. 9, 1465 (1976).
  * Peliti (1985) L. Peliti,  _Path integral approach to birth-death processes on a lattice,_ J. Phys.-Paris 46, 1469 (1985).
  * Täuber (2014) U. C. Täuber,  _Critical Dynamics_ (Cambridge University Press, 2014).
  * Pausch (2019) J. Pausch,  _Topics in Statistical Mechanics_ , Ph.D. thesis, Imperial College London (2019), http://hdl.handle.net/10044/1/73905.
  * Pausch and Pruessner (2019) J. Pausch and G. Pruessner,  _Is actin filament and microtubule growth reaction- or diffusion-limited?_J. Stat. Mech.: Theory Exp. 2019, 1 (2019).
  * Buzsáki and Draguhn (2004) G. Buzsáki and A. Draguhn,  _Neuronal oscillations in cortical networks,_Science 304, 1926 (2004).
  * Pausch  _et al._ (2020) J. Pausch, R. Garcia-Millan, and G. Pruessner, Time-dependent branching processes: a model of oscillating neuronal avalanches, (2020), arXiv:2003.03302.
  * Klimm  _et al._ (2014) F. Klimm, D. S. Basset, J. M. Carlson, and P. J. Mucha,  _Resolving structural variability in network models and the brain,_PLOS Comput. Biol. 10, 1 (2014).
  * Murray  _et al._ (2014) J. D. Murray, A. Bernacchia, D. J. Freedman, R. Romo, J. D. Wallis, X. Cai, C. Padoa-Schioppa, T. Pasternak, H. Seo, D. Lee, and X.-J. Wang,  _A hierarchy of intrinsic timescales across primate cortex,_Nat. Neursci. 17, 1661 (2014).
  * Pausch (2020) J. Pausch, Program for calculating moments of inter-spike intervals induced by pumped branching process, (2020), github.com/JohannesPausch/PumpedBranchingSpikeTimes.
  * Wagenaar  _et al._ (2006) D. A. Wagenaar, J. Pine, and S. M. Potter,  _An extremely rich repertoire of bursting patterns during the development of cortical cultures,_BMC Neurosci. 7 (2006), 10.1186/1471-2202-7-11, protocols were NIH approved.
  * Nakahama  _et al._ (1968) H. Nakahama, H. Suzuki, M. Yamamoto, S. Aikawa, and S. Nichioka,  _Statistcal analysis of spontaneous activity of central single neurons,_Physiol. Behav. 3, 745 (1968).
  * Bassant (1976) M. Bassant,  _Analyse statistique de l’activité des cellules pyramidales de l’hippocampe dorsal du lapin,_Electroencephalogr. Clin. Neurophysiol. 40, 585 (1976).
  * Levina and Priesemann (2017) A. Levina and V. Priesemann,  _Subsampling scaling,_Nat. Commun. 8, 1 (2017).
  * van Kampen (1992) N. van Kampen,  _Stochastic Processes in Physics and Chemistry_ (Elsevier, 1992).


# On Data-Driven Log-Optimal Portfolio:
A Sliding Window Approach

Pei-Ting Wang  Chung-Han Hsieh  Department of Quantitative Finance, National
Tsing Hua University, Hsinchu 30044, Taiwan R.O.C. (e-mail:
<EMAIL_ADDRESS>ch.hsieh@mx.nthu.edu.tw).

###### Abstract

In this paper, we propose a data-driven sliding window approach to solve a
log-optimal portfolio problem. In contrast to many of the existing papers,
this approach leads to a trading strategy with time-varying portfolio weights
rather than fixed constant weights. We show, by conducting various empirical
studies, that the approach possesses a superior trading performance to the
classical log-optimal portfolio in the sense of having a higher cumulative
rate of returns.

###### keywords:

Stochastic Systems, Data-Driven Approach, Financial Engineering, Empirical
Portfolio Optimization

††thanks: This paper is partially supported by the Ministry of Science and
Technology (MOST), Taiwan, R.O.C. under Grant: MOST–110–2222–E–007–005.

and

##  1 Introduction

The take-off point for this paper is the so-called log-optimal portfolio,
which is any portfolio that maximizes the Expected Logarithmic Growth (ELG) of
a trader’s wealth. This ELG maximization idea was originated from Kelly jr
(1956) to solve repeated coin-flipping gambling problems. Since then, various
ramifications and extensions along this line are studied extensively; e.g.,
see Algoet and Cover (1988); Rotando and Thorp (1992); Cover and Thomas
(2006); Thorp (2006). Specifically, the objective for solving the classical
ELG problem is to seek a portfolio with weight $K$ of its account value that
maximizes the ELG at the terminal stage. Some good and bad properties are
studied in MacLean et al. (2010). A rather comprehensive survey of the Kelly-
based approach can be found in MacLean et al. (2011) and the references
therein. While many of the existing papers contributed to Kelly’s problem and
its application to stock trading; e.g., see Rotando and Thorp (1992); Algoet
and Cover (1988); Cover and Thomas (2006); Lo et al. (2018); Kuhn and
Luenberger (2010); Thorp (2006); Hsieh et al. (2016), the resulting optimal
weights $K$ are typically time-invariant.

According to Cover and Thomas (2006); Cornuejols and Tütüncü (2006), it is
known that if the returns are independent and identically distributed (IID)
and known perfectly to the trader, then the constant weight $K$ is optimal.
However, in practice, the returns are typically neither IID nor perfectly
known to the trader; see Fama (2021) and Luenberger (2013). To this end,
various approaches are proposed to solve the log-optimal portfolio problem by
relaxing the IID assumptions; e.g., see Cover (1991) for proposing a universal
portfolio with unknown return distributions and Rujeerapaiboon et al. (2018,
2016); Sun and Boyd (2018) for solving various robust Kelly optimal portfolio
problems. In contrast, this paper, under somewhat weak assumptions on returns,
aimed at studying the data-driven log-optimal portfolio problem via a sliding
window approach. Unlike many of the papers mentioned above that considering a
portfolio with a fixed constant weight; see Algoet and Cover (1988); Thorp
(2006); Nekrasov (2014); Hsieh et al. (2018); Hsieh (2021, 2022), our approach
leads to a time-varying portfolio weight, which generalizes the case with a
fixed constant weight. To close this brief introduction, we also mention some
related work here; e.g., see Park and Irwin (2007) for a good survey of the
profitability of technical analysis, Wu et al. (2022) for an interesting
application of using Kelly-based approach in options trading, Hsieh (2020) for
analyzing log-optimal fraction using Taylor-based approximation approach.

###  1.1 Plans for Section to Follow

In Section 2, we provide some preliminaries and formulate a data-driven log-
optimal portfolio problem in a discrete-time setting. Then, in Section 3, a
sliding window algorithm is provided to solve the portfolio problem.
Subsequently, in Section 4, with the sliding window approach, we provide
various empirical studies using historical stock price data. Lastly, in
Section 5, some concluding remarks and future research directions are
discussed.

##  2 Problem Formulation

Fix an integer $N>1$. For stages $k=0,1,\dots,N-1$, we consider a trader who
is forming a portfolio consisting of $m\geq 2$ assets. If an asset is
riskless, its return is deterministic and is treated as a degenerate random
variable with value $X(k)=r\geq 0$ for all $k$ with probability one.
Alternatively, if Asset $i$ is a stock whose price at time $k$ is
$S_{i}(k)>0$, then its return is

$X_{i}(k)=\frac{S_{i}(k+1)-S_{i}(k)}{S_{i}(k)}.$

In the sequel, we assume that the return vectors

$X(k):=\left[X_{1}(k)\;X_{2}(k)\;\cdots\;X_{m}(k)\right]^{T}.$

The return vector $X(k)$ is drawn according to an unknown but identically
distributed distribution which is supported on $J$ distinct points. For
$k=0,1,\dots,N-1$, the corresponding joint probability mass function can be
estimated by

$P(X_{1}(k)=x_{1}^{j},\dots,X_{m}(k)=x_{m}^{j})=P(X(k)=x^{j})=p_{j}$

for $j=1,2,\dots,J$ with $p_{j}\geq 0$ and $\sum_{j=1}^{J}p_{j}=1.$ We also
assume that the returns satisfy $X_{\min,i}\leq X_{i}(k)\leq X_{\max,i}$ with
known bounds above and with $X_{\max,i}$ being finite and $X_{\min,i}>-1$.

###  2.1 Portfolio Weight and Account Value Dynamics

For $k=0,1,\dots$, let $V(k)$ be the trader’s account value at stage $k$. Now,
for $i=1,2,\dots,m$, take $0\leq K_{i}\leq 1,$ which represents a weight of
the portfolio allocated to the $i$th asset. In the sequel, we require that the
trade is long-only and cash-financed. Specifically, using a shorthand notation
$K:=[K_{1}\;K_{2}\;\cdots\;K_{m}]^{T},$ we require that the weight must
satisfy the classical unit simplex constraint; i.e.,

$K\in\mathcal{K}:=\left\\{K\in\mathbb{R}^{m}:K_{i}\geq 0\text{ for all }\ i,\
\sum_{i=1}^{m}K_{i}=1\right\\}.$

Now, at stage $k=0,$ we begin with the initial account value $V(0):=V_{0}>0$.
The associated account value dynamics can be described as the following
stochastic recursive equation

$\displaystyle V(k+1)$ $\displaystyle=(1+K^{T}X(k))V(k).$

The per-period log-return is defined as

$g(k):=\log\frac{V(k+1)}{V(k)}$

which is used as the objective of the maximization in the sections to follow.

###  2.2 Data-Driven Log-Optimal Portfolio Problem

To incorporate with available data, we modify the problem and consider the
data-driven log-optimal portfolio problem:

$\displaystyle\max_{K\in\mathcal{K}}\ \mathbb{E}\left[g(k)\right]$
$\displaystyle=\max_{K\in\mathcal{K}}\ \mathbb{E}[\ \log(1+K^{T}X(k))\ ]$
$\displaystyle=\max_{K\in\mathcal{K}}\ \sum_{j=1}^{J}p_{j}\log(1+K^{T}x^{j}).$
(1)

The problem above is readily verified as a concave program; i.e., a
maximization problem with a concave objective and a convex constraint set
$\mathcal{K}$; see Boyd and Vandenberghe (2004). In this paper, we will use a
modeling language called CVXPY (using Python language); see Diamond and Boyd
(2016). CVXPY enable us to solve the problem in a very efficient manner.

###### Remark 2.1

To close this section, it is worth mentioning that, unlike many of the
existing papers in solving log-optimal portfolio using fixed constant weight,
our approach seeks to find a portfolio weight $K\in\mathcal{K}$ by solving the
latest optimal portfolio problem with sliding window sizes $M\geq 1$. The
resulting optimal weight $K\in\mathcal{K}$ will be applied to the successive
stage $k+1$. Then we “slide” the window and resolve the problem. The procedure
is repeated until the terminal stage.

##  3 Sliding Window Algorithm

In this section, we are ready to introduce our sliding window approach aimed
at solving the data-driven log-optimal portfolio problem stated in Section 2.
Specifically, for each $k=0,1,\dots$, we rebalance our portfolio using a
sliding window with sizes $M$ days to obtain the latest optimal weight.

Using $s_{i}(k)$ to denote the realized prices for Asset $i$ at stage $k$, the
associated realized returns $x_{i}(k)$ satisfy

$\displaystyle x_{i}(k):=\frac{s_{i}(k+1)-s_{i}(k)}{s_{i}(k)}$ (2)

for $i=1,2,\dots,m$. Let

$x(k):=\left[x_{1}(k)\;x_{2}(k)\;\cdots\;x_{m}(k)\right]^{T}.$

Now sitting at $k=k_{0}$, the sliding window approach is aimed at obtaining an
optimal $K^{*}$ that will be applied to the next stage $k_{0}+1$. To this end,
we proceed as follows. We begin by selecting a collection of realized returns
$\\{x(j)\\}_{j=k_{0}-M}^{k_{0}-1}$, which is obtained in $M$ past stages.
Subsequently, we solve the maximization problem to obtain the $K^{*}$. The
detailed algorithm is summarized below.

Algorithm 1 Data-Driven Sliding Window Algorithm

1:Consider $m\geq 2$ assets, stock prices $\\{s_{i}(j)\\}_{k=0}^{N-1}$ for
$i=1,2,\dots,m$, and sliding window size $M$.

2:Optimal portfolio weight $K^{*}$ for each stage

3:Compute realized returns $\\{x_{i}(k)\\}_{k=0}^{N-1}$ with $x_{i}(k)$
defined in Equation (2).

4:if $k\geq M$ then

5: Choose $\\{x(j)\\}_{j=k-M}^{k-1}$ from $\\{x(k)\\}_{k=0}^{N-1}$ with joint
probability $p_{j}:=1/M$ for all $j=k-M,\dots,k-1$.

6: Solve the maximization problem  $\max_{K\in\mathcal{K}}\
\frac{1}{M}\sum\limits_{j=k-M}^{k-1}\ \log(1+K^{T}x(j)).$

7: Having obtained optimal $K^{*}(k):=K^{*}$, we apply it at stage $k+1$. Set
$k=k+1$ then back to Step 2.

8:end if

##  4 Empirical Studies

To illustrate our sliding window approach, two empirical studies are
considered: The first one is a three-asset portfolio using daily prices and
the second one is the same portfolio using the flipped upside down prices.
Consistent with the existing literature; e.g., Bodie et al. (2018), to compare
the trading performance, some typical metrics are used as follows.

###  4.1 Performance Metrics

For $k=0,1,\dots,N-1$, the first performance metric to be used in the analysis
to follow is the portfolio realized return in period $k$. That is,

$R^{p}(k):=\frac{V(k+1)-V(k)}{V(k)}.$

The (realized) cumulative return up to stage $k=N$ is given by
$(V(N)-V(0))/V(0)$ and the log-growth rate is the logarithm of the realized
cumulative return; i.e., $\log(V(N)/V(0))$. The excess return denoted

$R^{p}(k):=R^{p}(k)-r_{f}$

where $r_{f}$ is the risk-free rate. The realized (per-period) Sharpe ratio,
denoted by $SR$, of the portfolio is the average of the excess returns
$\overline{R^{p}}$ over the standard deviation of the excess returns $\sigma$;
i.e.,

$SR:=\frac{\overline{R^{p}}-r_{f}}{\sigma}$

and the N-period realized Sharpe ratio can be approximated by $\sqrt{N}\cdot
SR$; see Lo (2002). Lastly, other than standard deviation, to scrutinize the
downside risks over multi-period trading performance, we include the maximum
percentage drawdown as an alternative risk metric; i.e.,

$d^{*}:=\max_{0\leq\ell<k\leq N}\frac{V(\ell)-V(k)}{V(\ell)}.$

Lastly, we also report the running times of the sliding window approach. For
example, with $m=3$ assets, we generate $10,000$ of sample paths using the
historical daily price data with window sizes $M=10$. On a laptop with $2.5$
GHz with $4$GB RAM, the data-driven log-optimal problem can be solved about
$3.18$ seconds, which is in the same scale of solving classical log-optimal
portfolio.111The required running times justifies the sliding window approach
with daily price data since the time required to solve the problem is less
than a day.

######  Example 4.1 (Three Assets Portfolio)

Consider a portfolio consisting of three underlying assets: The first asset is
Vanguard Total World Stock (Ticker: VT), the second asset is Vanguard Total
Bond Market ETF (Ticker: BND), and the third asset is Vanguard Total
International Bond ETF (Ticker: BNDX).222It is worth mentioning that these
three assets (VT, BND, and BNDX) forms a well-diversified portfolio. Assuming
$r_{f}:=0$ and a two-year duration from February 14, 2018 to February 14,
2020, we solve the log-optimal portfolio problem to obtain the optimal weight
$K^{*}$. The corresponding prices are shown in Figure 1. Beginning with an
initial account value $V(0)=\$1$, we compare our sliding window approach with
the classical log-optimal portfolio.

Figure 1: Daily Stock Prices for VT, BND, and BNDX

Specifically, to obtain the classical log-optimal portfolio $K^{*}$, we solve
Problem (1) for the first year of the data from August 14, 2018 to February
13, 2019. The corresponding optimal $K^{*}$ is given by
$K^{*}=e_{3}:=[0\;0\;1]^{T},$ which suggests that one should invest all of
available funds on BNDX. The out-of-sample performance metrics under classical
log-optimal portfolio are summarized in Table 1. The corresponding account
value trajectory is shown in Figure 2 with a bold solid line in black color.
From the figure, we see that at the terminal stage, $V(N)\approx 1.0849$ with
a maximum percentage drawdown about $d^{*}=2.21\%$.

Table 1: Performance Metrics under Classical Log-Optimal Portfolio Classical
log-optimal portfolio $K^{*}$
---
Maximum percentage drawdown $d^{*}$ |  $2.21$%
Cumulative rate of return $\frac{V(N)-V_{0}}{V_{0}}$ |  $8.49$%
Realized Log-Growth $\log\frac{V(N)}{V(0)}$ |  $8.15$%
volatility (Annualized) $\sigma$ |  $2.85$%
Sharpe ratio $\sqrt{N}\cdot SR$ | $2.87$
Running Times (secs) | $3.17$

###  4.2 Log-Optimal Portfolio with Sliding Window Approach

With various window sizes $M=5,10,30,60,$ and $100$, we calculate the log-
optimal portfolio weight $K^{*}(k)$ via the sliding window approach mentioned
in Section 3. Having obtained $K^{*}(k)$, we then implement it on the next
trading stage $k+1$; see Figure 2 for an illustration. According to the
figure, we find that a relatively small $M$; e.g., $M=5$ or $M=10$, appears to
lead to a higher cumulative return. Instead of investing on BNDX only as
suggested by the classical log-optimal portfolio with a constant weight, the
sliding window approach yields a time-varying weight; e.g., see Figure 3.

Figure 2: Classical Log-Optimal Portfolio Versus Sliding Window Approach with
Sizes $M=5,10,30,60,$ and $100$ Figure 3: Time-Varying Portfolio Weights
$K(k)$ with size $M=10$ for the Three-Asset Portfolio

Various performance metrics are summarized in Table 2. It is worth noting that
one can view the window size $M$ as a new design variable in the following
sense: One seeks an optimal size $M^{*}$ that gives the largest return or
lowest volatility. In this example, $M=10$ appears to be the best choice in
terms of various performance metrics such as Sharpe ratio.

Table 2: Performance Metrics Under Sliding Window Approach Log-Optimal
Portfolio with Sliding Window Approach
---
Sliding window sizes $M$ | $100$ | $60$ | $30$ | $10$ | $5$
Maximum percentage drawdown $d^{*}$ | $6.39\%$ | $4.99\%$ | $6.44\%$ | $2.35\%$ | $6.27\%$
Cumulative rate of return $\frac{V(N)-V_{0}}{V_{0}}$ | $8.80\%$ | $14.53\%$ | $11.29\%$ | $20.55\%$ | $19.40\%$
Realized Log-Growth $\log\frac{V(N)}{V(0)}$ | $8.43\%$ | $13.56\%$ | $10.70\%$ | $18.69\%$ | $17.73\%$
Volatility (Annualized) $\sigma$ | $8.36\%$ | $7.70\%$ | $8.30\%$ | $6.30\%$ | $7.54\%$
Sharpe ratio (Annualized) $\sqrt{N}\cdot SR$ | $1.050$ | $1.801$ | $1.331$ | $2.998$ | $2.389$
Running Times (secs) | $3.614$ | $6.277$ | $3.325$ | $3.187$ | $2.880$

######  Example 4.2 (A Hypothetical Case)

Here we consider the same setting but with the hypothetical flipped upside-
down prices, which corresponds to a bear market situation.333A bear market
refers to as the market in which prices are falling or are expected to fall.
Using the sliding window algorithm in Section 3 to test for the flipped
upside-down case. Figure 4 shows the hypothetical prices of all these three
assets (VT, BND, BNDX) from February 14, 2018 to February 14, 2020. To obtain
a classical log-optimal portfolio $K^{*}$, we solve the log-optimal portfolio
problem for one year ago from February 14, 2018 to February 13, 2019, and
apply optimal $K^{*}$ into the next year from February 14, 2019 to February
14, 2020. In this case, corresponding log-optimal portfolio $K^{*}$ is given
by $K^{*}=e_{3}:=[0\;0\;1]^{T},$ which suggests one should invest all of
available money on asset BNDX. The corresponding account value trajectory is
shown in Figure 5 in a bold solid line with black color. Other performance
metrics are summarized in Table 3.

Table 3: Performance Metrics of Classical Log-Optimal Portfolio (Hypothetical
Case) Classical log-optimal portfolio $K^{*}$
---
Maximum percentage drawdown $d^{*}$ | $4.81\%$
Cumulative rate of return $\frac{V(N)-V_{0}}{V_{0}}$ | $-4.73\%$
Realized Log-Growth $\log\frac{V(N)}{V(0)}$ | $-4.84\%$
Volatility (Annualized) $\sigma$ | $1.89\%$
Sharpe ratio $\sqrt{N}\cdot SR$ | $-2.56$
Running Times (secs) | $4.36$

Similar to Example 4.1, we again implement Algorithm 1 with various window
sizes $M=5,10,30,60$, and $100$ days. The corresponding account value
trajectories are shown in Figure 5. Consistent with our theory, the resulting
strategy yields a time-varying portfolio weight; see Figure 6 for an example
of this fact with size $M=10$.

Figure 4: Hypothetical Daily Stock Prices for VT, BND, and BNDX Figure 5:
Classical Log-Optimal Portfolio Versus Sliding Window Approach with Sizes
$M=5,10,30,60,$ and $100$ (Hypothetical Case) Figure 6: Time-Varying Portfolio
Weights $K(k)$ with size $M=10$ (Hypothetical Case)

Even if in this (hypothetical) bear market, the sliding window approach
appears to maintain the account value without falling too much. In Table 4, we
see that the annualized rate of return in $M=10$ is even surprisingly
positive.

Table 4: Performance Metrics Under Sliding Window Approach (Hypothetical Case)
Log-Optimal Portfolio with Sliding Window Approach
---
Sliding window sizes $M$ | $100$ | $60$ | $30$ | $10$ | $5$
Maximum percentage drawdown $d^{*}$ | $7.96\%$ | $13.27\%$ | $10.60\%$ | $7.16\%$ | $7.56\%$
Cumulative rate of return $\frac{V(N)-V_{0}}{V_{0}}$ | $-7.08\%$ | $-11.75\%$ | $-5.20\%$ | $0.97\%$ | $-3.64\%$
Realized Log-Growth $\log\frac{V(N)}{V(0)}$ | $-7.34\%$ | $-12.49\%$ | $-5.34\%$ | $0.97\%$ | $-3.71\%$
Volatility (Annualized) $\sigma$ | $8.36\%$ | $8.81\%$ | $10.28\%$ | $9.62\%$ | $9.97\%$
Sharpe ratio (Annualized) $\sqrt{N}\cdot SR$ | $-0.858$ | $-1.375$ | $-0.468$ | $0.148$ | $-0.322$
Running Times (secs) | $4.826$ | $4.087$ | $3.909$ | $3.580$ | $3.431$

##  5 Concluding Remarks

In this paper, we studied a sliding window approach for solving a data-driven
log-optimal portfolio problem. In contrast to the classical log-optimal
portfolio with a constant weight, our approach leads to a time-varying weight,
which generalizes the classical case. We show, by example, that our approach
may be potentially superior to that with classical log-optimal portfolio in
terms of the cumulative account value. To close this brief conclusion section,
we provide two possible research directions.

Optimal Window Sizes $M$.  As mentioned previously in Section 4, the window
sizes $M$ can be viewed as a new design variable. Hence, determining an
“optimal” size would be of interest to pursue further. It is also worth
mentioning that the optimization using a sliding window approach with applying
the resulting optimal weights to the next trading day is indeed closely
related to the idea of model predictive control; e.g., see Mayne et al. (2000)
and Rawlings et al. (2017). Hence, pursuing this direction might be also
fruitful.

Computational Complexity Issues. While we solve a data-driven log-optimal
portfolio problem via a sliding window approach and report corresponding
running times for the algorithm, we did not address the computational
complexity of the proposed algorithm. Hence, it would be of great interest to
analyze the time computational complexities.

## References

  * Algoet and Cover (1988) Algoet, P.H. and Cover, T.M. (1988).  Asymptotic Optimality and Asymptotic Equipartition Properties of Log-Optimum Investment.  _The Annals of Probability_ , 16(2), 876–898.
  * Bodie et al. (2018) Bodie, Z., Kane, A., and Marcus, A. (2018).  _Investments_.  McGraw Hill.
  * Boyd and Vandenberghe (2004) Boyd, S. and Vandenberghe, L. (2004).  _Convex Optimization_.  Cambridge University Press.
  * Cornuejols and Tütüncü (2006) Cornuejols, G. and Tütüncü, R. (2006).  _Optimization Methods in Finance_.  Cambridge University Press.
  * Cover (1991) Cover, T.M. (1991).  Universal Portfolios.  _Mathematical Finance_ , 1(1), 1–29.
  * Cover and Thomas (2006) Cover, T.M. and Thomas, J.A. (2006).  _Elements of Information Theory_.  Wiley-Interscience.
  * Diamond and Boyd (2016) Diamond, S. and Boyd, S. (2016).  CVXPY: A Python-embedded modeling language for convex optimization.  _Journal of Machine Learning Research_ , 17(83), 1–5.
  * Fama (2021) Fama, E.F. (2021).  _Market Efficiency, Long-Term Returns, and Behavioral Finance_.  University of Chicago Press.
  * Hsieh (2020) Hsieh, C.H. (2020).  On Feedback Control in Kelly Betting: An Approximation Approach.  In _2020 IEEE Conference on Control Technology and Applications (CCTA)_ , 903–908.
  * Hsieh (2021) Hsieh, C.H. (2021).  Necessary and Sufficient Conditions for Frequency-Based Kelly Optimal Portfolio.  _IEEE Control Systems Letters_ , 5(1), 349–354.
  * Hsieh (2022) Hsieh, C.H. (2022).  Generalization of Affine Feedback Stock Trading Results to Include Stop-Loss Orders.  _Automatica_ , 136, 110051.
  * Hsieh et al. (2016) Hsieh, C.H., Barmish, B.R., and Gubner, J.A. (2016).  Kelly Betting Can Be Too Conservative.  In _Proceedings of the IEEE Conference on Decision and Control (CDC)_ , 3695–3701.
  * Hsieh et al. (2018) Hsieh, C.H., Gubner, J.A., and Barmish, B.R. (2018).  Rebalancing Frequency Considerations for Kelly-Optimal Stock Portfolios in a Control-Theoretic Framework.  In _Proceedings of the IEEE Conference on Decision and Control (CDC)_ , 5820–5825.
  * Kelly jr (1956) Kelly jr, J. (1956).  A New Interpretation of Information Rate.  _The Bell System Technical Journal_.
  * Kuhn and Luenberger (2010) Kuhn, D. and Luenberger, D.G. (2010).  Analysis of the Rebalancing Frequency in Log-Optimal Portfolio Selection.  _Quantitative Finance_ , 10(2), 221–234.
  * Lo (2002) Lo, A.W. (2002).  The Statistics of Sharpe Ratios.  _Financial Analysts Journal_ , 58(4), 36–52.
  * Lo et al. (2018) Lo, A.W., Orr, H.A., and Zhang, R. (2018).  The Growth of Relative Wealth and The Kelly Criterion.  _Journal of Bioeconomics_ , 20(1), 49–67.
  * Luenberger (2013) Luenberger, D.G. (2013).  _Investment Science_.  Oxford university press.
  * MacLean et al. (2010) MacLean, L.C., Thorp, E.O., and Ziemba, W.T. (2010).  Good and Bad Properties of the Kelly Criterion.  _Risk_ , 20(2), 1.
  * MacLean et al. (2011) MacLean, L.C., Thorp, E.O., and Ziemba, W.T. (2011).  _The Kelly Capital Growth Investment Criterion: Theory and Practice_.  World Scientific.
  * Mayne et al. (2000) Mayne, D.Q., Rawlings, J.B., Rao, C.V., and Scokaert, P.O. (2000).  Constrained model predictive control: Stability and optimality.  _Automatica_ , 36(6), 789–814.
  * Nekrasov (2014) Nekrasov, V. (2014).  Kelly Criterion for Multivariate Portfolios: A Model-Free Approach.  _Available at SSRN 2259133_.
  * Park and Irwin (2007) Park, C.H. and Irwin, S.H. (2007).  What do we know about the profitability of technical analysis?  _Journal of Economic surveys_ , 21(4), 786–826.
  * Rawlings et al. (2017) Rawlings, J.B., Mayne, D.Q., and Diehl, M. (2017).  _Model Predictive Control: Theory, Computation, and Design_ , volume 2.  Nob Hill Publishing.
  * Rotando and Thorp (1992) Rotando, L.M. and Thorp, E.O. (1992).  The Kelly Criterion and the Stock Market.  _The American Mathematical Monthly_ , 99(10), 922–931.
  * Rujeerapaiboon et al. (2018) Rujeerapaiboon, N., Barmish, B.R., and Kuhn, D. (2018).  On risk reduction in kelly betting using the conservative expected value.  In _Proceedings of the IEEE conference on decision and control (CDC)_ , 5801–5806.
  * Rujeerapaiboon et al. (2016) Rujeerapaiboon, N., Kuhn, D., and Wiesemann, W. (2016).  Robust Growth-Optimal Portfolios.  _Management Science_ , 62(7), 2090–2109.
  * Sun and Boyd (2018) Sun, Q. and Boyd, S. (2018).  Distributional Robust Kelly Gambling.  _arXiv preprint arXiv:1812.10371_.
  * Thorp (2006) Thorp, E.O. (2006).  The Kelly Criterion in Blackjack Sports Betting, and The Stock Market.  _Handbook of Asset and Liability Management: Theory and Methodology_ , 1, 385.
  * Wu et al. (2022) Wu, M.E., Syu, J.H., and Chen, C.M. (2022).  Kelly-Based Options Trading Strategies on Settlement Date via Supervised Learning Algorithms.  _Computational Economics_.