On the statistics of scaling exponents and the Multiscaling Value at Risk

Scaling and multiscaling financial time series have been widely studied in the literature. The research on this topic is vast and still flourishing. One way to analyse the scaling properties of time series is through the estimation of scaling exponents. These exponents are recognized as being valuable measures to discriminate between random, persistent, and anti-persistent behaviours in time series. In the literature, several methods have been proposed to study the multiscaling property and in this paper we use the generalized Hurst exponent (GHE). On the base of this methodology, we propose a novel statistical procedure to robustly estimate and test the multiscaling property and we name it RNSGHE. This methodology, together with a combination of t-tests and F-tests to discriminated between real and spurious scaling. Moreover, we also introduce a new methodology to estimate the optimal aggregation time used in our methodology. We numerically validate our procedure on simulated time series using the Multifractal Random Walk (MRW) and then apply it to real financial data. We also present results for times series with and without anomalies and we compute the bias that such anomalies introduce in the measurement of the scaling exponents. Finally, we show how the use of proper scaling and multiscaling can ameliorate the estimation of risk measures such as Value at Risk (VaR). We also propose a methodology based on Monte Carlo simulation, that we name Multiscaling Value at Risk (MSVaR), which takes into account the statical properties of multiscaling time series. We show that by using this statistical procedure in combination with the robustly estimated multiscaling exponents, the one year forecasted MSVaR mimics the VaR on the annual data for the majority of the stocks analysed.


Introduction
Nowadays, scaling and multiscaling are widely accepted as empirical stylized facts in financial time series and they need to be properly addressed and analysed since they provide important information to risk and asset managers.The (multi)scaling property of time series is particularly important in risk management, especially when the model used assumes independence of asset returns.If this assumption is not met risk measures might be severely biased, especially if there is long-range dependence and this is acting with a different degree across the time series statistical moments.In particular, multiscaling has been adopted as a formalism in two different branches of quantitative finance, i.e. econophysics and mathematical finance.The former devoted most of the attention to price and returns series in order to understand the source of multifractality from an empirical and theoretical point of view (Mantegna and Stanley, 1999;Dacorogna et al., 2001;Mantegna and Stanley, 1995;Di Matteo, 2007;Calvet and Fisher, 2002;Lux, 2004;Lux and Marchesi, 1999;Di Matteo, Aste, and Dacorogna, 2005;Buonocore et al., 2019) and has recently identified a new stylized fact which relates (nonlinearly) the strength of multiscaling and the dependence between stocks (Buonocore et al., 2019).The latter instead builds on the work of (Gatheral, Jaisson, and Rosenbaum, 2018) on rough volatility and has been used to construct stochastic models with anti-persistent volatility dynamics (Gatheral, Jaisson, and Rosenbaum, 2018;Takaishi, 2019;Fukasawa, Takabatake, and Westphal, 2019;Livieri et al., 2018).Even if the research question comes from different perspectives, it is important to recognize the relevance that its study has in finance.Multiscaling have been understood to originate from one or more phenomenon related to trading dynamics.In particular, it can be attributed to the fat tails, the autocorrelation of the absolute value of log-returns, liquidity dynamics, or (non-linear)correlation between high and low returns generated by the different time horizon of traders and the consequent volumes traded.It can also be caused by the endogeneity of markets for which a given order generates many other orders, this occurs especially in markets where algorithmic trading is prevalent.There are different methodologies used to extract the scaling exponent from time series.Among all, the Multifractal Detrended Fluctuation Analysis (MFDFA) proposed in (Kantelhardt et al., 2002), the Wavelet Transform Modulus Maxima (WTMM) introduced by (Muzy, Bacry, andArneodo, 1991, 1993) and the Structure function approach also known as the Generalized Hurst exponent (GHE) method (Di Matteo, Aste, and Dacorogna, 2003;Di Matteo, 2007;Sornette et al., 2018).In a recent paper, (Barunik and Kristoufek, 2010) tested different methodologies against some data specification and empirically showed that the GHE approach overperforms the other models.For this reason, throughout this work we will use the GHE method.Notwithstanding the importance of the correct estimation of the Hurst exponent, the analysis has been rarely addressed from a statistical point of view.In this paper, we propose a step-by-step procedure for a robust estimation and test of the multiscaling property.Application to simulated data and empirical data allows us also to demonstrate the impact of bias on scaling and multiscaling estimation.After the estimation and test of multiscaling, we also show how the use of proper scaling and multiscaling can ameliorate the estimation of risk measures such as Value at Risk (VaR).We finally propose a methodology based on Monte Carlo simulation, that we name Multiscaling VaR, which takes into account the statistical properties of multiscaling time series by using a multiscaling consistent data generating process.The paper is structured as follows.Sections 2 and 3 are devoted to a brief description of multiscaling in finance and to the statistical procedure proposed to consistently estimate and test the scaling spectrum.Section 4 shows the results of the proposed methodology applied to synthetic data while Section 5 reports the results of an empirical application to real financial time series.Section 6 is devoted to a practical application of scaling and multiscaling property to VaR while Section 7 concludes.

