Maximum Entropy approach to multivariate time series randomization

TL;DR

This paper introduces a data-driven, maximum entropy-based framework for hypothesis testing in complex multivariate time series, addressing challenges of non-stationarity and ergodicity in natural and social systems.

Contribution

It develops an unsupervised, statistical mechanical approach for creating ensembles of time series that preserve key empirical properties, enabling more robust hypothesis testing.

Findings

Framework successfully applied to financial portfolio analysis

Enables hypothesis testing without stationarity assumptions

Provides a new tool for analyzing complex interacting systems

Abstract

Natural and social multivariate systems are commonly studied through sets of simultaneous and time-spaced measurements of the observables that drive their dynamics, i.e., through sets of time series. Typically, this is done via hypothesis testing: the statistical properties of the empirical time series are tested against those expected under a suitable null hypothesis. This is a very challenging task in complex interacting systems, where statistical stability is often poor due to lack of stationarity and ergodicity. Here, we describe an unsupervised, data-driven framework to perform hypothesis testing in such situations. This consists of a statistical mechanical approach - analogous to the configuration model for networked systems - for ensembles of time series designed to preserve, on average, some of the statistical properties observed on an empirical set of time series. We showcase its possible applications with a case study on financial portfolio selection.

Figures

Tables

Table 1. Table 1: Fraction of empirical moments compatible with their corresponding ensemble distribution at different significance levels specified in terms of quantiles (e.g., 0.01-0.99 denotes that the 1st and 99th percentiles of the ensemble distribution are used as bounds to determine whether the null hypothesis of an empirical moment being compatible with the ensemble distribution can be rejected or not). Note that the confidence intervals used to obtain these results have not been adjusted for multiple hypothesis testing. Doing so (e.g., via False Coverage Rate  [ 31 ] ) would further suppress the number of true positives, resulting an even larger fraction of moments being compatible with the ensemble distribution. Moments are calculated both for each stock and each trading day. In the last column, we also report, for each moment, the median relative error between the empirical value and its ensemble average.

Returns		Significance null hypothesis			median rel. err.
Stat	Sample	0.01-0.99	0.05-0.95	0.1-0.9
Var	stock	0.95	0.76	0.59	0.2
	day	0.88	0.78	0.69	0.14
Skew	stock	1	0.98	0.95	0.13
	day	0.78	0.58	0.49	0.46
Kurt	stock	0.78	0.61	0.51	0.60
	day	0.85	0.68	0.55	0.1

Table 2. Table 2: Fraction of empirical return distributions (both for stocks and trading days) that are compatible with their ensemble counterparts based on Kolmogorov-Smirnov tests at different significance levels.

Ratio empirical aggregated pdfs not rejected by a K-S test	Aggregation level	K-S test significance
		0.01	0.05
	stocks	0.92	0.68
	days	0.82	0.75

Table 3. Table 3: Out-of-sample portfolio risk – quantified in terms of variance – with and without detrending the returns by subtracting their ensemble average. P 1 , 2 N superscript subscript 𝑃 1 2 𝑁 P_{1,2}^{N} (with N = 20 , 50 𝑁 20 50 N=20,50 ) refer to two different portfolios made of randomly selected S&P stocks, whereas q = N / T 𝑞 𝑁 𝑇 q=N/T denotes the portfolios’ “rectangularity ratio” (i.e., the ratio between the number of stocks and the length of the in-sample time window used to compute correlations and portfolio weights). The two top rows refer to portfolios whose weights are computed based on the raw returns, whereas the two bottom rows refer to portfolios whose weights are computed based on the detrended returns. In the latter case, the detrending is only performed in-sample to compute correlations and weights, and the out-of-sample risk is computed by retaining such weights on new raw returns. The numbers reported in each case refer to the average out-of-sample risk computed over a set of 30-days long non-overlapping time windows spanning the period September 2014 - November 2018.

	$P_1^{20}$	$P_2^{20}$	$P_1^{50}$	$P_2^{50}$
$q=2/3$	0.041 (0.027 , 0.091)	0.955 (0.026 , 0.815)	0.1029 (0.011 , 0.287)	0.1553 (0.015 , 0.461)
$q=1/4$	0.8488 (0.031 , 2.867)	1.001 (0.022, 3.011)	0.7846 (0.009 , 0.933)	0.0938 (0.009 , 0.136)
$q=2/3$	0.0093 (0.0053 , 0.0155)	0.0081 (0.0046 , 0.0124)	0.0034 (0.0021 , 0.0056)	0.0033 (0.0023 , 0.0053)
$q=1/4$	0.0113 (0.0055 , 0.0158)	0.0081 (0.0053 , 0.0111)	0.0041 (0.0022 - 0.0056)	0.0033 (0.0021 , 0.0054)

	$P_1^{20}$	$P_2^{20}$	$P_1^{50}$	$P_2^{50}$
$q=2/3$	0.032 (-0.19 , 0.31)	0.009 (-0.22 , 0.23)	-0.077 (-0.25 , 0.15)	-0.022 (-0.22 , 0.28)
$q=1/4$	-0.013 (-0.21 , 0.12)	-0.022 (-0.26, 0.23)	0.011 (-0.19 , 0.20)	-0.037 (-0.27 , 0.15)
$q=2/3$	0.042 (-0.14 , 0.22)	0.041 (-0.16 , 0.22)	0.081 (-0.17 , 0.34)	0.054 (-0.17, 0.25)
$q=1/4$	0.035 (-0.23 , 0.26)	0.035 (-0.17 , 0.24)	0.073 (-0.16 , 0.32)	0.062 (-0.13 , 0.34)

Table 5. Table 4: Number of tests passed by out-of-sample VaR estimates for the BNP stock at different significance levels α 𝛼 \alpha .

$\alpha$	Tests passed
	$M_3$	$M_2$	$M_1$
90%	6	4	4
95%	8	6	5
99%	8	6	6
99.99%	8	8	8

Table 6. Table 5: Number of tests passed by out-of-sample VaR estimates for the S&P Index at different significance levels α 𝛼 \alpha .

$\alpha$	Tests passed
	$M_3$	$M_2$	$M_1$
90%	5	4	4
95%	8	7	6
99%	8	7	7
99.99%	8	8	7

Equations

$$P(W)=\frac{e^{-H(W)}}{Z},$$

$$H(W)=\sum _{i=1}^d\sum _{t=1}^T\left[a_i+W_t{\alpha }_i+W_t^2{\beta }_i\right]\Theta ({\overline{W}}_t-{\overline{q}}_{{\xi }_{i-1}})\Theta (-{\overline{W}}_t+{\overline{q}}_{{\xi }_i}).$$

$$H_1=\sum _i\left[{\alpha }_i{N}_{{\xi }_i}+{\beta }_i{M}_{{\xi }_i}\right]$$

$$H_2=\sum _i\left[{\alpha }_i{N}_{{\xi }_i}+\beta M_{{\xi }_i}+\gamma M_{{\xi }_i}^2\right]$$

$$\sum _{W\in \mathcal{W}}\equiv \prod _{t=1}^T{\int }_{{\overline{q}}_{{\xi }_1}}^{{\overline{q}}_{{\xi }_d}}d{w}_t.$$

$$Z_1=\prod _{t=1}^T\sum _{i=1}^{d-1}{e}^{-{\alpha }_i}\frac{e^{-{\beta }_i{\overline{q}}_{{\xi }_i}}-e^{-{\beta }_i{\overline{q}}_{{\xi }_{i+1}}}}{{\beta }_i}$$

$$Z_2=\prod _{t=1}^T\sum _{i=1}^{d-1}\sqrt{\frac{\pi }{4\gamma }}{e}^{\frac{{\beta }^2}{4\gamma }-{\alpha }_i}\left(-\text{erf}\left[\frac{\beta +2\gamma {\overline{q}}_{{\xi }_i}}{2\sqrt{\gamma }}\right]+\text{erf}\left[\frac{\beta +2\gamma {\overline{q}}_{{\xi }_{i+1}}}{2\sqrt{\gamma }}\right]\right)$$

