The assumption of frictionless markets has been a longstanding topic of debate, attracting interest from both scholars and market practitioners. Broadly speaking, financial market frictions encompass any factors that impede natural trading and may include transaction costs, taxes, regulatory expenses, asset indivisibility, and liquidity constraints. Notably, the concept of market liquidity is a cornerstone in financial literature, with traditional liquidity measures generally categorized into four groups [1]: (a) transaction cost measures, such as the bid-ask spread and its extensions [2]; (b) volume-based measures, including transaction volume, turnover [3], and the Hui–Heubel Liquidity Ratio [4]; (c) price-based measures, for instance, the Amihud measure (ILLIQ) [5] and the Market Efficiency Coefficient (MEC) [6]; (d) market-impact measures, exemplified by the Market-Adjusted Liquidity Model developed by Hui and Heubel in 1984 [4].

The introduction of electronic trading systems in recent years have increased the attention toward another measure of liquidity: the order book. It captures the orders placed by traders to buy and sell stocks at different price points. In this direction, a very interesting work by Libman et al. [7] concludes that the deeper layers of the order book possess valuable information in the context of liquidity, a finding that is supported by other studies.

A radically different approach to modeling market liquidity is proposed by the Fractal Market Hypothesis (FMH) [8]. This posits that liquidity is driven by the heterogeneity of investment time horizons among market participants. Under normal conditions, investors operate on a variety of time horizons, from intraday (fractions of seconds or minutes) to several years (pension funds, institutional investors); liquidity originates from the different expectations (and probability assessments) of market participants looking at different time-scaled distributions of financial returns. Therefore, a lack of liquidity can occur when the distributions among time scales appear similar or when trading activity ceases on some relevant time scales and investors start trading on common (typically short-term) horizons. In these cases, usually occurring during periods of high volatility where rising uncertainty results in a shift toward quick buy-and-sell operations, trading can slow or even be freezed, simply because most of traders—who are supposedly rational—base their trading decisions on similar distributions which induce similar behaviors. Typically, short-term investors rely more heavily on technical information, whereas long-term investors focus on fundamental indicators [9]. In the event of a significant market crash that challenges fundamental assumptions, long-term investors may either exit the market or shift to short-term strategies [10]. The Fractal Market Hypothesis can also be used to assess the causes of financial (in)stability, as suggested by Anderson and Noss [11]. In particular, the analysis reveals some relevant implications to a number of ongoing debates regarding the regulation of financial markets and of their major participants. A review of the Fractal Market Hypothesis together with the basic principles of fractal geometry is provided by Blackledge and Lamphiere [12]. Specifically, they consider the intrinsic scaling properties that characterize a random walk and the Efficient Market Hypotheses and exploring the Lévy index, coupled with other metrics, such as the Lyapunov exponent and the volatility, a long-term forecasts analysis is developed. In a study of two developed market indices and one emerging market index, Karp and Van Vuuren [13] find a relationship between the change in a time series' fractal dimension (before breaching a threshold) and both the magnitude and direction of the subsequent change in the time series. This relationship was found to be prevalent during times of strong price persistence, characterized by elevated Hurst exponents, suggesting potential investment strategies.

From a mathematical perspective, while various methods—such as DFA, wavelet power spectra, variance scaling, and R/S analysis—have been proposed to quantify liquidity in terms of scale invariance (see [14] or [15] for a comprehensive discussion), most traditional approaches primarily examine single moments rather than full distributions, which may yield incomplete insights. The aim of this work is to present the implementation of a novel methodology that examines self-similarity and (il)liquidity by comparing the entire distributions of log-returns over pairs of time scales, rather than only specific moments. In particular, an algorithm and a ready-to-use Matlab software are discussed to evaluate liquidity through the analysis of a time indexed self-similarity matrix whose entries are the pairwise self-similarity parameters relating the log-returns distributions across different investment time horizons (from 1 day up to 6 months). Roughly speaking, the estimated parameters tell how diverse the paired distributions are.

The methodology we implement provides insights into both market liquidity and the evolution of return distributions across time scales. Thus, it can provide insights to market makers and regulators particularly during periods of market turbulence. Indeed, the 2007–2009 global financial crisis, as well as the recent pandemic and geopolitical crises, have underscored the importance of understanding liquidity and its underlying drivers in financial markets.

The remainder of this article is organized as follows: Section 2 outlines the theoretical framework and the self-similarity testing procedure employed in this study. Section 3 presents the proposed algorithm, and Section 4 concludes.

In subSection 2.1 the definition of self-similarity is introduced along with the justification of its usage in financial context. SubSection 2.2 shortly discusses how self-similarity can be evaluated using the Kolmogorov–Smirnov distribution. Both topics are introduced to clarify how the methodology has been implemented in the software introduced in Section 3.

The pseudocode Algorithm 1 provides a structured overview of the function SelfSimilarity, detailing the nested loop structure, the determination of rescaled sequences, and the two-sample Kolmogorov–Smirnov (KS) test implementation.

In detail, the algorithm iteratively estimates self-similarity on different time scales and intervals, aiming to find the parameter H that minimizes the KS statistic, thereby quantifying the self-similarity index between two paired lagged sequences. The function receives four input values: the time sequence X(t, m) (log-price for financial applications), which stores m paths of length t; the maximum number of time scales N and the indices of the startpoint and endpoint of the analysis. Within nested loops that iterate over time indices, time scales and a range of potential H values, the algorithm generates lagged sequences for each pair of scales and performs the KS test to assess distributional similarity. The KS statistic and its associated p-value are stored, enabling the identification of the optimal H for each pair of scales. Finally, the function outputs matrices containing the estimated self-similarity parameter, the KS statistic, and the p-value across all analyzed time scales and indices, providing a robust measure of self-similarity in the time series data.

