Recent Trends and Developments in Econophysics

Argyrakis, Panos (Ed.) Book Recent Trends and Developments in Econophysics Provided in Cooperation with: MDPI – Multidisciplinary Digital Publishing Institute, Basel Suggested Citation: Argyrakis, Panos (Ed.) (2024) : Recent Trends and Developments in Econophysics, ISBN 9783725804153, MDPI - Multidisciplinary Digital Publishing Institute, Basel, https://doi.org/10.3390/books978-3-7258-0416-0 This Version is available at: https://hdl.handle.net/10419/312690 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by-nc-nd/4.0/ mdpi.com/journal/entropy Special Issue Reprint Recent Trends and Developments in Econophysics Edited by Panos Argyrakis Recent Trends and Developments in Econophysics Recent Trends and Developments in Econophysics Editor Panos Argyrakis Basel •Beijing •Wuhan •Barcelona •Belgrade •Novi Sad •Cluj •Manchester Editor Panos Argyrakis Aristotle University of Thessaloniki Thessaloniki Greece Editorial Office MDPI St. Alban-Anlage 66 4052 Basel, Switzerland This is a reprint of articles from the Special Issue published online in the open access journal Entropy (ISSN 1099-4300) (available at: https://www.mdpi.com/journal/entropy/special issues/entropy econophys). For citation purposes, cite each article independently as indicated on the article page online and as indicated below: Lastname, A.A.; Lastname, B.B. Article Title. Journal Name Year,Volume Number, Page Range. ISBN 978-3-7258-0415-3 (Hbk) ISBN 978-3-7258-0416-0 (PDF) doi.org/10.3390/books978-3-7258-0416-0 © 2024 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license. The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons Attribution-NonCommercial-NoDerivs (CC BY-NC-ND) license. Contents Antonio Briola and Tomaso Aste Dependency Structures in Cryptocurrency Market from High to Low Frequency Reprinted from: Entropy 2022,24, 1548, doi:10.3390/e24111548 .................... 1 Charalampos M. Liapis, Aikaterini Karanikola and Sotiris Kotsiantis Investigating Deep Stock Market Forecasting with Sentiment Analysis Reprinted from: Entropy 2023,25, 219, doi:10.3390/e25020219 ..................... 19 Takuya Wada, Hideki Takayasu and Misako Takayasu Extraction of Important Factors in a High-Dimensional Data Space: An Application for High-Growth Firms Reprinted from: Entropy 2023,25, 488, doi:10.3390/e25030488 ..................... 49 Leticia P´erez-Sienes, Mar Grande, Juan Carlos Losada and Javier Borondo The Hurst Exponent as an Indicator to Anticipate Agricultural Commodity Prices Reprinted from: Entropy 2023,25, 579, doi:10.3390/e25040579 ..................... 72 Hongli Niu, Kunliang Xu and Mengyuan Xiong The Risk Contagion between Chinese and Mature Stock Markets: Evidence from a Markov-Switching Mixed-Clayton Copula Model Reprinted from: Entropy 2023,25, 619, doi:10.3390/e25040619 ..................... 83 Joe Scattergood and Steven Bishop A Network Model Approach to International Aid Reprinted from: Entropy 2023,25, 641, doi:10.3390/e25040641 .....................103 Tamara Kyrylych and Yuriy Povstenko Multi-Criteria Analysis of Startup Investment Alternatives Using the Hierarchy Method Reprinted from: Entropy 2023,25, 723, doi:10.3390/e25050723 .....................121 Frank Brennan Webb, Daniel Stimpson, Miesha Purcell and Eduardo L´opez Organizational Labor Flow Networks and Career Forecasting Reprinted from: Entropy 2023,25, 784, doi:10.3390/e25050784 .....................131 Ewa A. Drzazga-Szcz¸e´sniak, Piotr Szczepanik, Adam Zenon Kaczmarek and Dominik Szcz¸e´sniak Entropy of Financial Time Series Due to the Shock of War Reprinted from: Entropy 2023,25, 823, doi:10.3390/e25050823 .....................150 Nedim Bayrakdar, Valerio Gemmetto and Diego Garlaschelli Local Phase Transitions in a Model of Multiplex Networks with Heterogeneous Degrees and Inter-Layer Coupling Reprinted from: Entropy 2023,25, 828, doi:10.3390/e25050828 .....................162 Peter Tsung-Wen Yen, Kelin Xia and Siew Ann Cheong Laplacian Spectra of Persistent Structures in Taiwan, Singapore, and US Stock Markets Reprinted from: Entropy 2023,25, 846, doi:10.3390/e25060846 .....................197 Francisco Y´a˜nez Rodr´ıguez and Alberto P. Mu˜nuzuri A Goodwin Model Modification and Its Interactions in Complex Networks Reprinted from: Entropy 2023,25, 894, doi:10.3390/e25060894 .....................228 v Hua Zhong, Xiaohao Liang and Yougui Wang Transaction Entropy: An Alternative Metric of Market Performance Reprinted from: Entropy 2023,25, 1140, doi:10.3390/e25081140 ....................245 vi Citation: Briola, A.; Aste, T. Dependency Structures in Cryptocurrency Market from High to Low Frequency. Entropy 2022,24, 1548. https://doi.org/10.3390/ e24111548 Academic Editor: Panos Argyrakis Received: 21 September 2022 Accepted: 26 October 2022 Published: 28 October 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). entropy Article Dependency Structures in Cryptocurrency Market from High to Low Frequency Antonio Briola 1,2 and Tomaso Aste 1,2,3,* 1Department of Computer Science, University College London, London WC1E 6BT, UK 2Center for Blockchain Technologies, University College London, London WC1E 6BT, UK 3Systemic Risk Center, London School of Economics, London WC2A 2AE, UK *Correspondence: [email protected] Abstract: We investigate logarithmic price returns cross-correlations at different time horizons for a set of 25 liquid cryptocurrencies traded on the FTX digital currency exchange. We study how the structure of the Minimum Spanning Tree (MST) and the Triangulated Maximally Filtered Graph (TMFG) evolve from high (15 s) to low (1 day) frequency time resolutions. For each horizon, we test the stability, statistical significance and economic meaningfulness of the networks. Results give a deep insight into the evolutionary process of the time dependent hierarchical organization of the system under analysis. A decrease in correlation between pairs of cryptocurrencies is observed for finer time sampling resolutions. A growing structure emerges for coarser ones, highlighting multiple changes in the hierarchical reference role played by mainstream cryptocurrencies. This effect is studied both in its pairwise realizations and intra-sector ones. Keywords: complex systems; network science; econophysics; economics; financial markets; cryptocurrencies 1. Introduction Financial markets are complex systems [ 1 ]. The main source of complexity comes from the intricate interaction of heterogeneous actors following various strategies designed to impact at different time scales. They are highly stochastic environments with a low signal to noise ratio, dominated by strong non-stationary dynamics and characterized by feedback loops and non-linear effects [ 2 – 4 ]. Despite their complexity, financial systems are governed by a rather stable and partially identified framework of rules [ 5 ]. This last characteristic, jointly with the possibility to continuously monitor them across time, makes financial systems well suited for statistical characterization [ 6 ] and a good playground for the study of complex systems in general. In this paper, we analyse the behaviour of cryptocurrency market. A cryptocurrency is defined as a digital instrument for value transfer that exploits cryptography and distributed ledgers for security and decentralization [ 7 ]. As currencies, they have properties similar to fiat currencies [ 8 ]. The main differences being the exclusion of financial institutions as intermediaries [ 9 ] and not being controlled and regulated by any central authority [ 10 ]. Thanks to the above mentioned characteristics, cryptocurrency market is available 24 h a day, 7 days a week, and transactions take place between individuals with different physical locations across the globe [ 8 ]. Standard features of financial systems joined with peculiarities listed above, make cryptocurrencies highly volatile instruments. Finding assets with similar behaviours responding to endogenous or exogenous events is, hence, a challenging, but extremely valuable, exercise both from theoretical and applicative perspectives (e.g., risk management and investment). The ready access availability of large volumes of market data ease research on these instruments with respect to classical financial ones. Indeed, one of the main limits faced by research in the field of financial applications is the lack of easy access and share of high-quality data. In most cases they are sensible data, owned and managed by private financial institutions. Entropy 2022,24, 1548. https://doi.org/10.3390/e24111548 https://www.mdpi.com/journal/entropy 1 Entropy 2022,24, 1548 %7&86' (7+86' /7&86' (a) (7+86' /7&86' (b) Figure 2. Triangulated Maximally Filtered graphs representing log-returns time series’ dependency structure computed at ( a ) 15 s and ( b ) 1 day. Only hub nodes are labelled. The adopted colour mapping scheme follows the sectors’ taxonomy by [42] (see Appendix B). As a preliminary step into the study of the information level carried by the two networkbased information filtering approaches, Figure 3 shows how pairwise (see Figure 3a ) and the intra-sector (see Figure 3b) average Pearson’s correlation coefficient ρ evolves as a function of time horizon Δt . Figure 3a reports the mean Pearson’s correlation coefficient computed averaging over the n(n− 1 )/ 2 = 300 off-diagonal elements of the whole correlation matrix C at different time horizons. In order to give a more comprehensive view of the evolutionary dynamics of the mean pairwise correlation coefficient, we also report three meaningful percentile intervals. We observe that the average correlation coefficient ρ increases with time horizon Δt from a value equals to 0.19 at Δt= 15 to a value equals to 0.47 at Δt=86,400 . The value at Δt= 15 corresponds to the minimum average correlation coefficient across time horizons. On the other hand, the maximum average correlation coefficient does not coincide with the one computed at the maximum time horizon. It is instead detected at horizon Δt=14,400 , which corresponds to an intra-day resolution (i.e., 4 h). On average, the most prominent pairwise correlation weakenings are observed for most correlated pair of assets (i.e., those pairs of cryptocurrencies having a correlation coefficient included into highest percentiles). 15 60 900 3600 14,400 86,400 Horizon (seconds) 0.0 0.2 0.4 0.6 0.8 Correlation coefficient (ρ) ρ 1% - 99% 10% - 90% 25% - 75% (a) 15 60 900 3600 14,400 86,400 Horizon (seconds) 0.2 0.3 0.4 0.5 0.6 0.7 Correlation coefficient (ρ) Currencies Smart Contracts Centralised Exchanges (b) Figure 3. Evolutionary dynamics of the average correlation coefficient as a function of the time horizon Δt .( a ) reports the horizon-related mean Pearson’s correlation coefficient and three meaningful percentiles computed averaging over the n(n− 1 )/ 2 = 300 off-diagonal elements of the whole correlation matrix C .( b ) reports the horizon related mean Pearson’s correlation coefficients computed averaging over the ns(ns− 1 )/ 2 correlation coefficients of the ns assets belonging to one of three of the most relevant sectors defined by [ 42 ]: Currencies, Smart Contracts, Centralised Exchange sectors. 8 Entropy 2022,24, 1548 Figure 3b reports mean Pearson’s correlation coefficient computed averaging over the ns(ns− 1 )/ 2 correlation coefficients of the ns assets belonging to one specific sector [ 42 ]at different time horizons. Specifically, we report dynamics for Currencies, Smart Contracts, and Centralized Exchanges sectors. This choice is completed considering the relevance of the three sectors. The relevance of sectors is defined in relation to results discussed later in this section. An intra-sector scenario shows trends comparable to the ones observed in pairwise context. All the previously discussed dynamics are here more pronounced. In both cases, we observe the “Epps effect”, i.e., a decrease in pair correlations at finer time sampling resolutions. This effect has been extensively studied in equity markets by [ 22 , 30 ]. Results reported in Figure 3 show how, also in the cryptocurrency market, the intra-sector correlation increases faster than pairwise one. The “Epps effect” is, hence, more pronounced within each sector than outside it. Going deeper, in Appendix C, we compare the probability distribution of correlation coefficients in the empirical correlation matrix C with the probability distribution of correlation coefficients filtered, respectively, by the MST and by the TMFG at different time horizons. We also report the probability distribution of correlation coefficients for surrogate multivariate time series obtained by randomly shuffling log-returns time series of the 25 cryptocurrencies listed in Table 1. This step is performed in order to evaluate the null hypothesis of uncorrelated returns for the considered portfolio of cryptocurrencies. Results give us the possibility to asses the statistical significance of average correlation coefficients chosen both by MST and by TMFG networks. These findings are reported in a synthetic way in Table 2. The extended count and the corresponding statistical meaning of links having a value higher than the minimum and lower than the maximum correlation coefficient detected by shuffling log-returns time series at different time horizons for the three scenarios are reported in Appendix D. Table 2. Average absolute correlation coefficient |ρ| and quantiles (25–75%) computed on the empirical correlation matrix C , on the links filtered by MST and on the ones filtered by TMFG at different time horizons. Statistical significance of the average correlation coefficient is represented though asterisks. p-values > 0.05 are not marked. p-values ≤ 0.05 are marked as ∗ .p-values ≤ 0.01 are marked as ∗∗ .p-values ≤ 0.001 are marked as ∗∗∗ . The filtering power of the MST and TMFG is evident considering that the related mean correlation coefficients are always greater than the ones computed on the whole correlation coefficient matrix C . Results for both MST and TMFG are always robust across time horizons. ΔtC MST TMFG |ρ| 25% 75% |ρ| 25% 75% |ρ| 25% 75% 15 0.20 0.02 0.31 0.35 ∗∗∗ 0.26 0.49 0.31 ∗∗ 0.22 0.42 60 0.31 0.05 0.46 0.47 ∗∗∗ 0.42 0.63 0.44 ∗∗∗ 0.38 0.57 900 0.44 ∗∗ 0.22 0.62 0.60 ∗∗∗ 0.57 0.76 0.57 ∗∗∗ 0.55 0.69 3600 0.46 ∗∗ 0.28 0.63 0.62 ∗∗∗ 0.59 0.76 0.59 ∗∗∗ 0.56 0.70 14,400 0.49 ∗0.38 0.65 0.65 ∗∗∗ 0.64 0.77 0.62 ∗∗∗ 0.58 0.72 86,400 0.48 0.38 0.62 0.66 ∗∗ 0.64 0.77 0.61 ∗∗ 0.57 0.72 Average correlation coefficients for MSTs and TMFGs are always greater than the ones computed on the empirical correlation matrix C . The difference between cross-horizons mean of average correlation coefficients filtered by MSTs and cross-horizons mean of average correlation coefficients in C , is equal to 0.16. The difference between cross-horizons mean of average correlation coefficients filtered by TMFGs and cross-horizons mean of average correlation coefficients in C , is equal to 0.12. Correlation coefficients filtered by TMFGs are always lower than the ones filtered by MSTs. This depends on the fact that, as reported in Section 3.5, the TMFG contains, by construction, more information than the MST. The mean difference between average correlation coefficients filtered by MSTs and the ones filtered by TMFGs, is equal to 0.03. Results reported in Table 2 confirm that the two filtering approaches prune weakest correlations among considered cryptocurrencies keeping only the strongest ones. Differently from what happens for the empirical correlation matrix 9 Entropy 2022,24, 1548 C , results for both MST and TMFG are always statistical significant across time horizons. These results enforce the evidence that both MST and TMFG carry information about strongest interactions observed in the system, disregarding most of the links consistent with the null hypothesis of uncorrelated data. It is worth noting that such an analysis does not tell much about the statistical robustness of links selected by the two network-based information filtering approaches. In order to perform such an investigation, we adopt the technique proposed by [ 54 ]. For each time horizon Δt , we sample 1000 bootstrap replicas r= 1, ... ,1000 of the empirical log-returns time series data. The length of empirical data and the one of each replica is kept equal. We compute the MST*(r) and the TMFG*(r) for each replica r. For each time sampling resolution, we map each link of the original MST and TMFG to an integer number and we count the number of links present both in the MST and TMFG and in each of the MST*(r) and TMFG*(r). Table 3 reports, for each time sampling interval Δt , the number of links of the empirical MST and TMFG with a bootstrap value larger than 95%. Table 3. Percentage of links contained in empirical MST and TMFG at time horizon Δt with a bootstrap value larger than 95%. In the case of the MST, it is possible to notice how the robustness of the network structure decreases for coarser time sampling intervals. In the case of TMFG, on the contrary, the robustness is maintained across time horizons with low oscillations. ΔtMST TMFG 15 62.5% 28.9% 60 58.3% 37.7% 900 54.2% 36.2% 3600 58.3% 36.2% 14,400 41.6% 40.6% 86,400 25.0% 27.5% Results in Table 3 show how, in the case of the MST, the robustness of the underlying network structure decreases for coarser time sampling resolutions. A consistent result has been observed by [ 50 ] in equity markets. This finding can be explained in two different ways. The first and most straightforward explanation is the statistical one and can be resumed as follows: the higher number of samples at finer time sampling resolutions implies higher statistical significance, while the lower number of samples at coarser time sampling resolutions imply lower statistical significance. A second explanation can be given looking at the structure of the networks reported in Appendix A. At finer time sampling resolutions, we observe less structured networks where numerous small-degree nodes (spokes) coexist with few anchor ones (hubs) characterised by an exceptionally high number of links. At coarser time sampling resolutions we observe more structured networks with a less imbalanced degree distribution. Such a topological change directly implies a loss in the links’ statistical robustness. The case of TMFG is different. Statistical robustness of the network is maintained across horizons without significant draw-downs. Indeed, during the optimization phase of the objective function, TMFG tends to be marginally exposed to local minima, being robust to dramatic topological changes. These last findings can be formally characterised studying the evolution of the average shortest path in MST and in TMFG as a function of time sampling resolution. Figure 4 reports the significant different behaviour in compactness’ evolutionary dynamics of the two network-based information filtering approaches. In the case of MST, the minimum length of the average shortest path is equal to 2.46 and is detected at Δt= 15, while the maximum length is equal to 3.05 and is detected at Δt=86,400 . In the case of TMFG, we observe a strong compactness across time horizons. The minimum length of the average shortest path is equal to 1.83 and is detected at Δt= 3600, while the maximum length is equal to 1.9 at Δt= 60. In the case of MST, at the finest time sampling resolution (i.e., Δt= 15), we observe a structurally simple network with two cryptocurrencies (i.e., Ethereum and Bitcoin) acting as a hierarchical reference for the majority of other assets. 10 Entropy 2022,24, 1548 This topological structure persists switching to time horizon Δt= 60. Several changes in nodes’ reference roles can be observed for networks sampled at time horizons Δt= 900 and Δt=3600. 15 60 900 3600 14,400 86,400 Horizon ( seconds ) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 Average length of the shortest path Average shortest path MST Average shortest path TMFG Figure 4. Average length of the shortest path in MST and TMFG as function of the time horizon at which log-returns are computed. We observe a decreasing compactness of MST networks at coarser time sampling resolutions. Instead, the compactness of the TMFG turns out to be stable across time horizons with low oscillations. In both cases Ethereum maintains its reference role even reducing its centrality. Bitcoin, on the contrary, is gradually replaced in its role by Litecoin and FTX token (both part of the Bitcoin’s cluster at time horizon Δt= 60). This structural transition is evident at Δt=14,400 and fully realised at Δt=86,400. In the case of TMFG representations, it is harder to graphically detect similar dynamics. Figure 5 offers a comparative perspective between behaviours of the two networkbased information filtering approaches. It shows horizon dependent evolutionary dynamics of degree centrality (i.e., measurement of the number of connections owned by a node) [ 55 – 57 ] for Ethereum, Bitcoin, Litecoin, and FTX Token both for MST and for TMFG. Cross-assets similarities can be detected between the two types of graphs. In the case of MSTs, degree centrality is less sensitive to minor changes in reference roles played by mainstream cryptocurrencies across time horizons, amplifying only ‘extreme’ ones. In the case of TMFGs, on the contrary, the same centrality measure is able to capture even small variations in network structure. This observation can be easily explained considering the amount of information the two representations are able to express. This study can be further extended looking at sectors of cryptocurrencies instead of at singular assets. Figure 6 reports the evolution of degree centrality for the three sectors Ethereum, Bitcoin, Litecoin, and FTX Token belong to: the Currencies sector, the Smart Contracts sector, and the Centralized Exchanges sector. We remark that there is no consensus on a unique mapping between cryptocurrencies and sectors. The taxonomy adopted in the current paper is described in [42] and corresponds to the one used by Kraken [28]. Figure 6a shows how, in the case of MST, the average degree centrality for the Smart Contracts sector strongly decreases starting from time horizon Δt= 3600, following the trend of its leading representative: Ethereum cryptocurrency. The Currencies sector, on the other hand, does not experience a decreasing trend and tends to remain stable across time horizons with low level of oscillations. In this case the loss of centrality of Bitcoin after time horizon Δt= 900, is immediately compensated by Litecoin, which reaches a hierarchical reference role at coarser time sampling resolutions. The case of Centralized Exchanges sector is different. It is stable across time horizons, without experiencing any change in intra-sector reference role dynamics and always following the behaviour of its main representative, FTX token (see Figure 5a). This last finding can be explained considering the source of the data used in the current research work. As explained in Section 3.1, we fetch data from the FTX digital currency exchange. This can cause, on the one hand an over-estimation of the role played by the exchange specific token, FTX Token, in the whole 11 Entropy 2022,24, 1548 ecology of the system under investigation, and, on the other hand, can give a potentially biased stability to the sector the asset belongs to. 15 60 900 3600 14,400 86,400 Horizon ( seconds ) 0.0 0.2 0.4 0.6 0.8 1.0 Degree Centrality ETH BTC LTC FTT (a) 15 60 900 3600 14,400 86,400 Horizon ( seconds ) 0.0 0.2 0.4 0.6 0.8 1.0 Degree Centrality ETH BTC LTC FTT (b) Figure 5. Degree centrality computed on MST ( a ) and on TMFG ( b ) as a function of time sampling resolution. Results on the TMFG highlight the switch in the reference roles of mainstream cryptocurrencies. 15 60 900 3600 14,400 86,400 Horizon ( seconds ) 0.0 0.1 0.2 0.3 0.4 0.5 Average Degree Centrality Currencies Smart Contracts Centralized Exchanges (a) 15 60 900 3600 14,400 86,400 Horizon ( seconds ) 0.0 0.1 0.2 0.3 0.4 0.5 Average Degree Centrality Currencies Smart Contracts Centralized Exchanges (b) Figure 6. Group degree centrality computed on MST ( a ) and on TMFG ( b ) for Currencies sector, Smart Contracts sector, and Centralized Exchanges sector. Group degree centrality of a set of nodes is defined as the fraction of non-group members connected to group members. Sectors are defined following the taxonomy by [42]. 5. Conclusions We investigate how cryptocurrency market’s dependency structures evolve passing from high to low frequency time sampling resolutions. Starting from the log-returns of 25 liquid cryptocurrencies traded on the FTX digital currency exchange at 6 different time horizons spanning from 15 s to 1 day, we investigate pairwise correlations demonstrating that cryptocurrency market has an “Epps effect” which is comparable to the one widely studied in the equity market. Indeed, we show that the average correlation among assets increases moving from high to low frequency time horizons and we demonstrate how this dynamic is even more evident grouping cryptocurrencies into sectors. Using the concept of power dissimilarity measure, we review the building process of two network-based information filtering approaches: MST and TMFG. If, on the one hand, MST has been historically used in the description of dependency structures of different financial markets, on the other hand, this is the very first time TMFG is used to study interactions between digital assets at different time scales. Studying topologies of MSTs at finer time sampling resolutions, we observe structurally simpler networks characterised by an hub-and-spoke configuration with statistically robust links. We observe an increase in the complexity of the networks’ shape for coarser time sampling resolutions with a decrease in links’ statistical robustness. Such an horizon-dependent structural change is reflected by the average path length of the networks, characterised by an increasing trend moving from 12 Entropy 2022,24, 1548 high to low frequencies. TMFG offers a different perspective for the same problem. In this case, we do not observe dramatic changes in networks’ topologies across time horizons. Graphs are more compact and statistical robustness of links is maintained across time with negligible oscillations. As a consequence of this, the average path length is lower and almost constant across time horizons. Studying the relative position of assets in both MSTs and TMFGs through the usage of degree centrality measure, we outline the presence of multiple changes in the hierarchical reference role among the considered set of cryptocurrencies. These changes strongly characterise singular cryptocurrencies. We find that Ethereum acts as a hierarchical reference node for the majority of other assets and maintains this role across time, gradually losing its centrality at coarser time horizons. There is not a clear economic explanation for this result. We know that lots of other cryptocurrencies are based on the Ethereum’s blockchain technology but we do not think this represents a sufficient explanation to our finding. Other cryptocurrencies play a similar role with respect to smaller clusters of assets at specific time horizons. We refer specifically to Bitcoin, Litecoin, and FTX Token. Differently from Ethereum, their role does not emerge at finer time sampling resolutions and should be considered as the result of a structured evolutionary process across time horizons. We conclude stating that sectors’ dynamics captured by the chosen network-based information filtering approaches are poorly affected by the ones of their main representatives, efficiently absorbing horizon-dependent changes in cryptocurrencies dynamics. This is true especially for TMFG. Indeed, looking at the evolution of the degree centrality of the Smart Contracts and Currencies sectors, one can observe that dynamics captured by MST are strongly influenced by the ones of Ethereum and Bitcoin. This does not happen in the case of TMFG where sectors’ dynamics are typically detached from the ones of specific cryptocurrencies. Author Contributions: Conceptualization, A.B. and T.A.; methodology, A.B. and T.A.; software, A.B.; validation, A.B. and T.A.; formal analysis, A.B. and T.A.; investigation, A.B. and T.A.; writing— original draft preparation, A.B.; writing—review and editing, T.A.; visualization, A.B.; supervision, T.A.; project administration, T.A.; funding acquisition, T.A. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by ESRC (ES/K002309/1), EPSRC (EP/P031730/1) and EC (H2020-ICT-2018-2 825215). Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Data are accessible for free using the CCXT [41] Python Package. Acknowledgments: The authors acknowledge many members of the Financial Computing and Analytics Group at University College London. A special thank to Silvia Bartolucci, David Vidal- Tomás, and Yuanrong Wang. Additionally, thanks to Agne Kazakeviciute for fruitful discussions on foundational topics related to this work. Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results. Abbreviations The following abbreviations are used in this manuscript: DCE Digital Currency Exchange MST Minimum Spanning Tree PMFG Planar Maximally Filtered Graph TMFG Triangulated Maximally Filtered Graph 13 Entropy 2022,24, 1548 Appendix A %7&86' (7+86' (a) %7&86' (7+86' (b) %&+86' (7+86' )7786' /7&86' (c) %&+86' (7+86' )7786' /7&86' (d) (7+86' )7786' /7&86' (e) (7+86' )7786' /7&86' (f) Figure A1. Minimum Spanning Tree representing log-returns time series’ dependency structure computed at ( a )15s,( b ) 1 min, ( c ) 15 min, ( d )1h,( e )4h,and( f ) 1 day. The adopted colour mapping scheme follows the sectors’ taxonomy by [ 42 ]: red → currencies, green → smart contract platforms, blue → stablecoins, pink → centralized exchanges, orange → scaling, turquoise → decentralized exchanges, fuchsia → lending, and yellow → all the other sectors. Dashed, red edges represent negatives linear correlations among pairs of cryptocurrencies. Only hub nodes are labelled. 14 Entropy 2022,24, 1548 Appendix B %7&86' (7+86' /7&86' (a) %7&86' (7+86' /7&86' (b) $$9(86' %7&86' (7+86' (c) %7&86' (7+86' /7&86' (d) (7+86' )7786' /7&86' (e) (7+86' /7&86' (f) Figure A2. Triangulated Maximally Filtered Graphs representing log-returns time series’ dependency structure computed at ( a )15s,( b ) 1 min, ( c ) 15 min, ( d ) 1 h, ( e ) 4 h, and ( f ) 1 day. The adopted colour mapping scheme follows the sectors’ taxonomy by [ 42 ]: red → currencies, green → smart contract platforms, blue → stablecoins, pink → centralized exchanges, orange → scaling, turquoise → decentralized exchanges, fuchsia → lending, and yellow → all the other sectors. Dashed, red edges represent negatives linear correlations among pairs of cryptocurrencies. Only hub nodes are labelled. 15 Entropy 2022,24, 1548 Appendix C 0.00.10.20.30.40.50.60.70.8 Correlation Coefficient ( ρ) 0.00 0.05 0.10 0.15 0.20 0.25 0.30 Probability MST TMFG Min shuffled data Max shuffled data C Shuffled data −0.002 0.000 0.002 0.00 0.05 0.10 0.15 (a) 0.00.20.40.60.8 Correlation Coefficient ( ρ) 0.00 0.05 0.10 0.15 0.20 0.25 Probability MST TMFG Min shuffled data Max shuffled data C Shuffled data (b) 0.00.20.40.60.8 Correlation Coefficient ( ρ) 0.00 0.05 0.10 0.15 0.20 Probability MST TMFG Min shuffled data Max shuffled data C Shuffled data (c) 0.00.20.40.60.8 Correlation Coefficient ( ρ) 0.00 0.05 0.10 0.15 0.20 0.25 Probability MST TMFG Min shuffled data Max shuffled data C Shuffled data (d) 0.00.20.40.60.8 Correlation Coefficient ( ρ) 0.00 0.05 0.10 0.15 0.20 0.25 Probability MST TMFG Min shuffled data Max shuffled data C Shuffled data (e) −0.20.00.20.40.60.8 Correlation Coefficient ( ρ) 0.00 0.05 0.10 0.15 0.20 0.25 Probability MST TMFG Min shuffled data Max shuffled data C Shuffled data −0.10.00.1 0.000 0.005 0.010 0.015 0.020 (f) Figure A3. Probability distribution of correlation coefficients for the empirical correlation matrix C , MST, TMFG, and correlation matrix of shuffled log-returns time series computed at ( a )15s,( b ) 1 min, (c) 15 min, (d)1h,(e)4h,and(f) 1 day. 16 Entropy 2022,24, 1548 Appendix D Table A1. Number of links of the empirical correlation matrix C , of the MST and of the TMFG having a value higher than the minimum and lower than the maximum correlation coefficient detected by shuffling log-returns time series at different time horizons. Shuffling operation is repeated 100 times. Results can be interpreted as p-values of the average correlation coefficient computed for C , for the MST and for the TMFG. ΔtC MST TMFG 15 36 0 1 60 28 0 0 900 2 0 0 3600 2 0 0 14,400 16 0 0 86,400 34 1 3 References 1. Anderson, P.W. The Economy as an Evolving Complex System; CRC Press: Boca Raton, FL, USA, 2018. 2. Comerton-Forde, C.; Putnin ,š, T.J. Dark trading and price discovery. J. Financ. Econ. 2015,118, 70–92. [CrossRef] 3. Briola, A.; Turiel, J.; Marcaccioli, R.; Aste, T. Deep reinforcement learning for active high frequency trading. arXiv 2021 , arXiv:2101.07107. 4. Briola, A.; Turiel, J.; Aste, T. Deep learning modeling of limit order book: A comparative perspective. arXiv 2020 , arXiv:2007.07319. 5. Lux, T.; Marchesi, M. Scaling and criticality in a stochastic multi-agent model of a financial market. Nature 1999 ,397, 498–500. [CrossRef] 6. Aste, T.; Shaw, W.; Di Matteo, T. Correlation structure and dynamics in volatile markets. New J. Phys. 2010 ,12, 085009. [CrossRef] 7. Goodell, G. Tokens and Distributed Ledgers in Digital Payment Systems. arXiv 2022, arXiv:2207.07530. 8. Fang, F.; Ventre, C.; Basios, M.; Kanthan, L.; Martinez-Rego, D.; Wu, F.; Li, L. Cryptocurrency trading: A comprehensive survey. Financ. Innov. 2022,8, 1–59. [CrossRef] 9. Harwick, C. Cryptocurrency and the problem of intermediation. Indep. Rev. 2016,20, 569–588. 10. Rose, C. The evolution of digital currencies: Bitcoin, a cryptocurrency causing a monetary revolution. Int. Bus. Econ. Res. J. 2015 , 14, 617–622. [CrossRef] 11. Wasserman, S.; Faust, K. Social Network Analysis: Methods and Applications; Cambridge University Press: Cambridge, UK, 1994. 12. Newman, M.E.; Watts, D.J.; Strogatz, S.H. Random graph models of social networks. Proc. Natl. Acad. Sci. USA 2002 , 99, 2566–2572. [CrossRef] 13. Ronfeldt, D.F.; Arquilla, J. Networks and Netwars; Rand: Santa Monica, CA, USA, 2001. 14. Balcan, D.; Hu, H.; Goncalves, B.; Bajardi, P.; Poletto, C.; Ramasco, J.J.; Paolotti, D.; Perra, N.; Tizzoni, M.; Van den Broeck, W.; et al. Seasonal transmission potential and activity peaks of the new influenza A (H1N1): A Monte Carlo likelihood analysis based on human mobility. BMC Med. 2009,7, 45. [CrossRef] [PubMed] 15. Hufnagel, L.; Brockmann, D.; Geisel, T. Forecast and control of epidemics in a globalized world. Proc. Natl. Acad. Sci. USA 2004 , 101, 15124–15129. [CrossRef] [PubMed] 16. Sporns, O.; Tononi, G.; Kötter, R. The human connectome: A structural description of the human brain. PLoS Comput. Biol. 2005 , 1, e42. [CrossRef] [PubMed] 17. Hopkins, A.L. Network pharmacology. Nat. Biotechnol. 2007,25, 1110–1111. [CrossRef] [PubMed] 18. Wu, L.; Waber, B.N.; Aral, S.; Brynjolfsson, E.; Pentland, A. Mining Face-to-Face Interaction Networks Using Sociometric Badges: Predicting Productivity in an It Configuration Task. 2008. Available online: https://papers.ssrn.com/sol3/papers.cfm?abstract_id= 1130251 (accessed on 22 March 2022). 19. Mantegna, R.N. Hierarchical structure in financial markets. Eur. Phys. J. B Condens. Matter Complex Syst. 1999 ,11, 193–197. [CrossRef] 20. Bonanno, G.; Caldarelli, G.; Lillo, F.; Mantegna, R.N. Topology of correlation-based minimal spanning trees in real and model markets. Phys. Rev. E 2003,68, 046130. [CrossRef] 21. Bonanno, G.; Caldarelli, G.; Lillo, F.; Micciche, S.; Vandewalle, N.; Mantegna, R.N. Networks of equities in financial markets. Eur. Phys. J. B 2004,38, 363–371. [CrossRef] 22. Bonanno, G.; Lillo, F.; Mantegna, R.N. High-frequency cross-correlation in a set of stocks. Quant. Financ. 2001 ,1, 96–104. [CrossRef] 23. Wang, Y.; Aste, T. Dynamic Portfolio Optimization with Inverse Covariance Clustering. arXiv 2021, arXiv:2112.15499. 24. Procacci, P.F.; Aste, T. Portfolio Optimization with Sparse Multivariate Modelling. arXiv 2021, arXiv:2103.15232. 25. Wang, Y.; Aste, T. Sparsification and Filtering for Spatial-temporal GNN in Multivariate Time-series. arXiv 2022 , arXiv:2203.03991. 26. West, D.B. Introduction to Graph Theory; Prentice Hall: Upper Saddle River, NJ, USA, 2001; Volume 2. 17 Entropy 2023,25, 219 3. Experimental Procedure Information regarding the stages of the experimental procedure will now be presented. This presentation will be as detailed as possible given the necessary space constraints and content commitments in order not to disrupt the depictive nature of the paper. It has already been mentioned that to some extent, the “core” of the present work consists of an experimental procedure that aims, in its most abstract scope, to check the efficiency, on the one hand, of a number of state-of-the-art algorithms and, on the other, of incorporating sentiment analysis into predictive schemas. Thus, a total of 16 datasets × 67 combinations × 30 algorithms × 3time-shifts = 96,480 experiments were conducted. The dataset consisted of time series containing the daily closing values of various stocks along with a multitude of 67 different sentiment score setups. Specifically, 16 datasets of stocks containing such closing price values were used over a three-year period, beginning on 2 January 2018 and ending on 24 December 2020. Generated sentiment scores from relevant textual data extracted from the Twitter microblogging platform were used. Three different sentiment analysis methods were deployed. The sentiment score time series and the closing values were subjected to a 7-day and a 14-day rolling mean strategy, yielding a total of 12 distinct features. Various combinations of the created features resulted in a total of 67 distinct input setups per algorithm. The calculated sentiment scores along with the closing values were then tested under both univariate and multivariate forecasting schemes. Lastly, 30 state-of-the-art methods were investigated. Below, a more thorough presentation of the aforementioned experimental setting follows. 3.1. Datasets Starting with data, the process of collecting and creating the sets used will now be addressed. 3.1.1. Overview To begin with, Table 1 contains the names of the aforementioned datasets along with their corresponding abbreviations. These initial data included time series containing closing values for 16 well-known listed companies. All sets comprise three-year period data for dates ranging from 2 January 2018 to 24 December 2020. Table 1. Stock datasets. No Dataset Stocks 1 AAL American Airlines Group 2 AMD Advanced Micro Devices 3 AUY Yamana Gold Inc. 4 BABA Alibaba Group 5 BAC Bank of America Corporation 6 ET Energy Transfer L.P. 7 FCEL FuelCell Energy Inc. 8 GE General Electric 9 GM General Motors 10 INTC Intel Corporation 11 MRO Marathon Oil Corporation 12 MSFT Microsoft Corporation 13 OXY Occidental Petroleum Corporation 14 RYCEY Rolls-Royce Holdings 15 SQ Square 16 VZ Verizon Communications Essentially, the initial features were four: that is, the closing prices of each stock and three additional time series containing relative sentiment scores for the given period. Subsequently, and after applying 7- and 14-day rolling averages, a total of 14 features were extracted. Thus, for each share, the final input settings were composed by introducing 24 Entropy 2023,25, 219 altered features derived from stock values and a sentiment analysis process applied to an extended corpus of tweets. Figure 1 depicts a—rather abstractive—snapshot of the whole process from data collection to the creation of the final input setups. Figure 1. Feature setups: creation pipeline. 25 Entropy 2023,25, 219 3.1.2. Tweets and Preprocessing A large part of the process involved deriving sentiment scores related to stocks. Using the Twitter Intelligence Tool (TWINT) [ 53 ], a large number of stock-related posts written in English were downloaded from Twitter and grouped by day. TWINT is an easy-to-use yet sophisticated Python-based Twitter scraping tool. After a comprehensive search for stockrelated remarks that were either directly or indirectly linked to shares under consideration, a sizable amount of text data containing daily attitudes toward stocks were created. Then, the collected textual sets underwent the various preprocessing procedures necessary in order to be passed on to the classification modules for extracting their respective sentiment scores. Regarding preprocessing tweets, initially, irrelevant hyperlinks and URLs were removed using the Re Python library [ 54 ]. Each tweet was then converted to lowercase and split into words. Then, unwanted phrases from a manually produced list and various numerical strings were also dismissed. After performing the necessary joins to restore each text to its original structure, each tweet was tokenized in terms of its sentences using the NLTK [ 55 , 56 ] library. Lastly, using the String [ 57 ] module, punctuation removal was applied. The whole text-preprocessing step is schematically presented in Figure 2. Figure 2. Preprocessing. 3.1.3. Sentiment Analysis The subsequent process involved extracting sentiment scores from the gathered yet cleaned tweets. To perform the sentiment quantification step, three different sentiment analysis methods were utilized. Specifically, the procedure included extracting sentiment scores from TextBlob [ 58 ], using the Vader sentiment analysis tool [ 59 ], and incorporating FinBERT [ 60 ]. FinBERT is a financial-based fine-tuning of the BERT [ 61 ] language representation model. Using each of the above methods, daily sentiment scores were extracted for each stock. The daily mean was then extracted, forming the final collection, which constituted the sentimentvalued time series of every corresponding method. Then, 7- and 14-day moving averages were applied to the previously extracted sentiment score time series. This resulted in the extraction of nine sentiment time series, which, together with the application of the aforementioned procedure to the closing price time series, led to the final number of 12 generated time series used as features. Various combinations of the above features, along with the univariate case scenario, resulted in 67 different study cases. These data constituted the distinct experimental procedures that run for every algorithm. The use of three different methods of sentiment analysis has already been mentioned. Below, a rough description of these methods is given. For further information, the reader is advised to refer to the respective papers. • TextBlob: The TextBlob module is a Python-based library for performing a wide range of manipulations over text data. The specific TextBlob method used in this work is a rule-based sentiment-analysis scheme. That is, it works by simply applying manually created rules. This is how the value attributed to the corresponding sentiment score is calculated. An exemplified snapshot of the process would be counting the number of times a term of interest appears within a given section. This would modify the 26 Entropy 2023,25, 219 projected sentiment score values in line with the way the phrase is assessed. Here, within this experimental setup and by exploiting TextBlob’s sentiment property, a real number within the [− 1,1 ] interval representing the sentiment polarity score was generated for each tweet. The algorithm’s numerical output was then averaged using the individual scores of each tweet to obtain a single sentiment value representing the users’ daily attitudes; • Vader: Vader is also a straightforward rule-based approach for realizing general sentiment analysis. In the context of this work, the Vader sentiment analysis tool was used in order to extract a compound score produced by a normalization of sentiment values that the algorithm calculates. Specifically, given a string, the procedure outputs four values: negative, neutral, and positive sentiment values, as well as the aforementioned composite score used. A normalized average of all compound scores for each day was generated the usual way. The resulting time series contained daily sentiment scores that ranged within the [−1, 1]interval; • FinBERT: Regarding FinBERT, in this work, the implementation contained in [ 62 ] was utilized. Specifically, the model that was trained on PhraseBank presented in [63] was used. Again, first, the daily scores regarding sentiment attitudes were extracted to eventually form a daily average time series. Generally, the method is a pre-trained natural-language-processing (NLP) model for sentiment analysis. It is produced by simply fine-tuning the pre-trained BERT model over financial textual data. BERT, meaning bidirectional encoder representations from transformers, is an implementation of the transformers architecture used for natural language processing problems. The technique is basically a pre-trained representational model based on transfer learning principles. Given textual data, multi-layer deep representations are trained with a bidirectional attention strategy so that the various different contexts of each linguistic token constitute the content of the token’s embedding. Regardless of data references—here financial—the model can be fine-tuned in any domain by only using a single additional layer that addresses the specific tasks. 3.2. Algorithms In this section, the methods, algorithmic schemes, and architectures employed in the experiments are listed. Additional details are given on the implementation framework and the tools used. Regarding the algorithms used, a total of 30 different state-of-the-art methods and method variations were compared. The number of 30 methods used results from the supplementation of the set of well-known core methods with their variations. Further details can be found in the cited tsAI library [ 64 ], using which the implementation was carried out. However, it is this multitude of methods that apparently makes a detailed presentation practically impossible. Nevertheless, the reader is urged to track the cited papers. Table 2 contains the main algorithms utilized during the experimental procedure along with a corresponding citation. There, among others, one can notice that in addition to a multitude of state-of-the-art methods, implementations involving combinations of the individual architectures were also used. Note that in addition to the corresponding papers, information regarding the variations of the basic algorithms employed can be searched, inter alia, in notebook files taken from the library implementations. In order to carry out the experiments, the Python library tsAI [ 64 ] was used. The tsAI module is “an open-source deep learning package built on top of Pytorch and Fastai focused on state-of-the-art techniques for time series tasks like classification, regression, forecasting” [ 64 ], and others. Here, the forecasting procedure was essentially treated as a predictive regression problem. In the experiments, the initial parameters of the respective methods from the library were preserved with the implementation environment being kept fixed for all algorithmic schemes. Thus, all algorithms compared were utilized in the most basic configuration. That way, one can gain additional insight regarding implementing high-level yet low-code programming and data analysis in real-world tasks. Of the data, 27 Entropy 2023,25, 219 20% were used as the test set. Regarding prediction time horizons, three forecast scenarios were implemented: one single-step and two multi-step. In particular, with regard to multistep forecasts, and leaving aside the single-step predictions, estimates were provided for a seven-day window on the one hand and a fourteen-day window on the other. The results were evaluated according to the metrics presented in the following paragraph. Table 2. Algorithms. No. Abbreviation Algorithm 1 1 FCN Fully Convolutional Network [65] 2 FCNPlus Fully Convolutional Network Plus [66] 3 IT Inception Time [67] 4 ITPlus Inception Time Plus [68] 5 MLP Multilayer Perceptron[65] 6 RNN Recurrent Neural Network [69] 7 LSTM Long Short-Term Memory [70] 8 GRU Gated Recurrent Unit [71] 9 RNNPlus Recurrent Neural Network Plus [69] 10 LSTMPus Long Short-Term Memory Plus [69] 11 GRUPlus Gated Recurrent Unit Plus [69] 12 RNN_FCN Recurrent Neural—Fully Convolutional Network [72] 13 LSTM_FCN Long Short-Term Memory—Fully Convolutional Network [73] 14 GRU_FCN Gated Recurrent Unit—Fully Convolutional Network [74] 15 RNN_FCNPlus Recurrent Neural—Fully Convolutional Network Plus [75] 16 LSTM_FCNPlus Long Short-Term Memory—Fully Convolutional Network Plus [75] 17 GRU_FCNPlus Gated Recurrent Unit—Fully Convolutional Network Plus [75] 18 ResCNN Residual—Convolutional Neural Network [76] 19 ResNet Residual Network [65] 20 RestNetPlus Residual Network Plus [77] 21 TCN Temporal Convolutional Network [78] 22 TST Time Series Transformer [79] 23 TSTPlus Time Series Transformer Plus [80] 24 TSiTPlus Time Series Vision Transformer Plus [81] 25 Transformer Transformer Model [82] 26 XCM Explainable Convolutional Neural Network [83] 27 XCMPlus Explainable Convolutional Neural Network Plus [84] 28 XceptionTime Xception Time Model [85] 29 XceptionTimePlus Xception Time Plus [86] 30 OmniScaleCNN Omni-Scale 1D-Convolutional Neural Network [87] 1Methods and method variations used. 3.3. Metrics Regarding performance evaluation, six metrics were used. The use of the different metrics serves the necessity of having not only a presentation of the conclusions of a large comparison of methods and feature and sentiment setups but also a number of diverse extractions in terms of evaluation aspects that can be used in future research. This is exactly because each of the metrics exposes the results in different aspects, and therefore, an investigation would be incomplete if it focused on just one of them. Thus, regarding evaluating results, each one of the six performance indicators utilized has advantages and disadvantages. The metrics used are: • the Mean Absolute Error (MAE); • the Mean Absolute Percentage Error (MAPE); • the Mean Squared Error (MSE); • the Root Mean Squared Error (RMSE); • the Root Mean Squared Logarithmic Error (RMSLE); 28 Entropy 2023,25, 219 • the Coefficient of Determination R2. In what follows, a rather detailed description of aspects of the aforementioned wellknown evaluation metrics is given. The presentation aspires to provide details and some insight regarding the interpretation of the metrics. Below, the actual values are denoted by yaiand the forecasts are denoted by ypi. 3.3.1. MAE First is MAE: MAE =1 n n ∑ i=1ypi−yai(1) MAE stands for the arithmetic mean of the absolute errors, and it is a very straightforward metric and easy to calculate. By default, in terms of the difference between the prediction and the observation, the values share the same weights. The absence of exponents in the analytic form ensures good behavior, which is displayed even when outliers are present. The target variable’s unit of measurement is the one expressing the results. MAE is a scale-dependent error metric; that is, the scale of the observation is crucial. This means that it can only be used to compare methods in scenarios where every scheme incorporates the same specific target variable rather than different ones. 3.3.2. MAPE Next is MAPE: MAPE =1 n n ∑ i=1ypi−yai |yai|(2) MAPE is the mean absolute percentage error. It is a relative and not an absolute error measure. MAPE is common when evaluating the accuracy of forecasts. It is the average of the absolute differences between the prediction and the observations divided by the absolute value of the observation. A multiplication by 100 can afterwards convert this output to a percentage. This error cannot be calculated when the actual value is zero. Instead of being a percentage, in practice, it can take values in [0, ∞) . Specifically, when the predictions contain values much larger than the observations, then the MAPE output can exceed 100%. Conversely, in cases where both the prediction and the observation contain low values, the output of the metric may deviate greatly from 100%. This, in turn, can lead to a misjudgment of the model’s predictive capabilities, believing them to be limited when, in fact, the errors may be low. MAPE attributes more weight to cases where the predicted value is higher than the actual one. These cases produce larger errors. Hence, using this metric is best suitable for methods with low prediction values. Lastly, MAPE, being not scale-dependent, can be used to evaluate comparisons of a variety of different time series and variables. 3.3.3. MSE The next metric is MSE: MSE =1 n n ∑ i=1ypi−yai2(3) MSE stands for mean squared error. It constitutes a common forecast evaluation metric. The mean squared error is the average of the squares of the differences between the actual and predicted values. Its unit of measurement is the square of the unit of the variable of interest. Looking at the analytical form, first, the square of the differences ensures the nonnegativity of the error. At the same time, it makes information about minor errors usable. It is obvious, at the same time, that larger deviations entail larger penalties, i.e., a higher MSE. Thus, outliers have a big influence on the output of the error; that is, the existence of such extreme values has a significant impact on the measurements and, consequently, the evaluation. Furthermore, and in a sense the other way around, when differences are less than 1, there is a risk of overestimating the predictive capabilities of the model. Given the error’s differentiability, as one can observe, it can easily be optimized. 29 Entropy 2023,25, 219 3.3.4. RMSE Moving on to RMSE: RMSE =1 n n ∑ i=1ypi−yai2(4) RMSE stands for root mean squared error. It is a common metric for evaluating differences between estimated values and observations. To compute it, apparently, one just calculates the root of the mean squared error. From the numerical formulation, one can think of the metric as an abstraction that captures the representation of something of an average distance between the actual values and the predictions. That is, if one ignores the denominator, then one can observe the formula as being the Euclidean distance. The subsequent interpretation of the metric as a kind of normalized distance comes out of the act of division by the number of observations. Here also, the existence of outliers has a significant impact on the output. In terms of interpreting error values, the RMSE is expressed in the same units as the target variable and not in its square, as in the MSE, making its use straightforward. Finally, the metric is scale-dependent; hence, one can only use it to evaluate various models or model variations given a particular fixed variable. 3.3.5. RMSLE The next metric is also an error. The formula for RMSLE is as follows: RMSLE =1 n n ∑ i=1log(ypi+1)−log(yai+1)2(5) RMSLE stands for Root Mean Squared Logarithmic Error. The RMSLE metric seems as if it is a modified version of the MSE. Using this modification is preferred when predictions display significant deviations. RMSLE uses logarithms of both the observations and predicted values while ensuring non-zero values in the logarithms through the appropriate simple unit additions appearing in the formula. This modified version is resistant to the existence of outliers and noise, and it smooths the penalty that the MSE imposes in cases in which predictions deviate significantly from observations. The metric cannot be used when there are negative values. RMSLE can be interpreted as a relative error between observations and forecasts. This can be made evident by simply applying the following property to the radicand term of the square root: log(ypi+1)−log(yai+1)=logypi+1 yai+1(6) Since RMSLE gives more weight to cases where the predicted value is lower than the actual value, it is quite a useful metric for types of predictions where similar conditions require special care for the reliability of the application in real-world conditions, where lower forecasts may lead to specific problems. 3.3.6. R2 The last metric is the coefficient of determination R2: R2=1−SSRES SSTOT =1−∑n i=1ypi−yai2 ∑n i=1ypi−y2(7) The coefficient of determination R2 is not an error evaluation metric. It is the ratio depicted in the above equation. This metric is essentially not a measure of model reliability. R2 is a measure of how good a fit is: a quantification of how well a model fits the data. Its values typically range from 0 to 1. A rather simple interpretation would be this: the closer to 1 the value of the metric is, the better the model fits the observations, i.e., the predictions are closer, in terms of their values, to the observations. Thus, the value 0 corresponds to cases where the explanatory variables do not explain the variance of the dependent variable 30 Entropy 2023,25, 219 at all. Conversely, the value 1 corresponds to cases where the explanatory variables fully explain the dependent variable. However, this interval does not strictly constitute the set of values of the metric. There are conditions in which R2 could take negative values. Observing the formula, one can identify the above as permissible. In such cases, the model performs worse in fitting the data than a simple horizontal line, essentially being unable to follow the trend. Lastly, values outside the above range indicate either an inadequate model or other flaws in its implementation. 4. Results Returning to the dual objective of this work, the two case studies whose results will be presented in this chapter were: • On the one hand, the comparison of a large number of time series forecasting contemporary algorithms; • On the other hand, the investigation of whether knowledge of public opinion, as reflected in social networks and quantified using three different sentiment analysis methods, can improve the derived predictions. Accordingly, the presentation of the results of the experimental process is split into two distinct parts. In what follows, both various statistical analysis and visualization methods are incorporated. However, it should be noted that the number of comparisons performed yielded a quite large volume of results. Specifically, as already pointed out, in each case, the performance of the 30 predictive schemes and the 67 different feature setups was investigated over three different time frames (1, 7, and 14 day shifts). Note that these three time-shifting options have no—or at least no intended—financial consequences. Here, the primary goal in designing the framework was to forecast the stock market over short time frames, such as a few days. Then, an expansion was made to investigate the performance of both methods and feature setups over longer periods of time. Each of these schemas was evaluated with six different metrics, while the process was repeated for each of the datasets. Consequently, it becomes clear that the complete tables with the numerical results cannot contribute satisfactorily to the understanding of the conclusions drawn. Below, following a necessary brief reminder of the process, results are presented. As has already been mentioned, during the procedure, for each of the stocks, the following strategy was followed: each of the thirty algorithms to be compared was “ran” 67 times, each time accepting as input one of the different feature setups. This was repeated three times, once for each of the three forecast time frames. In each of the above runs, the six metrics used in the evaluation of the results were calculated. The comparison of the algorithms was performed by using Friedman’s statistical tests in terms of feature setups for each of the time shifts. Thus, given setups and stocks, the ranking of the methods per evaluation metric was extracted according to the use of the Friedman test [88] . Therefore, regarding this case study, a total of 67 × 6 × 3=1206 statistical tests were executed. In a similar way, the Friedman rankings of input feature setups were estimated in terms of metrics and time shifts, given the various algorithms and stocks. Here, a total of 30 ×6×3=540 statistical tests were performed. An additional abstraction of the results was derived as follows: For each of the 30 methods, the average rank achieved by each method in terms of feature setups and shares was calculated. So, for each metric and each of the three time frames, a more comprehensive display of the information was obtained based on the average value of the different setups. In an identical way, in the case of checking the effectiveness of features, the average value of the 30 algorithms for each of the 67 different input setups was calculated in each case. In both cases, the ranking was calculated based on the positions produced by the Friedman test, while at the same time, with the Nemenyi post hoc test [ 89 ] that followed, every schema was checked pair-wise for significant differences. The results of the Nemenyi post hoc tests are shown in the corresponding Critical Difference diagrams (CD-diagrams), in which methods that are not significantly different are joined by black horizontal lines. Two methods are considered not significantly different when the difference between their mean ranks is less than the CD value. 31 Entropy 2023,25, 219 Next, organized in both cases based on time frames, the results concerning the comparison of the forecast algorithms are presented, which are followed by those regarding the feature setups. 4.1. Method Comparison The presentation begins with results concerning the investigation of methods. The results are presented per forecast time shift. In each case, the Friedman Ranking results for all six metrics are listed. To save space, only methods that occupy the top ten positions of the ranking are listed. Full tables are available at: shorturl.at/FTU06 (accessed on 15 January 2023). The CD diagrams follow. There, we can visually observe which of the methods exhibit similar behavior and which differ significantly. Finally, box plots of results per metric are presented, again for the best 10 methods. The box plots present in a graphical and concise manner information concerning the distribution of the aforementioned data, that is, in our case, the average values of the sentiment setups per algorithm for all stocks. In particular, one can derive information about the maximum and minimum value of the data, the median, as well as the 1st and 3rd quartile values isolated by 25% and 75% of the observations, respectively. 4.1.1. Time Shift 1 With respect to the one-day forecasts, Table A1 lists the Friedman Ranking results for the top 10 scoring methods per metric. Although there is no single method that dominates all metrics and significant reorderings are also observed in the table positions, the TCN method achieves the best ranking in three out of six metrics (MAPE, R2, and RMSLE) and is always in the top four. Furthermore, from the box plots, it is evident that TCN has by far the smallest range of values. Apart from this, in all metrics, GRU_FCN is always in the top five. It is also observed that LSTM_FCN and LSTMPlus behave equally well. The latter shows a drop in the MAPE metric, but in all other cases, it is in the top three, while in two metricsm it ranks first. It should also be noted that the LSTMPlus method ranks first in two metrics, namely MAE and RMSE. In terms of R 2 and RMSLE, it occupies the second position of the ranking, while regarding MSE, LSTMPlus ranks third. However, at the same time, according to MAPE, the method is not even in the top ten. Thus, as will be seen in the following, TCN is the consistent choice. The results produced by Friedman’s statistical test, in terms of the six metrics, are presented in Table A1, while the corresponding CD diagrams and box plots are depicted in Figures 3 and 4. Figure 3. Box Plots: Methods—Shift 1. 32 Entropy 2023,25, 219 Figure 4. CD Diagrams: Methods—Shift 1. 4.1.2. Time Shift 7 At the one-week forecast time frame, the algorithms that occupy the top positions in the ranking produced by the statistical control appear to have stabilized. The corresponding ranking produced by the Friedman statistical test regarding the ten best methods with respect to the six metrics is presented in Table A2. In all metrics, the TCN method ranks first. From the CD diagrams, it can be seen that in all metrics—except for R2—this superiority is also validated by the fact that this method differs significantly from the others. Box plots show the method also having the smallest range around the median. Figures 5 and 6 contain the relevant results in the form of box plots and CD-diagrams. Figure 5. Box Plots: Methods—Shift 7. Other methods that clearly show some dominance over the rest in terms of given performance ratings are, on the one hand, TSTPlus, which ranks second in all metrics 33 Entropy 2023,25, 219 5.2. Feature and Sentiment Setups In relation to the second case study, the consideration of the results also points in some important directions. Of these, the main conclusion drawn seems to be that the use of information derived from both smoothed versions of the initial time series and sentiment analysis shows, in most cases, to have a beneficial effect on the derived forecasts. Not using sentiments in the feature setup of the inputs dominates the rest only in a small number of cases, and, as confirmed by the CD diagrams, only in two of them is this difference significant. Moreover, the answer to whether the use of sentiment setups specifically leads to the extraction of more accurate forecasts, as evidenced by the individual layouts of the weighted results, seems to be that, in general, sentiment analysis improves forecasts. Of course, it is also reasonable to investigate whether there is a specific sentiment setup that outperforms the rest. This would also lead to an assessment of the performance of the three sentiment analysis methods used. However, the answer to this question needs further investigation. However, even with the possibility of further inquiries within the framework of the experimental setup presented here, it is still not certain that firm conclusions will be drawn. Here, while such setups can be found for each time horizon, there is not one that dominates all three. In order, however, to illustrate a relative ranking of the three sentiment analysis methodologies used, regardless of the particular variation involved, an additional table was created. All variations of each method were placed under a corresponding class. The Friedman-aligned ranks [ 90 ] were then calculated. Hence, in order to draw a clearer picture of the way the three employed approaches to sentiment analysis performed, three sentiment classes were formed, one matching each of the previously described sentiment analysis methods. The arithmetic mean of all the sentiment setups that solely contain different variations of a particular sentiment analysis algorithm, that is, only one of the three incorporated, is used to represent the corresponding class concerning each metric. In other words, each class represents a sentiment analysis method, and each class corresponds to six sentiment setups that contain variations exclusively of the technique in question. Specifically, a representative value of a class, as it pertains to a particular method, is formed by the following setups: method, RM7method, RM14method,method + RM7method,method + RM14method, and RM7method + RM14method. The sum is then divided by six, which is apparently the number of setups, and this result is the output value to be depicted. This way, setups produced either by combining the various sentiment analysis methods or by using the target variable in variants containing rolling means are excluded in order to compare only the relative performances of the three individual techniques and their variations. Figure 16 illustrates these relative rankings of the three sentiment analysis methods per time shift. One can observe the relative performances in terms of individual wins with respect to each metric and time shift: the Blob and Vader classes top the ranking seven times each, while the Finbert class only has four wins. Again, a conclusion in terms of an obvious generality regarding a specific algorithm does not appear. Nevertheless, the identification of groups of such setups, even at the level of a specific time frame, can be particularly useful, with the methodology for the selection of individual setups needing more investigation. 40 Entropy 2023,25, 219 Figure 16. Sentiment rankings. Author Contributions: Conceptualization, C.M.L. and S.K.; methodology, C.M.L.; software, C.M.L.; validation, C.M.L., A.K. and S.K.; formal analysis, C.M.L. and A.K.; investigation, C.M.L. and A.K.; resources, S.K.; data curation, A.K.; writing—original draft preparation, C.M.L. and A.K.; writing—review and editing, C.M.L.; visualization, A.K.; supervision, S.K. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Institutional Review Board Statement: Not applicable. Data Availability Statement: URLs of the full Friedman Ranking results. (i) Methods rankings: shorturl.at/FTU06 (accessed on 15 January 2023). (ii) Feature setup rankings: shorturl.at/alqwx (accessed on 15 January 2023). Conflicts of Interest: The authors declare no conflict of interest. Appendix A Table A1. Friedman results: Algorithms—Shift 1. MAE MAPE R2 Method Friedman Score Method Friedman Score Method Friedman Score 1st LSTMPlus 8.266667 RNN_FCN 10.4 TCN 22.4 2nd LSTM 8.533333 GRU_FCN 10.53333 LSTMPlus 21.13333 3rd TCN 9.466667 TCN 11 LSTM 20.6 4th GRU_FCN 10 LSTM_FCN 11.2 GRU_FCN 20.53333 5th LSTM_FCN 10.2 GRU_FCNPlus 11.6 LSTM_FCN 19.73333 6th RNN 10.73333 RNN_FCNPlus 11.73333 LSTM_FCNPlus 19.13333 7th RNN_FCN 11.13333 RNN 11.93333 GRU_FCNPlus 18.93333 8th GRU_FCNPlus 11.33333 ResCNN 12.13333 RNN_FCN 18.93333 9th XCM 11.33333 LSTM_FCNPlus 12.13333 RNN 18.93333 10th LSTM_FCNPlus 11.4 FCNPlus 12.46667 XCMPlus 18.66667 41 Entropy 2023,25, 219 Table A1. Cont. MSE RMSE RMSLE Method Friedman Score Method Friedman Score Method Friedman Score 1st TCN 9 LSTMPlus 9.066667 TCN 7.733333 2nd GRU_FCN 9.266667 LSTM 9.6 LSTMPlus 9.333333 3rd LSTMPlus 9.6 GRU_FCN 9.6 LSTM 9.8 4th LSTM_FCN 9.8 TCN 9.733333 GRU_FCN 10 5th LSTM 9.933333 LSTM_FCN 10.06667 LSTM_FCN 10.2 6th RNN_FCN 10.33333 RNN_FCN 10.93333 RNN 11.13333 7th LSTM_FCNPlus 10.46667 LSTM_FCNPlus 11.06667 GRU_FCNPlus 11.26667 8th GRU_FCNPlus 10.8 RNN 11.13333 RNN_FCN 11.26667 9th RNN_FCNPlus 11.33333 GRU_FCNPlus 11.2 LSTM_FCNPlus 11.26667 10th FCNPlus 11.46667 RNN_FCNPlus 11.86667 GRU 12 Table A2. Friedman results: Algorithms—Shift 7. MAE MAPE R2 Method Friedman Score Method Friedman Score Method Friedman Score 1st TCN 3.733333 TCN 6.133333 TCN 25.86667 2nd TSTPlus 8.266667 XCMPlus 9.866667 TSTPlus 25.8 3rd XCMPlus 8.866667 RNNPlus 10.93333 XceptionTimePlus 19.7 4th XCM 10.53333 TSTPlus 11 XCMPlus 19.66667 5th RNN_FCNPlus 12.06667 RNN 11.06667 XceptionTime 19.53333 6th GRU_FCNPlus 12.13333 XCM 11.26667 XCM 18.53333 7th RNN_FCN 12.26667 LSTMPlus 13 RNN_FCN 16.9 8th GRU_FCN 13.2 GRU 13.66667 GRU_FCNPlus 16.66667 9th RNN 13.53333 ResCNN 13.86667 InceptionTime 16.5 10th LSTM_FCNPlus 13.53333 LSTM 14.06667 RNN_FCNPlus 16.16667 MSE RMSE RMSLE Method Friedman Score Method Friedman Score Method Friedman Score 1st TCN 3.666666667 TCN 3.933333333 TCN 3.8 2nd TSTPlus 8.733333333 TSTPlus 8.066666667 TSTPlus 8.066666667 3rd XCMPlus 8.933333333 XCMPlus 9.066666667 XCMPlus 9.4 4th XCM 11.93333333 XCM 10.8 XCM 10.06666667 5th RNN_FCNPlus 12.06666667 RNN_FCN 12.26666667 RNN 12.13333333 6th RNN_FCN 12.2 RNNPlus 12.8 RNNPlus 12.46666667 7th GRU_FCNPlus 12.6 RNN_FCNPlus 12.86666667 RNN_FCN 12.66666667 8th LSTM_FCNPlus 13 GRU_FCNPlus 13 GRU_FCNPlus 12.86666667 9th RNN 13.26666667 RNN 13.06666667 RNN_FCNPlus 13.06666667 10th FCN 13.4 LSTMPlus 13.66666667 GRU_FCN 14.06666667 42 Entropy 2023,25, 219 Table A3. TFriedman results: Algorithms—Shift 14. MAE MAPE R2 Method Friedman Score Method Friedman Score Method Friedman Score 1st TCN 6 TCN 8.2 TCN 25.33333333 2nd TSTPlus 8 TSTPlus 9.466666667 TST 21.6 3rd XCMPlus 10.4 RNN 9.533333333 TSTPlus 20.8 4th XCM 11.33333333 RNNPlus 9.666666667 XceptionTime 18.9 5th RNNPlus 11.86666667 XCMPlus 10.6 XCMPlus 17.96666667 6th LSTMPlus 11.93333333 LSTM 10.86666667 XceptionTimePlus 17.7 7th LSTM 12.13333333 XCM 11 RNNPlus 17.53333333 8th RNN 12.46666667 LSTMPlus 11.73333333 OmniScaleCNN 17.13333333 9th LSTM_FCNPlus 13.46666667 GRUPlus 13.13333333 LSTM 17.03333333 10th GRU_FCN 14.46666667 TST 13.8 RNN 16.93333333 MSE RMSE RMSLE Method Friedman Score Method Friedman Score Method Friedman Score 1st TCN 7.8 TCN 7.133333333 TCN 4.133333333 2nd TSTPlus 7.8 TSTPlus 7.6 TSTPlus 7.466666667 3rd XCM 10.13333333 XCM 10.46666667 XCMPlus 10 4th XCMPlus 10.6 XCMPlus 10.86666667 XCM 10.2 5th RNNPlus 10.86666667 RNNPlus 10.86666667 RNNPlus 10.73333333 6th LSTM 11.6 LSTMPlus 11.93333333 RNN 11.73333333 7th LSTMPlus 11.86666667 LSTM 12.06666667 LSTM 13.26666667 8th RNN 12.53333333 RNN 12.4 LSTMPlus 13.33333333 9th LSTM_FCNPlus 13.53333333 LSTM_FCNPlus 13.46666667 LSTM_FCNPlus 13.8 10th FCN 14.93333333 TST 14.73333333 InceptionTime 14.13333333 Appendix B Appendix B.1 Please use the abbreviation table below to read the corresponding results of the Friedman Ranks. Table A4. Feature Setups and Abbreviations. No. Abbreviation Feature Setup 1 U Univariate 2 B Blob 3 V Vader 4 F Finbert 5 RM7C Rolling Mean 7 Closing Value 6 RM14C Rolling Mean 14 Closing Value 7 RM7B Rolling Mean 7 Blob 8 RM14B Rolling Mean 14 Blob 9 RM7V Rolling Mean 7 Vader 10 RM14V Rolling Mean 14 Vader 11 RM7F Rolling Mean 7 Finbert 12 RM14F Rolling Mean 14 Finbert 43 Entropy 2023,25, 219 Appendix B.2 Table A5. Friedman results: feature setups—Shift 1. MAE MAPE R2 Feature Setup Friedman Score Feature Setup Friedman Score Feature Setup Friedman Score 1st B_RM7B 19.73333 V_F 19.2 U 54.93333 2nd RM7C_F 24.13333 B_RM7B 20.06667 RM7F 50.53333 3rd RM7F 24.53333 B_V 21.8 RM14C 47.8 4th V_F 25.33333 RM7F 24.8 RM7C_RM7F 47.73333 5th RM7C_B 26.8 RM7C_F 26.93333 RM7C 47.13333 6th B_V 27.26667 F_RM14V 27.2 RM14F 46.4 7th B 28.4 RM7B_RM14V 28.13333 RM7C_B 45.86667 8th RM7F_RM14F 28.4 RM7B_RM14F 28.2 RM7C_F 43.6 9th RM14F 30 RM7C_RM14B 28.26667 B 43.4 10th B_RM14V 30 B_RM14V 29.2 RM7C_RM14C 43.13333 MSE RMSE RMSLE Feature Setup Friedman Score Feature Setup Friedman Score Feature Setup Friedman Score 1st B_RM7B 21 B_RM7B 20.6 B_RM7B 20.93333 2nd V_F 21.06667 RM7F 24.13333 RM7C_F 21.86667 3rd B_V 22.66667 RM7C_F 24.2 B 24.8 4th RM7C_F 24.6 RM7C_B 26.13333 V_F 25.66667 5th RM7F 25.4 V_F 26.26667 RM7C_B 26.26667 6th B_RM14V 27.13333 B_V 28.26667 RM7C_RM14B 26.4 7th F_RM14V 27.73333 B 28.6 U 26.4 8th RM7B_RM14V 28.2 RM7F_RM14F 28.86667 RM7F 26.8 9th B_RM7F 28.4 RM7C_RM7B 29.06667 RM7C 27.53333 10th V_RM7V 28.6 RM7C_RM14B 29.13333 B_V 29 Table A6. Friedman results: feature setups—Shift 7. MAE MAPE R2 Feature Setup Friedman Score Feature Setup Friedman Score Feature Setup Friedman Score 1st B_RM7B 21.53333 V 22.33333 U 55.06667 2nd V 22.8 RM14F 24.8 RM14C 54.93333 3rd RM7B 24.46667 B_RM7B 25.2 RM7C 54 4th U 24.86667 V_RM7V 25.93333 RM7C_RM7V 49.33333 5th RM7V 25.33333 RM7B 26.33333 RM7C_RM14C 47.23333 6th RM14F 25.53333 RM7C_B 26.66667 RM14C_RM14F 45.46667 7th RM7C_B 25.86667 RM7C_RM14F 26.8 RM14C_B 45.33333 8th RM7C_RM7F 26.66667 RM7F_RM14F 27.53333 RM7C_RM14F 45.26667 9th RM7F 26.86667 U 27.73333 RM14F 45.13333 10th B 27.13333 V_RM14V 27.93333 RM14C_RM7V 44.93333 MSE RMSE RMSLE Feature Setup Friedman Score Feature Setup Friedman Score Feature Setup Friedman Score 1st V 21.93333 U 22.73333 U 17.6 2nd B_RM7B 23.13333 V 23.13333 RM7C_RM14F 21.13333 3rd RM7V 24 B_RM7B 23.66667 B_RM7B 22.53333 4th V_RM7V 24.8 RM7B 24.46667 RM14F 22.73333 5th RM7B 25.53333 RM14F 24.93333 V 23.2 6th RM14F 26.73333 RM7V 25.4 RM7F 24.46667 7th U 27.2 RM7C_RM7F 25.8 RM7C_RM7F 25.73333 8th RM7C_RM7F 27.73333 RM7C_B 25.93333 RM7C_B 26.4 9th B 27.86667 RM7C_RM14F 26.6 RM14C_RM7B 26.73333 10th RM7C_B 27.93333 RM7F 27.2 B 26.86667 44 Entropy 2023,25, 219 Table A7. Friedman results: feature setups—Shift 14. MAE MAPE R2 Feature Setup Friedman Score Feature Setup Friedman Score Feature Setup Friedman Score 1st RM7C_B 17.46667 RM7C_B 18.53333 RM14C 50.5 2nd RM7C_V 21.53333 RM14C_B 22.2 RM7C 48.9 3rd RM14C_B 22.26667 RM7C_V 22.93333 U 48.76667 4th RM7C 22.26667 RM7F_RM14F 23.73333 RM14C_RM7F 48.16667 5th U 23.73333 B_RM7V 24.2 RM14C_RM7B 46.83333 6th B_RM7V 23.8 V 25.6 RM7C_RM7B 46.06667 7th V 24.46667 RM7C_F 26.33333 RM7C_RM7F 45.76667 8th RM7C_F 24.86667 V_RM14F 27.4 RM7C_RM14C 45.56667 9th RM7B 25.86667 RM7C 27.66667 RM7F 43.83333 10th RM14B 27.33333 RM7B 28.13333 RM14C_F 43.36667 MSE RMSE RMSLE Feature Setup Friedman Score Feature Setup Friedman Score Feature Setup Friedman Score 1st RM7C_B 18.26667 RM7C_B 15.86667 RM7C 13.8 2nd RM7C_V 21.26667 RM7C 20.33333 RM7C_B 16 3rd B_RM7V 21.6 RM14C_B 21.26667 RM14C_B 21.06667 4th RM14C_B 23.86667 RM7C_V 21.4 U 22.33333 5th V 25.06667 U 22.8 RM7C_F 24.13333 6th RM7B 25.86667 RM7C_F 24 RM14B 24.33333 7th RM7C 26.2 V 24.06667 B_RM7V 25.06667 8th RM7C_F 26.26667 B_RM7V 24.46667 RM7C_V 26.06667 9th RM7B_RM7F 26.33333 RM7B 26 V 26.26667 10th RM7B_RM7V 26.66667 RM7C_RM14C 26.2 RM7B 26.46667 References 1. Basak, S.; Kar, S.; Saha, S.; Khaidem, L.; Dey, S.R. Predicting the direction of stock market prices using tree-based classifiers. N. Am. J. Econ. Financ. 2019,47, 552–567. [CrossRef] 2. Ren, R.; Wu, D.D.; Liu, T. Forecasting Stock Market Movement Direction Using Sentiment Analysis and Support Vector Machine. IEEE Syst. J. 2019,13, 760–770. [CrossRef] 3. Huang, W.; Nakamori, Y.; Wang, S.Y. Forecasting stock market movement direction with support vector machine. Comput. Oper. Res. 2005,32, 2513–2522. [CrossRef] 4. Zhong, X.; Enke, D. Predicting the daily return direction of the stock market using hybrid machine learning algorithms. Financ. Innov. 2019,5, 24. [CrossRef] 5. Abraham, B.; Ledolter, J. (Eds.) Statistical Methods for Forecasting; Wiley Series in Probability and Statistics, John Wiley & Sons, Inc.: Hoboken, NJ, USA, 1983. [CrossRef] 6. Armstrong, J.S.; Collopy, F.L. Integration of Statistical Methods and Judgment for Time Series Forecasting: Principles from Empirical Research. Forecast. Model. eJournal 1998, 269–293. 7. Bontempi, G.; Ben Taieb, S.; Le Borgne, Y.A. Machine Learning Strategies for Time Series Forecasting; Springer: Berlin/Heidelberg, Germany, 2013; Volume 138. [CrossRef] 8. Masini, R.P.; Medeiros, M.C.; Mendes, E.F. Machine learning advances for time series forecasting. J. Econ. Surv. 2021 ,37, 76–111. [CrossRef] 9. Cao, L.; Tay, F. Support vector machine with adaptive parameters in financial time series forecasting. IEEE Trans. Neural Netw. 2003,14, 1506–1518. [CrossRef] [PubMed] 10. Yang, A.; Li, W.; Yang, X. Short-term electricity load forecasting based on feature selection and Least Squares Support Vector Machines. Knowl.-Based Syst. 2019,163, 159–173. [CrossRef] 11. Sagheer, A.; Kotb, M. Time series forecasting of petroleum production using deep LSTM recurrent networks. Neurocomputing 2019,323, 203–213. [CrossRef] 12. Zhao, Z.; Chen, W.; Wu, X.; Chen, P.C.Y.; Liu, J. LSTM network: A deep learning approach for short-term traffic forecast. IET Intell. Transp. Syst. 2017,11, 68–75. 13. Graf, R.; Zhu, S.; Sivakumar, B. Forecasting river water temperature time series using a wavelet–neural network hybrid modelling approach. J. Hydrol. 2019,578, 124115. [CrossRef] 14. Kurumatani, K. Time series forecasting of agricultural product prices based on recurrent neural networks and its evaluation method. SN Appl. Sci. 2020,2, 1434. 45 Entropy 2023,25, 219 15. Khairalla, M.A.E.; Ning, X.; Al-Jallad, N.T.; El-Faroug, M.O. Short-Term Forecasting for Energy Consumption through Stacking Heterogeneous Ensemble Learning Model. Energies 2018,11, 1605. [CrossRef] 16. Alkandari, M.; Ahmad, I. Solar power generation forecasting using ensemble approach based on deep learning and statistical methods. Appl. Comput. Inform. 2020. [CrossRef] 17. Liapis, C.M.; Karanikola, A.; Kotsiantis, S.B. Energy Load Forecasting: Investigating Mid-Term Predictions with Ensemble Learners. In Proceedings of the AIAI, Crete, Greece, 17–20 June 2022. 18. Liapis, C.M.; Karanikola, A.C.; Kotsiantis, S.B. An ensemble forecasting method using univariate time series COVID-19 data. In Proceedings of the 24th Pan-Hellenic Conference on Informatics, Athens, Greece, 20–22 November 2020. 19. Liapis, C.M.; Karanikola, A.; Kotsiantis, S.B. A Multi-Method Survey on the Use of Sentiment Analysis in Multivariate Financial Time Series Forecasting. Entropy 2021,23, 1603. [CrossRef] [PubMed] 20. Siami-Namini, S.; Tavakoli, N.; Namin, A.S. A Comparison of ARIMA and LSTM in Forecasting Time Series. In Proceedings of the 2018 17th IEEE International Conference on Machine Learning and Applications (ICMLA), Orlando, FL, USA, 17–20 December 2018; pp. 1394–1401. 21. Çıbıkdiken, A.; Karakoyun, E. Comparison of ARIMA Time Series Model and LSTM Deep Learning Algorithm for Bitcoin Price Forecasting. In Proceedings of the 13th multidisciplinary academic conference in Prague, Hamburg, Germany, 27–30 August 2018 . 22. Yamak, P.T.; Yujian, L.; Gadosey, P.K. A Comparison between ARIMA, LSTM, and GRU for Time Series Forecasting. In Proceedings of the ACAI, Sanya, China, 20–22 December 2019. 23. Maleki, A.; Nasseri, S.; Aminabad, M.S.; Hadi, M. Comparison of ARIMA and NNAR Models for Forecasting Water Treatment Plant’s Influent Characteristics. KSCE J. Civ. Eng. 2018,22, 3233–3245. 24. Satrio, C.B.A.; Darmawan, W.; Nadia, B.U.; Hanafiah, N. Time series analysis and forecasting of coronavirus disease in Indonesia using ARIMA model and PROPHET. Procedia Comput. Sci. 2021,179, 524–532. 25. Paliari, I.; Karanikola, A.; Kotsiantis, S.B. A comparison of the optimized LSTM, XGBOOST and ARIMA in Time Series forecasting. In Proceedings of the 2021 12th International Conference on Information, Intelligence, Systems & Applications (IISA), Chania Crete, Greece, 12–14 July 2021; pp. 1–7. 26. Zhang, Y.; Yang, H.L.; Cui, H.; Chen, Q. Comparison of the Ability of ARIMA, WNN and SVM Models for Drought Forecasting in the Sanjiang Plain, China. Nat. Resour. Res. 2019,29, 1447–1464. 27. Tealab, A. Time series forecasting using artificial neural networks methodologies: A systematic review. Future Comput. Inform. J. 2018,3, 334–340. [CrossRef] 28. Sezer, O.B.; Gudelek, M.U.; Ozbayoglu, A.M. Financial Time Series Forecasting with Deep Learning: A Systematic Literature Review: 2005–2019. arXiv 2020, arXiv:abs/1911.13288. 29. Lara-Benítez, P.; Carranza-García, M.; Santos, J.C.R. An Experimental Review on Deep Learning Architectures for Time Series Forecasting. Int. J. Neural Syst. 2021,31, 2130001. [CrossRef] [PubMed] 30. Karanikola, A.; Liapis, C.M.; Kotsiantis, S. A Comparison of Contemporary Methods on Univariate Time Series Forecasting. In Advances in Machine Learning/Deep Learning-Based Technologies: Selected Papers in Honour of Professor Nikolaos G. Bourbakis—Volume 2; Tsihrintzis, G.A., Virvou, M., Jain, L.C., Eds.; Springer International Publishing: Cham, Switzerland, 2022; pp. 143–168. [CrossRef] 31. Wang, K.; Qi, X.; Liu, H. A comparison of day-ahead photovoltaic power forecasting models based on deep learning neural network. Appl. Energy 2019,251, 113315. 32. Rao, T.; Srivastava, S. Analyzing Stock Market Movements Using Twitter Sentiment Analysis. In Proceedings of the International Conference on Advances in Social Networks Analysis and Mining, Sydney, Australia, 31 July–3 August 2012. 33. Nguyen, T.H.; Shirai, K.; Velcin, J. Sentiment analysis on social media for stock movement prediction. Expert Syst. Appl. 2015 , 42, 9603–9611. 34. Kalyani, J.; Bharathi, H.N.; Jyothi, R. Stock trend prediction using news sentiment analysis. arXiv 2016, arXiv:abs/1607.01958. 35. Shah, D.; Isah, H.; Zulkernine, F.H. Predicting the Effects of News Sentiments on the Stock Market. In Proceedings of the 2018 IEEE International Conference on Big Data (Big Data), Seattle, WA, USA, 10–13 December 2018; pp. 4705–4708. 36. Souma, W.; Vodenska, I.; Aoyama, H. Enhanced news sentiment analysis using deep learning methods. J. Comput. Soc. Sci. 2019 , 2, 33–46. 37. Valle-Cruz, D.; Fernandez-Cortez, V.; Chau, A.L.; Sandoval-Almazán, R. Does Twitter Affect Stock Market Decisions? Financial Sentiment Analysis During Pandemics: A Comparative Study of the H1N1 and the COVID-19 Periods. Cogn. Comput. 2021 , 14, 372–387. 38. Sharma, V.; Khemnar, R.K.; Kumari, R.A.; Mohan, B.R. Time Series with Sentiment Analysis for Stock Price Prediction. In Proceedings of the 2019 2nd International Conference on Intelligent Communication and Computational Techniques (ICCT), Jaipur, India, 28–29 September 2019; pp. 178–181. 39. Pai, P.F.; Liu, C. Predicting Vehicle Sales by Sentiment Analysis of Twitter Data and Stock Market Values. IEEE Access 2018 , 6, 57655–57662. 40. Mohan, S.; Mullapudi, S.; Sammeta, S.; Vijayvergia, P.; Anastasiu, D. Stock Price Prediction Using News Sentiment Analysis. In Proceedings of the 2019 IEEE Fifth International Conference on Big Data Computing Service and Applications (BigDataService), Newark, CA, USA, 4–9 April 2019; pp. 205–208. 41. Mehta, P.; Pandya, S.; Kotecha, K. Harvesting social media sentiment analysis to enhance stock market prediction using deep learning. PeerJ Comput. Sci. 2021,7, e476. [CrossRef] 46 Entropy 2023,25, 219 42. Jin, Z.; Yang, Y.; Liu, Y. Stock closing price prediction based on sentiment analysis and LSTM. Neural Comput. Appl. 2019 , 32, 9713–9729. [CrossRef] 43. Wu, S.H.; Liu, Y.; Zou, Z.; Weng, T.H. S_I_LSTM: Stock price prediction based on multiple data sources and sentiment analysis. Connect. Sci. 2021,34, 44–62. 44. Jing, N.; Wu, Z.; Wang, H. A hybrid model integrating deep learning with investor sentiment analysis for stock price prediction. Expert Syst. Appl. 2021,178, 115019. 45. Smailovic, J.; Grcar, M.; Lavra, N.; Znidarsic, M. Stream-based active learning for sentiment analysis in the financial domain. Inf. Sci. 2014,285, 181–203. 46. Raju, S.M.; Tarif, A.M. Real-Time Prediction of BITCOIN Price using Machine Learning Techniques and Public Sentiment Analysis. arXiv 2020, arXiv:abs/2006.14473. 47. Abraham, J.; Higdon, D.W.; Nelson, J.; Ibarra, J. Cryptocurrency Price Prediction Using Tweet Volumes and Sentiment Analysis. SMU Data Sci. Rev. 2018,1,1. 48. Valencia, F.; Gómez-Espinosa, A.; Valdés-Aguirre, B. Price Movement Prediction of Cryptocurrencies Using Sentiment Analysis and Machine Learning. Entropy 2019,21, 589. [PubMed] 49. Deb, A.; Lerman, K.; Ferrara, E. Predicting Cyber Events by Leveraging Hacker Sentiment. Information 2018,9, 280. [CrossRef] 50. Masri, S.; Jia, J.; Li, C.; Zhou, G.; Lee, M.C.; Yan, G.; Wu, J. Use of Twitter data to improve Zika virus surveillance in the United States during the 2016 epidemic. BMC Public Health 2019,19, 761. 51. Chauhan, P.; Sharma, N.; Sikka, G. The emergence of social media data and sentiment analysis in election prediction. J. Ambient. Intell. Humaniz. Comput. 2021,12, 2601–2627. [CrossRef] 52. Tseng, K.K.; Lin, R.F.Y.; Zhou, H.; Kurniajaya, K.J.; Li, Q. Price prediction of e-commerce products through Internet sentiment analysis. Electron. Commer. Res. 2018,18, 65–88. [CrossRef] 53. Twintproject. Twintproject/Twint: An Advanced Twitter Scraping & OSINT Tool. Available online: https://github.com/ twintproject/twint (accessed on 7 October 2021). 54. Van Rossum, G. The Python Library Reference, Release 3.8.2; Python Software Foundation: Wolfeboro Falls, NH, USA, 2020. 55. Bird, S. NLTK: The Natural Language Toolkit. arXiv 2004, arXiv:cs.CL/0205028. 56. Bird, S.; Klein, E.; Loper, E. Natural Language Processing with Python; Packt Publishing Ltd.: Birmingham, UK, 2009. 57. String—Common String Operations. Available online: https://docs.python.org/3/library/string.html (accessed on 7 October 2021 ). 58. Simplified Text Processing. Available online: https://textblob.readthedocs.io/en/dev/ (accessed on 7 October 2021). 59. Hutto, C.J.; Gilbert, E. VADER: A Parsimonious Rule-Based Model for Sentiment Analysis of Social Media Text. In Proceedings of the International AAAI Conference on Web and Social Media, Ann Arbor, MI, USA, 1–4 June 2014. 60. Araci, D. FinBERT: Financial Sentiment Analysis with Pre-trained Language Models. arXiv 2019, arXiv:abs/1908.10063. 61. Devlin, J.; Chang, M.W.; Lee, K.; Toutanova, K. BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding. arXiv 2019, arXiv:abs/1810.04805. 62. ProsusAI. ProsusAI/finBERT: Financial Sentiment Analysis with Bert. Available online: https://github.com/ProsusAI/finBERT (accessed on 7 October 2021). 63. Malo, P.; Sinha, A.; Korhonen, P.J.; Wallenius, J.; Takala, P. Good debt or bad debt: Detecting semantic orientations in economic texts. J. Assoc. Inf. Sci. Technol. 2014,65, 782–796. [CrossRef] 64. timeseriesAI. Timeseriesai/Tsai: Time Series Timeseries Deep Learning Machine Learning Pytorch FASTAI: State-of-the-Art Deep Learning Library for Time Series and Sequences in Pytorch/Fastai. Available online: https://github.com/timeseriesAI/tsai (accessed on 7 October 2021). 65. Wang, Z.; Yan, W.; Oates, T. Time series classification from scratch with deep neural networks: A strong baseline. In Proceedings of the 2017 International Joint Conference on Neural Networks (IJCNN), Anchorage, AK, USA, 14–19 May 2017; pp. 1578–1585. 66. Oguiza, I. tsAI Models: FCNPlus. Available online: https://timeseriesai.github.io/tsai/models.fcnplus.html (accessed on 7 October 2021). 67. Fawaz, H.I.; Lucas, B.; Forestier, G.; Pelletier, C.; Schmidt, D.F.; Weber, J.; Webb, G.I.; Idoumghar, L.; Muller, P.A.; Petitjean, F. InceptionTime: Finding AlexNet for Time Series Classification. arXiv 2020, arXiv:abs/1909.04939. 68. Oguiza, I. tsAI Models: InceptionTimePlus. Available online: https://timeseriesai.github.io/tsai/models.inceptiontimeplus.html (accessed on 7 October 2021). 69. Oguiza, I. tsAI Models: RNNS. Available online: https://timeseriesai.github.io/tsai/models.rnn.html (accessed on 7 November 2022). 70. Hochreiter, S.; Schmidhuber, J. Long Short-Term Memory. Neural Comput. 1997,9, 1735–1780. [CrossRef] [PubMed] 71. Chung, J.; Gülçehre, C.; Cho, K.; Bengio, Y. Empirical Evaluation of Gated Recurrent Neural Networks on Sequence Modeling. arXiv 2014, arXiv:abs/1412.3555. 72. Oguiza, I. tsAI Models: RNN_FCN. Available online: https://timeseriesai.github.io/tsai/models.rnn_fcn.html (accessed on 7 November 2022). 73. Karim, F.; Majumdar, S.; Darabi, H.; Chen, S. LSTM Fully Convolutional Networks for Time Series Classification. IEEE Access 2018,6, 1662–1669. [CrossRef] 74. Elsayed, N.; Maida, A.; Bayoumi, M.A. Deep Gated Recurrent and Convolutional Network Hybrid Model for Univariate Time Series Classification. arXiv 2019, arXiv:abs/1812.07683. 47 Entropy 2023,25, 219 75. Oguiza, I. tsAI Models: RNN_FCNPlus. Available online: https://timeseriesai.github.io/tsai/models.rnn_fcnplus.html (accessed on 7 November 2022). 76. Zou, X.; Wang, Z.; Li, Q.; Sheng, W. Integration of residual network and convolutional neural network along with various activation functions and global pooling for time series classification. Neurocomputing 2019,367, 39–45. [CrossRef] 77. Oguiza, I. tsAI Models: ResNetPlus. Available online: https://timeseriesai.github.io/tsai/models.resnetplus.html (accessed on 7 November 2022). 78. Bai, S.; Kolter, J.Z.; Koltun, V. An Empirical Evaluation of Generic Convolutional and Recurrent Networks for Sequence Modeling. arXiv 2018, arXiv:abs/1803.01271. 79. Zerveas, G.; Jayaraman, S.; Patel, D.; Bhamidipaty, A.; Eickhoff, C. A Transformer-based Framework for Multivariate Time Series Representation Learning. In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, Singapore, 14–18 August 2021. 80. Oguiza, I. tsAI Models: TSTPlus. Available online: https://timeseriesai.github.io/tsai/models.tstplus.html (accessed on 7 November 2022). 81. Oguiza, I. tsAI Models: TSIT. Available online: https://timeseriesai.github.io/tsai/models.tsitplus.html (accessed on 7 November 2022). 82. Oguiza, I. tsAI Models: Transformermodel. Available online: https://timeseriesai.github.io/tsai/models.transformermodel.html (accessed on 7 November 2022). 83. Fauvel, K.; Lin, T.; Masson, V.; Fromont, E.; Termier, A. XCM: An Explainable Convolutional Neural Network for Multivariate Time Series Classification. arXiv 2021, arXiv:abs/2009.04796. 84. Oguiza, I. tsAI Models: XCMPlus. Available online: https://timeseriesai.github.io/tsai/models.xcmplus.html (accessed on 7 November 2022). 85. Rahimian, E.; Zabihi, S.; Atashzar, S.F.; Asif, A.; Mohammadi, A. XceptionTime: A Novel Deep Architecture based on Depthwise Separable Convolutions for Hand Gesture Classification. arXiv 2019, arXiv:abs/1911.03803. 86. Oguiza, I. tsAI Models: XceptionTimePlus. Available online: https://timeseriesai.github.io/tsai/models.xceptiontimeplus.html (accessed on 7 November 2022). 87. Tang, W.; Long, G.; Liu, L.; Zhou, T.; Blumenstein, M.; Jiang, J. Omni-Scale CNNs: A simple and effective kernel size configuration for time series classification. arXiv 2022, arXiv:2002.10061. 88. Friedman, M. The Use of Ranks to Avoid the Assumption of Normality Implicit in the Analysis of Variance. J. Am. Stat. Assoc. 1937,32, 675–701. [CrossRef] 89. Dunn, O.J. Multiple Comparisons among Means. J. Am. Stat. Assoc. 1961,56, 52–64. 90. Hodges, J.L.; Lehmann, E.L. Rank Methods for Combination of Independent Experiments in Analysis of Variance. Ann. Math. Stat. 1962,33, 403–418. [CrossRef] Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. 48 Citation: Wada, T.; Takayasu, H.; Takayasu, M. Extraction of Important Factors in a High-Dimensional Data Space: An Application for High-Growth Firms. Entropy 2023,25, 488. https://doi.org/10.3390/ e25030488 Academic Editor: Panos Argyrakis Received: 6 February 2023 Revised: 2 March 2023 Accepted: 8 March 2023 Published: 10 March 2023 Copyright: © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). entropy Article Extraction of Important Factors in a High-Dimensional Data Space: An Application for High-Growth Firms Takuya Wada 1, Hideki Takayasu 2,3 and Misako Takayasu 1,2,* 1Department of Mathematical and Computing Science, School of Computing, Tokyo Institute of Technology, Yokohama 226-8502, Japan 2Institute of Innovative Research, Tokyo Institute of Technology, Yokohama 226-8502, Japan; [email protected].co.jp 3Sony Computer Science Laboratories, Tokyo 141-0022, Japan *Correspondence: [email protected] Abstract: We introduce a new non-black-box method of extracting multiple areas in a high-dimensional big data space where data points that satisfy specific conditions are highly concentrated. First, we extract one-dimensional areas where the data that satisfy specific conditions are mostly gathered by using the Bayesian method. Second, we construct higher-dimensional areas where the densities of focused data points are higher than the simple combination of the results for one dimension, and then we verify the results through data validation. Third, we apply this method to estimate the set of significant factors shared in successful firms with growth rates in sales at the top 1% level using 156-dimensional data of corporate financial reports for 12 years containing about 320,000 firms. We also categorize high-growth firms into 15 groups of different sets of factors. Keywords: variable selection; feature selection; high-growth firms; Bayesian method; big data 1. Introduction We consider the general problem of extracting areas in a high-dimensional data space where points that satisfy specific conditions are concentrated. Generally, as factors associated with a specific condition are often unknown, we use the most available factors and examine their relevance to a particular condition [ 1 ]. However, the majority of the factors used are irrelevant or redundant, resulting in problems such as reduced accuracy of the analysis and increased analysis time [ 1 , 2 ]. Therefore, we are reducing the number of variables, a process called variable selection. Variable selection has various advantages, such as accuracy increase, analysis time reduction, and overfitting avoidance [ 2 – 4 ]. Many models have been proposed for this variable selection and used in various fields [ 4 – 6 ]. In recent years, machine learning models have been used to improve the accuracy of variable selection. For example, Genuer used random forests [ 7 ] to select significant variables in high-dimensional classification problems [ 8 ]. Grandvalet proposed a model that automatically performs relevance judgments and feature selection on support vector machines [ 9 ] and showed its effectiveness in facial expression recognition tasks [ 10 ]. However, machine learning models also have disadvantages; for example, generally their results are difficult to understand logically due to the complexity of these models and their black-box structure [ 11 , 12 ]. In addition, to the best of our knowledge, a general method for exhaustively extracting areas where the data that satisfy specific conditions are highly concentrated has not been established in the study of big data. In this paper, we propose a new method based on a non-black-box model to solve this general problem. We use indicators calculated using the Bayesian method and Szymkiewicz- Simpson coefficient as evaluation measures for variable selection and extraction of pairs of variables, respectively. The Bayesian method is a data analysis method that uses existing information [ 13 , 14 ]. This point differs from the likelihood method and gives the advantage Entropy 2023,25, 488. https://doi.org/10.3390/e25030488 https://www.mdpi.com/journal/entropy 49 Entropy 2023,25, 488 the 54 financial items. The distribution of the existence probability of high-growth firms and details of the areas extracted for one example of those financial items are presented in Figure 3 and Table 3. Figure 3. Existence probability of high-growth firms in each of the segmented areas, projected on the axis of the ratio of net income to sales (before amortization and after tax, %). The horizontal axis is the quantile from the beginning to the end of the segmented area, and the vertical axis is the existence probability of high-growth firms within the segmented area. The red dashed line represents 0.01, the percentage of high-growth firms in the overall area. For this financial item, the orange and green areas were extracted as the areas with densely populated high-growth firms, and the blue area was not extracted because it was not densely populated with high-growth firms. For the orange area, two areas were initially extracted: the 0–5.0% and 15.0–28.9% areas. These two areas and the areas in between where the existence probability of high-growth firms is low were merged into one area, as shown in Figure 2b. Table 3. Two areas extracted in the ratio of net income to sales (before amortization and after tax, %), orange and green, respectively, in Figure 3. The extracted areas are from lower to upper limits, which are denoted by percentage points within the financial item. The existence probability of high-growth firms of an area is calculated using the number of high-growth firms in the area, the number of all firms in the area, and Equation (5). The abbreviated names used in this table are defined in Table 1. Area Lower Limit Upper Limit NHF NF EPHF orange 0.0% 28.9% 7203 417805 0.017 green 94.6% 100.0% 1439 77645 0.017 These orange and green areas are where high-growth firms are about 1.7 times more dense than normal ones. These areas are the two edges of the financial items, and it is thought that firms grow high due to different factors. For validation, we performed the same one-dimensional extraction on five random data. We extracted 11, 11, 12, 13, and 13 areas, respectively. No multiple areas were extracted within a single financial item. The area with the highest existence of high-growth firms in these areas was about 1.08 times more dense than normal ones. These areas are used in Step2. 4.2. Reduction of Areas Containing Similar Data Points Similar areas were deleted in Step2 for the 197 areas of 143 financial items extracted in Step1. The result of calculating Equation (9) for all combinations of the 197 areas is presented in Figure 4. 56 Entropy 2023,25, 488 Figure 4. Cumulative distribution function of the values calculated for all combinations by using Equation (9). The horizontal axis is the value of the Szymkiewicz–Simpson coefficient, and the vertical axis is the cumulative distribution. The red dashed line represents 0.825, where the shape of the cumulative distribution function changes. This value was used as the threshold value. Figure 4 shows that the cumulative distribution function changes its slope around when the value of the Szymkiewicz–Simpson coefficient is 0.825. This value was used as the threshold value. In the combination of areas where the value of the Szymkiewicz–Simpson coefficient is greater than this value, the area with the smallest existence probability of high-growth firms was deleted. For example, the combination of an area with a turnover of current debt (months) of 7.44 or higher and an area with an increase/decrease in an investment of less than 0 (thousands of yen) resulted in a Szymkiewicz–Simpson coefficient value of 0.916. Therefore, we compared the existence probability of high-growth firms and removed the area with an investment volume of less than 0 (thousands of yen), which was a lower area. We finally extracted 67 areas of 51 financial items. For the five random data, the highest Szymkiewicz–Simpson coefficient was about 0.24 in the combination obtained from the areas of financial items obtained in each. Considering that this is smaller than the threshold value of 0.825 in the data for analysis and that no similarity exists among the financial items and among the areas as the data were randomly shuffled, none of the areas were removed. The 11, 11, 12, 13, and 13 areas obtained in Step1 were used in Step3–Step5. 4.3. Extraction of Two-Dimensional Areas The 67 areas of 51 financial items extracted in Step2 were used to extract the twodimensional areas. We checked all possible combinations, and the top five two-dimensional areas with the highest existence probability of high-growth firms are presented in Table 4. In the two-dimensional area where the existence of high-growth firms is in the first and second places, high-growth firms are about 20 times more dense than normal ones. Table 4 displays how many times the existence of high-growth firms is compared to when the two conditions are independent (Column Ratio), and these five areas are about five times as high. Therefore, some synergy must exist in the combination of these conditions. Figure 5 presents the extracted two-dimensional area of the first rank. 57 Entropy 2023,25, 488 Table 4. Top five two-dimensional areas in Step3. The existence probability of high-growth firms of an area is calculated using the number of high-growth firms in the area, the number of all firms in the area, and Equation (5). The ratio in this table is the existence probability of high-growth firms in two dimensions divided by the existence probability of high-growth firms calculated given that the two conditions are independent using Equation (8). The abbreviated names used in this table are defined in Table 1. Item Name EPHF(1D) Item Name EPHF(1D) EPHF(2D) Ratio CLR (M) 0.036 CACL (%) 0.013 0.199 5.142 CLR (M) 0.036 LACL (%) 0.012 0.196 5.232 CSGR (%) 0.025 GPE (T) 0.020 0.177 5.147 TCR (M) 0.022 FAR (M) 0.019 0.174 5.647 FAFL (%) 0.024 FAR (M) 0.023 0.173 4.702 Figure 5. Extracted two-dimensional area of the first rank. The vertical and horizontal axes are divided by the current liabilities to revenue ratio (months) and the current assets to current liabilities ratio (%), respectively. The size of the circle represents the number of firms in the area, and the radius is scaled in a logarithmic scale. The colors of the circles represent the proportion of high-growth firms in the area. It is drawn in the order of yellow, orange, red, brown, and black, starting from the lowest to the highest. The green box at the right top is the area extracted as the two-dimensional area with the highest concentration of high-growth firms. The blue box at the left bottom is the area that was not extracted because the existence probability of high-growth firms in this area is lower than that of high-growth firms using Equation (8) if the two conditions are independent. The green box area at the right top in Figure 5 is the area that satisfies the green areas in the turnover of current debt and the current ratio in the one-dimensional axes. It is 20 times more densely populated with high-growth firms than normal ones. It was also extracted as a two-dimensional area with the highest existence probability of high-growth firms. Meanwhile, the blue box area in Figure 5 is the area that satisfies the orange areas 58 Entropy 2023,25, 488 in the turnover of the current debt and the current ratio in the one-dimensional axes. The existence probability of high-growth firms in this area is 0.014. This value is lower than that of high-growth firms when the two conditions are independent, as calculated using Equation (8). Therefore, this area was not extracted as a two-dimensional area. We obtain 2211 two-dimensional areas using the 67 conditions used for the 67 areas extracted in Step2. Among them, we extracted 1036 areas that are more densely concentrated with high-growth firms than that when the conditions were independent. For the five random data, we check whether high-growth firms are densely populated in the two-dimensional areas using the conditions used for the areas extracted in Step2. The number of areas extracted as areas where the existence probability of high-growth firms is higher than that of high-growth firms calculated using Equation (8), under the condition that the two conditions are independent were 3, 4, 6, 7, and 9. Even in the area with the highest concentration of high-growing firms in any of the random data, the concentration of highest-growing firms is about 1.7 times the normal concentration. It was also about 1.5 times higher than when all conditions were independent, indicating no strong synergistic effect. These two-dimensional areas extracted as densely populated with high-growth firms in the random data are used in the analysis in step 4. 4.4. Extraction of Higher-Dimensional Areas For the 1036 two-dimensional areas extracted in Step3, we extract 1036 high-dimensional areas by repeatedly adding the 67 conditions used in the 67 areas extracted in Step2. The top two high-dimensional areas that are extracted are presented in Tables 5 and 6. Table 5. Eight-dimensional area with the first highest existence probability of high-growth firms among the extracted high-dimensional areas. The ratio in this table is the existence probability of high-growth firms in the n-dimensional area divided by that of high-growth firms calculated under conditions where the n -conditions are independent using Equation (8); n is the number of conditions in the row (Column NC). The abbreviated names used in this table are defined in Table 1. NC Item Name (Threshold) NHF EPHF Ratio 1 NOLR (%) (≤0) 3715 0.031 1.000 2 IR (M) (≤0) 2447 0.066 1.276 3 CAR (M) (≥10.2) 664 0.161 2.182 4 OITC (IC) (≤2) 217 0.408 4.194 5 GPE (T) (≤2727) 112 0.530 4.984 6 FAR (M) (≥10.19) 40 0.673 5.699 7 APR (M) (≤0) 33 0.748 5.382 8 CLR (M) (≥7.44) 26 0.779 4.837 9 OIR (CFS) (≤2) 26 0.779 4.133 Table 6. Seven-dimensional area with the second highest existence probability of high-growth firms among the extracted high-dimensional areas. The ratio in this table is the existence probability of high-growth firms in the n-dimensional area divided by that of high-growth firms calculated under conditions where the n -conditions are independent using Equation (8); n is the number of conditions in the row (Column NC). The abbreviated names used in this table are defined in Table 1. NC Item Name (Threshold) NHF EPHF Ratio 1 ARR (DT) (M) (≤0.25) 4336 0.028 1.000 2 Revenue to total capital ratio (IC) (≤3) 1790 0.058 1.604 3 OITC (IC) (≤2) 550 0.215 3.548 4 PPER (M) (≤0.16) 136 0.475 6.658 5 NCLR (M) (≤0) 92 0.606 7.117 6 Investment and financing returns (%) (≤0.02) 60 0.683 7.272 7 DR (%) (≤0) 37 0.771 7.322 8 OIR (CFS) (≤2) 36 0.765 5.690 59 Entropy 2023,25, 488 The existence probability of high-growth firms decreased when the 8th and 9th conditions were added to the areas in Tables 5 and 6. Therefore, the areas with the 7th and 8th dimensions in Tables 5 and 6 were extracted as areas with a high concentration of highgrowth firms. The existence probability of high-growth firms in these high-dimensional areas is about 0.77. This implies that high-growth firms in these areas are about 77 times more dense than normal ones. They are also about 5–7 times higher than that when all conditions were independent. Therefore, we can assume that some synergistic effects occur in the combinations of these conditions. As shown in Tables 5 and 6, we extract the high-dimensional areas from the 1036 two-dimensional areas obtained in Step3. The distribution of the existence probability of high-growth firms in the high-dimensional areas finally obtained is presented in Figure 6. Figure 6. Distribution of the existence probability of high-growth firms in the high-dimensional areas. The vertical axis and horizontal axes are the proportion of 1036 areas and the existence probability of high-growth firms, respectively. The red line represents 0.01, the percentage of high-growth firms in the overall area. As shown in Figure 6, 90% of the 1036 high-dimensional areas were able to extract areas where the high-growth firms are dense at 30 times or higher than the normal density. We have also extracted four areas where the high-growth firms are dense at less than three times the normal density, and all of these areas were two-dimensional ones. Subsequently, areas with a small number of data are called local ones. These areas became localized at the two-dimensional level, and no further high-dimensional areas could be extracted. Our method searched the entire area exhaustively, and the extracted areas include the local ones. For the 1036 high-dimensional areas obtained in these processes, we verified whether the existence probability of high-growth firms is also increased in the data for validation. The verification procedure is to add conditions in the same order as the conditions for the areas obtained in these processes until the existence probability of high-growth firms stops to increase. As specific examples, the results of the verification in the areas of Tables 5 and 6 are presented in the Tables 7 and 8, respectively. In the validation for both areas, the existence probability of high-growth firms decreased when the 5th condition was added. Thus, we confirmed the robustness of the results up to the four-dimensional area in these areas. In this validation, the existence probability of high-growth firms in the one-dimensional area in both validation results was almost the same as that when the data for analysis were used. The existence probability of high-growth firms in the four-dimensional area when the data for verification were used was about 0.33 and 0.21 for Tables 7 and 8, respectively. Although these values are lower than when using the data for analysis, we can assume that high-growth firms are concentrated at a high density, which cannot be considered coincidental. The reason for the lower existence probability of high-growth firms in the four-dimensional area, com- 60 Entropy 2023,25, 488 pared to that for analysis, and the failure of these areas to maintain robustness in the five-dimensional area can be attributed to the fact that the data for verification are one-fifth the number of data for analysis. That is the number of high-growth firms in the area at the four-dimensional area is about 15.7% and 14.0% in Tables 7 and 8 for validation compared to that for analysis. Thus, the number of high-growth firms in the area is reduced, and the results are no longer stable and robust in high dimensions. The same verification was conducted for the remaining 1034 high-dimensional areas. The distribution of the number of dimensions for which the existence probability of high-growth firms was maximized in the data for analysis and verification was checked (Figure 7). Table 7. Validation result for the high-dimensional area of Table 5 with the highest existence probability of high-growth firms. We add conditions in the same order as in Table 5 until the existence probability of high-growth firms stops to increase. The abbreviated names used in this table are defined in Table 1. NC Item Name Threshold NHF EPHF 1 NOLR (%) ≤0 551 0.031 2 IR (M) ≤0 379 0.064 3 CAR (M) ≥0.122 106 0.074 4 OITC (IC) ≤2 34 0.334 5 GPE (T) ≤2727 10 0.204 Table 8. Validation result for the high-dimensional area of Table 6 with the second-highest existence probability of high-growth firms. We add conditions in the same order as in Table 6 until the existence probability of high-growth firms stops to increase. The abbreviated names used in this table are defined in Table 1. NC Item Name Threshold NHF EPHF 1 ARR (DT) (M) ≤0.25 680 0.029 2 Revenue to total capital ratio (IC) ≤3 233 0.046 3 OITC (IC) ≤2 79 0.183 4 PPER (M) ≤0.16 19 0.208 5 NCLR (M) ≤0 11 0.207 Figure 7. Distribution of the number of dimensions for which the existence probability of high-growth firms was maximized in the data for analysis and verification. The vertical and horizontal axes are the number of dimensions in verification data and analysis data, respectively. The numbers represent the number of areas with each dimension in the analysis and validation data. The colors indicate that the darker the red color, the higher the value, i.e., the greater the number of areas. 61 Entropy 2023,25, 488 The numbers in Figure 7 represent the number of areas with each dimension in the analysis and validation data. For example, 77 with a vertical axis of 4 and a horizontal axis of 7 indicates that 77 areas have been extracted in seven-dimensional areas for analysis and verified to four-dimensional areas. Specifically, the area in Table 5 is contained in 72 with a vertical axis of 4 and a horizontal axis of 8, and that in Table 6 is contained in 77 with a vertical axis of 4 and a horizontal axis of 7 in Figure 7. Figure 7 presents that many highdimensional areas of more than three dimensions are robust for verification. In addition, we can observe a relationship whereby the areas with higher dimensionality for analysis also maintain a higher dimensionality for validation. There was also a 10-dimensional area for which robustness was confirmed up to nine dimensions for verification. The details of this area are provided in Table 9. Table 9. Ten-dimensional area for which robustness was confirmed in up to nine dimensions for verification. The abbreviated names used in this table are defined in Table 1. NC Item Name (Threshold) NHF (DA) EPHF (DA) NHF (DV) EPHF (DV) 1 Revenue (T) (≤108,917) 9577 0.032 1253 0.041 2 NOLR (%) (≤0) 3156 0.056 450 0.065 3 CAR (M) (≥10.2) 1046 0.134 143 0.123 4 OITC (IC) (≤2) 338 0.278 51 0.234 5 PPER (M) (≤0.16) 137 0.471 25 0.317 6 LR (M) (≥14.13) 78 0.562 17 0.358 7 NCLR (M) (≤0) 63 0.672 11 0.371 8 Revenue to total capital ratio (IC) (≤3) 59 0.691 11 0.404 9 IR (M) (≤0) 37 0.694 10 0.464 10 LACL (%)(≤41.45) 18 0.700 3 0.144 11 OIR (CFS) (≤2) 18 0.700 The area in Table 9 is the area where the high-growth firms are about 70 times more densely populated than usual for the analysis. This area maintains robustness up to nine dimensions. In the data for verification, the high-growth firms are about 46 times denser than usual in this nine-dimensional area. We also extracted high-dimensional areas that can retain such robustness. There are 165 areas where the increase in the existence probability of high-growth firms stops at one-dimensional areas for validation, despite that for analysis they are high-dimensional areas with six or more dimensions. In addition, in about half of the 1036 high-dimensional areas, an increase in the existence probability of high-growth firms stopped at three dimensions or less in the data for verification. Therefore, our method exhaustively searches the entire range and extracts local areas. In the following, we focus on somewhat larger areas wherein the number of highgrowth firms includes more than 1% (145 firms) of the total number of high-growth firms in the four-dimensional area in the data for analysis. There were 160 such high-dimensional areas. The areas in Tables 5 and 9 are included in these 160 areas, but the area in Table 6 is not. The distributions of the number of dimensions with the maximum existence probability of high-growth firms in the 1036 high-dimensional areas and the 160 non-local highdimensional areas for verification are presented in Figure 8a,b. 62 Entropy 2023,25, 488 E D Figure 8. Distribution of the number of dimensions with the maximum existence probability of highgrowth firms for verification. The vertical axis and horizontal axes are the proportion of 1036 areas in ( a ) and 160 areas in ( b ) and the number of dimensions, respectively. ( a ) In the 1036 high-dimensional areas. (b) In the 160 high-dimensional areas, which include more than 145 high-growth firms. As shown in Figure 8, the distribution of the number of dimensions that maximizes the existence probability of high-growth firms in the data for validation has changed significantly by narrowing down from 1036 high-dimensional areas to 160 high-dimensional areas, which include more than 145 high-growth firms. In most of the 160 areas, the number of dimensions in which the existence probability of high-growth firms is maximized in data for verification is four-dimensional or higher. Therefore, in these 160 areas, the robustness can be assumed to be up to four-dimensional. Focusing on these 160 areas, 1–3 in Section 3.2.4 of the method are performed on these areas. The first corresponds to 40 areas , the second to zero areas, and the third to two areas. We finally focused on the 118 four-dimensional areas. We extracted high-dimensional areas from each of the 29 two-dimensional areas extracted by the five random data. Consequently, we extracted seven three-dimensional areas and 22 two-dimensional areas. The results using the random data are presented in Table 10. Table 10. Results using random data. EPHF represents the value of the existence probability of high-growth firms in the area where the existence probability of high-growth firms is the highest among the extracted high-dimensional areas. Data NAE-1D NAE-2D NDEHA EPHF 1 11 3 2, 2, 3 0.0140 2 11 7 2, 2, 2, 2, 2, 2, 3 0.0226 3 12 9 2, 2, 2, 2, 2, 2, 2, 3, 3 0.0161 4 13 4 2, 2, 2, 2 0.0167 5 13 6 2, 2, 2, 3, 3, 3 0.0152 Table 10 shows that we did not extract any high-dimensional areas in any random data. The area with the highest existence probability of high-growth firms among all the random data was the area where high-growth firms were 2.3 times more densely populated than usual. A comparison of the results with the data for analysis indicates that the high-growth firms are much more densely populated than in the random data. Considering that the random data extracted a maximum of only nine areas, the data for analysis, which extracted 1036 high-dimensional areas, showed that the high-growth firms were densely concentrated in many areas. Therefore, we can assume that strong relations exist between high-growing factors of firms and financial items. 63 Entropy 2023,25, 488 4.5. Grouping We define groups of the 118 four-dimensional areas selected in Step4 via hierarchical clustering with the ward method, Step5. The result is presented in Figure 9. ⋕ ⋔ ⋓ ⋒ ⋑ ⋐ ⋏ ⋎⋍⋌ ⋋ ⋊⋉ ⋈⋇ Figure 9. Dendrogram of the result of hierarchical clustering for the 118 four-dimensional areas. The vertical and horizontal axes are the dissimilarity defined using Equation (10) and the result of grouping the 118 four-dimensional areas, respectively. We divided the 118 four-dimensional areas into 15 groups (Groups 1 to 15 ). Four-dimensional areas belonging to the same group have a common color. For example, Group 1  has green. Groups 2  to 15  are cyclically painted in six colors. We set the dissimilarity threshold used for grouping in Figure 9 to a value that has a condition common to most of the grouped four-dimensional areas. Thus, the threshold was set to 1, except for the one group on the left, which is grouped because 34 of the 36 fourdimensional areas have the same condition. Finally, we divided the 118 four-dimensional areas into 15 groups. The conditions common to each of the 15 groups are presented in Table 11. We focused on Groups 1  , 2  , 12  , and 14  , which are characteristic among the 15 groups. Here, 34 of the 36 four-dimensional areas in Group 1  have the common condition of small gross profit per capita (less than 2727). The small value indicates that the firms in these 34 areas have small sales and poor operating efficiency. The remaining two fourdimensional areas have the condition that the total capital (compared to all firms in the same industry) is small (smaller than three) and the turnover of total capital (month) is large (larger than 17.79). The total capital (compared to all firms in the same industry) is the value evaluated by TDB and takes the value 0–10. The small value indicates that the total capital is very small compared to other firms in the same industry. The turnover of total capital (month) is the value of total capital divided by sales. Specifically, a large value of that indicates that sales are smaller than the total capital, given that the total capital is very small. Therefore, these two areas extract firms with very small sales and low efficiency. Therefore, the 36 four-dimensional areas in Group 1  extract firms with small sales and low operating efficiency. These firms are considered to have improved their operations and increased their sales significantly after three years. 64 Entropy 2023,25, 488 Table 11. Conditions common to each of the 15 groups. The abbreviated names used in this table are defined in Table 1. If the variables common to a group include those with an alphabet in front of the variable name, all four-dimensional areas in the group have in common that one or more of them are satisfied. For example, all four-dimensional areas in Group 5 contain the condition of the turnover of current assets and one or more of either (a) or (b). If the variables common to a group include variables with an alphabet with tilde in front of the variable name, all four-dimensional areas in the group have in common that two or more conditions are satisfied in them. For example, all four-dimensional areas in Group 15 contain two or more of the conditions ã, ˜ b, or ˜ c. Group Item name Threshold 1 GPE (T) ≤2727 2 ARR (DT) (M) ≥3.86 3 PPER (M) (a) Financial account to revenue ratio (%) (b) DR (%) ≤0.00 ≤0.00 ≤0.00 4 Cash and deposits to revenue ratio (days) OITC (IC) ≥130.33 ≤2.00 5 CAR (M) (a) Interest coverage ratio (times) (b) Capital to revenue ratio (M) ≥10.20 ≤−8.49 ≤−0.81 6  OITC (IC) (a) CACL (%) (b) NCLR (M) (c) Capital to revenue ratio (M) ≤2.00 ≤78.68 ≤0.00 ≤−0.81 7 CAR (M) Non-operating income to revenue ratio (%) ≥10.20 ≥4.62 8 CAR (M) DR (%) ≥10.20 ≤0.00 9  PPER (M) (a) TCR (M) (b) Revenue to total capital ratio (IC) (c) OIR (CFS) ≤0.16 ≥17.79 ≤3.00 ≤2.00 10 IR (M) ≤0.00 11  CAR (M) (˜ a) Non-operating income to revenue ratio (%) (˜ b) OITC (IC) (˜ c) APR (M) ≥10.2 ≤0.00 ≤2.00 ≤0.00 12 ARR (DT) (M) ≤0.25 13  CAR (M) (a) Financial account to revenue ratio (%) (b) Investment and financing returns (%) ≥10.20 ≤0.03 ≤0.02 14  CAR (M) IR (M) Non-operating income to revenue ratio (%) ≥10.20 ≤0.00 ≤0.05 15  (˜ a) Total capital (CFS) (˜ b) Investment and financing returns (%) (˜ c) Capital to revenue ratio (M) ≤3 ≤0.02 ≥8.53 Next, we focus on Group 2  and Group 12  . These two groups are characterized by different areas of the single variable of the trade receivables (discounted and transferred) turnover periods (months) as shown in Figure 10. Therefore, there is no firm that belongs to both Group 2 and Group 12 . 65 Citation: Pérez-Sienes, L.; Grande, M.; Losada, J.C.; Borondo, J. The Hurst Exponent as an Indicator to Anticipate Agricultural Commodity Prices. Entropy 2023,25, 579. https://doi.org/10.3390/ e25040579 Academic Editor: Panos Argyrakis Received: 1 March 2023 Revised: 16 March 2023 Accepted: 23 March 2023 Published: 28 March 2023 Copyright: © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). entropy Article The Hurst Exponent as an Indicator to Anticipate Agricultural Commodity Prices Leticia Pérez-Sienes 1, Mar Grande 1,2, Juan Carlos Losada 1and Javier Borondo 1,2,3,* 1Grupo de Sistemas Complejos, ETS Ingeniería Agronómica, Alimentaria y de Biosistemas, 28040 Madrid, Spain 2AgrowingData, Navarro Rodrigo 2 AT, 04001 Almería, Spain 3Departamento de Gestión Empresarial, Universidad Pontificia de Comillas ICADE, Alberto Aguilera 23, 28015 Madrid, Spain *Correspondence: jbor[email protected] Abstract: Anticipating and understanding fluctuations in the agri-food market is very important in order to implement policies that can assure fair prices and food availability. In this paper, we contribute to the understanding of this market by exploring its efficiency and whether the local Hurst exponent can help to anticipate its trend or not. We have analyzed the time series of the price for different agri-commodities and classified each day into persistent, anti-persistent, or white-noise. Next, we have studied the probability and speed to mean reversion for several rolling windows. We found that in general mean reversion is more probable and occurs faster during anti-persistent periods. In contrast, for most of the rolling windows we could not find a significant effect of persistence in mean reversion. Hence, we conclude that the Hurst exponent can help to anticipate the future trend and range of the expected prices in this market. Keywords: efficient market; time series; agri-food; Hurst; market; prices; agriculture 1. Introduction Financial markets are extremely complex systems with a large number of interacting units, and anticipating their evolution is far from straightforward. Thus, their study has attracted the attention of researchers over the past decades. An active and relevant topic of discussion among researchers is whether or not the financial market prices display long memory properties. The importance of this question lies in its consequences for market theories and its predictability. The fact that a market presents a long time memory implies that prices do not follow a random walk, as there is autocorrelation, and they are therefore predictable. On the other hand, if there is no memory, the Efficient Market Hypothesis (EMH) [ 1 , 2 ] cannot be rejected. The EMH was introduced by Fama in 1970 and states that new information is immediately reflected in the asset prices and therefore show martingale behavior. According to this theory, price changes are not related to the historical behavior of price volatility, but represent a response to new information, and since this arrives randomly, the evolution of prices is unpredictable. In the current literature there are papers supporting both hypotheses. Several authors have shown evidence of markets that present long time memory [ 3 – 9 ] whereas other authors have found evidence supporting the EMH [10–12]. In nature we find several examples of physical systems that do present long time memory properties, such as radiation or rainfall. Thus, researchers from econophysics, inspired by the idea that the financial system may share the same properties, have also attempted to detect trends and patterns in financial time series that can help to anticipate its trend. This discipline is known as Technical Analysis [13,14]. Due to the high importance of long time memory for the predictability of time series, there is a clear need for a method that identifies the existence or not of this memory. Entropy 2023,25, 579. https://doi.org/10.3390/e25040579 https://www.mdpi.com/journal/entropy 72 Entropy 2023,25, 579 Currently, the Hurst exponent (H) is the most widely used and accepted test to measure long-term memory properties [ 10 , 15 ]. H is a measure for long term memory and fractality of a time series that quantifies the degree of persistence of similar change patterns. This analysis was originally introduced by Hurst in 1951 to study the storage capacity of reservoirs in the Nile River, taking into account the cyclical trends of the flow, drought periods, and floods. This work was popularized and extended to other disciplines in the 1960s by Benoit Mandelbrot [ 16 – 18 ], who claimed that this methodology was superior to the autocorrelation, the variance analysis, and to the spectral analysis. Since then, several other methods of Hurst calculation have been developed. The best known ones are Rescaled Range [ 15 ], Detrended Fluctuation Analysis [ 19 , 20 ], wavelet transforms [ 21 ], and Generalized Hurst Exponent [22]. H ranges between 0 and 1, and provides information on whether the series presents long-term or not. If H= 0.5, then each step is independent of the past values of the series. Thus, there is no memory and the series is equivalent to white noise. Under this setting, the time series is unpredictable and the EMH is fulfilled. When H≤ 0.5, the series is anti-persistent. In this scenario, the series is expected to display ‘mean-reversion’. This fact implies that increments are generally followed by a decrease, while drops are followed by an increment. Finally, when H≥ 0.5, the series is persistent. In a persistent regimen, the series is more likely to maintain the trend in a broader range than what is expected by pure random walk. Thus, a rise in the previous step will most likely be followed by another rise, while a fall will be followed by another fall. Long-term memory is an important feature of market dynamics with implications for its predictability. As a consequence, the Hurst exponent has been widely applied to study the stock, currencies markets, and, more recently, to cryptocurrencies [ 23 , 24 ]. For instance, Di Matteo et al. [ 12 ] show how H serves as an index to classify mature and emergent markets. In the same line of research, Bianchi et al. [ 25 ], used the Hurst–Hölder exponents to detect periods of efficiency and inefficiency in stock markets [ 26 ]. Other researchers have analyzed how to use H to find the most profitable trading pairs [ 27 ], concluding that H performs better when compared with the classical methods. Despite the wide use of H in the stock market, there is still limited research on its applicability to the agri-commodities market [28–30] . In this paper we will study the agri-commodities market, and more particularly the evolution of prices for four horticultural products. Understanding this market is becoming increasingly important to make the agri-food industry sustainable [31,32], as price crises result in a waste of food. The stock market and the agri-commodities market have some similarities, but at the same time also have some important differences. On the one hand, both of them represent a market that is driven by demand. Matia et al. [ 33 , 34 ] showed that the two markets share several properties, although they also found some differences. The cumulative distribution of returns can be adjusted to a power law for both markets. In addition, the returns for the stocks and commodities market exhibit a multifractal behavior. On the other hand, there are differences in the nature of these two markets. In contrast with the stock market, in agricommodities markets, commodities represent a physical product that has to be stored and transported, and in some cases it is even a fresh and perishable product. Moreover, for agri-commodities we can expect slower changes and response to the demand, since the market is very conditioned by the supply of each product. In this paper, we will explore the applicability of H to the agri-food market in order to anticipate the trend and range of the future price. In particular, we focus on fresh vegetables because price crises have a big impact on them, since they are perishable products that can not be stored. Thus, when co-ops fail to anticipate the price and can not market their production, it results in tons of wasted food. The effect of the long memory properties on the agri-commodities markets has still attracted little attention from researchers. Thus, there is a gap in the current scientific literature, which misses to fully understand the behavior of such markets. In Ref. [ 35 ] the authors analyze the auto-correlations and crosscorrelations of the volatility time series for the Brazilian stock and commodity markets. 73 Entropy 2023,25, 579 They found auto-correlations in the commodity market, which in fact are stronger than that observed for the stock market. In another study—see Ref. [ 30 ]—the authors computed the Hurst exponent for several commodity price series, and found that most commodity prices are consistent with the underlying assumption of a geometric Brownian motion. We will contribute to understanding the dynamics and properties of the market by analyzing the evolution of H over time, and evaluating whether the value of H can provide useful information to anticipate the future trend of the price. In addition, we will compare the results obtained when computing H for different time windows. The present paper is organized as follows. In Section 2, we will explain the methodology followed to compute the local Hurst exponent of the series and the mean reversion. Next, in Section 2.3, we describe our data. In Section 3, we expose our results. Finally, in Section 4, we present our conclusions and discuss the importance of our results. 2. Materials and Methods 2.1. Hurst Exponent The Hurst exponent (H) is used in time series analysis and fractal analysis as a measure of the long-term memory of a time series. In other words, H measures how chaotic or unpredictable a time series is. In the literature, we can find several methods to calculate H , such as re-scaled Range (RS) [ 15 ], Detrended Fluctuation Analysis (DFA) [ 19 , 20 ], wavelet transforms [21], and Generalized Hurst Exponent (GHE) [22]. In this work, we use the GHE algorithm in order to measure the long-term memory of the price time series of different agri-commodities. This method is based on the scaling behavior of the statistic: Kq(τ)=|X(t+τ)−X(t)|q |(X(t))q| , (1) which is given by Kq(τ)∝τqH, (2) where τ is the time scale and can vary between 1 and τmax , H is the Hurst exponent, <·> denotes the sample average on time t , and q represents the order of the moment considered. H is then calculated by taking logarithms in relation (2) for different values of τ . In this paper, we work with τ= 2 n(n= 0,1, ... , log2(N)− 2 ) , and q= 1, as H1 is the closest estimation to the classical Hurst exponent [12]. H ranges between 0 and 1, where H= 0.5 means that there is no memory and the series is equivalent to white noise. When H≤ 0.5, the series is considered anti-persistent and is expected to display ‘mean-reversion’. Finally, when H≥ 0.5, the series is considered persistent and is more likely to maintain the trend in a broader range than what is expected by a pure random walk. In order to prevent H from using future values of the time series, we calculate a local Hurst exponent with reference to a rolling window of 4, 8, 16, 32, and 52 weeks that ends the day of measurement. This method ensures that we use only past data to determine H. Note that in order to compute H, we have coded the described method using Python. 2.2. Days to Mean Reversion Mean reversion (MR), or reversion to the mean, is a theory used in finance that suggests that a measure of interest such as the price of a commodity or asset eventually reverts to its long-term average levels. Thus, this theory assumes that a variable that deviates far from its long-term trend will return, reverting to its average value. This concept has been used to define many investment strategies that seek to purchase or sell financial products whose recent market price differs greatly from their historical average [36]. In this work we are going to test whether this reversion is more likely to occur during anti-persistent periods rather than during persistent periods in the price time series. Thus, for each day i , we compute the number of days (d) that the price (p) of an agri-commodity lasts to revert to its average value (m)as 74 Entropy 2023,25, 579 d=j∗−i where j∗is the minimum jthat satisfies pj≤mi+,ifpi≥mi pi≥mi−,ifpi<mi with j≥i. The average value (m) of each price time series have been calculated with reference to a rolling window of 4, 8, 16, 32, and 52 weeks, thus ensuring that we use only past values of the series to determine m. 2.3. Dataset In this work, we analyze the price time series of four agri-commodities (aubergine, zucchini, green pepper, and cucumber) in the South region of Spain. We have focused on these four products because they are representative of the European vegetables market, as they represent a significant percentage of the total imports and exports. The units of the prices are measured in EUR per kilogram. Each time series consists of daily prices collected over five years, and we must note that a typical week consists of six observations, since the market is closed on Sundays. Figure 1 shows the evolution of the price of each commodity during the period 2015–2019 . As it can be seen, all products present high volatility and despite showing a seasonal component, the noise component is still very important. In addition, we have included in the figure the moving averages for different rolling windows. In Table 1, the Hurst exponent, average price, and standard deviation of each time series are shown. All of the products present a global Hurst below 0.45, i.e., antipersistent. Green pepper is the most antipersistent time series ( H= 0.18), and also the commodity with the highest and least volatile price. Figure 1. The figure shows the evolution of the price and its corresponding moving average for different rolling windows RW for (a) aubergine, (b) zucchini, (c) pepper, and (d) cucumber. 75 Entropy 2023,25, 579 Table 1. Total Hurst, H, and mean and standard deviation of price for the four time series shown in Figure 1. Aubergine Zucchini Green Pepper Cucumber H0.36 0.40 0.18 0.32 ¯ p±σp0.70 ±0.53 0.75 ±0.63 0.85 ±0.22 0.64 ±0.33 pmax −pmin 3.99 3.89 1.45 2.08 3. Results The main goal of our paper is to analyze if H can help to anticipate the future trend of the price of four different agri-commodities, and discuss the difference in performance when considering different rolling windows (RW). To this end, we will analyze whether the probability of MR in the short term is more likely during anti-persistent days than on persistent days or not. 3.1. Evolution of the Hurst Exponent To achieve this goal we begin by computing H as described in the methods section, over our time series for different values of RW (4, 8, 16, 32, and 52 weeks). Thus, for each day of our series we have computed the value of H for the mentioned windows. The evolution of H, for aubergine over time and its distribution, can be found in Figure 2. Panel B of the figure shows the distribution of H, which approximately follows a normal distribution. The mean value of H depends on the commodity and RW, but for most cases (except for green pepper) is close to 0.5. The exact values for each combination can be found in Tables 2 and 3. We found that for the four products, H varies over time, alternating periods of persistence, anti-persistency, and neutral regimens. Changes over days tend to be relatively smooth, and when the series enters one of the three regimens it keeps for a while. This effect is illustrated in panel B of Figure 2, which shows the case of aubergine. 3.2. Mean Reversion In the second step, we compute for each day and RW the days to MR—this is the number of days left before the price will return or cross the mean. The days to MR follow a heterogeneous distribution, where for all the RWs over 50% of the observations revert to the mean in less than 20 days, while a small fraction take over 100 days. This can be observed in Figure 3a), which shows the probability mass functions (PMF) of each RW for aubergine. Figure 2. ( a ) The evolution of the local Hurst exponent for aubergine and different rolling windows. ( b ) Histogram of aubergine local Hurst exponent values. In both panels the dashed lines reflect the selected thresholds of 0.45 and 0.55. 76 Entropy 2023,25, 579 Table 2. This table shows the median of days to MR and mean local Hurst ( ¯ H ) for the studied rolling windows RW = 4,8,15,32,52 weeks for aubergine and zucchini. The global value of H of each product is also shown in the top row. Aubergine (H=0.36) Zucchini (H=0.40) RW Days to MR ¯ H±σHDays to MR ¯ H±σH 4 10 0.57 ±0.18 10 0.57 ±0.20 8 11 0.56 ±0.14 12 0.57 ±0.15 16 15 0.54 ±0.12 22 0.55 ±0.12 32 20 0.58 ±0.08 24 0.52 ±0.10 52 14 0.46 ±0.06 17 0.49 ±0.07 Table 3. This table shows the median of days to MR and mean local Hurst ( ¯ H ) for the studied rolling windows RW = 4,8,16,32,52 weeks for green pepper and cucumber. The global Hurst of each product is also shown in the top row. Green Pepper (H=0.18) Cucumber (H=0.32) RW Days to MR ¯ H±σHDays to MR ¯ H±σH 4 3 0.32 ±0.16 11 0.52 ±0.19 8 4 0.28 ±0.10 11 0.53 ±0.12 16 5 0.26 ±0.07 14 0.50 ±0.07 32 5 0.25 ±0.05 15 0.46 ±0.05 52 4 0.22 ±0.04 14 0.41 ±0.06 Figure 3. ( a ) Probability mass distribution of days to MR for the aubergine time series and for various RW. (b) The corresponding cumulative probability distributions. When exploring the relation between H and days to MR, we find that the product (pepper) with a significantly smaller value of H, and in an anti-persistent regime (0.22–0.32), generally reverts to the mean more rapidly. For example, when conducting the analysis with RW = 16 weeks, the price of pepper returns to the mean in an average of 5 days, while for the other products this value ranges between 11 and 22 days. Next, we analyze if the probabilities that the price will revert to the mean are significantly different for persistent and anti-persistent periods for each of the rolling windows. To this end, we classify each observation of the time series into persistent ( H≥ 0.55), neutral (0.45 <H< 0.55), and anti-persistent ( h≤ 0.45). Thus, we classify each day into one of the three different groups. We calculate the cumulative probability distribution of the days to MR for the persistent group, the anti-persistent group, and the total population and compare them for all the RWs. Figure 4 shows these distributions for aubergine. 77 Entropy 2023,25, 579 We find that for all the RWs, except for 52 weeks, the cumulative distribution of days to MR significantly differs for the two regimens, the probability of MR being higher for antipersistent days. In agreement, with this observation, the curve for the total population lies in-between both. Thus, anti-persistent days revert to the mean more quickly than persistent ones, which at the same time revert with lower probability than the global population. For the mentioned RW of 52 weeks (1 year), we observe the contrary effect, where the probability of MR in d or less days is always smaller for the anti-persistent regimen, the persistent group and the total population exhibiting a very similar behavior. Figure 4. Cumulative probability of days to MR for the price time series of aubergine for different rolling windows RW . The blue curve represents the antipersistent regime, while the green curve represents the persistent regime. The horizontal line shows the 50% of probability, and the vertical one marks 6 days. The figure also shows the 95% level of confidence for the 500 random sub-samples. To test whether the observed behavior for the persistent and antipersistent groups could be random, we take random samples of the total population and compare their behavior to both groups. We do so, because the persistent and antipersistent groups are subsamples of the total population. Thus, there is a chance that by choosing a random subsample of similar size we find a similar effect. In such cases the effect we have observed would not be significant. Thus, we have randomly selected 500 subsamples and computed the cumulative distribution of days to MR for each one of them. In Figure 4, we have included the 2.5% and 97.5% percentiles so that they can be compared with the curves of the two groups. We find that for all the RWs, except 32 weeks, the effect of the persistent group is not likely to be significant as the curves lie in-between the two percentiles. However, for RW = 32 weeks the probability of MR in 10 or less days for persistent days is below 30% and for 20 days it is below 40%. This means that we are around 70% confident that the price will stay at the current side of the border of the mean during the following 10 days, and 60% that it will also stay in the next 20 days. Hence, when computed with a RW of 32 weeks, H≥ 0.55 seems to be a good and informative indicator to anticipate the range in which the price will move in the medium-long term. On the other hand, the anti-persistent group differs from the total population and the random samples for RWs of 4, 8, and 32 weeks, but not for 16 weeks. When analyzing our data, we can see that the most informative RWs are the shorter ones: 4 and 8, as the gap between the antipersistent group and the total population is the largest. This effect is 78 Entropy 2023,25, 579 especially relevant in the case of 8 weeks, where the probability of MR in one week is 58%, while the value for the total population is under 40%. Thus, an indicator detecting points where H≤ 0.45 would be informative for movements in the short term, as we will have 58% of probabilities that the price will return to its average. Hence, if the actual price is below the average we can anticipate an uptrend, and if the price is below we can expect a drop in the prices for the following week. 3.3. Paired t-Test To further test the effect of persistence and anti-persistence in the expected days to MR, we adopt a quasi-experimental design, through which we can compare persistent and antipersistent observations against a control group. To this end, we perform two dependent t-tests for paired samples. The first one is to measure the effect of persistence and the second one is to measure the effect of anti-persistence. We use a dependent paired t-test, because our observations are extracted from a time series of prices, and thus can not be considered independent. In particular, we control for the commodity and the week, since the time series of prices have a marked seasonal component. Thus, for the persistence experiment we match the number of days to MR of each persistent observation to the number of days to MR on a randomly selected observation of the same week and commodity in the control group (all non antiperspirant observations). The anti-persistence experiment is designed analogously. The results of both tests for the RWs of 4, 8, 16, and 32 weeks are summarized in Figure 5. As it can be seen for all the RWs except for RW = 16, the MR occurs significantly faster for antipersistent observations. The most informative RW is the one of 8 weeks, where for antipersistent days we can expect MR to happen on average 12 days faster. In contrast, the effect of persistence is not so evident in our data, and we only find a significant effect for RW = 16, where in persistent periods MR takes on average 12 days longer. Note that we have not included the RW = 52 case, as there were not enough paired observations for the results to be reliable. Figure 5. Results of the paired t-test for the persistence ( H≥ 0.55) and antipersistence ( H≤ 0.45) experiments. The figure shows the t-statistic, measured in days, for both experiments and the different RWs (RW = 4, 8, 16, 32 weeks). The significant differences have been marked with a red cross (p-value < 0.01). 4. Discussion In this paper we contribute to understand the agri-food market by exploring whether the market presents long term memory properties or not for several agri-commodities. 79 Entropy 2023,25, 579 To this end, we have computed the local Hurst exponent for several values of RW and measured the relation between its evolution and the probability that the price will revert to the mean. We found that in general for antipersistent days MR is more probable in the short term and occurs faster than for persistent and neutral days. The only value of RW for which we could not find a significant effect was RW = 16, while the most informative one was RW =8. The fact that MR is more probable and happens faster during antipersistent periods means that H can be a good indicator to anticipate the future motion and help actors operating in this market, such as co-ops or supermarkets. More particularly, it is important to discuss its implications to anticipate the price and operate in this market in the short-term. For RW = 8 weeks we found that for antipersistent days almost 60% of the times MR will happen in one week or less, a magnitude significantly larger than what is expected for the full population, where the chances of MR in one week are below 45%. This fact, shows that H can be helpful to operate in such a market as it provides information to anticipate the future trend of the price. If the price is below the average, the operator will know that there is a high chance that the price will go up, while when the price is above the average a downtrend is very probable. We focus on the one week resolution, because anticipating MR in the long term is not so useful in this market. For example, knowing that MR will happen during the following two months, but not knowing when, means that the operator has to trade daily tons of fresh products with a high uncertainty on when the price movement will happen. On the other hand, market indicators related to persistence can be related to the fact that MR is not very probable or will happen slowly. Thus, this kind of indicator is more useful to operate in the market in the long term. The fact that there is a low probability of MR for the following days is not too informative, as price time series present autocorrelation and the price from one day to another usually does not present big differences. In contrast, knowing that MR is not probable in the long term will help to anticipate the range in which the price will oscillate in the following months. When the price is above the average, knowing MR is not probable, and it is useful to know the lower barrier that the price will not cross in the following weeks and months. Likewise, when the price is below the average we have an upper frontier that the price is unlikely to cross. Thus, this information helps the actors of the market to negotiate long term contracts, which is very common between supermarkets and co-ops, where the second commits to provide a minimum quantity of tons during the following months to the second one for a fixed price. A relevant research topic that we plan to explore in future work is the development of a methodology to find the optimal rolling window to use when computing the local Hurst exponent for each price time series. Thus, we aim to analyze a wider variety of products that follow different dynamics, and find their corresponding best rolling window. Author Contributions: Conceptualization, J.C.L. and J.B.; Methodology, J.C.L. and J.B.; Validation, M.G., J.C.L. and J.B.; Formal analysis, L.P.-S. and M.G.; Investigation, L.P.-S.; Data curation, M.G.; Writing—original draft, J.B.; Writing—review & editing, L.P.-S., M.G., J.C.L. and J.B.; Visualization, L.P.-S. and M.G.; Funding acquisition, J.C.L. and J.B. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by Spanish Ministry of Science and Innovation under Contract No. PID2021-122711NB-C21 and DIN2018-010114, and by DG of Research and Technological Innovation of the Community of Madrid (Spain) under Contract No. IND2022/TIC-23716. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: The datasets used for this study are publicly available from Observatorio de Precios y Mercados, Junta Andalucia, Spain at https://www.juntadeandalucia.es/ agriculturaypesca/observatorio/servlet/FrontController?action=Static&url=introduccion.jsp and agroprecios.com at https://www.agroprecios.com/es/precios-subasta/. 80 Entropy 2023,25, 579 Conflicts of Interest: The authors declare no conflict of interest. References 1. Fama, E.F. Efficient capital markets: A review of theory and empirical work. J. Financ. 1970,25, 383–417. [CrossRef] 2. Fama, E.F. Efficient capital markets: II. J. Financ. 1991,46, 1575–1617. [CrossRef] 3. Greene, M.T.; Fielitz, B.D. Long-term dependence in common stock returns. J. Financ. Econ. 1977,4, 339–349. [CrossRef] 4. Hampton, J. Rescaled range analysis: Approaches for the financial practitioners, Part 3. Neuro Vest J. 1996,4, 27–30. 5. Lillo, F.; Farmer, J.D. The Long Memory of the Efficient Market. Stud. Nonlinear Dyn. Econom. 2004,8, 1–19. [CrossRef] 6. Barkoulas, J.T.; Baum, C.F. Long-term dependence in stock returns. Econ. Lett. 1996,53, 253–259. [CrossRef] 7. Wright, J.H. Long memory in emerging market stock returns. FRB Int. Financ. 2000 , pii: Discussion Paper No. 650. Available online : https://www.federalreserve.gov/econres/ifdp/long-memory-in-emerging-market-stock-returns.htm (accessed on 27 March 2023). 8. Kasman, S.; Turgutlu, E.; Ayhan, A.D. Long memory in stock returns: Evidence from the major emerging central European stock markets. Appl. Econ. Lett. 2009,16, 1763–1768. [CrossRef] 9. Cheong, C. Estimating the hurst parameter in financial time series via heuristic approaches. J. Appl. Stat. 2010 ,37, 201–214. [CrossRef] 10. Lo, A.W. Long-Term Memory in Stock Market Prices. Econometrica 1991,59, 1279–1313. [CrossRef] 11. Lo, A.W.; MacKinlay, A.C. Long-Term Memory in Stock Market Prices. A Non-Random Walk Down Wall Street, 1st ed.; Princeton University Press: Princeton, NJ, USA, 1999. 12. Di Matteo, T.; Aste, T.; Dacorogna, M.M. Long term memories of developed and emerging markets: Using the scaling analysis to characterize their stage of development. J. Bank. Financ. 2005,29, 827–851. [CrossRef] 13. Brown, D.P.; Jennings, R.H. On technical analysis. Rev. Financ. Stud. 1989,2, 527–551. [CrossRef] 14. Park, C.H.; Irwin, S.H. The Profitability of Technical Analysis: A Review. 2004. AgMAS Project Research Report No. 2004-04. Available online: http://dx.doi.org/10.2139/ssrn.603481 (accessed on 27 March 2023). 15. Hurst, H.E. Long Term storage capacity of reservoirs. Trans. Am. Soc. Civ. Eng. 1951,116, 770–799. [CrossRef] 16. Mandelbrot, B.B. When can price be arbitraged efficiently? A limit to the validity of the random walk and martingale models. Rev. Econ. Stat. 1971,53, 225–236. [CrossRef] 17. Mandelbrot, B. Statistical methodology for nonperiodic cycles from covariance to R/S analysis. Ann. Econ. Soc. Meas. 1972 ,1, 259–290. 18. Mandelbrot, B.; Wallis, J.R. Robustness of the rescaled range R/S in the measurement of noncyclic long-run statistical dependence. Water Resour. 1969,5, 967–988. [CrossRef] 19. Peng, C.K.; Buldyrev, S.V.; Havlin, S.; Simons, M.; Stanley, H.E.; Goldberger, A.L. Mosaic organization of DNA nucleotides. Phys. Rev. E 1994,49, 1685. [CrossRef] [PubMed] 20. Hu, K.; Ivanov, P.C.; Chen, Z.; Carpena, P.; Stanley, H.E. Effect of trends on detrended fluctuation analysis. Phys. Rev. E 2001 , 64, 011114. [CrossRef] 21. Simonsen, I.; Hansen, A.; Nes, O.M. Determination of the Hurst exponent by use of wavelet transforms. Phys. Rev. E 1998 , 58, 2779. [CrossRef] 22. Barabasi, A.L.; Vicsekt, T. Multifractality of self affine fractals. Phys. Rev. A 1991,44, 2730–2733. [CrossRef] 23. Caporale, G.M.; Gil-Ana, L.; Plastun, A. Persistence in the cryptocurrency market. Res. Int. Bus. Financ. 2018 ,46, 141–148. [CrossRef] 24. Dimitrova, V.; Fernández-Martínez, M.; Sánchez-Granero, M.A.; Trinidad Segovia, J.E. Some comments on Bitcoin market (in)efficiency. PLoS ONE 2019,14, e0219243. [CrossRef] [PubMed] 25. Bianchi, S.; Pianese, A. Time-varying Hurst-Hölder exponents and the dynamics of (in)efficiency in stock markets. Chaos Solitons Fractals 2018,109, 64–75. [CrossRef] 26. Cajueiro, D.O.; Tabak, B.M. Ranking efficiency for emerging markets. Chaos Solitons Fractals 2004,22, 349. [CrossRef] 27. Ramos-Requena, J.P.; Trinidad-Segovia, J.E.; Sánchez-Granero, M.A. Introducing Hurst exponent in pair trading. Physica A 2017 , 488, 39–45. [CrossRef] 28. Corazza, M.; Malliaris, A.G.; Nardelli, C. Searching for fractal structure in agricultural future markets. J. Future Mark. 1997 ,17, 433–473. [CrossRef] 29. Barkoulus, J.; Labys, W.C.; Onochie, J. Fractional dynamics in international commodity prices. J. Future Mark. 1997 ,17, 161–189. [CrossRef] 30. Turvey, C.G. A note on scaled variance ratio estimation of the Hurst exponent with application to agricultural commodities prices. Physica A A Stat. Mech. Its Appl. 2007,377, 155–165. [CrossRef] 31. Allen, P. Together at the Table: Sustainability and Sustenance in the American Agrifood System; Penn State Press: University Park, PA, USA, 2004. 32. Borsellino, V.; Schimmenti, E.; El Bilali, H. Agri-food markets towards sustainable patterns. Sustainability 2020 ,12, 2193. [CrossRef] 33. Matia, K.; Amaral, L.A.N.; Goodwin, S.; Stanley, H.E. Different scaling behaviors of commodity spot and future prices. Phys. Rev. E2002,66, 045103. [CrossRef] 81 Entropy 2023,25, 619 where St denotes the state variable and is assumed as the following Markov transition probability [48]: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ P11 =P(St=1|St−1=1)=exp(π1) 1+exp(π1) P12 =P(St=2|St−1=1)=1 1+exp(π1) P21 =P(St=1|St−1=2)=1 1+exp(π2) P22 =P(St=2|St−1=2)=exp(π2) 1+exp(π2) (12) 3.3. Markov-Switching Mixed-Clayton Copula Function In the Ms-M-Clayton copula, the correlations of lower–lower tail, lower–higher tail, higher–higher tail, and higher–lower tail are provided as follows [49]: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ λLL MS−M=lim α→0P(V≤α|U≤α)=0.5ωst2−1 α1 λLU MS−M=lim α→0P(V≥1−α|U≤α)=0.5(1−ωst)2−1 α2 λUU MS−M=lim α→1P(V≥α|U≥α)=0.5ωst2−1 α3 λUL MS−M=lim α→1P(V≤1−α|U≥α)=0.5(1−ωst)2−1 α4 (13) 3.4. Parameter Estimation Method We employ the maximum-likelihood (ML) function [ 50 ] as the basis for estimating parameters. Given that there are 12 parameters to be estimated, and a traditional approach, such as the interior-point method, easily falls into local optimum, we apply the genetic algorithm (GA) that performs well in global optimization of high-dimensional parameters to exact the solution of the model [51]. Referring to Equation (6), the joint probability density function of the Ms-M-Clayton copula model with variables xand yis given as: fXY(x,y)= 2 ∑ St=1 fX(x)fY(y)cu,v,θStP(St)(14) where P(St) is the prediction probability of St at time t− 1. P(St=1) and P(St=2) are defined as [52]: P(St=1)=P11 ∗c1 t−1P(St−1=1) c1 t−1P(St−1=1)+c2 t−1P(St−1=2)+P21 ∗c2 t−1P(St−1=2) c1 t−1P(St−1=1)+c2 t−1P(St−1=2)(15) P(St=2)=1−P(St=1)(16) where c1 t−1 and c2 t−1 represent the conditional probability density functions of the copula function in state 1 and state 2, respectively, at time t− 1. Then the logarithmic likelihood function of the copula model is expressed as: lnL=∑T t=1lncu,v;θStP(St)+∑T t=1lnfX(x)+∑T t=1lnfY(y)(17) 3.5. VaR and CoVaR This work employs the value-at-risk (VaR) to measure the downside and upside risks, which indicates the maximum loss that an investor may suffer within a certain time horizon and significant level by holding a long or a short position. For return series rt , we calculate the VaR based on its marginal distribution. With a given tail probability α , the VaRα,t D 88 Entropy 2023,25, 619 and VaRα,t U at time tis calculated by Prt≤VaRα,t D=α and Prt≥VaRα,t U= 1 −α respectively, which is formulated as: VaRα,t D=μt+σt·F−1 v(α) VaRα,t U=μt+σt·F−1 v(1−α)(18) where μt and σt represent the conditional mean and standard deviation determined by the marginal distribution model, and F−1 v(α)is the α-quantile of GED. The conditional VaR (CoVaR) is used to capture the risk spillover between markets [ 42 ]. The CoVaR is calculated based on the measurement of copula model, reflecting the VaR of a market conditional on the extreme volatility in another market. Let ri t and rj t denote the return series of market i and j , and the CoVaR in four different market statuses can be expressed as follows: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Pri t≤CoVaRβ,t iD|jDrj t≤VaRα,t jD =β Pri t≥CoVaRβ,t iU|jDrj t≤VaRα,t jD =β Pri t≤CoVaRβ,t iD|jUrj t≥VaRα,t jU =β Pri t≥CoVaRβ,t iU|jUrj t≥VaRα,t jU =β (19) where CoVaRβ,t iD|jD and CoVaRβ,t iU|jD represent the downside and upside VaRs of market i conditional on the extreme downside movement of market j given a confidence level β , while CoVaRβ,t iD|jU and CoVaRβ,t iU|jU , respectively, represent downside and upside VaR of market i conditional on the extreme upside movement of market j given a confidence level β. For example, the first row in Equation (19) can be written as: Fri trj tCoVaRβ,t iD|jD,VaRα,t jD Frj tVaRα,t jD=β(20) Therefore, the CoVaR requires the joint distribution function of ri t and rj t , and it can be represented by a copula function as Equation (4). Thus, Equation (19) can be written as: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ CFri tCoVaRβ,t iD|jD,α=αβ CFri tCoVaRβ,t iU|jD,α=α−αβ Fri tCoVaRβ,t iD|jU−CFri tCoVaRβ,t iD|jU,1−α=αβ Fri tCoVaRβ,t iU|jU−CFri tCoVaRβ,t iU|jU,1−α=α−αβ (21) Hence, the value of Fri tCoVaRβ,t iD|jD can be inferred by inverting the copula function for given values of α and β , which is denoted as ˆ Fri tCoVaRβ,t iD|jD , and the value of CoVaR can be inferred by inverting the marginal distribution function of ri t as CoVaRβ,t iD|jD = F−1 ri tˆ Fri tCoVaRβ,t iD|jD . Similarly, the other three types of CoVaR can be obtained. To validate the significance of the risk contagion, the Kolmogorov–Smirnov (K-S) test [ 20 ]is employed to implement the significance test. 4. Data and Descriptive Statistics This work adopts the China Securities Index 300 (CSI300), an important financial index jointly released by the Shanghai and Shenzhen Stock Exchanges on 8 April 2005 to represent 89 Entropy 2023,25, 619 the Chinese stock market. It consists of 300 stocks, accounting for approximately 70% of the total market capitalization of the Shanghai and Shenzhen stock markets. Compared with other stock indexes in China, the issuers of the constituent stocks in CSI300 are mostly mature companies that have the characteristics of strong resistance to manipulation, lower volatility, and strong liquidity. Therefore, it comprehensively reflects the performance of the Chinese stock market. According to [ 19 ], three risk areas, including Asia–Oceania, Europe, and the Americas, can be identified in risk contagion. Therefore, the S&P500 and GSPTSE indexes are selected to represent the Americas market, the DAX30 and FTSE100 indexes are selected to represent the European market, and the Nikkei225 and ASX200 indexes are selected to represent the Asia–Oceania market. The monthly price time series collected from Wind database are used for empirical analyses because: (1) it covers less noises than the daily and weekly prices and is widely employed in copula modeling, and (2) it contains more trend information than the yearly prices but does not suffer from manipulation [ 14 , 20 , 48 ]. The period is from July 2005, when CSI300 is officially released, to December 2020, with 186 data points containing multiple economic cycles and economic events. The logarithmic returns series rt reflecting the level of price changes are calculated as: rt=(lnPt−lnPt−1)×100%, where Ptdenotes the price at the end of month t. Figure 1 reports the prices and returns of the selected stock indexes. First, the stock market volatility in the same region is relatively similar, but those in different regions are quite different. Second, due to the global emergencies during the sample period, such as the global financial crisis, the European debt crisis, and the COVID-19 epidemic, the markets experienced several large fluctuations simultaneously, implying the potential risk contagion between Chinese and mature markets. Third, the volatility of Chinese market is significantly higher than mature markets, which may be caused by the large gap between Chinese and mature stock markets in terms of the completeness of risk supervision and the professionalism of market participants.         &6,      3ULFH 5HWXUQV         63      3ULFH 5HWXUQV           *6376(      3ULFH 5HWXUQV        '$;      3ULFH 5HWXUQV         )76(      3ULFH 5HWXUQV          1LNNHL      3ULFH 5HWXUQV         $6;      3ULFH 5HWXUQV Figure 1. Monthly prices and returns of the selected indexes. Table 1 reports the descriptive statistics of the return series, in which their average values are all positive. The CSI300 has the highest monthly average return with 0.0096, followed by the S&P500 and the DAX30, while the FTSE100 has the lowest monthly average return. The CSI300 has the highest volatility, with the standard deviation of 0.0858, followed by the Nikkei225. The lowest standard deviation 0.0402 is observed in the FTSE100. Moreover, the skewness statistics are all less than 0, suggesting that all the return series are featured as a long tail to the left, and there are more extreme negative returns. The skewness values of the Nikkei225 and the ASX200 are larger than others, and that of the CSI300 is closer to 0. Meanwhile, the Nikkei225 and the ASX200 have the highest kurtosis, implying the leptokurtosis feature in Asia–Oceania market is more prominent. 90 Entropy 2023,25, 619 The Jarque-Bera (J-B) test confirms that all return series are not normally distributed but featured as leptokurtosis. The Pearson correlation coefficients between CSI300 and other indexes proves a weak but positive correlation between Chinese and mature markets, and the correlations between Chinese market and the Americas, Asia–Oceania, and the European markets decreases in turn. Table 1. Descriptive statistics of monthly returns. CSI300 S&P500 GSPTSE DAX30 FTSE100 Nikkei225 ASX200 Mean 0.0096 0.0062 0.0030 0.0059 0.0013 0.0046 0.0023 Max. 0.2463 0.1194 0.0997 0.1550 0.1155 0.1401 0.0949 Min. −0.2991 −0.1856 −0.2168 −0.2131 −0.1413 −0.2722 −0.2380 Std. 0.0858 0.0436 0.0404 0.0543 0.0402 0.0570 0.0428 Skew. −0.524 −0.888 −1.649 −0.809 −0.668 −0.896 −1.470 Kurt. 4.691 5.245 9.930 5.010 4.236 5.347 7.963 J-B. 30.685 a63.540 a456.478 a51.620 a25.654 a67.595 a257.900 a Pearson. 1.000 0.402 a0.403 a0.374 a0.307 a0.361 a0.380 a Note: superscript a represent the significant levels at 1%. 5. Empirical Results This study uses Eviews 9 to perform a marginal distribution estimation and output the residual series and MATLAB 2018 to fit copula models. 5.1. Marginal Distribution Estimation A diagnostic test on stationarity, autocorrelation, and heteroscedasticity needs to be conducted before marginal distribution modeling. The results are reported in Table A1 (seen in Appendix A), showing that all the return series are stationary by ADF, PP, and KPSS tests. According to the Ljung–Box test, only the CSI300 have autocorrelation. The Q2(P) and ARCH(P) statistics ensure the presence of ARCH effects in all series except the DAX30. Thus, AR-GARCH is suitable to fit the marginal distribution. Considering the significance of parameters and the results of diagnostic test, the results of marginal distribution are provided in Panel A of Table A2 (see Appendix A). Most coefficients are significant at 5% level. Panel B of Table A2 reports the diagnostic results for the residuals, in which the autocorrelation and conditional heteroscedasticity are effectively overcome. Then, the standard residues are employed to conduct the risk dependence analyses with copula models. 5.2. Dynamic and Asymmetric Dependence Measured by MS-M-Clayton Copula The M-Clayton copula model is first employed to measure both positive and negative dependence structures (Wang et al., 2013; Ji et al., 2018), and Table 2 reports the results, in which all parameters are significant at the 1% level. It is worth noting that the weight parameter ω across different pairwise returns is various, indicating that the existence of negative dependence between Chinese and mature stock markets. Therefore, how to recognize the occurrence of different risk dependence structures and correlations has become an urgent problem to be clarified. Table 3 further reports the estimated results of the MS-M-Clayton copula model, where the model outperforms the invariant M-Clayton copula in terms of the logarithmic likelihood values. Most of the estimated parameters are significant at the 10% level, meaning that there are not only both positive and negative dependence structures but also dependence-switching between Chinese and mature stock markets. Overall, the risk dependence structures and correlations are different in each dependence state. Taking the CSI300-S&P500 as an example, the P22 of 0.864 is significant and higher than P11 , meaning that state 2 is the dominant dependence structure. Similarly, for the CSI300-GSPTSE, CSI300-FTSE100, and CSI300-Nikkei225 pairs, state 2 plays a dominant role, while state 1 is dominant in CSI300-DAX30 and CSI300-ASX200 pairs. 91 Entropy 2023,25, 619 Table 2. M-Clayton copula estimates of CSI300 with mature stock indexes. CSI300-S&P500 CSI300-GSPTSE CSI300-DAX30 CSI300-FTSE100 CSI300- Nikkei225 CSI300-ASX200 α10.682 a2.194 a0.889 a0.892 a0.488 a0.969 a α22.68 ×10−7a 8.67 ×10−8a 1.69 ×10−7a 3.24 ×10−8a 49.415 a5.23 ×10−7a α30.674 a0.513 a2.587 a7.08 ×10−9a 0.506 a4.705 a α43.0967 a5.58 ×10−8a 1.18 ×10−7a 9.049 a5.946 a2.45 ×10−7a ω0.924 a0.594 a0.571 a0.923 a0.947 a0.443 a Log-L−13.473 −8.808 −10.369 −10.453 −10.075 −9.776 Note: superscript a represent the significant levels at 1%. Table 3. MS-M-Clayton copula estimates of CSI300 with mature stock indexes. Copula CSI300- S&P500 CSI300-GSPTSE CSI300-DAX30 CSI300-FTSE100 CSI300- Nikkei225 CSI300-ASX200 αS1 13.176 0.262 a0.652 a1.622 0.527 a0.243 a αS1 21.82 ×10−10 a 9.12 ×10−8a 25.801 a4.11 ×10−9a 54.908 a2.68 ×10−8a αS1 30.601 c5.594 a0.998 a3.278 a2.66 ×10−9a 0.141 αS1 43.043 8.28 ×10−8a 3.95 ×10−9a 1.70 ×10−9a 3.62 ×10−9a 1.770 αS2 10.584 a1.637 a0.275 0.653 a0.461 a1.627 αS2 21.95 ×10−10 a 3.23 ×10−8a 1.11 ×10−10 a 3.28 ×10−9a 3.24 ×10−10 a 0.094 αS2 30.616 a1.39 ×10−7a 20.078 a1.93 ×10−10 a 0.776 a5.363 a αS2 43.077 8.03 ×10−9a 5.89 ×10−10 a 6.832 a6.844 a1.37 ×10−8a ωS10.656 b0.984 a0.978 a0.827 a0.676 a1.000 a ωS20.999 a0.676 a0.415 c0.817 a0.971 a0.932 c P11 0.517 0.943 a0.846 a0.941 a0.943 a0.881 a P22 0.864 a0.974 a0.710 a0.971 a0.993 a0.727 a Log-L−13.626 −11.276 −12.487 −12.369 −11.522 −11.777 Note: superscript a, b, and c represent the significant levels at 1%, 5%, and 10%, respectively. Table 4 reports the tail correlation coefficients based on the constructed copula. Specifically, the values of λUU are larger than that of λLL between CSI300 and S&P500, DAX30, and Nikkei225, meaning that the upside risk correlation triggered by positive factors is stronger than the downside risk correlation triggered by negative factors, while the opposite relationship occurs between CSI300 and GSPTSE, FTSE100, and ASX200. Moreover, compared with the Americas and European mature markets, the downside risk correlation between Chinese and Asia–Oceania markets manifesting in synchronized decline is the lowest, which is usually paid special attention in practice. Although the negative dependence is not in dominant in the dominant state, it is still asymmetric. Specifically, the upper–lower tail correlation between CSI300 and S&P500, FTSE100, and Nikkei225 is stronger than the lower–upper tail correlation, indicating the probability of extreme rises in Chinese market when extreme declines occur in the three mature markets. The opposite situation can be found between CSI300 and DAX30. As for the main dependence state between CSI300 and GSPTSE and ASX200 returns, the negative dependence correlation is not observed. Therefore, during the period of smooth economic operation denoted by the main state, except for monitoring the positive risk spillover, Chinese investors and managers should pay close attention to investment opportunities in the declines of S&P500, FTSE100, and Nikkei225 while managing exposure carefully in the rises of DAX30. 92 Entropy 2023,25, 619 Table 4. Tail correlation coefficients between CSI300 and mature stock indexes. State 1 State 2 λLL λLU λUU λUL λLL λLU λUU λUL CSI300-S&P500 0.264 0.000 0.103 0.137 0.152 0.000 0.162 0.001 CSI300-GSPTSE 0.035 0.000 0.435 0.000 0.221 0.000 0.000 0.000 CSI300-DAX30 0.169 0.011 0.244 0.000 0.017 0.000 0.200 0.000 CSI300-FTSE100 0.270 0.000 0.334 0.000 0.141 0.000 0.000 0.083 CSI300-Nikkei225 0.091 0.160 0.000 0.000 0.108 0.000 0.199 0.013 CSI300-ASX200 0.029 0.000 0.004 0.000 0.304 0.000 0.410 0.000 Figure 2 provides the trajectories of PS1 and PS2 , in which the state transitions are observed in the risk dependence between Chinese and most mature markets. For CSI300-S&P500, there is no state-switching, and state 2 is dominant during the entire sample period, implying the stable dependence and risk correlation between the two markets. For CSI300-GSPTSE, the state transitions occur concentrated in the periods from 2013 to 2015, corresponding to cyclical financial market bubbles and the post-COVID-19 [ 3 ], in which the secondary state should be paid more attention because more investment opportunities appear with a stronger upside tail correlation and a downside tail correlation close to 0. The state transitions of CSI300-DAX30 appear periodically around 2009 (may be affected by European debt crisis) with weak persistence [ 53 ]. In the secondary state, the upside tail correlation is significant, while the downside correlation decreases to near 0, increasing the investment motivation. For CSI300-FTSE100, state 2 with apparent downside risk correlation is dominant in most of the period. However, state 1 with both upside and downside risk correlations switches to be the main dependence structure temporarily around 2009 (European debt crisis) and since the COVID-19 epidemic [ 53 ]. For CSI300-Nikkei225, state 1 with reversal correlation was the main state before 2009 and in 2012, corresponding to the global financial crisis and the Asian financial turmoil led by the exchange rate system, respectively [ 54 ]. However, state 2 with positive dependence structure plays a dominant role in most of the period, especially in recent years. For CSI300-ASX200, state 1 with a relatively low tail correlation is dominant. The state-switching process occurs around 2012 and 2015 temporarily, which is accompanied by an increase in positive risk correlation caused by regional financial turmoil [ 54 ]. Moreover, in the comparison between markets in different regions, the Asia–Oceania markets have the relatively low risk association, especially the downside risk correlation that is paid much attention in practice, with the Chinese market.       3UREDELOLW\ &6,63       3UREDELOLW\ &6,*6376(        3UREDELOLW\ &6,'$;       3UREDELOLW\ &6,)76(       3UREDELOLW\ &6,1LNNHL       3UREDELOLW\ &6,$6; Figure 2. State transition probabilities between Chinese and mature markets (The blue line represents state 1, and the orange line represents state 2). 93 Entropy 2023,25, 619 5.3. Comparative Analysis 5.3.1. Static Dependence Measured by Invariant Copula Models To explain the similarity and differences between our findings and previous research, we first employ seven commonly used invariant copulas, including the Gaussian, Student’s t, Gumbel, 180 ◦ rotated Gumbel, Clayton, 180 ◦ rotated Clayton, and SJC copulas [ 53 ]to measure the risk dependence between Chinese and mature markets. The estimated results of invariant copulas are reported in Table 5. Table 5. Invariant copula estimates of the CSI300 with mature stock indexes. Copula CSI300- SP500 CSI300- GSPTSE CSI300- DAX30 CSI300- FTSE100 CSI300- Nikkei225 CSI300- ASX200 Gaussian ρ0.372 a0.348 a0.361 a0.266 a0.314 a0.349 a Log-L−12.496 −10.781 −11.660 −6.077 −8.637 −10.853 Student’s t ρ0.372 a0.347 a0.375 a0.285 a0.314 a0.356 a v99.899 a99.983 a7.966 a7.527 a99.320 a8.910 a Log-L−12.474 −10.596 −12.438 −7.069 −8.634 −11.710 Gumbel δ1.223 a1.183 a1.260 a1.141 a1.183 a1.228 a Log-L−6.541 −4.251 −8.307 −2.477 −4.729 −6.603 180◦rotated Gumbel δ1.326 a1.281 a1.329 a1.259 a1.256 a1.316 a Log-L−15.498 −11.794 −14.612 −10.533 −10.484 −14.249 Clayton ρ0.630 a0.543 a0.594 a0.523 b0.505 a0.593 a Log-L−16.645 −13.235 −14.419 −12.023 −12.126 −14.634 180◦rotated Clayton ρ0.307 c0.263 0.355 c0.130 0.237 0.310 c Log-L−4.233 −3.046 −5.244 −0.671 −2.478 −4.320 SJC λU2.83 ×10−74.77 ×10−75.57 ×10−81.85 ×10−74.21 ×10−74.96 ×10−7 λL0.380 0.404 0.366 a0.346 0.320 0.354 Log-L−16.508 −11.868 −14.474 −12.044 −11.464 −14.735 Note: superscript a, b, and c represent the significant levels at 1%, 5%, and 10%, respectively. According to the logarithmic likelihood values, it is found that the Clayton copula performs the best with significant estimated parameters, followed by the 180 ◦ rotated Gumbel copula, the Student’s t copula, and the Gaussian copula, successively, and the Gumbel copula and 180 ◦ rotated Clayton copula perform the worst. In the SJC copula measuring asymmetric positive dependence, the lower tail correlations are larger than the upper ones, but most parameters are not significant. The results suggest a positive but asymmetric risk dependence between Chinese and mature markets, and the downside correlation is stronger than the upside correlation. Overall, the results are generally consistent with the findings drawn from M-Clayton and MS-M-Clayton copulas but fail to capture the negative dependence structure and the upside correlations between CSI300 and S&P500, DAX30, and Nikkei225 effectively. Moreover, the static copulas are unable to capture the time-varying or dependence-switching characteristics of the correlations. 5.3.2. Dynamic Dependence Measured by Time-Varying Parameter Copula To assess the dynamic risk dependence correlation between Chinese and mature markets, Table 6 further reports the estimated results of four TVP copulas, in which most of the estimated parameters are significant at the 10% level. It can be found that TVP copulas perform better than the corresponding invariant copulas. Specifically, the TVP-180 ◦ rotated Gumbel copula describing the lower–lower tail correlation effectively captures the risk dependence between Chinese and mature markets, and the TVP-SJC copula also proves 94 Entropy 2023,25, 619 that the lower–lower correlation is more significant. The results confirm the positive risk dependence structure and the prominent downside risk correlation between Chinese and mature markets. The effectiveness of time-varying mechanism in depicting the dynamic risk correlation is also verified. Although TVP copulas provide an analytical view on dynamic risk correlation, a significant difference between them and the proposed MS-M-Clayton copula is that the potential negative dependence structure is not effectively depicted. Table 6. TVP copula estimates of the CSI300 with mature stock indexes. Copula CSI300-S&P500 CSI300-GSPTSE CSI300-DAX30 CSI300-FTSE100 CSI300- Nikkei225 CSI300-ASX200 TVP-Gaussian ψ00.255 a1.162 a0.091 a0.711 a0.321 a1.601 a ψ10.270 a0.258 a−0.169 a0.023 a−0.049 a−0.739 a ψ21.227 a−1.423 a2.042 a−0.635 a1.117 a−1.762 a Log-L−13.283 −10.836 −13.771 −6.078 −8.661 −11.657 TVP-180◦Rotated Gumbel ωL2.435 a1.144 a2.800 a1.557 a0.993 a−0.429 a αL−0.815 a−0.311 a−0.768 a−0.557 a−0.324 a0.683 a βL−2.851 a−0.766 a−5.030 a−1.291 a−0.263 a0.340 a Log-L−18.779 −11.955 −17.436 −10.738 −10.512 −14.862 TVP-Gumbel ωU2.338 a2.903 a3.366 a3.106 a−0.654 a−0.608 a αU−0.916 a−0.744 a−1.018 a−1.363 a0.942 a0.620 a βU−2.541 a−6.164 a−6.821 a−4.265 a−0.110 a1.100 a Log-L−9.867 −9.881 −15.758 −7.732 −5.015 −8.784 TVP-SJC ωU−14.830 a−14.363 a−15.343 a−15.242 a−14.593 a−14.490 a αU−0.012 a−0.002 b− 8.391 × 10 −4c −0.002 a−0.002 a−5.799 ×10−4b βU−0.003 a7.327 ×10−54.025 ×10−6−1.465 ×10−6−1.469 ×10−5−1.642 ×10−4 ωL2.792 a0.459 a5.150 a4.447 a−0.181 a−2.235 a αL−5.960 a−2.988 a−18.110 a−15.809 a−1.477 a1.071 a βL−4.505 a−1.209 a−4.230 a−4.203 a−0.858 a3.771 a Log-L−19.595 −12.379 −17.450 −12.804 −11.371 −14.976 Note: superscript a, b, and c represent the significant levels at 1%, 5%, and 10%, respectively. 5.4. Asymmetric Risk Spillover Measurement by VaR, CoVaR and Nomalized CoVaR To provide implications for risk supervision and portfolio risk management, we studied the extreme risk spillovers between Chinese and mature stock markets in different routes by VaR and CoVaR based on the information from marginal distribution and Ms- M-Clayton copula model. We set α and β equal to 0.05 for downside CoVaR and 0.95 for the upside CoVaR calculation. Table 7 reports the summary statistics of the VaR and the CoVaR, and Figure 3 shows the dynamic trajectories for intuitive observation. For stock index pairs except CSI300-FTSE100, the absolute values of upside VaR and CoVaR are larger than those of the downside, respectively, meaning that the upside risk is larger than the downside risk in Chinese market. Moreover, the VaR and CoVaR show phased extreme fluctuations, which may be related to the macroeconomic uncertainties, such as the periods around 2008, 2013, and 2015. For the positive risk contagion (3 and 6 rows in Table 7), the absolute values of CoVaR are all greater than that of VaR when measuring either upside or downside risks, indicating the synergistic risk spillover from mature markets to the Chinese market. In the measurement of negative risk contagion (4–5 rows in Table 7), the absolute values of CoVaR are generally smaller than that of VaR , implying the weak existence of reverse risk spillovers. Overall, the positive risk contagion from mature markets to the Chinese market are more significant than the negative contagion. It is noteworthy that the downside risk contagion between Chinese and Asia–Oceania markets is relatively weak, suggesting that the Asia–Oceania market can be considered as a potential choice for investors in the Chinese market to diversify their investment portfolios. 95 Entropy 2023,25, 619 Table 8 further reports the hypothesis testing results by K-S test, and the statistics are generally significant at 10% level, rejecting the null hypothesis that VaR is equal to CoVaR. Table 7. Summary statistics of the VaR and the CoVaR (The CoVaRβ,t CSI300(D)|Other(D) and CoVaRβ,t CSI300(U)|Other(D) denote the downside and upside VaRs of the CSI300 conditional on the extreme declines of mature markets, respectively; the CoVaRβ,t CSI300(D)|Other(U) and CoVaRβ,t CSI300(U)|Other(U) denote the downside and upside VaRs of the CSI300 conditional on the extreme rises of mature markets, respectively). CSI300- S&P500 CSI300- GSPTSE CSI300- DAX30 CSI300- FTSE100 CSI300- Nikkei225 CSI300- ASX200 VaRα,t CSI300,D−12.258 (4.483) VaRα,t CSI300,U14.720 (4.227) CoVaRβ,t CSI300(D)|Other(D)−19.060 (6.337) −18.840 (6.335) −18.791 (6.271) −19.021 (6.228) −18.275 (6.143) −18.369 (6.126) CoVaRβ,t CSI300(D)|Other(U)− 7.773 (3.309) −12.474 (5.122) − 8.736 (3.646) −11.459 (4.664) −10.170 (4.875) − 9.026 (3.706) CoVaRβ,t CSI300(U)|Other(D)13.260 (3.827) 14.060 (5.324) 10.639 (3.344) 16.491 (5.574) 12.942 (5.057) 11.247 (3.468) CoVaRβ,t CSI300(U)|Other(U)21.470 (6.056) 19.090 (4.497) 21.704 (6.116) 18.901 (4.802) 21.119 (5.988) 20.147 (5.650) Note: this table reports the means and the standard errors (in parentheses) of VaR and CoVaR.    ± ± ±    &6,63    ± ± ± ±      &6,*6376(    ± ± ±    &6,'$;    ± ± ± ±      &6,)76(    ± ± ± ±      &6,1LNNHL    ± ± ± ±      &6,$6; 9D5&6,' W &R9D5&6,'_2WKHU' W &R9D5&6,'_2WKHU8 W 9D5&6,8 W &R9D5&6,8_2WKHU' W &R9D5&6,8_2WKHU8 W Figure 3. The dynamic trajectories of the VaR and the CoVaR. 96 Entropy 2023,25, 619 Table 8. The hypothesis testing for equalities of CoVaR and VaR. Null Hypotheses CSI300- S&P500 CSI300- GSPTSE CSI300- DAX30 CSI300- FTSE100 CSI300- Nikkei225 CSI300- ASX200 CoVaRβ,t CSI300(D)|Other(D)= VaRα,t CSI300,D 0.593 a(0.000) 0.577 a(0.000) 0.577 a(0.000) 0.582 a(0.000) 0.550 a(0.000) 0.550 a(0.000) CoVaRβ,t CSI300(D)|Other(U)= VaRα,t CSI300,D 0.582 a(0.000) 0.077 (0.637) 0.456 a(0.000) 0.159 b(0.017) 0.330 a(0.000) 0.445 a(0.000) CoVaRβ,t CSI300(U)|Other(D)= VaRα,t CSI300,U 0.181 a(0.004) 0.220 a(0.000) 0.478 a(0.000) 0.132 (0.077) 0.324 b(0.047) 0.412 a(0.000) CoVaRβ,t CSI300(U)|Other(U)= VaRα,t CSI300,U 0.533 a(0.000) 0.456 a(0.000) 0.544 a(0.000) 0.418 a(0.000) 0.517 a(0.000) 0.473 a(0.000) Note: superscript a and b represent the significant levels at 1% and 5% respectively. To further evaluate the intensity of risk spillovers in different routes and analyze its asymmetry, Table 9 reports the summary statistics of the CoVaR normalized by VaR (CoVaR/VaR). It can be observed that the mean values of CoVaRβ,t CSI300(D)|Other(D) VaRα,t CSI300,D are greater than those of CoVaRβ,t CSI300(D)|Other(U) VaRα,t CSI300,D , and the mean values of CoVaRβ,t CSI300(U)|Other(U) VaRα,t CSI300,U are greater than those of CoVaRβ,t CSI300(U)|Other(D) VaRα,t CSI300,U , indicating that the positive and negative risk contagion effects are asymmetric, and the positive effect is stronger than the negative effect. Meanwhile, the mean values of CoVaRβ,t CSI300(D)|Other(D) VaRα,t CSI300,D are greater than those of CoVaRβ,t CSI300(U)|Other(U) VaRα,t CSI300,U , and the mean values of CoVaRβ,t CSI300(U)|Other(D) VaRα,t CSI300,U are greater than those of CoVaRβ,t CSI300(D)|Other(U) VaRα,t CSI300,D except in CSI300-GSPTSE pairwise returns, implying the asymmetry between upside and downside risk contagion effects, and the downside effect is generally stronger, while the opposite effect is in negative contagion. The analyses are statistically supported by K-S tests (see in Tables A3 and A4 of Appendix A). Table 9. Summary statistics of the CoVaR/VaR. CSI300- SP500 CSI300- GSPTSE CSI300- DAX30 CSI300- FTSE100 CSI300- Nikkei225 CSI300- ASX200 CoVaRβ,t CSI300(D)|Other(D) VaRα,t CSI300,D 1.571 (0.082) 1.551 (0.082) 1.548 (0.079) 1.569 (0.080) 1.504 (0.069) 1.513 (0.080) CoVaRβ,t CSI300(D)|Other(U) VaRα,t CSI300,D 0.625 (0.056) 1.014 (0.161) 0.707 (0.101) 0.926 (0.094) 0.828 (0.264) 0.729 (0.054) CoVaRβ,t CSI300(U)|Other(D) VaRα,t CSI300,U 0.903 (0.056) 0.943 (0.162) 0.721 (0.063) 1.110 (0.113) 0.871 (0.185) 0.762 (0.051) CoVaRβ,t CSI300(U)|Other(U) VaRα,t CSI300,U 1.461 (0.055) 1.317 (0.110) 1.478 (0.059) 1.299 (0.118) 1.441 (0.111) 1.372 (0.061) Note: this table presents the means and the standard errors (in parentheses) of the CoVaR/VaR. 6. Conclusions The risk contagion between Chinese and mature markets has attracted more and more attention from both scholars and market participants. In this work, we construct a novel Ms- M-Clayton copula model to identify both positive and negative dependences and revisit the risk contagion between Chinese market and six mature markets in the Americas, Europe, and Asia–Oceania. Four basic Clayton copulas with various rotations are weighted to capture different tail correlations, and a two-state transition mechanism following Markov 97 Entropy 2023,25, 641 research conducted on this topic, and sourcing and analysing additional data using statistical techniques, the results from this article will inform the adaptation and use of a weighted network model, first proposed by [ 4 ], that captures the properties, interactions and dynamics of the international aid system. Much of the literature and research into foreign aid dynamics focuses on the donor. Econometric techniques, primarily regression and ordinary least squares ([ 5 – 7 ]), are commonly employed in an attempt to reveal the relative importance of donor motivations and potential biases behind their aid allocation decisions. More recently, there has been research conducted on the growing field of network theory and the utilisation of related mathematical methods to model the allocation of aid that goes beyond regression ([4,8]). However, it is rarer to find research and analysis focusing on aid recipients. By understanding donor motivations and biases, aid recipients could exploit these ‘assets’ and potentially increase their aid receipts if they are viewed as a portfolio of investments. By identifying and quantifying significant donor motivations for allocating their aid budgets, inputting these variables into a general weighted network model [ 4 ] and then adapting it using financial mathematics (modern portfolio theory), aid recipients could use the model to optimise their aid income portfolio, treating donor variables similarly to assets in an investment portfolio. This is illustrated in this article using a simulation. The principal aim, then, of this article is to illustrate the power of network science and mathematical modelling when applied to the complex and dynamical system of international aid. The potential impact is an increase in transparency of the often-opaque motivations and biases of aid donors, which subsequently could be employed by recipients to increase their aid income. 2. Methods 2.1. Data and Data Analysis The first step to evaluating, and then adapting, the general weighted network model [ 4 ] is to identify the significant motivations and preferences shown by selected donors regarding the allocation of their aid budgets. These will be used as the model’s variables. Subsequently, the accuracy of the model’s mechanics and outputs can be tested against actual historical data for selected donors and recipients. Data from the OECD and World Bank were sourced and analysed using various statistical techniques to identify and understand the inter-connecting and inter-dependent variables that drive the data. The pertinent results are summarised here. 2.1.1. Economic and Foreign Aid Data To compare the economic fortunes of one country versus another, gross national income (GNI), a key measure of economic well-being and a superior metric for assessing the overall economic condition of a country, especially for countries that have large foreign receivables or outlays, will be used for identifying the level of need of an aid recipient (‘recipient need’). Furthermore, to assist with country comparisons, GNI per capita will be used rather than absolute GNI. Table 1 lists the top 10 aid recipients in 2019 by net official development assistance (ODA) receipts, classified as total net ODA flows from Development Assistance Committee (DAC) countries, multilateral organisations and non-DAC countries. When identifying the top donor countries, rather than looking at absolute aid donated, the affordability of a donor country to provide aid is assessed using aid donated as a percentage of country GNI. This is summarised in Figure 1, which lists the members of the DAC, a development committee of the OECD. 104 Entropy 2023,25, 641 Table 1. Top 10 ODA recipients, including significant regional aid donations, and figures for all developing counties for comparison [9]. All figures in US$m unless otherwise stated. Net ODA Receipts GNI/CAP (US$) GNI ODA/GNI (%) Country/Region 2015 2016 2017 2018 2019 2019 2019 2019 Syrian Arab Republic 4920 8900 10,428 9997 10,252 - - - Ethiopia 3239 4084 4125 4941 4810 850 95,641 5.03 Bangladesh 2593 2533 3782 3045 4518 1940 316,907 1.43 Yemen 1778 2301 3234 7985 4397 - - - Afghanistan 4274 4069 3812 3792 4285 540 19,402 22.08 Nigeria 2432 2498 3359 3305 3531 2030 433,449 0.81 Kenya 2464 2188 2480 2491 3251 1750 93,578 3.47 Democratic Republic of the Congo 2599 2102 2293 2514 3026 520 45,879 6.59 Jordan 2141 2728 2980 2526 2797 4300 43,429 6.44 India 3174 2679 3198 2462 2611 2130 2,843,902 0.09 Regional (not specific to any country) South of Sahara 2435 2635 2759 3137 3410 Africa region 2184 2777 3017 3241 3201 All developing countries 146,742 158,811 165,090 166,540 168,588 511,750 292,854,611 0.58 Figure 1. ODA grant equivalent as percentage of GNI in 2020 for DAC donors in the OECD. The grey bars identify those countries that contribute less than the United Nations (UN) target of 0.7%, the blue bar shows total of the DAC countries as a percentage of their GNI and the green bars highlight those countries who contribute over the UN target of 0.7%. Figure 1 is sourced from the OECD website [ 10 ] and arranged in descending order based on the percentage of DAC-country GNI donated in 2020, with Sweden donating the highest percentage of its GNI at 1.15%. Donor affordability is epitomised by the 0.7% target agreed by the United Nations (UN) in 1970 for aid contributions by DAC countries to developing countries. It is reasonable then to assume that since the UN members agreed to 0.7%, they can therefore afford to 105 Entropy 2023,25, 641 donate 0.7% of their GNI. However, as shown in Figure 1, this target is not being met by most UN countries, including the USA. Bilateral aid flows—aid given directly from a country donor to a recipient donor— comprised circa. 67% of total ODA donated in 2019, with the remaining third being flows from multilateral institutions and international financial institutions (for example, the World Bank). However, the proportion of bilateral aid reduced significantly in 2020 by 36% on 2019 levels to 42% of the total aid donated, with multilateral institutions taking up the slack, due mainly to the impact of COVID [11]. 2.1.2. Aid Flows from Donors to Former Colonies There are robust conclusions in the research performed by [ 3 , 5 , 6 ], among many others, that a strong motivation behind aid allocation decisions by donors lies in whether the recipient is an ex-colony or not. Former colonies receive proportionately more aid from their former colonial masters than other recipients. Indeed, according to [ 5 ], between 1970 and 1994, France gave 57% of its total bilateral aid to its former colonies, the UK gave 78% and Portugal 99.6%. Moreover, according to the OECD [ 12 ], in 2009, the largest recipient of UK aid was India and, by 2019, this was Pakistan, both former UK colonies. Thus, colonial history is positively correlated with aid, as identified by [5] and confirmed by own analysis performed. 2.1.3. Trade Activity Before correlation techniques were applied to detect any interdependencies between trade activity and aid donations, the raw data were analysed. Trade data are sourced from the World Integrated Trade Solution (WITS) website, a sister site of the World Bank specifically focused on trade [ 13 ], for the period 1993 to 2019. By charting this trade data with aid data sourced from the World Bank [ 9 ], a pattern of aid versus trade can be viewed over time. This suggested a positive correlation, confirmed by calculating correlations between the two data sets over many periods. This result is also backed by research performed by [5,6,14]. 2.1.4. Recipient Need The literature is mixed regarding the relative importance of recipient need as a variable in a donor’s aid allocation decisions. In [ 5 ], the authors are clear on donor motivations being based mainly on self-interest and political and strategic considerations over aid recipient needs. However, later studies, such as [ 14 ], dispute this conclusion stating that self-interest, while still a significant input into aid allocation decisions, is not as important as recipient need. Moreover, [ 6 ] conclude that the USA behaves very differently from all other aid donors, except Japan, by putting much less emphasis on recipient need and much more emphasis on donor self-interest. 2.1.5. The Herding Phenomenon (the Bandwagon Effect) Another variable to consider for inclusion in a mathematical model of foreign aid is herding behaviour often exhibited by donors, also termed the ‘bandwagon effect’. This refers to the actions and impulses of a group of agents, countries, politicians, or financial traders to follow the actions of the ‘crowd’ rather than trust their own individual judgment. The phenomenon has similar attributes to ‘groupthink’. It is an emergent behaviour of a dynamical system due to the many interactions taking place within that system. Herding is commonly associated with financial market behaviour, for example asset bubbles [ 15 ]. Grounded in behavioural finance, herd mentality refers to investors’ bias to follow what other investors are doing, being largely influenced by emotion and intuition, rather than by their own evaluations of potential investments. In terms of aid allocation, the bandwagon effect manifests itself when a recipient receives more aid from one donor, leading to an increase in aid from many more donors. In 106 Entropy 2023,25, 641 other words, the more aid a recipient receives, the more it attracts. It likely depends on the relative influence of the lead donor rather than characteristics of the recipient. Research conducted by [ 6 ] attempted to measure the effect using regression and incorporating aid from other sources, not only ODA. They find that there is some support for the herding argument, but it is far from conclusive. Ref. [ 16 ] gives the phenomenon a more thorough review, concluding that there is around an 11% impact on aid donations through donor herding, which is relatively significant. 2.2. A Network Model for Foreign Aid The principal outcome of the research conducted, using data analysis and statistical techniques, is the identification of the following significant aid determinants: •Past colonial relationships; •Trade activity and commercial interests; •Poverty alleviation (recipient need); and •Bandwagon impacts (‘herding’). These variables will now be incorporated into a weighted network model to demonstrate how such a model can be modified and used by an aid recipient to treat their various donor-sourced aid receipts as an investment portfolio and maximise their aid income using modern portfolio theory. The weighted network model introduced allows for additional variables to be incorporated and, indeed, further variables were considered for inclusion. Variables such as the occurrence of war, migration and recipient corruption could be reflected in the model; however, the focus here is on long-term and relatively stable determinants of aid. Furthermore, the bandwagon impact may partially and indirectly incorporate these variables; for example, the war in Afghanistan in the early 2000s led to significant amounts of aid donated by the USA to Afghanistan, swiftly followed by aid donated by other donors. 2.2.1. The General Weighted Network Model The general model proposed by [ 4 ] follows a weighted network model approach utilising donor-specific preference functions to measure donor motivations and biases when deciding aid allocations. The preference functions quantify the relative contributions of aid determinants used in aid allocation input decisions, such as poverty, trade activity, past colonial relationships and bandwagon impacts (‘herding’), into ‘weights’ which are applied to a network model, revealing donor behaviours and the relative importance placed on these aid determinants. Figure 2 is an archetypal bipartite network model which has nation donors on the left representing the set of nodes D and the recipients on the right representing the set of nodes R . D and R are disjoint sets of nodes in which links can exist only between the two sets and not within each set, thus illustrating the flow of aid which is directed only from elements of D to elements of R . In total, there are six nodes split into two disjoint sets of three donors, di∈D , and three recipients, rj∈R , where i , j= 1,2,3 represent donors and recipients, respectively. The n -vector node-specific information uk β represents the n quantities (or vector elements) associated with each country, k=i , j , in the network, where β= 1, ... , n is used to denote the numbered element of the vector uk . For example, recipient Ethiopia’s node in Figure 2 is labelled r3 . In the case that this node’s specific information contains poverty levels, u3 1 , colonial history, u3 2 , and trade activity, u3 3 , then the vector u3 has 3 elements n= 3, denoted as u3 β, where β=1, 2,3. Further, the links between each set holds m -vector link weights lij α , representing the m relationships between donors i and recipients j , where α= 1, ... , m denotes the numbered element of the vector lij . For example, the vector representing trade activity and colonial history between the UK, d2 , and Bangladesh, r1 , are denoted as l21 1and l21 2 , respectively, where m=2 in this example, and is denoted as l21 α, where α=1,2. 107 Entropy 2023,25, 641 Figure 2. A complete bipartite graph of a weighted network containing notation which underpins the model. The node and link vector information, as defined, can be quantified into weights and input into a preference Function (1), supplemented by input Functions (2) and (3), outputting a percentage of aid allocated by a donor, di, to a recipient, rj. Pijlij,uj:=∏ α fi αlij α,li• α.∏ β gi βuj β,u• β(1) The superscript • in preference Function (1) denotes all recipients in the set of nodes R , and its role is seen in the denominators of (2) and (3). The input functions, fi α and gi β , quantify donor preferences towards specific determinants of aid, such as trade activity and recipient need, into proportioned weights before input into (1): fi αlij α,li• α=lij α ∑k∈V2lik αμi α (2) gi βuj β,u• β=⎡ ⎣ uj β ∑k∈V2uk β⎤ ⎦ ηi β (3) where the exponent parameters, μi α and ηi β , hold only non-negative real values. These are referred to as ‘power’ parameters. The terms in the square brackets in (2) and (3) are functional inputs holding information on chosen aid determinants positively correlated with aid allocation and expressed as a proportion. The greater the proportion, the higher the value generated by the function and therefore the greater the weight due to a particular aid determinant for input into (1). There is a difference in usage between the input Functions (2) and (3). The function fi αlij α,li• α in (2) is used for link-specific weights, lij α , and quantifies behaviours and relationships that exist between a donor di and a recipient rj , for example the levels of trade activity. The function gi βuj β,u• β in (3) is used for node-specific weights, uj β , and quantifies a specific recipient metric, such as the poverty ratio among recipients, which is a determinant of aid that bears no direct relationship to a particular donor. Note uj β in 108 Entropy 2023,25, 641 (3) is specific to a recipient rj ; however, it can also be specific to a donor, ui β , to quantify determinants specific to a donor, di, such as donor affordability. The form of (2) and (3) assumes that a positive correlation exists between the determinant in question and aid allocated. For negative correlations, these terms are modified by subtracting the functions from (1), resulting in a recipient with a lesser proportion receiving a higher weight of preference and therefore aid allocated relative to the other recipients in the model. As an example, assume ‘recipient need’ was chosen to be an aid determinant by a particular donor. This variable can be measured in several ways. First, say it is measured by poverty levels per capita. This measure is assumed positively correlated with aid: the more people in poverty, then the more aid the recipient should attract, resulting in a relatively higher proportional output by Equation (3)—the recipient-specific weight—assuming the power parameter ηi β is unity. The output from (3) is then an input into the preference Function (1), resulting in a higher percentage of aid to that recipient. On the other hand, if recipient need was instead measured using recipient GNI, the output of Function (3) would need to be subtracted from (1) because recipient GNI is assumed negatively correlated with aid received. The result of these two approaches in terms of the output by preference Function (1) should be roughly equal. To allow for biases when allocating aid to certain recipients based on aid determinants, the power parameters μi α and ηi β on Functions (2) and (3), respectively, allow for a choice to be made by a donor with regards to the relative contribution and, thus, importance of particular determinants on the final aid allocations output by (1), which is in the form of percentages of the total aid budget. Donors can dial up or dial down the level of influence that their selected determinants have on the outcome by changing the values of the power parameters. For example, if a donor wanted to allocate more aid to recipients with which it experiences large amounts of trade activity over those recipients with higher poverty levels, then the donor will choose a higher value for the power parameter applicable to the relevant functional equation that is quantifying trade activity. These parameters then, also provide a means for deducing historical donor behaviours and biases in simulations. If a mathematical model is to be used by politicians, countries and organisations, then it needs to be simple, effective and able to be communicated. An important feature of this weighted network model is indeed its simplicity and transparency with the inputs into Equations (2) and (3) and, in turn, into preference Function (1), determined by verifiable, properly sourced, factual data. The weighted network model’s initial purpose was to reflect the decisions made by donors with regards to their motivations towards aid allocation based on certain factors, such as trade and recipient poverty. Donors can decide how much emphasis these factors have on the final allocation of their aid budgets. With historical data, sourced from the World Bank and OECD, input into Functions (2) and (3), and with the historical aid allocation figures which are the outputs from preference Function (1) also known, this leaves the power parameter values as the only unknowns. These values can be determined by playing the role of balancing figures and, from these estimated values, donor motivations and the relative importance placed on certain aid determinants can be studied. 2.2.2. Adaptation of the General Model Before adaptation of the general weighted network model for use by an aid recipient, the model Functions (1)–(3) need to first reflect the analysis conducted in Section 2.1 and the network model in Figure 2. Therefore, three donors, three recipients and the four identified significant aid determinants are to be incorporated into the general model. The set of donors, D , are Germany, the UK and the USA, respectively; d1,d2,d3∈D , each having its own preference function. The set, R , contains the three recipients: Bangladesh, Afghanistan and Ethiopia, respectively; r1,r2,r3∈R. 109 Entropy 2023,25, 641 From analysis performed, trade relationships were found to be positively correlated with aid and bilateral trade activity was identified as one of the four significant aid determinants that a donor considers when allocating their aid budgets. For incorporation into the general model, trade activity is to be represented by the variable tij . For example, the superscript i= 1 identifies the donor as Germany, d1 , and superscript j= 1 represents the recipient Bangladesh, r1 . The trade ratio between a donor and its aid recipients, measured in terms of exports to the recipient in US$, was calculated from [ 13 ]. For example, for Germany, t11 :t12 :t13 ≡13:1:5 for the year 2019. Recipient poverty was another significant aid determinant identified and is to be represented in the adapted model by pj , which denotes the poverty levels in recipient rj . The ratio of poverty levels is determined by using GNI per capita to quantify recipient need [ 9 ]. Correlation analysis indicated that this metric is negatively correlated with aid; therefore, there will be a subtraction from 1 in the poverty-specific functional input equation. Another of the significant aid determinants identified was colonial relationships, represented by the variable cij between donor di and recipient rj . This variable will be quantified using a binary zero-one integer programming variable defined as cij =1+1, colonial relationship existed between diand rj 0, no colonial relationship (4) The only colonial relationships of relevance to the simulation are Afghanistan and Bangladesh, both ex-UK colonies and protectorates, and thus c21 = 2 and c22 = 2; with cij = 1 for all other combinations of iand j, where i,j=1,2, 3. The final significant aid determinant identified was the bandwagon effect, or herding, discussed in Section 2.1.5. Simply, it refers to the tendency of aid donors to follow other donors in allocating aid to certain recipients, who then gain ‘star’ status in the network. This can reveal itself when a recipient attracts a larger proportion of the total aid donated for no discernible reason, controlling for other factors; see [ 16 ]. The weighted network model framework can quantify the bandwagon effect, to be denoted bj , by capturing the phenomenon using aid received by a recipient, rj , in the previous period as a proportion of total aid donated by all donors in the entire network. This captures the herding effect by measuring recipients’ previous success in receiving aid relative to other recipients, thus becoming a ‘star’ node in the model network. The four determinants have now been allocated specific variables and associated data to be input into an adapted model. The aid allocated by the three donors to the three recipients in Figure 2 is to be used as the output values of the adapted model’s preference function, sourced from [ 9 ]. Thus, the only remaining unknowns are the values of the four power parameters in the four input equations, each representing one of the four aid determinants. These parameter values can be estimated by running the model using the known inputs (the aid determinants) and known outputs (actual historical data) to provide important insights into donors’ individual and relative motivations and behaviours with regards to allocating their aid budgets. The higher the power parameter value, the more bias has been baked into the aid allocation output from that aid determinant. To use the model over multiple consecutive time periods, it needs to be made temporal. Starting at time t= 1, the weighted network model can be iterated forward in time with the outputs of the equations changing at each t due to the recipients’ economic response to aid receipts, which feed back into donors’ decisions on the allocation of aid at t+ 1, acting as a feedback mechanism. For example, assume aid donated at time t led to a reduction in poverty in a recipient. By rolling the model forward to the next time-period, t+ 1, this reduced level of poverty will be fed into the model at t+ 1, producing a different aid allocation output percentage for the donor at t+1 compared to t. By denoting recipient poverty as pj t , trade relationships as tij t , and the bandwagon effect as bj t where the subscript t represents the time-period, the model becomes dynamic 110 Entropy 2023,25, 641 with respect to time. Note that colonial history, cij , quantified in (4) is a static measure: it does not change with time and, consequently, has no subscript. After the alterations discussed, the preference Function (1) and input Functions (2) and (3) in the general model have been adapted to create Equations (5)–(9), with the four determinant functions in the preference Function (5) and four input Equations (6)–(9). Pijcij,tij t,pj t,bj t:=fi 1cij,ci•·fi 2tij t,ti• t·gi 1pj t,p• t·gi 2bj t,b• t(5) fi 1cij,ci•=cij ∑k∈V2cik μi 1 (6) fi 2tij t,ti• t=tij t ∑k∈V2tik tμi 2 (7) gi 1pj t,p• t=1−pj t ∑k∈V2pk tηi 1 (8) gi 2bj t,b• t=bj t ∑k∈V2bk tηi 2 =∑3 m=1Amj t−1 ∑3 k=1∑3 n=1Ank t−1ηi 2 (9) where Amj t−1 in (9) is the amount of aid donated by donor m to recipient j in the previous period t− 1, and the denominator of (9) quantifies the total aid donated within the network at period t−1. The model is iterated forward starting from t= 1 (year 2015) to t= 5 (year 2019). Note that t− 1 is the year 2014, for which actual aid data is input into Function (9). By iterating forward, the power parameters can be backward calculated for each year from 2015 to 2019. The values are shown in Table 2. Entries marked N/A for not applicable are included where the relationship was not relevant for the year in question. Table 2. Power parameters required by the model to recreate the actual aid allocation results for each year 2015 to 2019 for donors Germany, UK and USA and recipients Afghanistan, Bangladesh and Ethiopia. Data that is not relevant is labelled N/A for not applicable. Donor Aid Determinant 2015 2016 2017 2018 2019 Germany Colonial history N/A N/A N/A N/A N/A Trade relationship 1.5 1.4 2.0 0.8 - Poverty 1.1 1.6 1.1 2.0 1.0 Bandwagon 4.0 6.0 9.0 5.0 2.2 UK Colonial history - - 0.5 - - Trade relationship 0.4 0.4 0.4 0.1 - Poverty - - 0.6 0.2 0.1 Bandwagon 1.0 1.0 0.6 0.3 0.1 USA Colonial history N/A N/A N/A N/A N/A Trade relationship - - - 1.1 - Poverty 1.4 1.6 1.5 1.6 0.5 Bandwagon 1.0 0.6 0.4 0.3 0.9 111 Entropy 2023,25, 641 The values in Table 2 can be put into matrix form. For each year the model was iterated, the matrix of power parameter values was input into (6)–(9) and fed into (5) to create the next period’s aid allocations. For example, for the year 2015, the matrix was the following: Φi α=⎛ ⎜ ⎝ μ1 1μ1 2η1 1η1 2 μ2 1μ2 2η2 1η2 2 μ3 1μ3 2η3 1η3 2 ⎞ ⎟ ⎠=⎛ ⎜ ⎝ 0.0 1.5 1.1 4.0 0.0 0.4 0.0 1.0 0.0 0.0 0.4 1.0 ⎞ ⎟ ⎠(10) The general weighted network model represented by Equations (1) to (3) has been tested, modified to Equations (5) to (9), and the relevant power parameters now calculated. All significant inputs and outputs for the years 2015 to 2019 are now known. Next, the model is adapted for use by an aid recipient by incorporating modern portfolio theory before performing a simulation to illustrate how that recipient could optimise their aid receipts. 2.2.3. An Aid Recipient’s Investment Strategy Using Network Theory If a recipient invests in increasing its trade activity with donors, then that recipient should expect an increase in aid receipts from those donors who place a relatively high value on commercial trade in their aid allocation decisions. This increase would then be compounded by the herding effect, leading to additional receipts that can be re-invested back into trade activity or other similar investment ‘assets’, creating a virtuous cycle of investment and increasing returns. Paradoxically, recipients may not have an incentive to reduce poverty since it may lead to a fall in aid receipts. Instead, if recipients focus primarily on increasing trade activities, then their GNI should naturally increase and poverty should be reduced. This argument is limited however as it depends on other limiting factors such as the quality of governance and institutions in the recipient country. The fruits of increased trade activity may also fuel corruption rather than being devoted to alleviating poverty. Often, increased trade activity is performed by state-owned companies with the recipient’s President as the main shareholder. Donors may wish to accommodate this in their aid decisions, which the weighted network model can do. Despite these complications, the main interest here is regarding the ability of the weighted network model, illustrated by Figure 2 and Equations (5)–(9), to be used by a recipient as an investment tool to maximise their aid receipts. By creating a foreign aid network model, a recipient would initially discover how influential it is in the network using centrality measures, the links it holds with donors and those that it does not. Specific weights can be added to links and nodes containing proportions of aid received, trade activity and other recipient–donor dyad information. This network model could also indicate if the recipient should seek out new donors, invest in current donors or a combination of the two. Recipients can treat their aid network model much like a company seeking to attract funding. They could view the aid determinants used by donors as an ‘asset portfolio’, safeguarding and maximising the value of those assets by treating them as investments. Recipients can invest their aid income into the asset portfolio, for example by investing in trade relationships with donors. The recipient may also need to invest in other sub-activities such as governance quality and public relations activities, which the network model and portfolio can identify. An investment plan for a typical aid recipient is outlined as follows: Step (1): Create a weighted network model, providing insights into links, level of influence and current donors in the recipient’s foreign aid network. Analyse each donor’s aid determinant preferences, motivations and biases. Step (2): Produce an asset portfolio representing the donor preferences identified, e.g., trade activity and poverty alleviation, with the USA being highlighted as a highly influential donor. 112 Entropy 2023,25, 641 Step (3): Identify those assets in the portfolio that provide the highest returns, then invest in these. For example, the recipient could invest to increase trade activity with the USA, and in the related governance quality and infrastructure. Step (4): The investment should lead to higher returns in the form of increased aid income, which is re-invested into the asset investment portfolio; e.g., increased trade activity with the USA should lead to further aid receipts donated by the USA, which then feeds back into the donor’s aid allocation model for the following years. The herding phenomenon then compounds the effect. Treating aid determinants like assets in a portfolio implies the existence of an optimal mix of such variables which provides maximum return for minimal risk. There are in fact two main models that can be used for asset portfolio analysis: Modern Portfolio Theory (MPT) and the Capital Asset Pricing Model (CAPM). The CAPM model is more robust with fewer inputs; whereas the MPT model, though elegant, loses some practicality from the attempt to find asset returns, volatilities and correlations. Unfortunately, the CAPM model’s principal purpose is for modelling and pricing equity market assets and their equivalents, where risk and returns are measured against some trade index such as the FTSE100. There is no equivalent transparently priced market for aid determinants and, therefore, the CAPM model cannot be used here. Instead, MPT is used to illustrate the concept using a simulation. Let us assume that the portfolio for aid recipient rj contains two controllable assets, N= 2, ‘owned’ by the recipient: trade activity, tj t , and poverty, pj t , at time subscript t , denoted in a set by Pj t=#tj t,pj t$(11) These assets are ‘investable’ with varying risk-reward ratios and could be correlated or uncorrelated since increasing trade volumes do not always translate into reducing poverty, dependent on the recipient country and its regime as discussed earlier. For simplicity in this simulation, it is assumed that the assets are uncorrelated (ρ=0) ; however, equations can be adapted for the case when the assets are correlated and the correlation coefficient ρ=0, discussed in Section 3.1. The mean and variance of the two-asset portfolio (11) can be written as μPj t =Wμtj+(1−W)μpj(12) σ2 Pj t =W2σ2 tj+2W(1−W)ρtj t,pj t σtjσpj+(1−W)2σ2 pj(13) with the correlation between the assets subject to the constraint −1≤ρtj t,pj t≤1. In (12) and (13), W∈ [0,1] is a parameter that determines the proportion of aid receipts invested in trade activity, i.e., W is the weight of the trade activity ‘asset’, tj t , in the portfolio. The weight on the poverty alleviation ‘asset’, pj t, must be 1 −W, because N ∑ i=1 Wi=1 (14) Further, if the two assets are uncorrelated, then ρtj t,pj t = 0 and the variance (13) becomes σ2 Pj t =W2σ2 tj+(1−W)2σ2 pj(15) The value of parameter W is important since the mean and standard deviation of returns of each asset should technically be known. As W is varied, the risk and reward dynamics of the portfolio change in response. 113 Entropy 2023,25, 641 investments they want to make. An obvious one is trade activity: donors can use the model to ensure that more aid is allocated to countries that provide higher levels of trade activity with the donor, such as that seen between Germany and China [17]. Finally, by providing a framework to explore the properties of foreign aid networks and the impact that decision variables will have on those properties, including the final aid allocations, the weighted network model can help donors and recipients, and potentially multilateral organisations, with one of the issues associated with foreign aid: aid spillage, by reducing the costs arising from inappropriate use of foreign aid budgets. In conclusion, the weighted network model, underpinned by network theory, has been demonstrated to successfully model the international aid system and is able to shed new light on the complexity and interactions inherent in foreign aid networks. Author Contributions: Writing—review and editing, S.B. and J.S. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Informed Consent Statement: Not applicable. Data Availability Statement: Data regarding aid donations are obtained from reference [1]. Conflicts of Interest: The authors declare no conflict of interest. References 1. OECD. Development Finance Data. 10 June 2021. Available online: https://www.oecd.org/dac/financing-sustainable- development/development-finance-data/ (accessed on 27 February 2023). 2. Ramalingam, B. Aid on the Edge of Chaos; Oxford University Press: Oxford, UK, 2011. 3. Schraeder, P.J.; Hook, S.W.; Taylor, B. Clarifying the Foreign Aid Puzzle: A Comparison of American, Japanese, French, and Swedish Aid Flows. World Politics 1998,50, 294–323. [CrossRef] 4. Downes, R.J.; Bishop, S.R. Aid Allocation: A Complex Perspective. In Global Dynamics: Approaches from Complexity Science; Wiley: New York, NY, USA, 2016; pp. 271–290. 5. Alesina, A.; Dollar, D. Who gives foreign aid to whom and why? J. Econ. Growth 2000,5, 33–63. [CrossRef] 6. Harrigan, J.; Wang, C. A New Approach to the Allocation of Aid Among Developing Countries: Is the USA different from the Rest? World Dev. 2011,39, 1281–1293. [CrossRef] 7. McGillivray, M. Modelling Aid Allocation: Issues, Approaches and Results; WIDER Discussion Papers/World Institute for Development Economics: Helsinki, Finland, 2003. 8. Swiss, L. Foreign Aid Allocation from a Network Perspective: The Effect of Global Ties. Soc. Sci. Res. 2016 ,63, 111–123. [CrossRef] [PubMed] 9. World Bank. June 2021. Available online: https://data.worldbank.org/indicator (accessed on 27 February 2023). 10. OECD. Data Visualisations. 6 June 2021. Available online: https://www.oecd.org/dac/financing-sustainable-development/ datavisualisations/ (accessed on 27 February 2023). 11. Development Initiatives. Aid Data 2019–2020: Analysis of Trends before and during COVID. 8 February 2021. Available online: https://devinit.org/resources/aid-data-2019-2020-analysis-trends-before-during-covid/#section-1-3 (accessed on 27 February 2023). 12. OECD. DAC List of ODA Recipients. 5 June 2021. Available online: https://www.oecd.org/dac/financing-sustainable- development/development-finance-standards/daclist.htm (accessed on 27 February 2023). 13. WITS. World Integrated Trade Solution. 2021. Available online: https://wits.worldbank.org/ (accessed on 27 February 2023). 14. Hoeffler, A.; Outram, V. Need, Merit or Self-Interest-What Determines the Allocation of Aid? University of Oxford: Oxford, UK, 2008. 15. Bikhchandani, S.; Sharma, S. IMF. 1 March 2000. Available online: https://www.imf.org/external/pubs/ft/wp/2000/wp0048.pdf (accessed on 27 February 2023). 16. Frot, E.; Santiso, J. Herding in Aid Allocation. KYKLOS 2011,64, 54–74. [CrossRef] 17. Karnitschnig, M. Facing China. 10 September 2020. Available online: https://www.politico.eu/article/germany-china-economy- business-technology-industry-trade-security/ (accessed on 27 February 2023). Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. 120 Citation: Kyrylych, T.; Povstenko, Y. Multi-Criteria Analysis of Startup Investment Alternatives Using the Hierarchy Method. Entropy 2023,25, 723. https://doi.org/10.3390/ e25050723 Academic Editors: Stanisław Dro˙ zd˙ z and Panos Argyrakis Received: 24 February 2023 Revised: 7 April 2023 Accepted: 24 April 2023 Published: 27 April 2023 Copyright: © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). entropy Article Multi-Criteria Analysis of Startup Investment Alternatives Using the Hierarchy Method Tamara Kyrylych * and Yuriy Povstenko Department of Mathematics and Computer Sciences, Faculty of Science and Technology, Jan Dlugosz University in Czestochowa, al. Armii Krajowej 13/15, 42-200 Czestochowa, Poland *Correspondence: [email protected] Abstract: In this paper, we discuss the use of multi-criteria analysis for investment alternatives as a rational, transparent, and systematic approach that reveals the decision-making process during a study of influences and relationships in complex organizational systems. It is shown that this approach considers not only quantitative but also qualitative influences, statistical and individual properties of the object, and expert objective evaluation. We define the criteria for evaluating startup investment prerogatives, which are organized in thematic clusters (types of potential). To compare the investment alternatives, Saaty’s hierarchy method is used. As an example, the analysis of three startups is carried out based on the phase mechanism and Saaty’s analytic hierarchy process to identify investment appeal of startups according to their specific features. As a result, it is possible to diversify the risks of an investor through the allocation of resources between several projects, in accordance with the received vector of global priorities. Keywords: multi-criteria analysis; criteria composition; investment; startup; Saaty’s method; global priority vector; choosing alternatives 1. Introduction The positive tendencies towards economic development require updated business entities according to the current market conditions and the emergence of new structural units, all of which form a competitive economic system. The active development of any economy is not possible without the constant emergence of new economic enterprises. This process stimulates the formation of the market environment with healthy competition and ensures scientific and reproducible functioning. Currently, we observe the positive tendency towards building potential for realizing business ideas through the creation of startups, whose business concepts have been dictated by the needs of the modern society and industries. A startup is a strategic economic unit with innovative concepts with the potential to enter the market. First, we outline the essential features of startups: (1) the innovation of an idea; (2) the necessity of capital investment; (3) reproducibility (possibility to sell the inventive solution multiple times); (4) business expansion; (5) the existence of a detailed and structured business plan; (6) generally, a startup is a project in initial stages of implementation; (7) the possibility of significant growth of the project; (8) often, startups propose new technologies; (9) uniqueness; (10) the potential team of professionals; (11) the riskiness of the investments; (12) the concentration of management decisions by the startup founders; (13) the flexibility as well as quick and efficient adaptation to changes in the environment; Entropy 2023,25, 723. https://doi.org/10.3390/e25050723 https://www.mdpi.com/journal/entropy 121 Entropy 2023,25, 723 (14) the possibility to individualize the products, according to the demands of consumers; (15) the dependence on credit resources; (16) the close relations between the founder and the employees, etc. Currently, one of the biggest problems is finding investors for startups, the qualified and objective evaluation of the concepts in terms of costs and benefits for future investment, and the successful presentation of the project to investors. Often, this work is entrusted to consulting agencies that professionally evaluate innovative ideas. For the investor, it is important to have a final estimation containing not only a list of factors justifying the appropriateness of investments in the suggested startups but also the method used for comparing several investment alternatives. For an objective and comprehensive assessment of a startup, a large number of criteria should be taken into account; however, this complicates the evaluation process and prolongs its execution. To assess investment alternatives, many methods, mechanisms, techniques, and tools enable the investigation of investments from different points of view. Research has been concentrated in several directions: the economic basis of startups, the mechanisms of their initiation, and the behaviors of investors. The most substantiated and successful in practice are mathematical models that predict the best investment alternatives. Based on the startup founder’s viewpoint, a comprehensive analysis of investment alternatives should involve the requirements from the idea to launch, from the gathering and successful use of information to the potential of the startup’s innovation in a functioning market. This step-by-step mechanism for building a business was precisely outlined in [1]. The behaviors of investors (especially, business “angels”) towards newly created enterprises in the early stages of their development, the ways of evaluating those enterprises, and the interactions of investors and entrepreneurs were described in [ 2 , 3 ]. The basics of practical venture capital management and the details of the cooperation of venture capitalists and entrepreneurs were presented in [ 4 ]. Practical advice and the confirmation of the importance of a correct, accurate assessment of the business opportunities of startups were given in [ 5 ]. An analysis of venture capital from the viewpoint of current and future investing in an uncertain environment and the high level of competition confirms complexity of the investment choice [6]. An important step towards identifying the most attractive startup for investment involves not only formulating the list of criteria but also establishing their importance (weights). Today, many consulting companies use expert assignment methods to identify the weights of the criteria, but sometimes, the methods are too subjective and dependent on the composition of the expert team, the expert engagement, and lobbying interests. In this area, special attention is paid to the decision-making theory and the Saaty hierarchy method. The multi-criteria decision-making analysis, known as the analytic hierarchy process, was elaborated by Saaty [ 7 – 13 ]. This approach has been applied to many areas, such as economics, management, engineering, mathematics, information systems, cybernetics, mechanics, design, chemistry, health service, etc. The literature on this subject is considerable, including the following books [ 14 – 19 ] and review articles [ 20 – 33 ], where additional references can be found. The choice and the comparison of the criteria are important parts of decision-making. As the criteria and their weights can significantly influence decision-making, several approaches to solve this problem have been elaborated. In the analytic hierarchy process (AHP), several prioritization methods have been used for deriving weights, such as the eigenvalue (EV) method [ 8 , 10 , 34 , 35 ], the logarithmic least squares (LLS) method [ 36 , 37 ], the weighted least squares (WLS) method [ 38 , 39 ], the fuzzy preference programming (FPP) method [ 40 – 43 ], and the cosine maximization method (CMM) developed in [ 44 ]. A good description of several of the most-used methods was given by Srdjevic [ 45 ]. The main feature of the step-wise weight assessment ratio analysis (SWARA) [ 46 , 47 ] is the possibility to estimate the opinions of experts and interested groups according to the significance ratio of the criteria in the process of their weight determination. In the best–worst method (BMW) [ 48 – 50 ], two vectors of pair-wise comparison were used to determine the weights 122 Entropy 2023,25, 723 of the criteria. The full consistency method (FUCOM) [ 51 – 53 ] is based on the pairwise comparison of the criteria and the satisfaction of the mathematical transitivity conditions. The level based weight assessment (LBWA) model [54,55] is suitable for use in complex multi-criteria models with a large number of criteria, and it allows for the additional corrections of the values of the weight coefficients, depending on the preferences of the decision-makers. The main purpose of this article is to provide a comparison of several startups from the investor viewpoint. In this paper, we discuss the use of a multi-criteria analysis for investment alternatives as a rational, transparent, and systematic approach that reveals the decision-making process during the study of influences and relationships in complex organizational systems. The proposed methods can be useful for consulting agencies, investors, and also for startups founders, who can then assess their competitive position against offers from other competitors in the selected economic branch or industrial sector. The procedure of the startup assessment, especially during the initial stages of implementation (development, operation, execution phases, etc.) is often subjective and challenging, as it requires the determination and account of many indexes as well as extended expert consultation, the formation of criteria, and so on. We propose new criteria and a new criteria composition for evaluating the investment appeal of startups. As an example, we consider three alternative investments in startups: the production of LED traffic lights, the manufacture of information–reference electronic terminals, and the manufacture of rotor-reactive turbo-rotational heaters of liquids. The analysis of the three startups is carried out based on the phase mechanisms and Saaty’s analytic hierarchy process to identify the investment appeal of the startups accounting for their specific features. The consistency index, the consistency ratio, and the global priority vector are calculated. As a result, it is possible to diversify the risks of an investor through the allocation of resources among several projects, in accordance with the calculated vector of global priorities. 2. Criteria Composition for Evaluating Investment Attractiveness of Startups Based upon the review of the literature, the study of the practice of founding and launching startups, successful experiences of investing in startup enterprises, and the results of our previous research, we suggest the following criteria composition, which are consolidated into 12 blocks (Table 1). Similar grouping of sub-criteria into blocks was considered, for example, in [ 41 ]. We used several of the block-criteria discussed in [ 56 – 59 ], and then we supplemented and extended these according to our own criteria. This criteria could be adjusted according to the economic branch or industrial sector, according to the special features of the business plans presented to the investor. The criteria allow us to analyze the characteristics of startups in a variety of ways, and grouping the proposed criteria could enable potential investors to predetermine the priority groups of the criteria and use the proposed “sketch” of the influential factors to focus attention on the current trends. This criteria-composition model aims to draw the attention of the researcher (investor, consultant) not only on the “classical” list of basic investment indicators (such as payback period and the value of investments) but also to the governmental support of the industry, the innovation and autonomy of startups, time, and resources, as well as the social, scientific, technical, informational, and environmental characteristics. 123 Entropy 2023,25, 723 Table 1. Criteria composition for evaluating investment appeal of a startup. No Type of Potential Criteria 1. Value of investment 2. Payback period (PBP) 3. Expected profitability 1. Financial strength 4. Risk level 5. Full or partial investor control of the startup 6. Possibility of reverse repurchase (RRP) 7. Possibility of tranche-funding, depending on the stage of the project 2. Product/service 8. Availability of samples or models of the product potential 9. Startup position in the market 10. Forecasted level of demand for the product/service 11. Level of competition in the economic branch or industrial sector 3. Marketing potential 12. Evaluation of startup competitiveness 13. Significant target audience 14. Availability of marketing strategy 15. Requirements to attract and interact with customers within the startup initial stage 4. Organizational 16. Availability of organizational plan potential 17. Innovation of idea 5. Scientific and 18. Innovation of technology technical potential 19. Availability of project plan for technical realization 20. Availability of intellectual property rights 6. Staff potential 21. Availability of potential specialists 22. Uniqueness of specialists Potential of the governmental, 23. The level of development of economic branch or sector in which the startup 7. international, economic, will operate and political situation 24. The level of governmental support of industry branch 25. Period of project completion 8. Time potential 26. Stage of project development 27. Duration of product introductory period/start of retail service 9. Autonomy 28. Dependence of the startup on other economic branches potential or industrial sectors 29. Dependence of the startup on other similar projects 10. Ecological potential 30. Level of negative impact on the environment 11. Social potential 31. Accessibility of project’s social utility 12. Information potential 32. Availability, reliability, and quality of information in economic branch or industrial sector in which the startup will operate Figure 1 presents the structure of the Saaty method as the operational algorithm, indicating the priority of investments in startups. 124 Entropy 2023,25, 723  Creation of criteria composition and category groups (k – number of criteria; m – number of alternatives).  Selection and appointment of an expert council from potential investors.  Construction of a decision-making graph of the investment priority as a model of the dominant hierarchy.  Formation of a matrice of pairwise comparisons = (×) for establishing weights of criteria, using Saaty’s scale («scale 1-9») on the basis of consistency which provides that   =1,  =  , ,=1,      , as well as k matrices of pairwise comparisons  () =󰇡  () 󰇢× for m alternatives on the consistency basis  ()=1,  ()=  (), ,=1,;=1,. 5 Determination of priority vectors (normalized vectors of geometric means) in each row of the matrix of criteria weights and in each row of matrices of alternatives by the following formulas: =∏    ∑∏      , = 1, ,  () =∏ ()    ∑∏ ()      , = 1, k ,=1,. 6 Determination of the consistency index (CI) in the weight matrix. The following steps should be performed: - multiplication of the matrix = (×) on the right by the vector  ; ,=1,; - calculation of the maximal eigenvalue   of the matrix  ; - determination of the consistency index CI=( −)/(−1) . 7 Determination of the consistency index (CI)()) of matrices of alternatives. The following steps should be performed: - multiplication of the matrix () =󰇡 ()󰇢× on the right by the vector  (); ,=1,; - calculation of the maximal eigenvalues  () of the matrices  (); - determination of the consistency index (CI)() =( () −)/(−1) , =1,. 8 Calculation of the consistency ratio (CR) as a ratio of the consistency index (CI) and the random consistency index (RI): CR = CI/RI; (CR)() =(CI)()/(RI)() . It is considered as optimal when the consistency ratio does not exceed 0.1; when it is more than 0.2, the opinions of experts should be reviewed and, if necessary, experts should be changed. 9 Calculation of vector of global priorities: =∑    (),=1,. Figure 1. Saaty’s analytic hierarchy process for the identification of the investment appeal of startups based on their specific features. 3. Implementation of the Saaty Method for Identified Criteria Composition To illustrate practically the Saaty method, we analyze three investment alternatives of startups: the production of LED traffic lights, the manufacture of information–reference electronic terminals, and the manufacture of rotor-reactive turbo-rotational heaters of liquids. The structure of the method is first presented as the dominant hierarchy model in an oriented graph (Figure 2). After considering the business plans of three investment alternatives and establishing the criteria for assessing the prerogatives of investing in the compared startups, we identified the investment priorities. First, we determined the weights of the criteria according to the sequence of the algorithm; this was the fourth step of the hierarchical procedure, as shown in Figure 1. Table 2 presents the results of the criteria comparison for evaluating the startups using Saaty’s scale (“scale 1–9”) [ 60 , 61 ]. Therefore, we obtained the matrix of pairwise comparisons for establishing the weights of the criteria. The numbers 1–12 in 125 Entropy 2023,25, 723 the top row and the first column correspond to the name of the criteria in Table 1. The priority vector ( μi ) is calculated as the normalized geometric means in accordance with step 5 (see Figure 1). The column RM presents the results of the multiplication of the paired comparison matrix Bij on the right by the vector μj . The column DV is obtained by dividing the component of the vector in the column RM by the corresponding component of the vector μj . The approximation of the maximal eigenvalue is calculated as the arithmetic mean of the components of the vector in the column DV and equals λmax = 12.72. The consistency index CI =( 12.72 − 12 )/ 11 = 0.06545. According to [ 7 ], for k= 12, the random consistency index RI = 1.48; therefore, the consistency ratio is CR =CI/RI = 0.06545 and does not exceed 0.1. Scientific and technical potential Decision on investment priority Financial strength Type of potential Product / service potential Organizational potential Potential of the governmental, international, economic, and political situation Time potential Autonomy p otential Ecological potential Social potential Information potential Staff potential Marketing potential Production of LED traffic lights Manufacture of informationreference electronic terminals Manufacture of rotor-reactive turbo-rotational heaters of liquid Alternatives Figure 2. The dominant hierarchical representation of the problem of choosing investment alternatives in startups. Table 2. The matrix of pairwise comparisons to determine the validity of 12 groups of criteria. Groups of 1 2 345 6789101112 μiRM DV Criteria 113 342 12245 4 30.1893 2.3312 12.31 2 1/3 11/2 1 1/2 1/2 1 1/2 3 3 3 4 0.0800 1.1099 13.87 3 1/3 2 12 1 1/2 1 1 2 2 2 2 0.0911 1.1345 12.45 4 1/4 1 1/2 11/2 1/5 1 1/2 1 1 1 1 0.0490 0.6099 12.45 5 1/2 2 1 2 112122 2 20.1058 1.2968 12.26 6 1 2 251 11 2 2 2 2 2 0.1282 1.6571 12.93 7 1/2 1 1 1 1/2 1 11 1 1 1 1 0.0666 0.8525 12.80 8 1/2 2 1 2 1 1/2 1 12 2 2 2 0.0942 1.1661 12.38 9 1/4 1/3 1/2 1 1/2 1/2 1 1/2 11 1 1 0.0483 0.5950 12.32 10 1/5 1/3 1/2 1 1/2 1/2 1 1/2 1 11/2 1/2 0.0422 0.5329 12.63 11 1/4 1/3 1/2 1 1/2 1 1 1/2 1 2 11/2 0.0511 0.6742 13.19 12 1/3 1/4 1/2 1 1/2 1/2 1 1/2 1 2 2 10.0542 0.6975 12.87 A similar analysis was performed for the 12 matrices with 3 alternatives. The results for the group of criteria “Financial strength” are shown in Table 3. In this case, we obtain λ(1) max = 3.0183. The consistency index (CI)(1)=( 3.0183 − 3 )/ 2 = 0.0092. The random consistency index (RI)(1)= 0.52 for m= 3[ 7 ]. The consistency ratio (CR)(1)= 0.0158 and does not exceed 0.1. Taking into account the 12 criteria groups, the final results are shown 126 Entropy 2023,25, 723 in Table 4. The conducted research allows us to assert that the startup for manufacturing information–reference electronic terminals is most attractive for investment, as its global priority of 0.3855 is the highest among the analyzed investment proposals. At the same time, the values of the global priorities for the startups producing LED traffic lights and manufacturing rotor-reactive turbo-rotational heaters of liquids are equal to 0.2547 and 0.3599, respectively. Table 3. The matrix of pairwise comparisons for the group of criteria “Financial strength”. Production Manufacture of Manufacture of Priority Vector of LED Information– Rotor-Reactive (The Normalized Startup Traffic Reference Turbo-Rotational Vector of Geometric RM DV Lights Electronic Heaters Means) Terminals of Liquids ν(1) r Production of LED traffic lights 11/3 1 0.20984 0.63337 3.01835 Manufacture of information–reference 3 12 0.54994 1.65990 3.01833 electronic terminals Manufacture of rotor-reactive 1 1/2 10.24021 0.72503 3.01832 turbo-rotational heaters of liquids Table 4. The optimal choice of startups according to the investment alternatives, based on the groups of criteria. Manufacture of Manufacture of Production of LED Information–Rotor-Reactive Turbo- Investing Alternatives in Startups Traffic Lights Reference Electronic Rotational Heaters Terminals of Liquids No Groups of Criteria Priority Vectors 1. Financial strength 0.2098 0.5499 0.2402 2. Product/service potential 0.2000 0.4000 0.4000 3. Marketing potential 0.2000 0.4000 0.4000 4. Organizational potential 0.2500 0.5000 0.2500 5. Scientific and technical potential 0.1634 0.5396 0.2970 6. Staff potential 0.1958 0.3108 0.4934 7. Potential of governmental, international, economic, and political situation 0.2500 0.2500 0.5000 8. Time potential 0.5936 0.1571 0.2493 9. Autonomy potential 0.1634 0.2970 0.5396 10. Ecological potential 0.2500 0.5000 0.2500 11. Social potential 0.3333 0.3333 0.3333 12. Information potential 0.3325 0.1396 0.5278 13. Vector of global priorities 0.2547 0.3855 0.3599 4. Concluding Remarks New criteria and new criteria composition for the comparison of investment alternatives were proposed. Considering the sub-criteria could aid establishing weights of groups of criteria. The criteria and alternatives are mutually independent. A multi-criteria approach based on the analytic hierarchy method was used providing a gradual, clear, and logically structured assessment of the parameters of the given alternatives to ensure a successful solution. The proposed approach also has some limitations. For a large number of criteria and alternatives, Saaty’s scale 1–9 could not be enough. The decision-making process could also be time consuming for a large number of criteria and alternatives. For 127 Entropy 2023,25, 723 example, in the case of 32 sub-criteria, there appears a large matrix of pairwise comparisons, and for k> 15 in the literature there is no value of the random index (RI) and only an approximate estimation of the consistency ratio (CR) can be obtained. Therefore, we grouped the 32 new sub-criteria proposed in this study into 12 blocks (potentials). Despite these limitations, the AHP approach is one of the most popular and objective methods for multi-criteria decision-making. The proposed use of the Saaty method for optimal decision-making has a number of advantages, as well. It does not require the unification of the units of measurement for different criteria. It ensures the accuracy of the evaluation by increasing the possibility of intra-matching within the selected criteria. In addition, the presence of a numeric scale allows the relations between the factors to be clearly identified. Finally, this method is adaptable, enabling the criteria composition to be modified by adding or eliminating factors. We compared the maximal eigenvalues λmax obtained as the arithmetic mean of the vector in column DV and the value of λmax obtained using the available mathematical package. With a precision of four digits, the results were the same. It should be emphasized that the consistency ratio (CR) of the pairwise comparison of the 12 groups of criteria, as well as all the 12 consistency ratios CR(p) , p= 1,2, ... ,12, did not exceed 0.1; therefore, the evaluation was consistent. In the future, we are planning to extend our research to compare our results with results obtained by other techniques, in particular, using the Bellman–Zadeh fuzzy set approach. Author Contributions: Conceptualization, T.K. and Y.P.; methodology, T.K.; validation, Y.P.; formal analysis, T.K.; investigation, T.K. and Y.P.; writing—original draft preparation, T.K. and Y.P; writing—review and editing, T.K. and Y.P.; supervision, Y.P. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Institutional Review Board Statement: Not applicable. Data Availability Statement: Not applicable. Acknowledgments: The authors would like to thank the reviewers for the helpful comments that allowed us to improve the final version of the paper. Conflicts of Interest: The authors declare no conflict of interest. References 1. Blank, S.; Dorf, B. The Startup Owner’s Manual: The Step-by-Step Guide for Building a Great Company; K&S Ranch Press: Pescadero, CA, USA, 2012. 2. Benjamin, G.A.; Margulis, J.B. The Angel Investor’s Handbook: How to Profit from Early-Stage Investing; Bloomberg Press: Princeton, NJ, USA, 2001. 3. Belton, V.; Stewart, T. Multiple Criteria Decision Analysis: An Integrated Approach; Springer: New York, NY, USA, 2002. 4. Campbell, K. Smarter Ventures: A Survivor’s Guide to Venture Capital through the New Cycle; Prentice Hall: Harlow, UK, 2003. 5. Kessler, A. Eat People: Furthermore, Other Unapologetic Rules for Game-Changing Entrepreneurs; Penguin Group: New York, NY, USA, 2011. 6. Li, Y. Duration analysis of venture capital staging: A real options perspective. J. Bus. Ventur. 2008,23, 497–512. [CrossRef] 7. Saaty, T.L. The Analytical Hierarchy Process: Planning, Priority Setting, Resource Allocation; McGraw-Hill: New York, NY, USA, 1980. 8. Saaty, T.L. Multicriteria Decision Making: The Analytical Hierarchy Process; RWS Publications: Pittsburgh, PA, USA, 1988. 9. Saaty, T.L. Decision Making with Dependence and Feedback: The Analytical Network Process, 2nd ed.; RWS Publications: Pittsburgh, PA, USA, 2001. 10. Saaty, T.L. Fundamentals of Decision Making and Priority Theory with the Analytical Hierarchy Process, 2nd ed.; RWS Publications: Pittsburgh, PA, USA, 2006. 11. Saaty, T.L.; Vargas, L.G. Decision Making with the Analytical Network Process: Economic, Political, Social and Technological Applications with Benefits, Opportunities, Costs and Risks; Springer: New York, NY, USA, 2006. 12. Saaty, T.L. Decision Making for Leaders: The Analytical Hierarchy Process for Decisions in a Complex World, 3rd ed.; RWS Publications: Pittsburgh, PA, USA, 2012. 13. Saaty, T.L.; Vargas, L.G. Models, Methods, Concepts & Applications of the Analytical Hierarchy Process, 2nd ed.; Springer: New York, NY, USA, 2012. 14. Brunelli, M. Introduction to the Analytical Hierarchy Process; Springer: Cham, Switzerland, 2015. 128 Entropy 2023,25, 723 15. Roy, U.; Majumder, M. Vulnerability of Watersheds to Climate Change Assessed by Neural Network and Analytical Hierarchy Process; Springer: Singapore, 2016. 16. De Felice, F.; Saaty, T.L.; Petrillo, A. (Eds.) Applications and Theory of Analytical Hierarchy Process—Decision Making for Strategic Decisions; IntechOpen: London, UK, 2016. 17. Ozsahin, D.U.; Hüseyin Gökçeku¸s, H.; Uzun, B.; LaMoreaux, J. (Eds.) Application of Multi-Criteria Decision Analysis in Environmental and Civil Engineering; Springer: Cham, Switzerland, 2021. 18. Thakkar, J.J. Multi-Criteria Decision Making; Springer: Singapore, 2021. 19. Kulakowski, K. Understanding the Analytical Hierarchy Process; Chapman and Hall/CRC: Boca Raton, FL, USA, 2022. 20. Pohekar, S.D.; Ramachandran, M. Application of multi-criteria decision making to sustainable energy planning—A review. Renew. Sustain. Energy Rev. 2004,8, 365–381. [CrossRef] 21. Vaidya, O.S.; Kumar, S. Analytical hierarchy process: An overview of applications. Eur. J. Oper. Res. 2006,169, 1–29. [CrossRef] 22. Ho, W. Integrated analytical hierarchy process and its applications—A literature review. Eur. J. Oper. Res. 2008 ,186, 211–228. [CrossRef] 23. Liberatore, M.J.; Nydick, R.L. The analytical hierarchy process in medical and health care decision making: A literature review. Eur. J. Oper. Res. 2008,189, 294–307. [CrossRef] 24. Ishizaka, A.; Labib, A. Review of the main developments in the Analytical Hierarchy Process. Expert Syst. Appl. 2011 ,38, 14336–14345. 25. Subramanian, N.; Ramanathan, R. A review of applications of Analytical Hierarchy Process in operations management. Int. J. Prod. Econ. 2012,138, 215–241. [CrossRef] 26. Schmidt, K.; Aumann, I.; Hollander, I.; Damm, K.; von der Schulenburg, J.M.G. Applying the Analytical Hierarchy Process in healthcare research: A systematic literature review and evaluation of reporting. BMC Med. Inform. Decis. Mak. 2015 ,15, 112. [CrossRef] 27. Russo, R.F.S.M.; Camanho, R. Criteria in AHP: A systematic review of literature. Procedia Comput. Sci. 2015 ,55, 1123–1132. [CrossRef] 28. Nisel, S.; Özdemir, M. AHP/ANP in sports: A comprehensive literature review. Int. J. Anal. Hierarchy Process 2016,8, 405–429. 29. Emrouznejad, A.; Marra, M. The state of the art development of AHP (1979–2017): A literature review with a social network analysis. Int. J. Prod. Res. 2017,55, 6653–6675. [CrossRef] 30. Rajput, V.; Kumar, D.; Sharma, A.; Singh, S.; Rambhagat. A literature review on AHP (Analytical Hierarchy Process). J. Adv. Res. Appl. Sci. 2018,5, 349–355. 31. Darko, A.; Chan, A.P.C.; Ameyaw, E.E.; Owusu, E.K.; Pärn, E.; Edwards, D.J. Review of application of analytical hierarchy process (AHP) in construction. Int. J. Constr. Manag. 2019,19, 436–452. 32. Goyal, P.; Kumar, D.; Kumar, V. Application of multi-criteria decision analysis in the area of sustainability: A literature review. Int. J. Anal. Hierarchy Process 2020,12, 512–545. 33. Madzík, P.; Falát, L. State-of-the-art on analytical hierarchy process in the last 40 years: Literature review based on Latent Dirichlet Allocation topic modelling. PLoS ONE 2022,17, e0268777. [CrossRef] 34. Saaty, T.L.; Hu, G. Ranking by eigenvector versus other methods in the Analytical Hierarchy Process. Appl. Math. Lett. 1998 , 11, 121–125. [CrossRef] 35. Saaty, T.L. Decision-making with the AHP: Why is the principal eigenvector necessary. Eur. J. Oper. Res. 2003 ,145, 85–91. [CrossRef] 36. Crawford, G.B. The geometric mean procedure for estimating the scale of a judgment matrix. Math. Model. 1987 ,9, 327–334. [CrossRef] 37. Csató, L. A characterization of the Logarithmic Least Squares Method. Eur. J. Oper. Res. 2019,276, 212–216. [CrossRef] 38. Wang, L.; Xu, L.; Feng, S.; Meng, M.Q.-H.; Wang, K. Multi–Gaussian fitting for pulse waveform using Weighted Least Squares and multi-criteria decision making method. Comput. Biol. Med. 2013,43, 1661–1672. [CrossRef] 39. Wu, S.; Fu, Y.; Lai, K.K.; Leung, W.K.J. A Weighted Least-Square Dissimilarity Approach for multiple criteria ABC inventory classification. Asia-Pac. J. Oper. Res. 2018,35, 1850025. [CrossRef] 40. Mikhailov, L. Fuzzy programming method for deriving priorities in the Analytical Hierarchy Process. J. Oper. Res. Soc. 2000 , 51, 341–349. [CrossRef] 41. Wang, J.; Fan, K.; Wang, W. Integration of fuzzy AHP and FPP with TOPSIS methodology for aeroengine health assessment. Expert Syst. Appl. 2010,37, 8516–8526. [CrossRef] 42. Almulhim, T.; Mikhailov, L.; Xu, D.-L. A fuzzy group prioritization method for deriving weights and its software implementation. Int. J. Artif. Intell. Interact. Multimed. 2013,2, 7–14. [CrossRef] 43. Fallahpour, A.; Wong, K.Y.; Rajoo, S.; Olugu, E.U.; Nilashi, M.; Turskis, Z. A fuzzy decision support system for sustainable construction project selection: An integrated FPP-FIS model. J. Civ. Eng. Manag. 2020,26, 247–258. [CrossRef] 44. Kou, G.; Lin, C. A cosine maximization method for the priority vector derivation in AHP. Eur. J. Oper. Res. 2014 ,235, 225–235. [CrossRef] 45. Srdjevic, B. Combining different prioritization methods in the analytical hierarchy process synthesis. Comput. Oper. Res. 2005 , 32, 1897–1919. [CrossRef] 129 Entropy 2023,25, 784 From the standpoint of observation, we want to know among the C< , what fraction of them were observed to have transitions during T> . Labeling the observed transitions as C(o) >T, the fraction of transitions matching candidate links is given by ℘(o)=|C<∩C(o) >| |C<|(3) The statistic of interest is finally defined as the ratio between the two quantities, which we call the excess probability xW, given by xW=℘(o) ℘(s). (4) If xW is above 1 with a large degree of certainty (has a very small p -value), we conclude that the threshold W leads to an OLFN that is useful for prediction. To provide intuition for this statement, note that xW measures the network averaged increase in probability with respect to the random model that a transition during T> occurs along a pair of nodes with W or more transitions during T< . To illustrate, if xW is 2, the transitions actually observed during T> are twice as likely to occur along pairs of nodes that had W transitions during T< than what would be expected from the random model. Therefore, a value of xW> 1 (and the greater the better) means that transitions during T> are predictable on the basis of transitions during T< because they prefer to occur along candidate links by a factor of xW than along random node pairs. Finally, as a technical point, the p -values can be determined semi-analytically (or analytically in the case of the uncorrelated random model, where only the total number of system transitions is preserved) by using the methodology in [10]. 2.2.3. Career Sequences and Their Probability Distributions To study careers, we are interested in the non-degenerate version of the sequences encoded in Equation (1). To illustrate what this means, consider an employment sequence cα in which α spends from t to t+Δt working at location i ,or cα(t)=···=cα(t+Δt)=i but cα(t− 1 )=cα(t) and cα(t+Δt)=cα(t+Δt+ 1 ) . We will refer to such a time period of uninterrupted work at a given location as a spell. Since our primary interest is in the locations (or career steps) individuals take, we create a non-degenerate version of cα called uα such that only the location of a spell is recorded but not the number of time steps spent in a location. Thus, for example, if α ’s career is spent only in two locations, i and j and cα={i , i , i , ... , i , j , j , ... , j} , the corresponding career sequence is uα={i , j} . We should note that uα preserves temporal ordering so that if α first worked in location i and then in j , these appear in that same order in uα . Our sequences also possess the feature that if an individual were to return to a previous location, this would be captured in the sequence. Thus, an individual with an employment sequence of the form {i,i,j,j,i,i}would have career sequence {i,j,i}. The frequencies with which career sequences occur is very useful information because they offer insights on the sorts of choices individuals make under the constraints of the opportunities that become available within the organization (an individual cannot change into a job that is not offered, an important observation from the perspective of modeling made by the seminal work of White [ 24 ]). In order to understand how common or rare specific career sequences are, we define the distribution of observed sequences ˆ φ(u) , where uis the random variable of career sequences. For a given time period of observation, ˆ φ(u=u)=∑αδuα,u ∑αuα(5) where uα corresponds to the career sequence of α , δuα,u is the Kronecker delta equal to 1 when α ’s career matches the desired sequence u and 0 otherwise, and the denominator is 136 Entropy 2023,25, 784 the total number of distinct careers observed. Described intuitively, Equation (5) precisely defines how we count careers to determine their probability of occurring. Because careers can be sensitive to the initial location, we further specialize our analysis to distinguish careers on the basis of their initial location. Let us label the first location of career u as uo (or uα,o when the career refers to that of individual α ). Then, we are interested in the set of conditional distributions ˆ φ(u=u|uo=i)=∑αδuα,uδuα,o,i ∑αδuα,o,i . (6) 2.2.4. Temporal Statistics: Length of Service Research on manpower identified early on some important features about the study of workforces inside organizations. When the emphasis is not on specific individuals, manpower studies are very similar to population studies with one critical difference: in the latter, survival times of segments of the population can be known quite well and vary slowly over time (the number of people of a certain ethnicity of a given age) whereas in the former the population of employees is much more changeable [ 22 ]. Thus, the concept of the completed length of service emerged [6,25]. The key conceptual point still carries over in terms of career forecasting: as an individual enters an organization, it is important to anticipate how long that individual is likely to stay in the organization. For simplicity, we approach this question here in a similar way to the manpower literature. In fact, we hinted at this point already in our definition of employment sequences (Equation (1)), where we introduced the quantity τα to represent α ’s job tenure in the organization. This quantity corresponds to the length of service random variable τ . Given the e individuals in the data, to determine the length of service distribution ψ(τ) we exclude from E all those employment sequences for which the last location recorded occurs in the last time unit in the data. This is because at this point, we are not capable to tell if any of those individuals exit the organization in that very last time unit, or if they continue in the organization. Due to the sensitive nature of the data, we do not report the specific distribution of length of service of individuals in the organization, but use it in order to model careers in the ways we explain next (Section 2.2.5). 2.2.5. Markov Models of Career Sequences To test the usefulness of OLFNs in modeling the movement of personnel across an organization, we construct two Markov chains, one which relies solely on the network structure (based on [ 9 , 10 ]) and another that uses the network structure plus memory (when applicable) about the prior transition [ 18 ]. At the most basic level, Markov chains require that one defines states of the system and probabilities to go between states. Our method based solely on network structure uses as states the current job (node) held by an employee, and the probability to transition between jobs is estimated on the basis of the transitions made by all workers over some selected period of time of the data (for example, the first half of the years in the data). On the other hand, our method to include memory generally defines as a state the tuple made of the current and previous job a worker has held (with exceptions needed to handle the first job of the worker), and the transition probabilities are estimated from other workers and the last two jobs they held. We now describe these details. Let us start by clarifying that both models simply lead to the creation of simulated employment sequences and their associated career sequences. Since we mostly focus on career sequences, we introduce r(o) and r(1) to represent random career sequences respectively created from the Markov network model or the Markov model with one-step memory. These random variables are characterized by the distributions φ(o)(r(o)) and φ(1)(r(1)) . These distributions are created from a large number of model realizations. There are two kinds of such realizations. On the one hand, a single random walker can only 137 Entropy 2023,25, 784 generate a single career, not enough to generate useful distributions φ(o)(r(o)) or φ(1)(r(1)) . Therefore, to generate these distributions, we use Mw walkers which correspondingly generate Mw careers from which to create the distributions. A second way to introduce multiple realizations is to generate Md distributions φ(o)(r(o)) and φ(1)(r(1)) so that no one single realization of Mw walkers dominates the results. Our ultimate goal is to determine the quality of the models, which we do by defining below a set of metrics that compare each of these distributions to ˆ φ(u). A common feature to both models is the fact that individuals can begin to work at the organization at any time in any one of its locations. For the purposes of modelling their career sequences, one could ignore the specific point in time unless there were reasons to assume that temporal interactions play an important role. The initial location, on the other hand, is always relevant in terms of the number of either employment of career sequences generated. Thus, we should keep in mind that all the distributions we study are in reference to careers that start at each specific location (node) in the network. One last feature shared by both models is that the number of time steps an individual travels is drawn from the length of service distribution ψ(τ) . The effect is that each individual has a randomly drawn, fixed lifetime in the organization so that after time τ , the individual’s career sequence (either r(o) or r(1) ) is completed and counted toward the appropriate distribution. The model based on [ 10 ] makes use of the network structure but, deviating from that article, also includes weights to construct the transition rates between nodes. In the model, a simulated individual located at i at time step t has a probability pij to choose j as their next location, and this probability is constant in time. To determine pij , we make use of all the employment sequences in Equation (1). Such sequences can be used from the entire data (all the time points) or limited to parts of the time (e.g., T< which would require some small adjustments like redefining work spells). Assuming we are using the entire data, we first count the number of moves fij from node i to j on the basis of the number of sequences (and the number of times in that sequence) where a transition occurs from i to j . Concretely, fij =∑ α tα,o+τα−1 ∑ t=tα,o δi,cα(t)δj,cα(t+1). (7) This equation states that fij is given by the number of times any individual makes a transition from i to j . For the Markov process, the probability of the transition i to j is then given by the proportion of all transition out of i that go to j with respect to all transitions out of i,or pij =fij ∑jfij [i,j∈N]. (8) Note that the definitions of fij and pij include diagonal terms. Thus, the diagonal of the transition matrix of the Markov chain accounts for the very frequent occurrence of individuals remaining in their locations. In contrast to the pure network model, the model that keeps track of the previous step (if the career has visited at least one other node) makes use of a slightly more complicated transition matrix. Note that when an individual enters the network at a node and has not yet made transitions to other nodes, the model is applied as if it was the pure network model described above; only after one transition can memory begin to play a role. To make use of memory, let us focus on a node j . The probability that an individual transitions from j to h given that it had previously transitioned from i to j is based on the number of careers that have previously made the same sequence of moves. Therefore, if f(i,j),(j,h) is given by f(i,j),(j,h)=∑ α tα,o+τα−2 ∑ t=tα,o δi,cα(t)δj,cα(t+1)δh,cα(t+2), (9) 138 Entropy 2023,25, 784 the probability for an individual to go from jto hgiven that they came from iis given by p(i,j),(j,h)=f(i,j),(j,h) ∑hf(i,j),(j,h) [i,j,h∈N]. (10) In both types of models, it is possible that the probabilities are 0 for an individual to move beyond their current location. If that is the case, the individual merely remains in the node until either the simulation finishes or the number of time units τ assigned to the individual are complete. We should note that a single realization for a walker can last up to the length of time we choose to model. 2.2.6. Evaluating Predicted Career Sequences Next, we describe the metrics we use to assess the quality of the models. Essentially, we are interested in knowing whether the models tend to produce with high probability the careers actually observed, along with their observed frequencies. Symbolically, this is equivalent to testing for the similarity of the numerical values between ˆ φ(u=u|uo=i) and φ(m)(r(m)=r|ro=i) when u=r over the space of possibilities of u (the sample space), where m= 0,1 for the memoryless Markov model or the one-step memory model, respectively. As a practical matter, we note that because all careers are distinguished by their initial location i , all the quantities we define are computed according to their initial location. Stated in plain English, the data show certain career paths and the models try to imitate these. Therefore, evaluating the models is done by checking how “similar” the imitation created by the models is to the observed careers. In an ideal scenario, two distributions are similar if their sample spaces are similar and the probabilities of events (the elements of the sample space) are also similar. To be precise about what similar means, we now proceed to introduce several different quantitative measures of that similarity and highlight how each focuses on a particular aspect of that similarity. Let us first concentrate on the similarity between probabilities ˆ φ(u=u|uo=i) and φ(m)(r(m)=r|ro=i) . In this case, similarity means that the observed and modeled probabilities of the same career u starting at node i have similar values, i.e., φ(m)(r(m)= u|ro=i)≈ˆ φ(u=u|uo=i) . But this comparison has to be done carefully because for any given initial node i , u is not independent of other careers starting from i . Let us denote all the observed careers starting from i as U(i)={uα}α∈E;uo=i . Then, they are related by the fact that ∑u∈U(i)ˆ φ(u=u|uo=i)= 1 which is the normalization condition for ˆ φ . Modeled careers also satisfy a similar relation; calling the set of these careers R(m)(i)= {r(m) θ}{θ};ro=i for model m , they satisfy ∑u∈R(m)(i)φ(m)(r(m)=u|ro=i)= 1. Note that U(i)={uα}α∈E;uo=i and R(m)(i)={r(m) θ}{θ};ro=i are, respectively, the sample spaces of the observed and modeled careers starting at i . The relation between the probabilities of all careers starting at a single node means that it is not enough to know that one particular career u is such that φ(m)(r(m)=u|ro=i)≈ˆ φ(u=u|uo=i) . Instead, we need to know that the entire collection of careers starting from i have approximately equal values of probability between observation and model. An effective way to study this is through information theoretic methods. Here we apply the Jensen-Shannon divergence (JSD) for this purpose [ 19 ]. This quantity measures information divergence between distributions in such a way that, unlike the Kullback-Liebler divergence, is efficient in handling possible mismatches in the sample spaces of the distributions. Defining the entropy of a random 139 Entropy 2023,25, 784 variable X with distribution P(X) as H(P)=−∑XP(X)log P(X) , the JSD applied to ˆ φ(u=u|uo=i)and φ(m)(r(m)=r|ro=i)takes the form JSD(m)(i)=H1 2ˆ φ(u|uo=i)+1 2φ(m)(r(m)|ro=i) −1 2)H(ˆ φ(u|uo=i)) + H(φ(m)(r(m)|ro=i))*. (11) Intuitively, the Jensen-Shannon divergence measures how much information two distributions share, with a value of 0 if they share all information (the distributions are identical), and a maximum possible value of log(2)when one distribution has no information about the other. Since the distributions φ(m)(r(m)=r|ro=i)are generally different between different Monte Carlo realizations, we generate one JSD(m)(i) for each of the Md realizations. To perform a complete test in terms of JSD, we create two versions of it, one that computes the JSD between pairs of distributions φ(m)(r(m)=r|ro=i) emerging from the Monte Carlo realizations (providing Md(Md− 1 )/ 2 distinct values of JSD) and another comparing the real distribution ˆ φ(u|uo=i) of careers against the simulated distributions (providing Md values of JSD). To explain this strategy further (using Md realizations), note that the random distribution φ(m)(r(m)=r|ro=i) and ˆ φ(u|uo=i) are both sample distributions. First, the modeled distribution φ(m)(r(m)=r|ro=i) emerges from generating Mw walks that begin at i and generate a set of walks R(m)(i) . Second, the distribution ˆ φ(u|uo=i) is formed by all the observed careers beginning at i . Because both distributions emerge from a finite number of samples, even if either of the models m= 0 or 1 was perfectly correct, one cannot expect the two distributions to overlap perfectly. Thus, a more realistic evaluation of their similarity comes from observing how much ˆ φ(u|uo=i) typically differs from φ(m)(r(m)=r|ro=i) . This leads us to the need for creating Md versions of φ(m)(r(m)=r|ro=i) to compare against ˆ φ(u|uo=i) . When needed, we label each such realization by the index q= 1, ... , Md . Finally, note that the comparison between simulated career distributions allows us to develop a baseline for how well the observed career distribution is expected to match simulations. As a practical matter regarding numerical estimation of entropy, our situation is dominated by careers out of virtually all starting nodes where the most common career is to stay at that node; this means that we are able to estimate entropy via simple naive methods as in our case these are not particularly affected by problems such as those highlighted in the literature on entropy estimation [26–28]. Shifting to sample space testing, we introduce the Jaccard index which determines how similar two sets are by checking for the proportion of elements that are common between the sets; when both sets have the same elements the Jaccard index is 1, and when they share no elements it is 0. Thus, for a given location i , we define the Jaccard index J(m)(i)of node idue to model mas J(m)(i)=|U(i)∩R(m)(i)| |U(i)∪R(m)(i)|(12) which quantifies how much the sets U(i) and R(m)(i) resemble each other. Since R(m)(i) is a product of simulations, one does not expect J(m)(i) to be the same for every realization. One simple approach (that we adopt here) to deal with this is to create a union of the simulated careers, ∏Md q∪R(m) q(i) and compare this set with U(i) . Note that the choice to check against the union over R(m) q(i) is well justified on the basis that we are not after a test of probability, only sample space. As a final check, we introduce a ratio test for careers. This check is useful for several purposes. For one, it can identify particular career sequences that are especially rare compared to random expectation. Another advantage is that it can be put to use in 140 Entropy 2023,25, 784 generating career profiles for each starting node that provide a sense for how well the collection of modeled careers match the collection of observed careers. A final use comes as an alternative to the measurements from JSD and can be readily applied to obtaining full descriptions of a model over the entire network. All these depend on the definition d(m)(u,i)=logˆ φ(u=u|uo=i) φ(m)(r(m)=u|ro=i), (13) which compares the observed probability of career u with initial location i against its simulated probability. The quantity approaches 0 as the simulated and observed probabilities of a career become more similar (i.e., φ(m)(r(m)=u|ro=i)≈ˆ φ(u=u|uo=i) ). On the other hand, if a model overestimates the frequency of u , d(m)(u , i)> 0; if it is underestimated, d(m)(u,i)<0. Using d(m)(u,i)over all observed careers beginning at iprovides another way to test the models. This can be done, for a given node, by measuring the average d(m)(u , i) over observed career paths, or d(m)(i)=∑u∈U(i)d(m)(u,i) |U(i)|. (14) As indicated, this quantity can also serve as a measure of the quality of a model at the level of each individual starting point for careers. A related quantity that can be derived is the variance of d(m)(u,i), defined as var(d(m)(i)) = ∑u∈U(i))d(m)(u,i)−d(m)(i)*2 |U(i)|(15) which provides a measure of how well models capture the totality of the careers predicted to start at i. A final use for d(m)(u , i) is introduce is the creation of a profile for the effectiveness of each model to recover individual observed careers. Let us create a rank-ordered list of careers u∈U(i) so that u0 is the most probable career departing i (that is ˆ φ(u=u0|uo=i) >ˆ φ(u=u|uo=i) for u=u0 ). Similarly, u1 is the second most probable career from i , which means that ˆ φ(u=u0|uo=i)>ˆ φ(u=u1|uo=i)>ˆ φ(u=u|uo=i) for u=u0 , u1 . After ordering all careers, we can construct the curve π(m) i(c)=(c ,10 d(m)(uc,i)) where c= 0,1, ... , |U(i)|− 1. This profile for node i shows in decreasing order of importance how closely model m is able to reproduce careers in i . A perfect model will tend to produce a flat curve of the form (c ,1 ) . On the other hand, if some careers deviate strongly, there will be noticeable jumps. 3. Results 3.1. The Validity of Organizational Labor Flow Networks To verify that OLFNs are in fact informative, we apply the method in Section 2.2.2 where the time steps are monthly periods and the T< and T< are quarterly periods (3 months). To test that the information of previous transitions is strong enough, we simply impose W= 1 and measure the time series of x1 over the years of data we possess. Given our ability to choose the definition of locations, we explore the three versions mentioned above, operating units, occupational series code, and geographic location (in this case, at the state level). The results are shown in Figure 1. The model used corresponds to fixed strength of nodes based on candidate links, the most demanding test based on results from [ 10 ]. For all choices of the definition of location, the excess probabilities xW are considerably above 1 which means that defining and OLFN on the basis of any of these locations produces networks on which a walker (representing an employee) can travel along careers that are likely to be found in the real data. However, the value of xW is larger 141 Entropy 2023,25, 784 for units than other definitions of location (solid blue line). This result in interesting in that it reinforces the value of work done in [8–10] where nodes are defined on the basis of firms in the economy. The similarity is that, just like firms, operating units are the actual administrative units within which people work. Figure 1. Quarterly excess probabilities xW over the time frame of the data. We make the measurements with three different definitions of locations, operating units (dotted line), occupational series codes (circle), and US state where the employee is located (dashed line). The model corresponds to fixed strength of nodes based on candidate links, the most demanding test based on results from [ 10 ]. Even in this case, it is clear that xWis markedly above 1. Given the effectiveness of using operating units for predicting job change, we further explore this definition of network. In Figure 2, we study the effect of the threshold W on the excess probability xW . The temporal tracking is the same as in Figure 1. In this case, we see that increasing W leads to modest gains in predictive ability of the network, yet remaining within the same order of magnitude as W=1. Figure 2. Quarterly excess probabilities xW and the values of ℘(o) and ℘(r) across the time frame of the data for units in the AAW, tested across increasing W . The dotted lines correspond to W= 1, the circles to W= 2, and dashed lines to W= 3. The bundle of curves in the middle of the plot correspond to ℘(o) . The lower bundle of curves represent ℘(s)(W) due to random models. Finally, the excess probabilities xWare represented by the upper bundle of curves. 142 Entropy 2023,25, 784 Based on this analysis, we conclude that even a single observed transition ( W= 1) between a node pair has considerable predictive power regarding future transitions and therefore, in the absence of some pre-established tolerance level, we adopt even a single job transition to be an acceptable link in an OLFN. Clearly, our result confirms that the idea of OLFNs is not just theoretical, but one that actually captures real employment affinity and can help predict future job changes. The results of this analysis also dictate how we define the probabilities of transitions in our Markov models (see Section 2.2.5). 3.2. Structure of Organizational Labor Flow Networks Once networks are generated, we check their general topological characteristics. As indicated above, three possible definitions of nodes can be used, operational unit, occupational series code, or geography. However, given that geographic location appears to provide the smallest values of xW above 1, we concentrate on the topological features of the OLFNs generated with locations defined as operational units and occupational series. In Figure 3, we present the degree distributions Pr(k) of the OLFNs defined with W= 1, where k represented the degree of a node. Given the small number of nodes present in the network built on occupational codes, the degree distribution (left) does not seem to provide a clear structure. The effect of the number of nodes on this lack of structure is another reason why our expectations for obtaining systematic results based on state locations as nodes are low, further justifying our obviating this analysis (while there are about 100 distinct occupations, there are only 50 states in the US; this small number of nodes is unlikely to show much connectivity structure). Figure 3. Degree distributions for OLFNs defined by occupational series ( left ) or operational units ( right ). The plots are shown in log-log scale. For units, we add a reference solid line that decays as k−1.2. On the other hand, when the network is defined in terms of operational units, much more topological information can be seen. First, the right panel of Figure 3 exhibits a long tail distribution of degree, with close to two decades of steady, near-linear decay in double logarithmic scale which is consistent with a power-law. Assuming this shape of the degree distribution (a power-law), we find by inspection a decaying slope of a value of ≈− 1.2, or Pr(k)∼k−1.2 . This slope is close to the value observed in much on the literature on the firm-size distribution, known to show exponents in a range near − 1 but with considerable variation that includes the value − 1.2 (see [ 29 , 30 ]). Although here we are reporting the probability of a node to have degree k , the degree is a consequence of job transitions which are proportional to the size of units. Consequently, the exponent we measure can be directly compared to that of the firm size distribution. In addition, no prior empirical work has addressed the internal structure of firms, and the simulation studies that have been performed [ 31 ] have predicted that the distribution of unit sizes inside a firm should grow, which is the opposite prediction to our observations. The agreement between the exponent value found here and exponents in the literature on the firm-size distribution suggests that large organizations, even if they have highly 143 Entropy 2023,25, 784 controlled structures, somehow organize themselves in a way that mimics the organization of entire economies. After the seminal paper by Simon and Bonani recognizing this phenomenon [ 12 ], and given the abundant literature on this topic (see e.g. [ 13 , 17 ]), we do not attempt to explain this phenomenon here. However, we do note the importance of this finding in the context of this debate because it suggest that the phenomenon is a truly emergent feature of the functioning of economic entities. 3.3. Jensen-Shannon Divergence Moving beyond the macro-structure of the system, we now focus on the probabilistic structure of careers. For this purpose, we apply the JSD explained in Section 2.2.6. Given the limited value shown in defining careers in terms of geographic locations, we narrow our focus to operating units and occupational series only. Evaluating the numerical values of JSD requires establishing a baseline, as explained above, that compares careers among the random distributions versus the comparison of careers between a random and the observed distribution. In Figure 4, we illustrate the nature of the results of our analysis. The left panel contains JSD distributions for one illustrative occupational series code. The model without memory is represented by the red and green distributions. The red distribution corresponds to Md distinct values of the JSD between the distribution of observed careers ˆ φ commencing in the occupational code of interest and the Md modeled distributions φ(o) of careers starting at the same occupation code. In contrast, the green distribution is constructed from the Md(Md− 1 )/ 2 distinct JSD values that emerge from comparing all the pairs of distributions Md distributions φ(o) with each other. From the figure, we see that the green distribution among the random career realizations is characterized by lower values of JSD. This should be expected from the fact that the careers generated by the model are fundamentally similar to each other. The red distribution, in contrast, has larger values of JSD because observed and random careers need not be as similar. It is notable, however, that the JSD values are small indicating that both random models perform well. Figure 4. Distributions of values of JSD when comparing careers generated by random modeling and observed careers. Locations defined by operating unit on panel ( A ) and by occupational series are shown on panel ( B ). The model without memory can be seen on both panels, with the green distributions showing the pairwise comparisons between the random distributions of careers and the red showing the comparisons between the observed distribution against each of the random distributions. Similarly, the one-step memory model can also be seen on both panels (in blue), with the distributions showing the pairwise comparisons between the random distributions of careers and the orange showing the comparisons between the observed distribution against each of the random distributions. When memory is introduced, generating random career distributions φ(1) , the orange and blue JSD distributions emerge. Once again, the JSD values that emerge from comparing Md(Md− 1 )/ 2 random distributions in pairs are lower (blue) than the Md distinct JSD 144 Entropy 2023,25, 784 values that compare ˆ φ with φ(1) . We can also observe that memory lowers the JSD values of these distributions in comparison to the ones from the model without memory. Similar results can be gleaned when locations are redefined to operating units (Figure 4, right panel). While the examples presented in Figure 4 correspond to a particular unit and occupational code, the qualitative characteristics observed are consistent for the remaining nodes and definitions of locations. 3.4. Jaccard Index Having tested the similarity of the distributions, we are now in a position to determine if the structure of careers predicted by the models is similar to real observed careers. As explained above in Section 2.2.6, the Jaccard index eliminates the advantage that comes to popular careers when evaluated through the JSD. Instead, all careers are compared on equal footing, providing much more clarity about the difference between the models and the real-world. Although it would be perfectly informative to generate distributions of values of the Jaccard index, it is very useful to compare the two models we use directly on the basis of their ability to achieve large values of Jaccard index approaching 1. Each point in Figure 5 corresponds to a starting career location i (left are occupational series, right are units) where the horizontal coordinate represents the Jaccard index of U(i)∩)∏Md∪R(1)(i)* and the vertical coordinate to the Jaccard index of U(i)∩)∏Md∪R(0)(i)*. Figure 5. Comparison of Jaccard indices calculated from models with memory and without memory for units ( B ) and occupational codes ( A ) for locations in the OLFN . Each point is an initial node for careers. The horizontal coordinate captures the Jaccard index of the collected Md careers created in the one-step memory model, and the vertical the Jaccard index of the collected Md careers created with the memoryless model. The solid line highlights the diagonal of the plot. The results clearly illustrate the situation. The memoryless model is hardly ever able to approach the value 1, generating values that are almost exclusively confined in the range between 0 and 0.1 (with some exceptions). On the other hand, the model with one-step memory is partially successful at achieving Jaccard indices of 1, as well as generating other less optimal, yet better performing values between 0 and 1 in comparison to the memoryless model. 3.5. Career Profiles and Overall Evaluation of Career Forecasts In order to develop better intuition about the ability of our models to replicate observation, we also study the career profiles generated. In Figure 6 we present the profiles π(m) i(c) for the same operating unit and occupational code as those in Figure 4 with both the memoryless and one-step memory models. In both panels, it is clear that generally the one-step memory model performs better than the memoryless model. Deviations tend to be more attenuated. In both examples, the quality of the forecast of the most likely careers starting from each of the nodes (the points to the 145 Entropy 2023,25, 1140 limitations of traditional market equilibrium analyses, while emphasizing the importance of transaction uncertainty. The remaining sections of this paper are organized as follows. Section 2 formulates the functions of supply and demand based on the concept of willingness price. In Section 3, we analyze market performance, including transaction quantity, market surplus, and transaction entropy, using the rationing rate. Additionally, we compare the total entropy in centralized and decentralized markets and discuss policy implications based on the comparison results. Section 4 presents the simulation settings and results, demonstrating the generating process of each variable that characterizes market performance. In Section 5, we discuss the importance of market transaction uncertainty by highlighting the shortcomings of the Walrasian general equilibrium and Marshall partial equilibrium approaches. We also discuss the plausible applications of transaction uncertainty analyses in real-world scenarios. Section 6 draws the conclusions. 2. The Expression of Demand and Supply with Willingness Price A partial equilibrium analysis (PEA) is a widely used tool for understanding market performance. It argues that supply and demand collectively represent two sides of traders in a market, making it simple to analyze the consequence of their interaction by tracing the equilibrium point and social welfare implications [ 27 ]. However, the PEA also needs to be improved, since it fails to clearly identify how sellers and buyers constitute supply and demand curves correspondingly. To solve this problem, Wang and Stanley introduced the concept of willingness price and formulated supply and demand functions to restate the PEA in a goods market [ 28 ]. The major advantage of this approach is that the laws of supply and demand can be derived directly, and the efficiency of market equilibrium can be strictly proved. In this paper, we follow their approach to describe the supply and demand in a goods market. We assume that each trader is willing to make a trade of one unit of goods and has a willingness price before participating in the trade. For one seller, their willingness price is defined as the minimum price that they are willing to sell one unit of goods. On the other side, the willingness price of a buyer is defined as the maximum price that they are willing to spend for one unit of goods. Supposing that a seller with a willingness price vs meets a buyer with a willingness price vb , their deal can be made only if vb≥vs is valid. Although we cannot identify all traders’ willingness prices in real markets, we know that they exist there and govern whether a deal can be made or not. As all participants’ willingness prices are exogenously given, the willingness prices of sellers and buyers must have a distribution correspondingly. It is reasonable to assume that willingness prices spread over the domain of (0, + ∞ ). This spread can be characterized by probability density functions, fs(v) and fB(v) for sellers and buyers, respectively. Supposing that the numbers of the sellers and buyers are given exogenously, denoted as NS and NB , respectively, then we can use Fs(v)=Ns×fs(v) and FB(v)=NB×fB(v) to characterize such distributions. From the normalization condition, we have the integrals of Fs(v) and FB(v)over the whole region of willingness prices, which are Nsand NB, respectively, ∞ 0Fs(v)dv =Ns, (1) ∞ 0FB(v)dv =NB. (2) For any one seller, given a market price of p , they will make their choice by comparing the willingness price and market price, that is to say, the necessary condition for the seller to sell one unit of goods can be expressed as p≥vS. (3) Otherwise, the seller will withdraw their offer. 248 Entropy 2023,25, 1140 Equation (3) implies that only the sellers whose willingness price is not greater than the actual market price are willing to sell their goods. Combining (1) and (3), we can obtain the supply function with a given market price QS(p), which can be written as QS(p)=p 0Fs(v)dv. (4) The above rationale can also be applied to derive the demand function. For a buyer, only if his willingness price vB is higher than or equal to the market price p , will he buy one unit of goods in a market, i.e., p≤vB. (5) Otherwise, he will give up on his purchase. Combining (2) and (5), we can obtain the demand function with a given market price QD(p)of the market, which is given by, QD(p)=∞ pFB(v)dv. (6) As is well known, there are many factors that can affect the supply and demand in a market. From the expressions of supply and demand given by Equations (4) and (6), the implicit governing factor of supply and/or demand is the willingness prices of the market participants. Thus, we can infer that most relevant factors take their effects through the willingness prices of sellers and buyers. As a result, any change in any variable that impacts these willingness prices will have an impact on the supply and demand of the goods. In addition, the extent of a market determines the total quantities of the goods demanded and supplied, which also has an impact on the supply and demand functions. Another important inference of supply and demand functions is that we can prove the laws of supply and demand by taking a derivative of these two formulas. The first derivatives of the supply and demand functions can be expressed, respectively, as the following, dQS dp =Fs(p)>0, (7) dQD dp =−FB(p)<0. (8) The results show that the relationship between the quantity supplied and the market price is positive. In other words, the higher market price, the more goods supplied in the market. On the contrary, the relationship between the quantity demanded and the market price is negative. Fewer goods are demanded as the price rises. The interaction between supply and demand determines the equilibrium price level and quantity of transactions. Combining Equations (4) and (6), we can obtain the equilibrium price p=p∗ .The equilibrium transaction quantity T∗ can be derived directly, which can be expressed as, T∗=p∗ 0Fs(v)dv =∞ p∗FB(v)dv. (9) Figure 1 illustrates the supply and demand curves in a commodity market. The supply curve is upward sloping, and the demand curve is downward sloping. The cross-point of these two curves specifies the market equilibrium, which corresponds to the equilibrium quantity and market-clearing price of the market. 249 Entropy 2023,25, 1140 Figure 1. A simplified diagram of supply and demand curves in a market. The shortage region is marked in yellow color, while the surplus region is in green color. 3. The Market Performance with Formulated Supply and Demand Functions In this section, our primary focus is on evaluating various aspects of market performance using the newly formulated supply and demand functions. Specifically, we analyze three key dimensions: transaction quantity, market surplus, and market uncertainty caused by a quantity mismatch of the supply and demand in a disequilibrium market. To quantify this uncertainty, we propose the concept of transaction entropy, which is derived from information entropy. 3.1. The Quantity of Transactions Supply and demand represent two parties of a goods market, and their interaction determines not only the market price, but also the quantity of transactions. In this section, we set the market price as being given exogenously, and investigate how transaction quantity is determined by supply and demand as the price varies. The state of a market depends on the level of given price. The market is in equilibrium when the price makes the market clear. Otherwise, the market is in disequilibrium. This disequilibrium can be divided into two cases, one is shortage and the other is surplus. When the price is lower than the equilibrium level, it corresponds to a state of shortage, where there is more quantity demanded than the quantity supplied in the market. When the price is higher than the equilibrium level, it corresponds to a state of surplus, where there is more quantity supplied than the quantity demanded in the market. As shown in Figure 1, the regions of shortage and surplus are marked in yellow color and green color, respectively. According to the short-side principle, the realized quantity of transactions is determined by the short side. The short side refers to the trading party with fewer willing exchanges, and those with more are at the long side. At equilibrium, the quantity supplied is equal to the quantity demanded. In this case, the quantity of realized transactions T∗ given by Equation (9) is equal to the quantity supplied and demanded. In a shortage market, the quantity demanded exceeds the quantity supplied. Therefore, the quantity of realized transactions is determined by the quantity supplied. The expression of the realized quantity of transactions in a shortage market TST(p)can be expressed as, TST(p)=p 0Fs(v)dv p <p∗. (10) For a surplus market, the quantity demanded is less than the quantity supplied. In contrast, the quantity of realized transactions in a surplus market TSP(p) can be written as follows, TSP(p)=∞ pFB(v)dv p >p∗. (11) 250 Entropy 2023,25, 1140 Based on the preceding analyses, the transaction quantity in the various states of a market can be given by, T(p)=⎧ ⎪ ⎨ ⎪ ⎩ p 0Fs(v)dv p <p∗, p∗ 0Fs(v)dv =∞ p∗FB(v)dv p =p∗ ∞ pFB(v)dv p >p∗. , (12) Figure 2 shows the computational results of the relationship between transaction quantity and market price based on Expression (12), represented by the blue line. Obviously, the quantity of transactions increases with an increase in market price when p<p∗ , and decreases when p>p∗ . The quantity of transactions reaches its maximum when the market price attains its equilibrium level. Figure 2. The relationship between transaction quantity and market price. The blue line represents the computational results, while red dots denote the simulation results. The participants’ willingness prices and market prices are in the range of [ 2 , 18 ], and the market reaches equilibrium at a price of p∗=10. For further details about the simulation settings, see the section of Simulation Results. 3.2. Market Surplus 3.2.1. The Rationing Rates According to the short-side principle, we know that all the participants at the long side are willing to make transactions, nevertheless, some of them cannot achieve their desired outcome. Thus, we define the rationing rate as the ratio of the quantity of actual transactions to the quantity of desired exchanges. The sellers’ and buyers’ rationing rates can be used in the following analysis of market surplus and transaction entropy. Their expressions (Gsand GB) are given as follows, respectively, Gs=T QS , (13) GB=T QD . (14) It is obvious that Gs and GBare in the range of [0, 1]. The quantities supplied and demanded will change with a variation in the market price. Therefore, the level of rationing rate will be altered as the market price varies. When the market price equals the equilibrium one, the rationing rates of either sellers or buyers equal one. Thus, we obtain, Gs(p∗)=GB(p∗)=1. (15) In the shortage region, i.e., p<p∗ , all sellers can fulfill their willing exchanges, where only a portion of buyers can successfully match with the sellers and achieve their desired 251 Entropy 2023,25, 1140 transactions. As a result, the sellers’ rationing rate is one, while the buyers’ rationing rate would be less than 1. Thus, we obtain, Gs(p)=1, (16) GB(p)<1. (17) Meanwhile, with an increasing market price, there are more commodities supplied and less demanded. The rationing rate of sellers remains constant with the increase in price, while the rationing rate of buyers increases. We then obtain, dGs(p) dp =0, (18) dGB(p) dp >0. (19) In contrast, the above rationale can also be applied to the surplus region, where p>p∗ . The rationing rate of sellers is lower than 1, and the buyers’ rationing rate is one. Then, we obtain, GS(p)<1, (20) GB(p)=1. (21) The relationship between the rationing rate and market price in a surplus market can also be derived. In this case, as the market price increases, sellers are less likely to obtain their rations, because the quantity supplied increases while the quantity demanded decreases. Meanwhile, the rationing rate of the buyers will not change. The derivatives of the rationing rates of sellers and buyers have the following properties, dGs(p) dp <0, (22) dGB(p) dp =0. (23) Figure 3 depicts the dependence of these rationing rates on market price. As shown in this figure, when a market is in a shortage, the rationing rate of buyers is less than one, whereas the rationing rate of sellers is equal to one. In contrast, the rationing rate of sellers is smaller than 1, while the buyers’ rationing rate equals one when a market is in surplus. When a market is in equilibrium, the rationing rates of either the sellers or buyers are 1. 3.2.2. The Formulation of Market Surplus Market surplus, used to measure market efficiency, is another essential component of traditional market performance analyses. The surplus of one seller (buyer) can be defined as the difference between the actual (willingness) price and the willingness (actual) price. In the transactions of a goods market, only a portion of participants will be able to realize their willing exchanges, and a surplus will be generated. Therefore, it is reasonable to take rationing rates into account when formulizing the surplus of a market. For sellers, given a market price p, the total realized surplus of these sellers (Zsr) in the market could be calculated as follows, Zsr(p)=p 0Fs(v)(p−v)Gs(p)dv. (24) 252 Entropy 2023,25, 1140 On the other side, given a market price p , the total realized surplus of the buyers (ZBr) in the market could be given by, ZBr(p)=∞ pFB(v)(v−p)GB(p)dv. (25) The total realized market surplus for a price Zr(p)is the sum of them, i.e., Zr(p)=p 0Fs(v)(p−v)Gs(p)dv +∞ pFB(v)(v−p)GB(p)dv. (26) Taking the first derivatives of Equation (26), the expression of the relationship between the derivation of surplus and market price can be expressed as, ∂Zr(p) ∂p=p 0Fs(v)Gs(p)dv −∞ pFB(v)GB(p)dv +p 0Fs(v)(p−v)∂Gs(p) ∂pdv +∞ pFB(v)(v−p)∂GB(p) ∂pdv.(27) Combining Equations (4)–(6), (13) and (14), Equation (27) can be rewritten as, ∂Zr(p) ∂p=p 0Fs(v)(p−v)∂Gs(p) ∂pdv +∞ pFB(v)(v−p)∂GB(p) ∂pdv. (28) When the market is in a shortage, we can obtain the following expression by combining Equations (18), (19) and (28), ∂Zr(p) ∂p>0. (29) When the market is in surplus, we can obtain the following expression by combining Equations (22), (23) and (28), ∂Zr(p) ∂p<0. (30) Figure 4 depicts the relationship between market surplus and market price. From this figure, we can find that the market surplus increases when p<p∗ and decreases when p>p∗. When the market is at equilibrium, the market surplus attains its maximum. Figure 3. The relationship between rationing rates and market price. The green line and red dots represent the computational and simulation results of Gs(p) , respectively, while the orange line and blue dots are computational and simulation results of GB(p) , respectively. For details about the simulation settings, see the section of Simulation Results. 253 [Document text truncated for crawler view.]