$$H(W)=\sum_{i=1}^{N}\sum_{t=1}^{T}\bigg{[}\left(\alpha^{N}_{i}+\alpha^{T}_{t}\right)A^{+}_{it}+\left(\beta^{N}_{i}+\beta^{T}_{t}\right)A^{-}_{it}+\left(\gamma_{i}^{N}+\gamma_{t}^{T}\right)w^{+}_{it}+\left(\sigma_{i}^{N}+\sigma_{t}^{T}\right)w^{-}_{it}\bigg{]}\ ,$$

$$\sum_{W\in\mathcal{W}}\equiv\prod_{i=1}^{N}\prod_{t=1}^{T}\int_{-\infty}^{+\infty}dw_{it}=\prod_{i=1}^{N}\prod_{t=1}^{T}\sum_{\vbox{\Let@\restore@math@cr\default@tag\halign{\hfil$\m@th\scriptstyle#$&$\m@th\scriptstyle{}#$\cr&(0,1)\\ (A_{it}^{+},A_{it}^{-})=&(1,0)\\ &(0,0)\crcr}}}\;\int_{0}^{+\infty}dw^{+}_{it}\;\int_{0}^{+\infty}dw^{-}_{it}\ ,$$

$$Z=\sum _{W\in \mathcal{W}}{e}^{-H(W)}=\prod _{i=1}^N\prod _{t=1}^T{Z}_{it}=\prod _{i=1}^N\prod _{t=1}^T\left[1+\frac{e^{-({\alpha }_i^N+{\alpha }_t^T)}}{{\gamma }_i^N+{\gamma }_t^T}+\frac{e^{-({\beta }_i^N+{\beta }_t^T)}}{{\sigma }_i^N+{\sigma }_t^T}\right]=\prod _{i=1}^N\prod _{t=1}^T\left(1+e^{\frac{{\mu }_{it}^1-{ϵ}_{it}}{T_{it}}}+e^{\frac{{\mu }_{it}^2-{ϵ}_{it}}{T_{it}}}\right),$$

$$P(W)=\prod _{i=1}^N\prod _{t=1}^T{\left[P_{it}^{+}\right]}^{A_{it}^{+}}{\left[P_{it}^{-}\right]}^{A_{it}^{-}}{\left[1-P_{it}^{+}-P_{it}^{-}\right]}^{1-A_{it}^{+}-A_{it}^{-}}{\left[Q_{it}^{+}(w_{it}^{+})\right]}^{A_{it}^{+}}{\left[Q_{it}^{-}(w_{it}^{-})\right]}^{A_{it}^{-}},$$

$$P(W_{it}=x)=(1-P_{it}^{+}){\lambda }_{it}^{-}{e}^{{\lambda }_{it}^{-}x}\Theta (-x)+P_{it}^{+}{\lambda }_{it}^{+}{e}^{-{\lambda }_{it}^{+}x}\Theta (x),$$

$$S(W)=-\sum _{W\in \mathcal{W}}P(W)\ln P(W),$$

$$⟨{\mathcal{O}}_{\ell }(W)⟩=\sum _{W\in \mathcal{W}}{\mathcal{O}}_{\ell }(W)P(W)={\mathcal{O}}_{\ell }(\overline{W})={\overline{O}}_{\ell },$$

$$\sum _{W\in \mathcal{W}}P(W)=1.$$

$$\frac{\partial }{\partial P}\left[S+\alpha \left(1-\sum _{W\in \mathcal{W}}P(W)\right)+\sum _{\ell =1}^L{\beta }_{\ell }\left(O_{\ell }-\sum _{W\in \mathcal{W}}{\mathcal{O}}_{\ell }(W)P(W)\right)\right]=0,$$

$$P(W)=\frac{e^{-H(W)}}{Z}.$$

$$\sum_{W\in\mathcal{W}}\equiv\prod_{i=1}^{N}\prod_{t=1}^{T}\sum_{\vbox{\Let@\restore@math@cr\default@tag\halign{\hfil$\m@th\scriptstyle#$&$\m@th\scriptstyle{}#$\cr&(0,1)\\ (A_{it}^{+},A_{it}^{-})=&(1,0)\\ &(0,0)\crcr}}}\;\int_{0}^{+\infty}dw^{+}_{it}\;\int_{0}^{+\infty}dw^{-}_{it}\ .$$

$$Z=$$

$$=$$

$$T_{ij}$$

$${ϵ}_{ij}$$

$${\mu }_{ij}^2$$

$$Z=\prod _{i,t=1}^{N,T}{Z}_{it}=\prod _{i,t=1}^{N,T}\left[\frac{e^{-({\alpha }_i^N+{\alpha }_t^T)}}{{\gamma }_i^N+{\gamma }_t^T}+\frac{1}{{\sigma }_i^N+{\sigma }_t^T}\right].$$

$$P$$

$$H(W)=\sum _{i=1}^N\sum _{t=1}^T\left[\left({\alpha }_i^N+{\alpha }_t^T\right)A_{it}^{+}+\left({\gamma }_i^N+{\gamma }_t^T\right)w_{it}^{+}+\left({\sigma }_i^N+{\sigma }_t^T\right)w_{it}^{-}\right],$$

$$Z=\prod _{i=1}^N\prod _{t=1}^T{Z}_{it}=\prod _{i=1}^N\prod _{t=1}^T\left[\frac{e^{-({\alpha }_i^N+{\alpha }_t^T)}}{{\gamma }_i^N+{\gamma }_t^T}+\frac{1}{{\sigma }_i^N+{\sigma }_t^T}\right].$$

$$\underset{\mathbf{\pi }}{\min }\sum _{i,j=1}^N{C}_{ij}{\pi }_i{\pi }_j$$

$$\sum _{i=1}^N{\pi }_i{\mu }_i=\mu ;\sum _{i=1}^N{\pi }_i=1,$$

$${\pi }_i(\mu )=\sum _{j=1}^N{C}_{ij}^{-1}(\ell (\mu )+\mu g(\mu )),$$

$$\ell (\mu )=\frac{c-b\mu }{ac-b^2},g(\mu )=\frac{a\mu -b}{ac-b^2}$$

Full text

1 Introduction

Hypothesis testing lies at the very core of the scientific method. In its general formulation, it hinges upon contrasting the observed statistical properties of a system with those expected under a null hypothesis. In particular, hypothesis testing allows to discard potential models of a system when empirical measurement that would be exceedingly unlikely under them are made.

However, there is often no theory to guide the investigation of a system’s dynamics. What is worse, in many practical situations one may be given a single - and possibly unreproducible - set of experimental data. This is indeed the case when dealing with most complex systems, whose collective dynamics often are markedly non-stationary, ranging from climate [1, 2] to brain activity [3] and financial markets [4, 5, 6]. This, in turn, makes hypothesis testing in complex systems a very challenging task, that potentially prevents from assessing which properties observed in a given data sample are “untypical”, i.e, unlikely to be observed again in a sample collected at a different point in time.

This issue is usually tackled by constructing ensembles of artificial time series sharing some characteristics with those generated by the dynamics of the system under study. This can be done either via modelling or in a purely data-driven way. In the latter case, the technique most frequently used by both researchers and practitioners is bootstrapping [7, 8], which amounts to generating partially randomised versions of the available data via resampling, that can then be used as a null benchmark to perform hypothesis testing. Depending on its specificities, bootstrapping can account for autocorrelations and cross-correlations in time series sampled from multivariate systems. However, it relies on assumptions, such as sample independence and some form of stationarity [9], which limit its power when dealing with complex systems.

As far as model-driven approaches are concerned, the literature is extremely vast [10]. Broadly speaking, modelling approaches rely on a priori structural assumptions for the system’s dynamics, and on identifying the parameter values that best explain the available set of observations within a certain class of models (e.g., via Maximum Likelihood [10]). A widely used class for multivariate time series is that of autoregressive models, such as VAR [11], ARMA [12], and GARCH [13], which indeed were originally introduced, among other reasons, to perform hypothesis testing [12]. In such models, the future values of each time series are given by a linear combination of past values of one or more time series, each characterised by their own idiosyncratic noise to capture the fluctuations of individual variables. Such a structure is most often dictated by its simplicity rather than by first principles. As a consequence, once calibrated, autoregressive models produce rather constrained ensembles of time series that do not allow to explore scenarios that differ substantially from those observed empirically.

Another modelling philosophy places more emphasis on capturing the structural collective properties of multivariate systems rather than their dynamical ones. Random Matrix models are a prime example in this direction, which usually rely on ansatzes on the correlation structure of the system under study trying to strike a balance between the resulting models’ analytical tractability and their adherence to empirical observations. One of the first - and still most widely used - Random Matrix models is the Wishart Ensemble [14, 15], which in its simplest form leads to the much celebrated Marčenko-Pastur distribution [16] for uncorrelated systems, up to rather recent developments to tackle non-stationarities in financial data [17].

Here we propose a Maximum Entropy approach - inspired by Statistical Mechanics - to perform hypothesis testing on sets of time series. Starting from the Maximum Entropy principle, we will introduce a (gran)canonical ensemble of correlated time series subject to constraints based on the properties of an empirically observed set of measurements. This, in turn, will result in a multivariate probability distribution which allows to unbiasedly sample values centred on such measurements, which represents the main contribution of this paper. The theory we propose in the following shares some similarities with the canonical ensemble of complex networks [18, 19, 20, 21], and, as we will show, inherits its powerful calibration method based on Likelihood maximization [22].

The paper is organized as follows. In the next Section, we outline the general formalism of our approach. Then, as a formative example, we show how the methodology introduced can be used to reconstruct an unknown probability density function from repeated measurements over time. After that, we proceed to study the most general case of multivariate time series, showing how our approach recovers collective statistical properties of interacting systems without directly accounting for such interactions in the set of constraints imposed on the ensemble. Before concluding with some final remarks, we briefly mention an interesting analogy between our approach and Jaynes’ Maximum Caliber principle [23].

General Framework Description

Let $\mathcal{W}$ be the set of all real-valued sets of $N$ time series of length $T$ (i.e., the set of real-valued matrices of size $N\times T$), and let $\overline{W}\in\mathcal{W}$ be the empirical set of data we want to perform hypothesis testing on (i.e., $\overline{W}_{it}$ stores the value of the $i$-th variable in the system sampled at time $t$, so that $\overline{W}_{it}$ for $t=1,\ldots,T$ represents the sampled time series of variable $i$). Our aim is to define a probability density function $P(W)$ on $\mathcal{W}$ such that the expectation values $\langle\mathcal{O}_{\ell}(W)\rangle$ of a set of observables ($\ell=1,\ldots,L$) coincide with the value $\overline{O}_{\ell}$ of the corresponding quantity as empirically measured from $\overline{W}$.

Following Boltzmann and Gibbs, we can achieve the above by invoking the Maximum Entropy principle, i.e., by identifying $P(W)$ as the distribution that maximises the entropy functional $S(W)=-\sum_{W\in\mathcal{W}}P(W)\,\ln P(W)$, while satisfying the $L$ constraints $\langle\mathcal{O}_{\ell}(W)\rangle=\sum_{W}\mathcal{O}_{\ell}(W)P(W)=\overline{O}_{\ell}$ and the normalisation condition $\sum_{W}P(W)=1$. As is well known [24, 25], this reads

[TABLE]

where $H(W)=\sum_{\ell}\beta_{\ell}\ \mathcal{O}_{\ell}(W)$ is the Hamiltonian of the system, $\beta_{\ell}$ ($\ell=1,\ldots,L$) are Lagrange multipliers introduced to enforce the constraints, and $Z=\sum_{W}e^{-H(W)}$ is the partition function of the ensemble, which verifies $\langle\mathcal{O}_{\ell}(W)\rangle=\partial\ln Z/\partial\beta_{\ell}\;,\forall\,\ell$.

Figure 1 provides a sketch representation of the ensemble theory just introduced. The rationale of enforcing the aforementioned constraints is that of finding a distribution $P(W)$ that assigns low probability to regions of the phase space $\mathcal{W}$ where the observables associated to the Lagrange multipliers $\beta_{\ell}$ take values that are exceedingly different from those measured in the empirical set $\overline{W}$, and high probability to regions where some degree of similarity with $\overline{W}$ is retained (it should be noted that in some cases this does not necessarily leads to the distribution $P(W)$ being peaked around the values $\overline{O}_{\ell}$). This, in turn, allows to test whether other properties (not encoded in any of the constraints) of $\overline{W}$ are statistically significant by measuring how often they appear in instances drawn from the ensemble. The existence and uniqueness of the Lagrange multipliers ensuring the ensemble’s ability to preserve the constraints $\overline{O}_{\ell}$ is a well known result, and they are equivalent to those that would be obtained from the maximization of the Likelihood of drawing the empirical matrix $\overline{W}$ from the ensemble [26].

In the two following Sections, we shall illustrate our general framework on two examples – one devoted to the single time series case, one to the multivariate case. In both examples, we shall make a selection of possible constraints that can be analytically captured by the approach, i.e., constraints for which the resulting ensemble’s partition function can be computed in closed form. In particular, since we will apply our approach to financial data in a later section, we will choose constraints that have a clear interpretation in the analysis of financial time series. However, it should be kept in mind that such constraints are by no means to be interpreted as general prescriptions, and all results presented in the following could be reobtained – depending on the applications of interest – with any other set of constraints allowing for analytical solutions.

Single time series

As a warm up example to showcase our approach, we shall consider a simple case of a univariate and stationary system with no correlations over time. This amounts to a time series made of independent and identically distributed random draws from a probability density function, which we aim at reconstructing. In the following section we will then proceed to consider multivariate and correlated cases.

Let us then consider an $1\times T$ empirical data matrix $\overline{W}$ coming from $T$ repeated sampling of an observable of the system under consideration. If the processes is stationary and time-independent, this is equivalent to sampling $T$ times a random variable from its given, unknown, distribution and therefore the task of the model can be translated into reconstructing the unknown distribution given the data. Let us consider a vector $\xi\in[0,1]^{d}$ and the associated empirical $\xi$-quantiles $\overline{q}_{\xi}$ calculated on $\overline{W}$. In order to fully capture the information present in the data $\overline{W}$, we are going to constrain our ensemble to preserve, as averages, one or more quantities derived from $\overline{q}_{\xi}$. Possible choices may be:

• •

The number of data points falling within each pair of empirically observed adjacent quantiles:

$\overline{N}_{\xi_{i}}=\sum_{t}\Theta(\overline{W}_{t}-\overline{q}_{\xi_{i-1}})\;\Theta(-\overline{W}_{t}+\overline{q}_{\xi_{i}})$

• •

The cumulative values of the data points falling within each pair of adjacent quantiles:

$\overline{M}_{\xi_{i}}=\sum_{t}\overline{W}_{t}\;\Theta(\overline{W}_{t}-\overline{q}_{\xi_{i-1}})\;\Theta(-\overline{W}_{t}+\overline{q}_{\xi_{i}})$

• •

The cumulative squared values of the data points falling within each pair of adjacent quantiles:

$\overline{M}_{\xi_{i}}^{2}=\sum_{t}\overline{W}_{t}^{2}\;\Theta(\overline{W}_{t}-\overline{q}_{\xi_{i-1}})\;\Theta(-\overline{W}_{t}+\overline{q}_{\xi_{i}})$

In each of the above constraints we assumed $i=2,\ldots,d$, and we have used $\Theta(\cdot)$ to indicate Heaviside’s step function (i.e., $\Theta(x)=1$ for $x>0$, and $\Theta(x)=0$ otherwise). In general we are not required to use the same $\alpha$ for all constraints, for example we can freely choose to impose on the ensemble the ability to preserve $\overline{N}_{\xi_{i}}\forall i\in[1,d]$ together with the total cumulative values $\overline{M}=\sum_{i}\overline{M}_{\xi_{i}}$, and total cumulative squared values $\overline{M}^{2}=\sum_{i}\overline{M}_{\xi_{i}}^{2}$, as well as each $\overline{M}_{\xi_{i}}$ and $\overline{M}_{\xi_{i}}^{2}$ separately. Note that the first constraint in the above list effectively amounts to constraining the ensemble’s quantiles.

As we will discuss more extensively later on, the above set of constraints is discretionary. A possible strategy to get rid of such discretionality, would be to partition the data based on a binning procedure aimed at compromising between resolution and relevance, according to its definition introduced in Ref. [27].

A defined set of constraints will lead to a different Hamiltonian, to a different number of Lagrange multipliers and therefore to a different statistical model. If we choose, for example, to adopt all the constraints specified above, the Hamiltonian $H$ of the ensemble will depend on a total of $3(d-1)$ parameters:

[TABLE]

The freedom to choose the amount of constraints of course comes with a cost. First of all, it must be noted that the Likelihood of the data matrix $\overline{W}$ will be in general a non linear function of the Lagrange multipliers and therefore of the constraints. These latter can vary both in magnitude (by choosing different values for the entries of $\alpha$) and in size (by choosing a different $d$). In general, tackling the issue of finding the optimal positions for the constraints, given their number $d$, can become highly not trivial and goes out of the scope of the present work. However, loosely speaking, the Likelihood of finding $\overline{W}$ after a random draw from the defined ensemble is an increasing function of the number constraints, coherently with the idea that a larger number of parameter leads to better statistics on the data used to train the model. As a result, in order to avoid overfitting, given a set of constraints, we can compare different values of $d$ by using standard model selection techniques such as the Bayesian [28] or Akaike information criteria [29].

In the following we are going to show how to apply the methodology just outlined to a synthetic dataset. For this example, let us assume that the data generating process follows a balanced mixture of a truncated standard Normal distribution and a truncated Student’s t-distribution with $\nu=5$ degree of freedom. The two models we are going to use to build the respective ensembles are specified by the following Hamiltonians:

[TABLE]

The model resulting from $H_{1}$ will have a total of $2(d-1)$ parameter and will preserve the average number of data point contained within each pair of adjacent quantiles together with their cumulative values, while the model coming from $H_{2}$ will be characterised by $d+1$ parameters and will preserve the average number of data point contained within each pair of adjacent quantiles together with the overall mean and variance calculated across all data points. In order to find the Lagrange multipliers able to preserve the chosen constraints, we first need to the find the partitions functions of the two ensembles $Z_{1,2}=\sum_{W}e^{-H_{1,2}(W)}$. In order to do that, we first need to specify the sum over the phase space:

[TABLE]

The above expression leads to the following partition functions (in Appendix A we present a detailed derivation of the partition function shown in the following Section; the following result can be obtained with similar steps):

[TABLE]

where with erf we indicate the Gaussian error function $\text{erf}(z)=\frac{2}{\pi}\int_{0}^{z}e^{-t^{2}}dt$.

In Figure 2 we show how the models resulting from the partition functions $Z_{1}$ and $Z_{2}$ are able to reconstruct the underlying true distribution starting from different amount of information (i.e different sample sizes) and quantiles vector set to $\overline{q}=[-\infty,\overline{q}_{0.25},\overline{q}_{0.5},\overline{q}_{0.75},\infty]$ (the steps to obtain the ensemble’s probability density function from its partition function are outlined in Appendix A for the case discussed in the next Section). First of all we note that, as expected, estimating the unknown distribution from more data gives estimates that are closer to the real underlying distribution. Moreover, looking at Figure 2, we can qualitatively see that the model described by $Z_{1}$ does a better job than $Z_{2}$ at inferring the unknown pdf. We can verify both statements more quantitatively by calculating the the Kullback–Leibler divergence of the estimated distributions from the true one: for the case with 40 data points we observe $D_{KL}(P_{Z_{1}}|P_{T})=0.10$ and $D_{KL}(P_{Z_{2}}|P_{T})=0.19$ while for the case with 4000 samples we have $D_{KL}(P_{Z_{1}}|P_{T})=0.01$ and $D_{KL}(P_{Z_{2}}|P_{T})=0.08$. Of course, we cannot conclude yet that $Z_{1}$ gives overall a better model for our reconstruction task than $Z_{2}$ since they are described by a different number of parameters. As mentioned above, in order to complete our model comparison exercise, we need to rely on a test able to assess the relative quality of the models for a given set of data. We choose the Akaike information criterion which uses as its score function the following $\text{AIC}=2k-2\log\hat{L}$, where $k$ is the number of estimated parameters and $\hat{L}$ is the maximum value of the likelihood function for the model. We end up with $\text{AIC}_{Z_{1}}=130$ and $\text{AIC}_{Z_{2}}=950$ for the 40 data points case and $\text{AIC}_{Z_{1}}=1,15\times 10^{4}$ and $\text{AIC}_{Z_{2}}=2.11\times 10^{5}$ when 4000 data points are available. Of course it is worth repeating that all the steps mentioned above are performed given a fixed vector $\overline{q}$ common to the two models.

\floatsetup

[figure]style=plain,subcapbesideposition=top

Multiple time series

In this Section we proceed to present the application of the approach introduced above to the multivariate case.

Let us consider an $N\times T$ empirical data matrix $\overline{W}$ whose rows have been rescaled to have zero mean, so that $\overline{W}_{it}>0$ ($\overline{W}_{it}<0$) will indicate that the time $t$ value of the $i$-th variable is higher (lower) than its empirical mean. Also, without loss of generality, let us assume that $\overline{W}_{it}\in\mathbb{R}_{\neq 0}$, and that $\overline{W}_{it}=0$ indicates missing data. For later convenience, let us define $A^{\pm}=\Theta(\pm W)$ and $w^{\pm}=\pm W\Theta(\pm W)$ (we shall denote the corresponding quantities measured on the empirical set as $\overline{A}^{\pm}$ and $\overline{w}^{\pm}$), and let us constrain the ensemble to preserve the values of the following observables:

• •

The number of positive (above-average) and negative (below-average) values $\overline{N}_{i}^{\pm}=\sum_{t}\overline{A}^{\pm}_{it}$, and the number of missing values $\overline{N}_{i}^{0}=T-\overline{N}_{i}^{+}-\overline{N}_{i}^{-}$ recorded for each time series ($i=1,\ldots,N$).

• •

The cumulative positive and negative values recorded for each time series: $\overline{S}_{i}^{\pm}=\sum_{t}\overline{w}_{it}^{\pm}$ ($i=1,\ldots,N$).

• •

The number of positive, negative, and missing values recorded at each sampling time: $\overline{M}_{t}^{\pm}=\sum_{i}\overline{A}^{\pm}_{it}$, $\overline{M}_{t}^{0}=N-\overline{M}_{t}^{+}-\overline{M}_{t}^{-}$ ($t=1,\ldots,T$).

• •

The cumulative positive and negative value recorded at each sampling time: $\overline{R}_{t}^{\pm}=\sum_{i}\overline{w}^{\pm}_{it}$ ($t=1,\ldots,T$).

Note that the second constraint in the above list indirectly constrains the mean of each time series. As mentioned in the General Framework section, we selected the above constraints inspired by potential financial applications (and indeed we will assess the ensemble’s ability to capture time series behaviour on financial data). In such context, the above four constraints respectively correspond to: the number of positive and negative returns of a given financial stock, the total positive and negative return of a stock, the number of stocks with a positive or negative return on a given trading day, the total positive and negative return across all stocks on a given trading day. Such constraints amount to some of the most fundamental “observables” associated with financial returns. As we will detail in the following, forcing the ensemble to preserve them on average also amounts to effectively preserving other quantities that are of paramount importance in financial analysis, such as, e.g., the skewness and kurtosis of return distributions, and some of the correlation properties of a set of financial stocks (which are central to financial portfolio analysis and selection.)

The above list amounts to $8(N+T)$ constraints, and the Hamiltonian $H$ depends on the very same number of parameters:

[TABLE]

where we have introduced the Lagrange multipliers associated to all constraints. This choice for the Hamiltonian naturally generalizes the framework introduced in [30].

Let us remark that none of the above constraints explicitly accounts for either cross-correlations between variables or for correlations in time. Accounting for these would amount to constraining products of the type $\sum_{t=1}^{T}w_{it}w_{jt}$ (in the case of cross-correlations) and $\sum_{i=1}^{N}w_{it}w_{it^{\prime}}$ (in the case of temporal correlations), which introduce a direct coupling between the entries of $W$, resulting in a considerable loss in terms of analytical tractability of the approach. However, as we shall see in a moment, the combination of the above constraints is enough to indirectly capture some of the correlation properties in the data of interest.

In order to calculate the partition function $Z=\sum_{W}e^{-H(W)}$, we first need to properly specify the sum over the phase space. Given the matrix representation we have chosen for the system, and the fact that $w_{it}^{\pm}=w_{it}A^{\pm}_{it}$, this reads:

[TABLE]

where the sum specifies whether the entry $A_{it}$ stores a positive, negative or missing value, respectively. This signifies that negative and positive events (this in general holds for any discretization of the distribution of the entries of $W$), i.e., values above and below the empirical mean of each variable, obviously cannot coexist in an entry $W_{it}$, which, once occupied, cannot hold any other event. In this respect, we can anticipate that negative and positive events will effectively be treated as different fermionic species populating the energy levels of a physical system. Following this line of reasoning, the role of $w_{it}^{\pm}$ is that of general coordinates for each of the two fermionic species. In principle, the integrals in Eq. (8) could have as upper limits some quantities $U^{\pm}_{it}$ to incorporate any possible prior knowledge on the bounds of the variables of interest.

The above expression leads to the following partition function (fully derived in Appendix A):

[TABLE]

where the quantities $\mu_{it}^{1,2}$, $\epsilon_{it}$, and $T_{it}$ are functions of the Lagrange multipliers (specified in Appendix A).

Some considerations about Eq. (9) are now in order. First of all, the partition function factorises into the product of independent factors $Z_{it}$, and therefore into a collection of $N\times T$ statistically independent sub-systems. However, it is crucial to notice that their parameters (i.e., the Lagrange multipliers) are coupled through the system of equations specifying the constraints ($\langle\mathcal{O}_{\ell}(W)\rangle=\partial\ln Z/\partial\beta_{\ell}\;,\forall\,\ell$). As we shall demonstrate later, this ensures that part of the original system’s correlation structure is retained within the ensemble. Moreover, with the above positions, the aforementioned physical analogy becomes clear: the system described by Eq. (9) can be interpreted as a system of $N\times T$ orbitals with energies $\epsilon_{it}$ and local temperatures $T_{it}$ that can be populated by fermions belonging to two different species characterised by local chemical potentials $\mu_{it}^{1}$ and $\mu_{it}^{2}$, respectively.

From the partition function in Eq. (9) we can finally calculate the probability distribution $P(W)$:

[TABLE]

where $P_{it}^{\pm}$ and $Q_{it}^{\pm}(w^{\pm}_{it})$ are functions of the Lagrange multipliers (specified in Appendix A) and correspond, respectively, to the probability of drawing a positive (negative) value for the $i$-th variable at time $t$ and to its probability distribution.

\floatsetup

[figure]style=plain,subcapbesideposition=top

\floatsetup

[figure]style=plain,subcapbesideposition=top

As an example application of the ensemble defined above, let us consider the daily returns of the $N=100$ most capitalized NYSE stocks over $T=560$ trading days (spanning October 2016 - November 2018). In this example, the aforementioned constraints force the ensemble to preserve, on average, the number of positive and negative returns and the overall positive and negative return for each time series and for each trading day, leading to $6(N\times T)$ constraints. When these constraints are enforced, an explicit expression for the marginal distributions can be obtained (see Appendix A):

[TABLE]

where $\lambda_{it}^{\pm}$ are also functions of the Lagrange multipliers (specified in Appendix A). The above distribution allows both to efficiently sample the ensemble numerically and to obtain analytical results for several observables. Remarkably, it has been shown [32] that sampling from a mixture-like density such as the one in Eq. (11) can result in heavy tailed distribution, which is of crucial importance when dealing with financial data.

Figure 3 and Tables 1 and 2 illustrate how the above first-moment constraints translate into explanatory power of higher-order statistical properties. Indeed, in the large majority of cases, the empirical return distributions of individual stocks and trading days and their higher-order moments (variance, skewness, and kurtosis) are statistically compatible with the corresponding ensemble distributions, i.e., with the distributions of such quantities computed over large numbers ($10^{6}$ in all cases shown) of time series independently generated from the ensemble. Notably, this is the case without constraints explicitly aimed at enforcing such level of agreement. This, in turn, further confirms that the ensemble can indeed be exploited to perform reliable hypothesis testing by sampling random scenarios that are however closely based on the empirically available data.

In that spirit, in Figure 4 we show an example of ex-post anomaly detection, where the original time series of a stock is plotted against the 95% confidence intervals obtained from the ensemble for each data point $W_{it}$. As it can be seen, the results are non-trivial, since the returns flagged as anomalous are not necessarily the largest ones in absolute value. This is because the constraints imposed on the ensemble reflect the collective nature of financial market movements, thus resulting in the statistical validation of events that are anomalous with respect to the overall heterogeneity present in the market.

Following the above line of reasoning, in Figure 4 we show a comparison between the eigenvalue spectrum of the empirical correlation matrix of the data, and the average eigenvalue spectrum of the ensemble. As is well known, the correlation matrix spectrum of most complex interacting systems normally features a large bulk of small eigenvalues which is often approximated by the Marchenko-Pastur (MP) distribution [33] of Random Matrix Theory (i.e., the average eigenvalue spectrum of the correlation matrix of a large system of uncorrelated variables with finite second moments) [34, 35, 15], plus a few large and isolated eigenvalues that contain information about the relevant correlation structure of the system (e.g., they can be associated to clusters of strongly correlated variables [36]). As it can be seen in the Figure, the ensemble’s average eigenvalue spectrum qualitatively captures the same range of the empirical spectral bulk (for reference, we also plot the MP distribution), and the ensemble distribution for the largest eigenvalue is very close to the one empirically observed, demonstrating that the main source of correlation in the market is well captured by the ensemble. Conversely, the average distance between the empirically observed largest eigenvalue and its ensemble distribution can be interpreted as the portion of the market’s collective movement which cannot be explained by the constraints imposed on the ensemble.

In Appendix B, we also apply the above ensemble approach to a dataset of weekly and hourly temperature time series recorded in North-American cities. We do this to showcase the approach’s ability to capture inherent time periodicities in empirical data – which would be very hard to capture directly – by means of the constraints already considered in the examples above.

Applications to financial risk management

In this section we push the examples of the previous section – where we demonstrated the ensemble’s ability to partially capture the collective nature of fluctuations in multivariate systems – towards real-world applications. Namely, we will illustrate a case study devoted to financial portfolio selection.

Financial portfolio selection is an optimization problem which entails allocating a fixed amount of capital across $N$ financial stocks. Typically, the goal of an investor is to allocate their capital in order to maximise the portfolio’s expected return while minimising the portfolio’s expected risk. When the latter is quantified in terms of portfolio variance, the solution to the optimization problem amounts to computing portfolio weights $\pi_{i}$ ($i=1,\ldots,N$), where $\pi_{i}$ is the amount of capital to be invested in stock $i$ (note that $\pi_{i}$ can be negative when short selling is allowed). As is well known, these are functions of the portfolio’s correlation matrix [37], which reflects the intuitive notion that a well balanced portfolio should be well diversified, avoiding similar allocations of capital in stocks that are strongly correlated. The mathematical details of the problem and explicit expressions for the portfolio weights are provided in Appendix C.

The fundamental challenge posed by portfolio optimization is that portfolio weights have to be first computed “in-sample” and then retained “out-of-sample”. In practice, this means that portfolio weights are always computed based on the correlations between stocks observed over a certain period of time, after which one observes the portfolio’s realized risk based on such weights. This poses a problem, as financial correlations are known to be noisy[35, 34] and heavily non-stationary[6], so there is no guarantee that portfolio weights that are optimal in-sample will perform well in terms of out-of-sample risk.

A number of solutions have been put forward in the literature to mitigate the above problem. Most of these amount to methods to “clean” portfolio correlation matrices[38], i.e., procedures aimed at subtracting noise and unearthing the “true” correlations between stocks (at least over time windows where they can be reasonably assumed to be constant). Here, we propose to exploit our ensemble approach in order to apply the same philosophy directly on financial returns rather than on their correlation structure.

Using the same notation as in the previous section, let us assume that $\overline{W}_{it}$ represents the time-$t$ return of stock $i$ ($i=1,\ldots,N$; $t=1,\ldots,T$). Let us then define detrended returns $\tilde{W}_{it}=\overline{W}_{it}-\langle W_{it}\rangle$, where $\langle W_{it}\rangle$ denotes the ensemble average of the return computed from Eq. (11). The rationale for this is to mitigate the impact of returns that may be anomalously large (in absolute value), i.e., returns whose values are markedly distant from their typical values observed in the ensemble (as in the example shown in Fig. 4).

In Table 3 we show the results obtained by performing portfolio selection on the aforementioned detrended returns, and compare them to those obtained without detrending. Namely, we form four portfolios (two of size $N=20$ and two of size $N=50$) with the returns of randomly selected S&P500 stocks in the period from September 2014 to October 2018, and compute their weights based on the correlations computed over non-overlapping period of lengths $T=N/q$, where $q\in(0,1)$ is the portfolio’s “rectangularity ratio”, which provides a reasonable proxy for the noisiness of the portfolio’s correlation matrix (which indeed becomes singular for $q\rightarrow 1$). The numbers in the Table represent the average and $90\%$ confidence level intervals of the out-of-sample portfolio risk – quantified in terms of variance – computed over a number of non-overlapping time windows of $30$ days within in the aforementioned period (see caption for more details), with the first two rows corresponding to the raw returns and the two bottom rows corresponding to the detrended returns (note that in both cases out-of-sample risk is still computed on the raw returns, i.e., detrending is only performed in-sample to compute the weights). As it can be seen, detrending by locally removing the ensemble average from each return reduces out-of-sample risk dramatically, despite the examples considered here being plagued by a number of well-known potential downsides, such as small portfolio size, small time windows, and exposure to outliers (the returns used here are well fitted by power law distributions – using the method in [39] – whose median tail exponent across all stocks is $\alpha=3.9$). In Appendix D we report the equivalent of Table 3 for the portfolios’ Sharpe ratios (i.e., the ratio between portfolio returns and portfolio variance over a time window), with qualitatively very similar results.

In Appendix E, we provide details of an additional application of our ensemble approach to financial risk management, where we compute and test the out-of-sample performance of estimates of Value-at-Risk (VaR) – the most widely used financial risk measure [40] – based on our ensemble approach. This application is specifically aimed at demonstrating that – despite the large number of Lagrange multipliers necessary to calibrate the more constrained versions of our ensembles – the approach does not suffer from overfitting issues. Quite to the contrary – in line with the literature on configuration models for networked systems [19, 20], which are a fairly close relative of the approach proposed here – we find the out-of-sample performance of our method to improve even when increasing the number of Lagrange multipliers quite substantially. The reason for this lies in the fact that in classic cases of overfitting one completely suppresses any in-sample variance of the model being used (e.g., when fitting $n$ points with a polynomial of order $n-1$). This is not the case, instead, with the model at hand. Indeed, being based on maximum entropy, our approach still allows for substantial in-sample variance even when building highly constrained ensembles. We illustrate this in the aforementioned example in Appendix E by obtaining progressively better out-of-sample VaR performance when increasing the number of constraints (and therefore or Lagrange multipliers) the ensemble is subjected to.

Relationship with Maximum Caliber principle

Before concluding, let us point out an interesting connection between our approach and Jaynes’ Maximum Caliber principle [23]. It has recently been shown [41] that the time-dependent probability distribution that maximizes the caliber of a two-state system evolving in discrete time can be calculated by mapping the time domain of the system as a spatial dimension of an Ising-like model. This is exactly equivalent to our mapping of a time-dependent system onto a data matrix, where the system’s time dimension is mapped onto a discrete spatial dimension of the lattice representing the matrix.

From this perspective, our ensemble approach represents a novel way to calculate and maximize the caliber of systems sampled in discrete time with a continuous number of states. This also allows to interpret some recently published results on correlation matrices in a different light. Indeed, in [21] the authors obtain a probability distribution on the data matrix of sampled multivariate systems starting from a Maximum Entropy ensemble on their corresponding correlation matrices. Following the steps outlined in our paper, the same results could be achieved via the Maximum Caliber principle by first mapping the time dimension of the system onto a spatial dimension of a corresponding lattice, and by then imposing the proper constraints on it.

Discussion

In this paper we have put forward a novel formalism - grounded in the ensemble theory of statistical mechanics - to perform hypothesis testing on time series data. Whereas in physics and in the natural sciences, hypothesis testing is carried out through repeated controlled experiments, this is rarely the case in complex interacting systems, where the lack of statistical stability and controllability often hamper the reproducibility of experimental results. This, in turn, prevents from assessing whether the observations made are consistent with a given hypothesis on the dynamics of the system under study.

The framework introduced here tackles the above issues by means of a data-driven and assumptions-free paradigm, which entails the generation of ensembles of randomized counterparts of a given data sample. The only guiding principle underpinning such a paradigm is that of entropy maximization, which allows to interpret the ensemble’s partition function in terms of a precise physical analogy. Indeed, as we have shown, in our framework events in a data sample correspond to fermionic particles in a physical system with multiple energy levels. In this respect, our approach markedly differs from other known methods to generate synthetic data, such as bootstrap. Notably, even though the Hamiltonians used throughout the paper correspond to non-interacting systems, and therefore the correlations in the original data are not captured in terms of interactions between particles (as is instead the case in Ising-like models), the ensemble introduced here is still capable of partially capturing properties typical of interacting systems through the ‘environment’ the particles are embedded in, i.e., a system of coupled local temperatures and chemical potentials.

All in all, our framework is rather flexible, and can be easily adapted to the data at hand by removing or adding constraints from the ensemble’s Hamiltonian. From this perspective, the number and type of constraints can be used to “interpolate” between very different applications. In fact, loosely constrained ensembles can serve as highly randomized counterparts of an empirical dataset of interest, and therefore can be used for statistical validation purposes, i.e., to determine which statistical properties of the empirical data can be explained away with a few basic constraints. The opposite situation is instead represented by a heavily constrained ensemble, designed to capture most of the statistical properties of the empirical data (the examples detailed in previous sections go in this direction). From a practical standpoint, this application of our approach can be particularly useful in situations where sensitive data cannot be shared between parties (e.g., due to privacy restrictions), and where sharing synthetic data whose statistical properties closely match those of the empirical data can be a very valuable alternative.

For example, the constraints in the applications showcased here (i.e., on the sums of above and below average values) result in two fermionic species of particles. More stringent constraints (e.g., on the data belonging to certain percentiles of the empirical distribution) would result in other species being added to the ensemble.

As we have shown, our framework is capable of capturing several non-trivial statistical properties of empirical data that are not necessarily associated with the constraints imposed on the ensemble. As such, it can provide valuable insight on a variety of complex systems by both allowing to test theoretical models for their structure and by allowing to uncover new information in the statistical properties that are not fully captured by the ensemble. We have illustrated some of these aspects with a financial case study, where we demonstrated that a detrending of stock returns based on our ensemble approach yields a dramatic reduction in out-of-sample portfolio risk.

Once more, let us stress that the main limitation of our approach is that of not accounting explicitly for cross-correlations or temporal correlations. As mentioned above, these could in principle be tackled by following the same analytical framework presented here, but would result in a considerably less tractable model. We aim to explicitly deal with the case of temporal correlations in future work.

Data availability

The financial data employed in this paper are freely available to download form Yahoo Finance.

Acknowledgements

G.L. acknowledges support from an EPSRC Early Career Fellowship in Digital Economy (Grant No. EP/N006062/1).

Author contributions statement

Authors contribution: R.M. and G.L. designed research; R.M. performed research and analyzed data; R.M. and G.L. wrote the paper.

Additional information

Competing Interests: The authors declare no conflict of interest.

Appendix A Explicit calculation of the partition function

We want to find a probability density function $P(W)$ on $\mathcal{W}$ such that the expectation values of a set of observables coincide with their empirical value, i.e $\langle\mathcal{O}_{\ell}(W)\rangle=\overline{O}_{\ell}$ ($\ell=1,\ldots,L$), where $\overline{W}\in\mathcal{W}$ is the empirical set of measurements. At first, this problem may appear almost impossible to solve, given that $P(W)$ may be determined by a number of degrees of freedom way larger than the number of constraints we are imposing. However, as introduced in the main text, this can be done by using the Maximum Entropy principle, or, in other words, by adding another (functional) constraint on the probability distribution, which requires that $P(W)$ should also maximise the Gibbs entropy:

[TABLE]

while preserving the constraints:

[TABLE]

and the normalization:

[TABLE]

Eqs. (12)-(14) define a constrained optimization problem, whose solution is found by solving the following equation:

[TABLE]

where, as usual in such scenarios, each constraint has been coupled with a Lagrange multiplier $\alpha,\beta_{\ell}$ ($\ell=1,\ldots,L$). Defining $H(W)=\sum_{\ell}\beta_{\ell}\mathcal{O}_{\ell}(W)$ as the Hamiltonian of the ensemble and $Z=e^{\alpha+1}=\sum_{W}e^{-H(W)}$ its partition function, the solution of Eq. (15) reads:

[TABLE]

This is the general probability density function ruling the ensemble theory we are proposing. Of course, the sum $\sum_{W}$ on the phase space of the system used in the above equations still needs to be properly specified.

We are going to do so while considering the Hamiltonian specified in the main paper. As pointed out above, in order to find the partition function $Z$ of the system, we just need to sum $e^{-H(W)}$ over all possible configurations, i.e., over the set of all the $N\times T$ real valued matrices $\mathcal{W}$. Recalling the notations introduced in the main text $A^{\pm}=\Theta(\pm W)$ and $w^{\pm}=\pm W\Theta(\pm W)$, we can write the sum over the phase space as follows:

[TABLE]

We can now calculate the partition function $Z$ of the ensemble:

[TABLE]

where all the Lagrange multipliers must be positive and we have defined the following quantities in order to make apparent the analogy with the two species fermionic gas introduced in the main text:

[TABLE]

From the above partition function, via Eq. (16) we can derive the probability density function in Eq. (5) of the main text, which quantifies the probability of drawing a specific instance $W$ from the ensemble. The quantities defining such probability distribution have a well defined physical meaning, and read as follows:

$P_{it}^{+}=\frac{e^{-\left(\alpha^{N}_{i}+\alpha^{T}_{t}\right)}}{\left(\gamma_{i}^{N}+\gamma_{t}^{T}\right)Z_{it}}$

Probability of observing a positive value in the $i$-th time series at time $t$

$P_{it}^{-}=\frac{e^{-\left(\beta^{N}_{i}+\beta^{T}_{t}\right)}}{\left(\sigma_{i}^{N}+\sigma_{t}^{T}\right)Z_{it}}$

Probability of observing a negative value in the $i$-th time series at time $t$

$1-P_{it}^{+}-P_{it}^{-}$

Probability of observing a missing value in the $i$-th time series at time $t$

$Q^{+}_{it}(w)=(\gamma_{i}^{N}+\gamma_{t}^{T})e^{-(\gamma_{i}^{N}+\gamma_{t}^{T})w}$

Probability distribution of a positive value $w$ for the $i$-th time series at time $t$

$Q^{-}_{it}(w)=(\sigma_{i}^{N}+\sigma_{t}^{T})e^{-(\sigma_{i}^{N}+\sigma_{t}^{T})w}$

Probability distribution of a negative value $w$ for the $i$-th time series at time $t$

When no data are missing, i.e. $(A_{it}^{+},A_{it}^{-})\neq(0,0)$, the sum defined in Eq. (17) changes and, as a result, the partition function (18) becomes :

[TABLE]

After noticing that $A_{it}^{+}=0\Rightarrow w_{it}^{+}=0\,\wedge\,w_{it}^{-}>0$, the probability of drawing from the ensemble an instance $W$ can be easily found:

[TABLE]

where the quantities in the above expression are defined as those above.

Looking at Eq. (19), we can understand how we have obtained Eq. (6) in the main text. In order to simulate a drawing of a set of time series $W$ from the ensemble, we first need to construct a “topology” of positive events by placing a positive event in entry $W_{it}$ with probability $P_{it}^{+}$ and a negative event otherwise. Then we need to place a weight $W_{it}=x$ using one of the two exponential distributions $Q_{it}^{\pm}$ defined above, depending on the type of event that was assigned to $W_{it}$. This procedure is encompassed by the hyperexponential distribution in Eq. (6) of the main text, which can be obtained via the standard generating function approach, and whose parameters read $\lambda_{it}^{+}=(\gamma_{i}^{N}+\gamma_{i}^{T})^{-1}$, and $\lambda_{it}^{-}=(\sigma_{i}^{N}+\sigma_{i}^{T})^{-1}$.

Appendix B Application to a set of temperature time series

\floatsetup

[figure]style=plain,subcapbesideposition=top

We now apply the framework introduced in the main text to sets of time series featuring temperatures recorded at different frequencies (week/day/8 hours) in $N=30$ different North American cities 111Vancouver, Portland, San Francisco, Seattle, Los Angeles, San Diego, Las Vegas, Phoenix, Albuquerque, Denver, San Antonio, Dallas, Houston, Kansas City, Minneapolis, Saint Louis, Chicago, Nashville, Indianapolis, Atlanta, Detroit, Jacksonville, Charlotte, Miami, Pittsburgh, Toronto, Philadelphia, New York, Montreal, Boston (weekly data range from July 2013 to July 2018, daily data range from July 2016 to July 2018, 8 hour data range from January 2017 to July 2018). We do so in order to test the ability of our ensemble approach to capture the main features of time series whose most relevant statistical properties are markedly different from those of financial returns, which we instead studied in the main paper. In particular, our main focus will be on the ability of the ensemble to capture the periodicities that characterize temperature data at different time scales.

As done in the main text, we indicate as $\overline{W}$ the $N\times T$ data matrix (with $T=264,730,2321$ in the case of temperatures recorded at the weekly, daily, and 8 hour frequency, respectively) with values rescaled to have zero mean, and we indicate as $W$ any generic instance drawn from the corresponding ensemble. We also redefine here for convenience the matrices $A^{\pm}=\Theta(\pm W)\,,\,w^{\pm}=\pm W\Theta(\pm W)$. The ensemble we are going to use is fully specified by the $6(N+T)$ constraints enforced in the following Hamiltonian (there are no missing data, which leads to $2(N+T)$ fewer constraints with respect to the general formulation outlined in the main paper):

[TABLE]

leading to the partition function:

[TABLE]

In Figure 5 we show that, independently from the frequency at which temperatures are sampled, the average ensemble power spectral density captures well the relevant frequencies that characterize the empirical time series of each city. Indeed, as can be seen from panels a and b, the ensemble power spectra based on the data recorded at the weekly and daily frequency perfectly capture the six-months periodicity associated with the seasons’ cycle. Panel c shows that the same frequency is also captured in the data recorded every 8 hours, and that when calibrating the ensemble on such data, the power spectrum also perfectly captures the daily frequency associated with the day-night cycle (see inset).

In Fig. 6 we expand the above analysis to the periodicities of moments. Panel a shows the empirical daily variance of temperatures recorded across the 30 cities mentioned above against the corresponding ensemble average. At first sight, the latter seems to be largely uncorrelated from the former. Yet, the corresponding power spectrum shown in panel b highlights that the relevant frequencies in the data (six months and one day) are captured very well, although the ensemble places additional power on such frequencies.

\floatsetup

[figure]style=plain,subcapbesideposition=top

A somewhat similar phenomenon is shown in panels c and d, which show the daily skewness computed across all cities and its corresponding power spectra. Once again, the average ensemble spectrum places more power on the six-months and daily frequencies with respect to the empirical one. This results in a clearly discernible oscillating pattern, which significantly deviates from the empirical behavior. Nevertheless, these results are interesting. Indeed, as it can be seen in panel c positive (negative) skewness values take place during the summer (winter) months, as a reflection of higher (lower) average temperatures. Although this is a fairly trivial example, it highlights how the ensemble approach can reveal stylized trends that are genuinely informative about the dynamics of the system under study.

Appendix C Optimal portfolio selection

Let us consider a matrix $r_{it}$ ($i=1,\ldots,N$, $t=1,\ldots,T$) of daily financial returns, and let $C$ be the corresponding correlation matrix (i.e., $C_{ij}$ denotes the Pearson correlation coefficient between stocks $i$ and $j$, computed over the time period $[1,T]$). The optimal portfolio problem then amounts to solving the following optimization problem for a vector $\bm{\pi}=(\pi_{1},\ldots,\pi_{N})\in\mathbb{R}^{N}$:

[TABLE]

subject to

[TABLE]

where Eq. (22) expresses the minimization of the portfolio variance, while the equations in (23) express constraints on the expected returns ($\mu_{i}$ denotes the expected out-of-sample return of stock $i$, while $\mu$ denotes the portfolio’s desired expected return) and on the available capital (conventionally set to one unit). The expected returns $\mu_{i}$ are computed – as is often customary – based on mean reversion, i.e., assuming that the return on day $t+1$ will be minus the return on day $t$.

With the above positions, the solution to the optimization problem reads (see Ref. [37] in the main paper)

[TABLE]

where

[TABLE]

Appendix D Sharpe ratios of portfolios in Table 3

In the Table below we report the Sharpe ratios (defined as the ration between the portfolio return and portfolio variance over a time interval) for the same portfolios considered in Table 1 of the main paper.

As it can be seen, applying the “cleaning” procedure outlined in the main paper has beneficial effects also in this case, leading to increased average Sharpe ratios and reduced uncertainty around them across board.

Appendix E Testing for overfitting with an application on the estimation of Value-at-Risk

In this section we expand on the financial application of our ensemble approach discussed in the main paper, with a specific focus on potential overfitting issues. Namely, the number of Lagrange multipliers the ensemble depends on increase linearly with the number of constraints one wants to enforce. For instance, the multivariate case detailed in the main paper (Eq. (7)) depends on $8(N+T)$, which, for small numbers of variables $N$ and small sample sizes $T$, can be of the same order of magnitude of the number of data points ($N\times T$) used to calibrate the ensemble. This, in turn, may raise concerns about potential overfitting issues.

We tackle the above issue by showcasing the approach’s performance when using it to compute out-of-sample Value-at-Risk (VaR) estimates based on increasingly constrained (and therefore parametrized) versions of the ensemble. VaR is the most widely adopted measure of financial risk. For a given significance level $\alpha$ it simply refers to the $1-\alpha$ quantile of the return distribution of a certain financial stock or portfolio. The simplest procedure to estimate VaR is via historical estimation, which amounts to computing the in-sample $1-\alpha$ quantile of a financial time series of interest. However, due to non-stationarities, historical estimates are known to typically be unreliable out-of-sample, and there is a vast literature devoted to enhancing historical estimates with Monte Carlo simulations and other techniques to generate synthetic scenarios. In the following, we will apply a similar line of reasoning by computing VaR estimates as quantiles of the distributions generated by our ensemble starting from an empirical time series of interest, and we will assess their out-of-sample performance based on a number of standard statistical tests.

In the following, we consider two financial time series of length $T=1000$ and $T=1500$ days corresponding, respectively, to BNP returns from May 7, 2008 to March 13, 2013, and to S&P Index returns from May 7, 2008 to July 26, 2014. For each time series we proceed to compute VaR estimates with a rolling window approach. Namely, we compute an in-sample VaR estimate over a time window $[t_{0},t_{0}+\tau]$, with $\tau=150$ days, and assess its out of sample performance on day $t_{0}+\tau+1$. We do this based on the following three versions of our ensemble approach:

• Model $M1$

– This corresponds to a loosely constrained ensemble based on the single time series case discussed in the paper, where we only constrain the ensemble to preserve the empirical time series’ variance and the cumulative values of the data falling within each pair of adjacent quartiles (denoted as $\overline{M}_{\xi_{i}}$ in the main paper, with $\xi_{i}=0.25,0.5,0.75$). Overall, these correspond to 4 constraints and Lagrange multipliers.

• Model $M2$

– This corresponds to a deliberately highly parametrized model based on an adaptation of the multiple time series case. Namely, let us consider the 150 returns of interest to compute a new risk estimate and let us denote them as $r_{1},\ldots,r_{150}$. We then form a $25\times 126$ (which roughly amount to the length of a trading month and half of a trading year, respectively) matrix with such returns with the following circulant structure

[TABLE]

The quantity $\epsilon$ in the upper-right entry of the matrix denotes the unknown out-of-sample return on day 151. We assume as possible values for it $\epsilon=\pm\min|r_{t}|$, and then generate the corresponding ensemble constraining it to preserve the cumulative positive and negative values for each row and column (denoted respectively as $\overline{S}_{i}^{\pm}$ and $\overline{R}_{t}^{\pm}$ in the main paper), which correspond to $2(25+126)=302$ constraints and Lagrange multipliers222It can be shown that as long as the matrix $R$’s sizes $L_{1}$ and $L_{2}$ are not multiple of each other, then such constraints are all linearly independent. In the case of linear dependence, the effective number of constraints decreases by at most $\max(L_{1},L_{2})$, which still amounts to an over-parametrized model.. We generate the ensembles for both aforementioned values of $\epsilon$ and combine the two resulting distributions for them in order to compute a VaR estimate for the return on day 151.

• Model $M3$

– The same as model $M2$ with additional constraints on the number of positive and negative returns recorded in each column (the equivalent of the quantity denoted as $\overline{M}_{t}^{\pm}$ in the main paper). In addition to the constraints mentioned above, this gives a total of $3(25+126)=453$ constraints and Lagrange multipliers.

Models $M2$ and $M3$ are highly constrained (and therefore highly parametrized) ones, as they force the corresponding ensemble to preserve a very large number of local properties of the time series.

We calibrate the above models on all time windows of length $\tau=150$ days starting on days $t_{0}=1,2,\ldots,T-151$ and compute out-of-sample VaR estimates for each of them, resulting in $849$ estimates for BNP and $1349$ estimates for the S&P Index, respectively. We then pool such estimates for both time series and assess their out-of-sample performance by means of 8 standard tests widely adopted in the financial literature. These are the traffic light, binomial, proportion of failures, time until first failure, conditional coverage, conditional coverage independence, time between failures, and time between failures independence tests (see, e.g., Nieppola, O., Backtesting Value-at-Risk Models for their definitions).

The results – reported as the number of tests passed – are shown in Table 4 for BNP and in Table 5 for the S&P Index, for varying significance levels $\alpha$. As it can be seen, the out-of-sample performance systematically improves when increasing the number of constraints, regardless of the significance level, even when pushing these to numbers close to the number of available data points. Remarkably, all tests are passed when using model $M3$ at significance $95\%$ or higher.