Multiscaling in finance
In this Section, we explain the importance of the multifractal (multiscaling) formalism in financial markets.Let us first fix the notation by defining the prices time series as P t and the log-prices p t = ln(P t ).From this, the log-returns over a time aggregation τ are r τ (t) = p (t+τ ) −p t , where τ is expressed in days.Financial models are usually based on the assumption that log-prices follow a Brownian Motion and the for this model, the rescaled second moment of the log-returns over time aggregation τ follows where σ τ is the standard deviation at aggregation horizon τ while σ is the standard deviation at daily aggregation.This Equation is usually referred to as the square root of time rule and it is widely applied in quantitative finance Danielsson and Zigrand (2006); Wang, Yeh, and Cheng (2011).Examples are the Black and Scholes model in which the volatility evolves as στ 1 2 , or the VaR which under Basel regulatory framework can be computed for higher time aggregation, e.g. the τ days VaR can be computed as the daily VaR multiplied by τ 1 2 .In the analysis of the Nile river, Hurst found that the scaling behaviour described by a Brownian Motion was not in line with the empirical data (Hurst, 1956).This formalism has been later introduced in finance (Mandelbrot, 1963(Mandelbrot, , 2013;;Sornette et al., 2018).To detect multiscaling, it is necessary to study the non-linearity of the scaling exponents of the q-order moments of the absolute value of log-returns (Di Matteo, 2007).In particular, the process (p t ) with stationary increments is said to be multiscaling if where q = {q 1 , q 2 , . . ., q M } is the set of evaluated moments, τ = {τ 1 , τ 2 , . . ., τ N } is the set of lags used to compute the log-returns, K q is the q-moment for τ = 1, H q is the so called generalized Hurst exponent which is a function of q.Finally, the function qH q is concave (Mandelbrot, Fisher, and Calvet, 1997) and codifies the scaling exponents of the process.A process is uniscaling when the function H q does not depend on q, i.e.H q = H (Di Matteo, 2007), while it is multiscaling otherwise.If H = 0.5, the process does not behave as a standard Brownian Motion and neglecting this feature, would significantly bias the estimation of the true risk.In particular, if H < 0.5(H > 0.5) the process is said to be anti-persistent (persistent) while if H = 0.5 the process has independent increments.Given Equation 2, a possible way to define a multiscaling proxy is by quantifying the degree of non-linearity of the function qH q .The standard procedure used in order to extract qH q consists in running a linear regression in log-log scale of Equation 2, which reads as where τ is defined in the range τ = [τ min , τ max ] and q = [q min , q max ] (Di Matteo, 2007).A multiscaling proxy can be obtained by fitting the measured scaling exponent with a second degree polynomial (Buonocore, Aste, and Di Matteo, 2016;Buonocore et al., 2019) of the form where A and B are two constants.In this mathematical setting, the measured B, B, represents the curvature of qH q .If B = 0, the process is uniscaling, while if B = 0, the process is multiscaling (Buonocore, Aste, and Di Matteo, 2016;Buonocore et al., 2019).What is of vital importance in oder to widely apply the multiscaling formalism in finance is the ability to correctly estimate the value of qH q and consequently, of A and B.

Methodology
As highlighted in the previous Sections, estimating the Hurst exponent from empirical data is a difficult task.The difficulties can be categorized in two different classes: • Those due to the statistical procedure adopted; • Those linked to the financial data themselves.
In the first class, we identify two main issues related to the following two points: • The statistical model used to compute the scaling exponents; • The input variables used in the statistical procedure.
In the second class, issues arise mainly from the following question: • If the data contain an anomaly, how this is impacting the estimation of the scaling and multiscaling exponents?
In this work we first focus on the statistical procedure and the implication on financial times series with and without anomalies.Then, we we discuss practical implications related to finance with a special attention to Value at Risk.

Statistical procedure
Multiscaling properties of financial time series have been understood to come from one or more phenomena related to trading dynamics.From the point of view of the financial microstructure, scaling can be attributed to liquidity dynamics, endogeneity of markets or any other trading dynamic existing in the market.In particular, the superimposition of different strategies and investment horizons generates the long-range dependence and its variation when evaluated at different order moments.In this Section we propose a methodology to estimate the Hurst exponent qH q and the multiscaling depth (curvature) coefficient B in a robust manner.As specified in Equation 3, the estimation of scaling laws is generally performed through a linear regression in log-log scales.The statistical problem which might arise in this context is that the regression is performed minimizing the squared log-errors instead of the true errors.This procedure might, in case of strong deviation from the assumed statistical model for the errors, severely impact the results.The solution to this problem consists of applying a nonlinear regression to the original (i.e.not transformed) data, comparing the fit of the two specifications to the original data and using the one which performs better.Another issue related to the statistical model is the variation of the intercept for the q regressions.In particular, we can exactly compute the value of K q rather than estimating it, thus eliminating possible errors and bias.We can define the standardized Ξ(τ, q) as from which it is possible to notice that Equation 2 becomes Equation 6 eliminates the possible bias introduced by the estimation of K q via regression.
To easily exploit and model the multiscaling behaviour, we define the q-order normalized moment as ... Ξ (τ, q) = Ξ(τ, q) 1 q (7) which transforms Equation 6 in ... Ξ (τ, q) ∼ τ Hq .(8) Within this new formulation, the analysis is much easier since now, all the q regressions have a 0 intercept and the multiscaling is present only if the regression coefficients H q differ for different values of q.In fact, for uniscaling time series all regression lines are overlapped while, for multiscaling time series they diverge.Given the formalism introduced by Equation 8, it is easy to check with the naked eye whether a process is multiscaling or not.In addition, we can now rewrite Equation 4 for the normalized and standardized structure function of Equation 8as This Equation, even if mathematically equivalent to Equation 4, has a statistical advantage.
Eliminating the multiplication by q from both sides of the Equation, we reduce the possibility of spurious results in case q is a dominant factor in the multiplication.Indeed, the interpretation is equivalent, i.e.A is the linear scaling index while B is the multiscaling proxy.Let us now define the relative structure function between two consecutive moments q i and q j , with q j > q i as ... Ξ (τ, q i , q j ) = ... Ξ (τ, q j ) ...
This formalization helps in the statistical analysis since we can now test if a process is statistically multiscaling using a significance test on the estimated H(q i , q j ).In fact, for uniscaling time series we have that H q = H, which implies that the difference between different order moments is always 0.2 On the contrary, for multiscaling time series it should be different from 0 for all q.This reduces to a t-test on the regression coefficients estimated using Equation 10.Besides the multi-regression approach, it is possible to perform a multivariate regression by rewriting Equation 10as where M is the maximum number of moments used.This is a multivariate nonlinear regression which can be easily solved via nonlinear optimization algorithms.Such methodological approach implies a possible relationship between the q-moments used in the regression.Depending on the model assumptions, one can use Equation 11or perform M separate regressions for each exponent.In the first case, it is then possible to use an F-test to test if all the coefficient except for the first one (H(0, q 1 )) are jointly equal to 0 against the alternative that some coefficients are different from 0. This is a less restrictive multiscaling test compared to the multiple t-tests.We call strongly multiscaling processes those processes which reject both the null hypothesis for all the t-tests and the null of the F-test.Conversely, we call weakly multiscaling processes those processes for which the null hypothesis of all the t-tests is rejected but not the null of the F-test.3This is quite intuitive since if a process is multiscaling, all the relative increments are statistically significant.However, if the process reconstructed with a single exponent is statistically equivalent to the one reconstructed with the full multiscaling spectrum, it means that such multiscaling behaviour is weak.As already mentioned, it is recommended to perform the model in the log-log scale and in the original coordinate system, and base the choice of the model on the goodness of fit.

The choice of q and τ
An important point in the statistical evaluation of the multiscaling exponent is the choice of q and τ .They must be selected using specific statistical criteria.The wrong values of q and τ can severely bias the result.Regarding q, many research papers about multiscaling systems propose the use of a vast spectrum of qs.This approach has two fallacies.The first one lies in the the fact that multiscaling processes are such even for small values of q.Secondly and most importantly, given a distribution of returns with tail exponent α, for q ≥ α we have that E[r q ] diverges.Hence, we need to have q < α.Any multiscaling behaviour found by neglecting or ignoring this fact, is severely biased and possibly false.The method used to set q can come from two different approaches, i.e. established research results or direct tail computation.Since it is empirically proven that for financial time series α ∈ (1, 2) (fat tails), a conservative approach would be to use q ≤ 1.Alternatively, it is possible to estimate α on the empirical distribution through a tail estimator (e.g.(Clauset, Shalizi, and Newman, 2009;Virkar, Clauset et al., 2014)) and use it as threshold for q.In this paper, we will follow the conservative approach and use q ≤ 1.In fact, as already stated if the multiscaling phenomenon is present, it can be extrapolated from this range of moments.
Regarding the time aggregation τ , a general rule would be to use the minimum possible value of τ ,τ * , such that the autocorrelation information of the series is preserved.The autocorrelation ρ of the return series at lag τ is defined as: where µ and σ 2 τ are respectively the mean and variance of r(t).It is a well known stylized fact that returns are expected to be uncorrelated at daily frequency while the absolute and squared returns exhibit long-persistence (Cont, 2001;Chakraborti et al., 2011).Among the different procedures used to estimate τ * , we mention: • Segmented regression on the structure function (Yue et al., 2017); • Autocorrelation significance test (Buonocore, Aste, and Di Matteo, 2017).
The first one computes the structure function for each q-moment and then, fits a segmented regression in log-log coordinates between τ and Ξ(τ, q), and finds two slopes -one for the scaling component and one for the non-scaling component (Yue et al., 2017).The second approach instead, chooses the value of τ prior to compute the structure function, setting the value of τ * as the minimum value of τ for which the autocorrelation is not statistically significant (Buonocore, Aste, and Di Matteo, 2017).In this paper we propose a new approach which takes the advantages of both methods.We name it Autocorrelation Segmented Regression.The rationale behind this approach is to perform a segmented regression on the autocovariance function computed on the absolute returns and take τ * as the value which minimizes the sum of squared residuals for the high autocorrelation state and the random noise state (plateau). 4This approach has the advantage of setting the value of τ in advance avoiding ad-hock solutions and reducing computations.Nevertheless, the method is less sensitive to a unique non-significant lag.In fact, in noisy data it is possible to have a lag for which the autocorrelation is not significant unless for a considerable amount of consequent lags.The Equation for the proposed Autocorrelation Segmented Regression (ACSR) takes the form where α is the intercept of the regression and can be fixed to be equal to ρ τ with τ = 1, β is a slope parameter for the autocorrelation function, τ is the lag at which the autocorrelation is computed, and τ * is the value of aggregation which maximizes the autocorrelation information. 5We use ln(τ ) instead of τ for better estimation of τ * .Figure 1 shows how this method works.We generated a process with known τ max = 250 and run the ACSR on the autocorrelation function.As shown in the Figure, we get a value of ln(τ * ) = 5.51 which corresponds to τ * = 247.

Multiscaling estimation and testing procedure
Before turning the attention to the simulation experiments, let us here recall the full procedure to follow in order to robustly extract the scaling exponents: 1. Compute τ * with the Autocorrelation Segmented Regression method; 2. Compute q = α or rely on the empirical evidence available in the literature; 3. Perform the linear and nonlinear regressions with the above parameters (equation 11); 5. Compute the multiscaling curvature using Equation 9and test for statistical significance.
Concerning the last point, in this paper we propose a full procedure to run what we call the multiscaling test.The testing procedure is divided in four steps.In the first step, we test if each scaling increment H(q i , q j ) is statistically significant through a t-test.The second step is devoted to the F-test.In particular, we perform the F-test using the predicted relative moments from the regression with the full estimated scaling spectrum and the one in which only the first scaling is different from 0, i.e.
... Ξ (τ, q) and ... Ξ (τ, q), where q = {q 1 , 0, . . ., 0}.If the null hypothesis is rejected, it means that the full spectrum is necessary to recover all the relative moments. 6The third step of the procedure consists of a random walk (RW) hypothesis.Assuming the multiscaling parameter B = 0, we perform the regression of Equation 9 with only the constant A and test if A = 0.5 with a t-test.In case the null hypothesis is rejected, this means that the RW scaling is incorrect and the use of the squareroot of time rule severity creates a bias in the risk measures.The last step involves a confirmatory test of the results deriving from the first and second steps of the test procedure.In particular, we perform the full regression of Equation 9 and test for A = 0.5 and B = 0 using a t-test.If this test gives a conflicting result with respect to the first and second steps, we cannot assert anything precise on the process and a deeper analysis needs to be performed by controlling for different input specifications.

Simulation experiment
In the simulation experiment, we focus on one of the most used models to generate multifractal time series: the Multifractal Random Walk (MFRW), proposed by (Bacry, Delour, and Muzy, 2001b,a).This model is capable to generate multifractal time series with a known multiscaling spectrum.In addition, this model is able to generate time series with the empirically observed stylized facts.In the discrete version of the MFRW, the process r τ (t) = p(t + τ ) − p(t) is defined as (Bacry, Delour, and Muzy, 2001b): whith ∆t ∼ N (0, σ 2 ∆t), ω ∼ N (−λ 2 ln(L/∆t), λ 2 ln(L/∆t)) where λ is called intermittency parameter and determines the strength of the multifractality, L is the autocorrelation length, σ 2 is the variance of the process and ∆t is the discretization step.The distinctive feature of the MFRW is that, even if the ∆t (k) are independent, the ω ∆t (k) are not, having autocovariance In the continuous limit, the scaling exponents of this model are The power of this model is that it encompasses all the major stylized facts using only three parameters (λ, L, σ).In fact, this model is able to reproduce fat tails, volatility clustering and a multiscaling spectrum.For the purpose of simulation, we generated 100 paths each of dimension T = 10000 and we set the model parameter to L = 250,σ = 1 and λ = 0.05, 0.1, 0.3 in accordance to empirical findings (Bacry, Kozhemyak, and Muzy, 2008a,b;Løvsletten and Rypdal, 2012).As explained in previous Section, q max = 1.In particular, we use q ∈ [0.02, 1] with pace 0.02 which converts to 50 evaluated moments.Since the process is generated to be endogenously correlated and with volatility clustering, the anomaly procedure does not find any anomaly.To select τ max , we use the τ * estimated by the ACSR.Table 1 shows the result of the method for the different specifications of λ.As it is possible to observe, the procedure is quite accurate and the 95% confidence intervals (C.I.) always contain the value of L = 250 which is the truncation parameter.Once the parameters are estimated, we compute the multiscaling exponents and evaluated their statistical significance.Since Equation 14gives us the true multiscaling spectrum, we can easily test the performance of the GHE approach and compare it with the new proposed methodology of this paper based on the normalized and standardized structure function (NSSF) proposed in Equation 8 and the relative normalized and standardized structure function (RNSSF) proposed in Equation 11, which we name Normalized and Standardized Generalized Hurst Exponent (NSGHE) and the Relative Normalized and Standardized Generalized Hurst Exponent (RNSGHE), respectively.We will then use the latter methodology to test the multiscaling spectrum.Tables 2 and 3 present the root mean squared errors (RMSE) of the different methodologies computed over the 100 realizations for both A and B of Equation 4. 7 As it is possible to notice, the RNSGHE generally overperformes with respect to the other specifications.It is important to highlight that the elimination of the slope ambiguity has improved considerably the results.In fact, the standard GHE approach has the highest RMSE among all the specifications.This result was expected since the new methodology helps to remove uncertainty and thus the estimation performance.We finally report the results for the multiscaling test.
Given the superiority of the RNSGHE L and RNSGHE N L models, from this stage on we will use these models in the paper.Now, we show the nature of a process by performing the multiscaling test.Figure 2 shows the p-values of all the 50 coefficients related to the q moments Equations for a realization of the MFRW model with the same parameters specification as above and λ = 0.05.What we can observe is that if we choose a confidence level of GHE RNSGHE L NSGHE NL RNSGHE NL λ = 0.05 0.0112 0.0039 0.0047 0.0039 λ = 0.10 0.0098 0.0043 0.0046 0.0043 λ = 0.30 0.0167 0.0144 0.0137 0.0144 Table 3: RMSE for the parameter B for the different methodologies.Subscript L refers to the linear regression while N L to the non-linear regression.
5%, the null hypothesis of increments scaling increments equal to 0 for all the moments will be rejected.However, if we set a more stringent confidence level, for example 1%, the null is not rejected for some coefficients, resulting in a uniscaling process.
Figure 2: P-values of all the 50 coefficients related to q moments Equations for a multifractal random walk with T = 10000, λ = 0.05, L = 250 and σ = 1.
For λ = 0.1, 0.3 all the p-values are almost 0.This is not unexpected since we generated processes with a non negligible amount of multiscaling.Once the t-test has been carried out, we perform an F-test on the overall scaling spectrum (equation 11).Results are reported in Table 4.It is possible to infer that the process generated with λ = 0.05 does not reject the null for which all the scaling increments are equal to 0 while, for λ = 0.10, 0.30 the null is rejected and the full scaling spectrum is necessary to reconstruct the relative moments.By combining this result with the outcome of the t-test, we can state that the process generated with λ = 0.05 is weakly multiscaling at 5% confidence level but it is not multiscaling at 1% confidence level.The other two specifications are strongly multiscaling at any reasonable confidence level.
Since these processes have been generated such that they have specific multiscaling spec-RNSGHE L RNSGHE NL λ = 0.05 0.99999 0.99999 λ = 0.10 0.00003 0.00003 λ = 0.30 0.00000 0.00000 Table 4: P-values of the F-test for the null that only H(0, q 1 ) is different from 0. Subscript L refers to linear regression while N L to non-linear regression.
trum, the other two tests of the multiscaling test procedure have trivial results.

Empirical application
In this Section we perform the empirical application of our statistical methodology.

Data
The dataset used for the analyses is composed by stocks listed in the Dow Jones (DJ).
In particular, close prices of stocks are recorded on a daily basis from 03/05/1999 up to 20/11/2019, i.e. 5363 trading days.We use 27 over the 30 listed stocks as they are the ones for which the entire time series is available.For the purpose of our analysis, we use log-prices and log-returns.Table 5 reports the summary statistics of the stocks under analysis.
As it is possible to notice from Table 5, all the log-returns are centered at 0, the skewness is (in most of the cases) different from 0 while the Kurtosis clearly shows on of the stylized facts of most financial time series: fat tails.

Multiscaling test
In this Section, we report the results of the multiscaling test. 8We report the results of all the steps of the testing procedure explained at the end of Section 3.1.2.Results are presented using the RNSGHE L because as explained in the previous Section, it has the best performance in the correct estimation of the scaling spectrum.9Results are summarized in Table 6.The second column of the Table presents the τ * calculated using the ACSR methodology, which is presented in Section 3.1.1.For its estimation, we fix a maximum value for the choice of τ equal to T 5 = 1072 in order not to bias the scaling estimation with too few values.Hence, we notice that several stocks reach the boundary value, suggesting a very high rate of persistency in the time series.The third column of the Table reports the response to the weak multiscaling (weak M-S) process hypothesis, i.e. hypothesis that a single is not rejected, namely Cisco and Pfizer.However, this is a first order approximation of the process and do not check if the process is multiscaling.The sixth column summarizes results of the confirmation test, which is equivalent to test A = 0.5 and the statistical significance of B in the full regression model of Equation 9. Finally, the last three columns of the Table reports the estimated Hurst exponent H computed for the RW test, the linear scaling index A and the multiscaling proxy B. These results point out that multiscaling is a stylized fact and can be statistically tested by rewriting the structure function in a convenient way.In the next subsection, these results are used to compare the standard VaR calculation using square root of time rule and using the estimated scaling.

Effect of anomalies in the multiscaling estimation
Mltiscaling time series are generated from trading dynamics.One of the fundamental aspects of multiscaling systems is the strong endogeneity of the sample paths, aspect which is considered to be originated by financial trading dynamics.For this reason, transient exogenous shocks only distort the analysis and consequently, the estimation procedure.Hence, the statistical procedures used to analyse multiscaling systems are highly sensitive to exogenous shocks.In this context, we refer to an exogenous shock as an unexpected and transient behaviour of the stock price, not explainable by the market conditions or by the price path.In addition, it is possible to have anomalies in the time series due to errors or algorithmic trading crashes.The anomalies in financial time series can be grouped in 3 main categories: spikes, jumps and contamination errors.Figure 3 shows these possible anomalies.The top left panel is dedicated to the original log-price time series for Verizon.The time series is quite volatile and in fact, the log-returns has a Kurtosis index equal to 9.However, although the distribution of log-returns is fat tailed, there are not clear anomalies.The top right panel of Figure 3 depicts the same log-price time series to which a strong fat tailed series (Kurtosis larger than 1000) is added.This is the case of contamination error.This is generally due to machine errors in the data transmission process.The bottom left panel reflects the Verizon log-price time series with a random spike added.The spike can arise from multiple sources, among all algorithmic trading errors or contamination errors due to data manipulations.The last panel in the bottom right corner, represents the log-price series with an added jump.Jumps per se can arise from endogenous or exogenous shocks, but if they come from an endogenous driving force, they persist in the jump direction.If they come from an exogenous source instead, they tend to be transient.In a relatively recent paper (D.Sornette and Muzy., 2003), the authors explain that huge financial crashes can be originated from endogenous shocks which have a huge persistence behaviour.This kind of shocks are inherent in the price process so they are not transient anomalies.
In a mostly technical paper, (Katsev and L'Heureux, 2003) showed both theoretically and experimentally, that such data anomalies can strongly bias the results, especially for short datasets.In particular, the paper shows that under certain circumstances, these irregularities can generate spurious scaling.For these reason, it is suggested to analyse the time series and eliminate such anomalies before proceeding with the scaling estimation.In order to do so, we propose a methodology based on financial stylized facts.More precisely, we use volatility clustering and long-range dependence of asset returns (Cont, 2001;Chakraborti et al., 2011).In this empirical context, the quantities that we name cumulative variance should be very similar, except when an exogenous (unexpected) anomaly exists.The volatility clustering drives the similarity in the short period since |r i+1 | and |r i | are expected to be very similar (same cluster), while the long-range dependence drives the similarity of the two measures over the long-run.Figure 4 represents the two quantities for the Verizon stock.These two quantities are approximately equal, confirming that even with high volatility and many tail events, the time series does not contain exogenous shocks.
The difference between the two cumulative series is given by D(t) = CV (t) − CAV (t).By running a change-point detection in the mean level of D(t), it is possible to detect the anomalies in the price series and replace the corresponding values on the original series according to a specific rule, e.g. the mean of previous and subsequent datapoints.The top Given the fact that such events (spikes or jumps) are rare and have unconventional magnitude, their removal can only benefit the analysis.Let us now estimate the multiscaling exponent for the Verizon stock when the anomalies reported in Figure 3 (bottom panels) are not removed by the time series.Table 7 reports the results.As it can be noticed, the estimated values changes considerably, especially in the scenario in which a spike is added.

H
A B Spike 0.5290 0.6013 -0.1417 Jump 0.4942 0.5220 -0.0543Table 7: Results of the multiscaling estimation on the times series reported in Figure 3 (bottom panels) with anomalies not removed.
For completeness, we also performed a t-test with the null hypothesis of no difference between the estimates with anomalies and the one reported in Table 6.The null hypothesis for all the coefficient is strongly rejected at any confidence level.
To show how these anomalies can generate spurious multiscaling, we generate 100 fractional Brownian motions (uniscaling process) of length 1000 with Hurst exponent H = 0.47 (the one estimated for Verizon).At these simulated time series we add a spike and a jump and estimate H, A and B for both the series with and without the anomalies.The results are reported in Table 8.As we can see, when the anomalies are not present in the time series, the average values for H, A and B are in line with the true values and not statistically different from them.In the case of a spike, we can see that H, A and B are severely biased.In particular, the spikes make the times series looking multiscaling, while we know it is not.Finally, for the case of the added exogenous jump, we have that the scaling exponent curves and the parameter B is not equal to 0. Also the other two estimated coefficients are upward biased and statistically different from the theoretical ones.These results clearly show that the scaling exponents are sensitive to such anomalies and estimates can be biased if such anomalies are not taken with care.Therefore, it is fundamental to analyse the data before performing the multiscaling analysis.

Practical application of scaling and multiscaling to VaR
In this section we show that a simple VaR configuration without a scaling or multiscaling consideration might bias the VaR estimation at higher aggregation scales.We use daily stocks data from the Dow Jones index to estimate the multiscaling spectrum and to perform the multiscaling test described in Section 5.After the procedure is carried out, we estimate the Historical and Gaussian VaR at 1 day and we successively use it to compute the yearly VaR using the square root of time rule.We then compare it with the fractional VaR with proper scaling and highlight eventual biases.To conclude, we propose a methodology to compute a multifractal VaR.

Value at Risk
VaR is an easy and intuitive way to quantify risk for assets and portfolios.Let V aR(τ, 1 − α) be the Value at Risk at frequency τ for a confidence level equal to 1 − α which satisfies where r τ (t) are the log-returns at frequency τ .Several methodologies are used to compute the VaR.Among all, we recall the Historical VaR (HVaR) and the Gaussian VaR(GVaR).The former is a non-parametric approach which uses historical data to compute the VaR, while the latter assumes a Gaussian distribution of stock returns and applies the Gaussian formula for the percentiles computation to extract the VaR at a given confidence level.The issue faced in applied finance is that the square root of time rule works only under the assumption of iid Gaussian returns.However, this technique is widely adopted regardless of its assumptions.
In our analysis, VaR is computed using the two aforementioned approaches at τ = 1 day at 95% confidence level (V ar(1, 95%)).Annual VaR (V ar(250, 95%)) is calculated with the scaling exponent equal to 0.5, i.e.V ar(250, 95%) = V ar(1, 95%) × 250 0.5 and to the estimated H, i.e.V ar H (250, 95%) = V ar(1, 95%) × 250 H .12 We further estimate the true V ar(250, 95%) using annual returns (τ = 250 days) and compare them.Results are shown in Figures 7 and 8.These Figures show that when we compare the VaR with H scaling time rule and the VaR with square root of time rule, the deviation from the true VaR is lower when the former approach is used.In fact, the bias (with respect to the true VaR) of the VaR with H scaling time is considerably lower for most of the stocks in both the HVaR and GVaR settings.This is due to the fact that over the long-run, even a small divergence from the assumption of scaling exponent equal to 0.5 can make a substantial impact.
To conclude the analysis, we also report the relative error, i.e.where K is equal to 0.5 for the VaR computed with the square root of time rule or equal to the estimated H, H, reported in Table 6.This helps us to identify the magnitude of the bias from the true VaR and to compare the two scaling approaches.Figure 9 shows the results.As the Figure shows, using the correct scaling results in a smaller relative error.This confirms the fact the choice of a proper scaling exponent should not be neglected by the financial community considering that its estimation and testing is relatively simple.
Figure 9: Relative error between the true VaR calculated using annual data, HVaR and GVaR computed using daily data scaled by the factor 0.5 and by the estimated factor H.

Multiscaling consistent VaR
In the previous subSection, we showed that using the correct scaling contributes to reduce the computation error for VaR at smaller frequencies.However, as explained in Section 5.2, all the time series analysed present a strong multiscaling.To deal with such situations, we report here a possible solution.While VaR is related to the log-returns, multiscaling is a property of the moments of the log-returns.For this reason, there is not a straightforward formula to compute VaR which takes into account multiscaling.An exception to this is the Multifractal VaR proposed in (Lee, Song, and Chang, 2016), where the author introduces a VaR consistent with the multifractality of financial time series using the Multifractal Model of Asset Returns (MMAR).In a previous paper, (Batten, Kinateder, and Wagner, 2014) performs a similar analysis but relies on the MMAR Monte Carlo simulations and then, computes VaR on the simulated time series.The Monte Carlo approach has the advantage of letting the researcher use the model that best depicts the data.In fact, one can calibrate the MRW or the MMAR and generate a large number of sample paths which can be used to compute VaR.In case of moderate multiscaling, the difference can be low but for multiscaling processes with a |B| > 0.05, neglecting such feature can strongly distort the VaR.In this work, we use the MFRW to simulate 250 trading days (i.e. 1 year) of log-returns and compute the VaR of the simulated paths.For this purpose, three parameters need to be estimated: the variance σ 2 , the autocorrelation scale parameter L and the intermittency parameter λ.The variance can be estimated from the log-returns time series as σ 2 = V ar(r τ (t)), with τ = 1, the parameter L is set to be equal to τ * , while the intermittency parameter λ can be extracted by equating the estimated coefficients of Equation 4 to the parameters in Equation 14getting two (possibly different) estimates of λ, i.e. λ A and λ B . 13 For each stock we estimate the three parameters, σ 2 , L and λ, and we generate 100000 independent paths of daily returns for a year (i.e.250 days).Finally, we compute the VaR which we name Multiscaling VaR (MSVaR) by using a 95% confidence level on these simulations.Results are depicted in Figure 10.It is possible to appreciate that the MSVaR computed on the simulated paths has comparable size to the Historical VaR computed on annual data.It is also important to notice that the values predicted using λ A and λ B are very similar, suggesting that the stocks log-returns can be adequately approximated by the MFRW model.Nevertheless, we remark the importance of the full multiscaling estimation and testing procedure which leads to the MSVaR.In fact, if the previous analysis is bypassed the estimated risk metrics can be severely biased and inconsistent.

Conclusions
In this paper we proposed a step by step procedure to robustly estimate and test multiscaling in financial time series.By rewriting the structure function in a convenient way we were able to make multiple tests on the scaling spectrum and asses the statistical significance of multiscaling, discriminating between weak and strong multiscaling.We have shown the effect of anomalies in financial time series and studied the impact on the estimated scaling exponents.Moreover, we have showed how the use of proper scaling can help reducing the error in the forecasting of VaR at smaller frequency with respect to the commonly used square root of time rule.Finally, we have proposed a Multiscaling consistent VaR using a Monte Carlo two estimated coefficients, i.e. λ A and λ B , are very similar.
MRW simulation calibrated to the data and on which the VaR is then computed.Results are encouraging and confirms the goodness of such methodological approach.Multiscaling is a stylized fact which can make the difference in the assessment of risk measures and in building quantitative models.It can be easily extrapolated from data and should not be overlooked by risk managers and authorities.

Figure 1 :
Figure 1: ACSR methodology computed on the autocorrelation function of the absolute values of the log-returns time series.The red circle shows the breakpoint where the regression line has a break.

Figure 3 :CV
Figure 3: Verizon log-prices time series and some anomalies of financial time series.

Figure 4 :
Figure 4: Cumulative variance (CV) and cumulative auto-covariance (CAV) for the Verizon time series.

Figure 5 :
Figure 5: Cumulative variance (CV (t)) and cumulative auto-covariance (CAV (t)) for the Verizon time series (top panel) and the measure D(t) (bottom panel) in case of the added spike.
Figure 7: HVaR using annual data and using daily data rescaled by the factor 250 0.5 and by 250 H .The stocks are sorted in order of the magnitude of the true VaR.

Figure 8 :
Figure 8: GVaR using annual data and using daily data rescaled by the factor 250 0.5 and by 250 H .The stocks are sorted in order of the magnitude of the true VaR.

Figure 10 :
Figure 10: Annual Historical VaR (HVaR) and Multiscaling VaR (MSVaR).Confidence level 95%.Stocks are sorted according to the magnitude of the Historical VaR.

Table 1 :
Results of the ACSR for the estimation of τ max .95% C.I. computed over 200000 bootstrapped samples.

Table 2 :
RMSE for the parameter A for the different methodologies.Subscript L refers to the linear regression while N L to the non-linear regression.

Table 6 :
Results of the multiscaling estimation and testing procedure.