Regarding the two-sample Kolmogorov–Smirnov test used by the algorithm, it should be noted that this test is readily available in the standard function libraries of most major programming languages. Implementations in languages such as Python, MATLAB, and R provide accessible and optimized functions that allow users to efficiently perform the KS test to assess whether two samples come from the same distribution. These library functions typically offer flexibility in terms of input parameters, enabling control over test options such as significance levels and alternative hypotheses.

To demonstrate the effectiveness of the proposed algorithm with real data, it was implemented in the MATLAB R2024a environment and tested on a computer with a 64-bit Windows operating system, an Intel Core i7-10510U CPU (1.80 GHz–2.30 GHz), and 16 GB of RAM. The input matrix X was constructed using the logarithmic prices of 183 assets continuously listed in the S&P500 index (Standard and Poor's, USA) from January 2000 to December 2023, totaling 6,037 observations. Even if 183 out of 500 individual stocks were consider as representative of the index, constructing the distribution from the all stocks that make up the market index makes it possible to eliminate the inertial effect that would be observed if—instead of the individual stocks—only the time-marginal distribution of the index were considered.

The average execution time for each iteration was ~30 s.

Figure 1 illustrates three contrasting market conditions in terms of liquidity: the top row represents a period of near-total liquidity shortage, revealed by the self-similarity surface (left-hand plot) and the rescaled log-return densities, reported for the pair a = 58 and b = 72 (right-hand plot) on March 16, 2020. The self-similarity surface, dominated by a uniform blue color across all investment pairs, indicates very low values of the self-similarity parameter, signifying an extremely low liquidity level in the market due to the fact that the distributions almost overlap as they are, the scaling exponent being close to zero. This outcome aligns with expectations, as this date was chosen during the U.S. quarantine or lockdown period, where markets were affected by heightened concerns regarding the COVID-19 pandemic's progression. Additional evidence appears in the right panel, where both rescaled densities exhibit a pronounced left-skewed tail.

In contrast, the plots in the mid row illustrate market liquidity conditions on June 8, 2022, when the S&P index rose by 1% or 39.25 points, with positive performance across most sectors. The self-similarity surface (left-hand plot) displays high values (yellow shading) for nearly all the considered investment pairs, indicating a high-liquidity market. In the corresponding right-hand plot, the rescaled log-return densities for the pair a = 83 and b = 100 are shown, featuring right-skewed tails. Finally, the plots in the bottom row illustrate a scenario in which self-similarity is rejected at a 95% significance level on August 29, 2008. In this case, “holes" appear in the similarity surface (blank areas), as the two densities (right-hand plot) are significantly different. At H = 0.4, the distance between the distributions is minimized but remains statistically significant; therefore, self-similarity is rejected (i.e., no value of the parameter H satisfies the distributional equality for the pair a = 3 and b = 91).

Figure 2 illustrates two summary measures for the daily distributions of Ĥ(t) and the estimated diameters δ^(ΨH^(t)). To improve the readability of the graphs, the entire dataset was downsampled using a 20-day step, corresponding to an average trading month. The median value and the InterQuartile Range (IQR = P₇₅(·) − P₂₅(·)) were then computed for the distributions of both Ĥ(t) and δ^(ΨH^(t)). The data are particularly interesting because the dynamics of P₅₀Ĥ(t) are fully consistent with the known estimates of the Hurst parameter reported in the literature (see, e.g. [21–24]). Furthermore, the fact that the values oscillate around 1/2 is of particular financial significance, as this is the only value consistent with the absence of arbitrage opportunities in the market—a sort of benchmark to which the market always returns after deviations that can be more or less prolonged depending on the shocks it experiences.

In this paper, we have implemented a new methodology for assessing market liquidity. Starting from the Fractal Market Hypothesis's (FMH) assumption that liquidity emerges from the heterogeneity of investment time scales, the scale invariance or self-similarity between pairs of investment horizons is examined with regard to the whole distribution rather than, as usual, its individual moments. Specifically, the self-similarity exponent is obtained as the unknown argument that minimizes the diameter of the space of rescaled distributions, which is observed to reduce to the Kolmogorov–Smirnov statistic. A numerical procedure is built to determine this value with respect to a number of time scales flowing through time and an algorithm provides as output a dynamical estimation of the pairwise self-similarity parameter of the considered scaled distributions. Implemented in MatLab, the algorithm uses nested cycles that, for each trading day, iterates over both time scales and a range of potential values of H, and performs the KS test to assess distributional similarity. The KS statistic and its associated value are stored, allowing the optimal H for each pair of scales to be identified. Low values of the self-similarity parameter, typically occurring during volatile periods where long-term investors may either exit the market or shift to short-term scales, indicate potential liquidity shortage. In contrast, high values of the self-similarity parameter describe liquid market phases in which each time horizon is populated by a sufficient share of market participants to match trades. Empirical evidence that proves the effectiveness of the proposed algorithm with real data is offered by analyzing 183 assets continuously listed in the in the S&P500 index (Standard and Poor's, USA) from January 2000 to December 2023.