Fluctuation analysis of electric power loads in Europe: Correlation multifractality vs. Distribution function multifractality
aa r X i v : . [ q -f i n . S T ] J un Fluctuation analysis of electric power loads in Europe: Correlation multifractality vs.Distribution function multifractality
Hynek Lavička
1, 2, 3, 4, ∗ and Jiří Kracík † Alten Belgium N.V., Chaussée de Charleroi 112, B-1060 Bruxelles, Belgium DYMO BVBA, Industriepark-Noord 30, 9100 Sint-Niklaas, Belgium Department of Institutional, Environmental and Experimental Economics,University of Economics in Prague, W. Churchilla 4, CZ-130 67 Praha 3, Czech Republic Institute for Theoretical Physics, Celestijnenlaan 200D, KU Leuven, B-3001 Leuven, Belgium Charles University in Prague, Faculty of Social Sciences,Institute of Economic Studies, Opletalova 26. CZ-11000 Prague 1, Czech Republic
We analyze the time series of the power loads of the 35 separated countries publicly sharinghourly data through ENTSO-E platform for more than 5 years. We apply the Multifractal De-trended Fluctuation Analysis for the demonstration of the multifractal nature, autocorrelation andthe distribution function fundamentals. Additionally, we improved the basic method described byKanterhardt, et al using uniform shuffling and surrogate the datasets to prove the robustness of theresults with respect to the non-linear effects of the processes. All the datasets exhibit multifrac-tality in the distribution function as well as in the autocorrelation function. The basic differencesbetween individual states are manifested in the width of the multifractal spectra and in the locationof the maximum. We present the hypothesis about the production portfolio and the export/importdependences.
PACS numbers: 02.50.Fz, 05.10.Gg, 05.40.Fb, 05.45.Tp, 88.05.-b
I. INTRODUCTION
Recently unprecedented processes influenced the globaleconomy; “globalization” as an integration of the globalproduction system on one side and a slow but stable ban-ishment of the national overlook of the markets - “liberal-ization” on the other side. The energy sector despite of itscrucial importance and appropriate level of the govern-ment control was also influenced. In Europe, two jointmarkets, Nordpool and CE, have been established andthe production of the electric energy is not necessarilyproduced and consumed within one country but local dis-tributors are free to buy electric energy from abroad.On the other hand, from the ecological point of view,the Kyoto protocol and the Doha amendment prescribethe governments to focus on renewable sources to assuretransition to sustainable development. A general fea-ture of the renewable sources in contrast to the classicalsources is low yield as well as the fact that installationsare placed far from the centers of consumption that cre-ates problems of additional currents on the internationalscale, see e.g., [1]. We also note that each country hasspecific portfolio of the power plants and they operate tofollow needs of consumers.Electric grid is operated on national level by the Trans-mission System Operator (TSO) that controls a swift andeffective flow of energy from the powerplants to the con-sumers. The process is based on the large quantities ofdata measured in real-time and instantaneous action is ∗ hynek.lavicka@fjfi.cvut.cz; Thanks to I. Jex, P. Exner, J. Tolarfor support of the work. † [email protected] taken if necessary to balance the electric network to beproportional to consumption. Naturally, the process isinfluenced by deterministic trends as well as random ef-fects, see [2, 3].In the framework of physics there was developed theLangevin equation describing a path of a particle in arandom environment d X ( t ) = µ ( t, X ( t )) · d t + d W ( t, X ( t )) , (1)where µ ( t, X ( t )) is the deterministic drift that standsfor deterministic effects in electric power system, while dW ( t, X ( t )) is a term representing noise depending onactual state X ( t ) that reflects random effects within theelectricity power grid. Analysis of Eq. 1 in the frameworkof Brownian motion has been done in Ref. [4]. Usual as-sumption on the properties of the noise is the Gaussiandistribution and memorylessness h X ( t ) · X ( t + ∆ t ) i ∼ δ (∆ t ) or short range memory h X ( t ) · X ( t + ∆ t ) i ∼ exp (cid:16) − ∆ tL c (cid:17) , where L C is a correlation length.However, in the real systems the properties may notbe always satisfied, e.g., the probability distribution is anon-Gaussian one. Even more, it may be fat-tailed and(or) the autocorrelation function posses long-range mem-ory h X ( t ) · X ( t + ∆ t ) i ∼ ∆ t − γ , where γ is the exponentof the power law, see [5]. A typical field of the studywhere the theory was used represents a stock exchange,see [6], where large datasets have been stored and pre-pared for analysis.In recent decades computers were widely used for thetime-series analysis and among the methods that areused belong the Detrended Fluctuation Analysis (DFA)and its derivation the Multifractal Detrended Fluctua-tion Analysis (MFDFA), see [7–11]. Historically it iserived from the R/S method used in the field of hy-drology, see Ref. [12]. Recently the time series analysishas focused on the problem of the trends and its impacton the analysis, see Refs. [13–15]. The literature alsocontain a non-orthodox modification of the MFDFA, seeRef. [16, 17] as well as modification usable for analyz-ing the cross-correlations among the time series, see Ref.[15, 18]. We note that the literature also contains modi-fication that goes beyond one dimension, see Ref. [19].Versatility of the method has been proved in a num-ber of applications. Originally it was used to analyzethe datasets in biophysics, see Refs. [9, 20, 21] but re-cently the method is employed in the number of the stud-ies spanning from surface roughness via human gain andeconophysics to earthquakes and clouds, see [6, 22–30].Generally speaking, the method demonstrated extractseffectively the information about the scaling propertieswithin the dataset and then the properties of the auto-correlation function.Focusing on electroenergetics there is a number of thestudies focusing either on electric power loads and pricesof electricity, see Ref. [2, 31–34]. Recently a competitionof the teams that predict price or electric power loads hasgrown around the globe and they use various methods,see Refs. [35–37]. General weakness of number of stud-ies is natural presence of oscillations that spoils correctestimation of the Hurst exponent. The gap was filled in[38].In this study we would like to examine the datasets ofthe electric power loads in the countries of Europe. Weintend to perform the MFDFA and we also plan to exploitvarious enhancements of the method to reveal propertiesof both the distribution and autocorrelation functions.Moreover we test how the properties of the functionsare condensed among parameters – we intend to ques-tion multifractality of the datasets. Finally we comparefunctionality of the states by the multifractal properties.The structure of the paper is as follows. Firstly, wedescribe the method used in the study. Then we exploitpower of computers and show the results of the analysis.Finally, we come with the conclusions focusing on thestochastic properties of the datasets and then we comparepower loads of the European countries. II. METHODS
To analyze a timeserie { X ( i ) } Ni =1 with N datapoints,usually as a result of a measurement of a physical vari-able, we perform modification of the Multifractal De-trended Fluctuation Analysis (MFDFA). If the timeserieis non-stationary where the driving forces of the systemare stronger than the fluctuations within them we pro-pose a method that is based on the following 10 steps: a. Fourier transform: Firstly we preform theFourier transform to detect presence of the carrier sig-nal in the sample of the timeserie L o a d [ G W ] Evolution of electric power load
Jan 1 2010 Jan 1 2011 Jan 1 2012 Jan 1 2013 Oct 31 2013 − 10− 50510 L o a d [ G W ] Figure 1. The electric power load of three countries of Europe(upper figure) and the signals modulated X car ( i ) on the car-rying signal X car ( i ) after performing Step 1 (bottom figure)spanning period since Jan 1 2010 to Oct 31 2013. Red, greenand blue curves are Germany, Italy and Norway respectively.We show the typical examples of the the electric power load ofcountries with respect to seasonality and inter-day and inter-week changes. Germany is an example of the strong season-ality with the strong inter-day and inter-week changes. Italyshows the strong inter-day and inter-week changes with thelimited seasonality. Finally, Norway exhibits the strong sea-sonality with the limited inter-day and inter-week changes. b X ( ω ) = 1 √ N N X j =1 X ( j ) exp ( − π · i · j · ω ) . (2)Generalizations of the Fourier transforms are described in[39–46] and they can be used to strengthen abilities of themethod to the wavelets of various types. We note thatcorresponding inversion formula must be used below.For the colored noise (fluctuations) the Fourier trans-form follows ˆ X ( ω ) ∼ | ω | β with a number of the addi-tional peaks caused by the carrier signal upon the fluctu-ations are modulated. We decompose the signal into twocomponents b X ( ω ) = b X fluc ( ω ) + b X car ( ω ) , (3)where b X fluc ( ω ) ∼ | ω | β is the Fourier transform of themodulated signal and the Fourier transform of the carriersignal is b X car ( ω ) and it defines the profile average X car ( j ) = 1 √ N N X ω =1 b X car ( ω ) · exp (2 π · i · j · ω ) . (4)2 m10 Q A C ( m ) Test statistics
GermanyItalyNorwaySwitzerlandCyprus(χ ) −1 (0.99) Figure 2. The test statistics for the selected countries of Eu-rope. Its comparison with χ m -distribution at the signifi-cance level with the m degrees of freedom to the test hypoth-esis of none autocorrelations (region (cid:2) , χ m (0 . (cid:3) ). Quantile of Gauss distribution Q u a n t i l e o f E x p e r i m e n t a l d i s t r i b u t i o n Comparison of quantiles
Ideal fitNorth IrelandEstoniaPolandAustria SpainFranceIcelandGermanySweden
Figure 3. Comparison of quantiles of stochastic part of thetime series for selected countries of the Europe. b. Test of presence of autocorrelations
The first stepof the procedure is to test of the presence of autocorre-lation within the sample. We define the test statisticsmade of the modulated signal X fluc ( i ) as Q AC ( m ) = N m X i =1 C i N − i , (5)where the sample autocorrelation is defined by formula C i = P Nk = i +1 X fluc ( k ) · X fluc ( k − i ) P Nk =1 X fluc ( k ) . (6)The test statistics Q AC ( m ) follows χ m with m degrees of freedom if none autocorrelations are present in thesample. We fix the level of significance α and varying m .If Q AC ( m ) < (cid:0) χ m (cid:1) − ( α ) , where (cid:0) χ m (cid:1) − ( α ) stands for α -th quantile of χ m distribution [47], at the significancelevel α we conclude independence of x k s for the segmentof length m . In the opposite case Q AC ( m ) > (cid:0) χ m (cid:1) − ( α ) the hypothesis of none autocorrelations is rejected.We admit that the step is optional and it serves as anindicator of the presence of autocorrelations for the Gaus-sian distribution. The indicator for the non-Gaussiandistribution would need to derive appropriate χ m distri-bution. c. Construction of profile Construction of a profile X prof ( i ) of the timeserie from the modulated signal X fluc ( i ) is performed by formula X prof ( i ) = i X j =1 X fluc ( i ) . (7)The step is auxiliary and it helps to perform the anal-ysis accurately even for the anti-persistent processes.Please, see Ref. [7] for the discussion. d. Separation into segments We separate the profileof the timeserie X prof ( i ) into the N s windows X seg,w ( i ) covering the whole dataset with the length s where eachwindow is denoted by a number w . The minimal numberof the windows is ⌊ Ns ⌋ , where ⌊ x ⌋ is the largest integersmaller than or equal to a number x . Thus the number ofthe timesteps N is not generally multiple of s . In orderto obtain better statistics we use as small overlap of theconsecutive windows where possible. e. Construction of trend We establish the localpolynomial trend X prof,w ( i ) within each window w ofthe size s using the least square fit of the dataset. Us-ing the trend we detrend the data and we calculate thesample variances in the window F ( s, w ) = 1 s P si =1 (cid:0) X seg,w ( i ) − X prof,w ( i ) (cid:1) , (8)where the calculation is performed for the windows w =1 , . . . , N s . f. Calculation of fluctuation function The fluctua-tion function of the q -th order is defined by the formula: F q ( s ) = q = 0 N s P N s w =1 (cid:16) F ( s, w ) q (cid:17) q q = 0 exp (cid:16) N s P N s w =1 ln (cid:0) F ( s, w ) (cid:1)(cid:17) . (9) g. Calculation of generalized Hurst exponent Thefluctuation function behaves as F q ( s ) ∼ s h ( q )+1 , (10)where +1 correction stands for the correction to the in-tegrated time series in Step . Generally, for positive3 , h ( q ) describes scaling behavior of large fluctuationswithin a segment while for negative q , h ( q ) describes scal-ing behavior of small fluctuations. Independence h ( q ) of q means monofractal behavior of the dataset. If we as-sume that the fluctuations in the dataset are stable andfollow the Gaussian distribution then the Hurst exponent H = h ( q = 2) is then < H < for the Gaussian distri-bution of the noise. For H = , long range autocorrela-tion is not present and the dataset may be independent. H > means long-range autocorrelations (persistence)and H < stands for long range anti-autocorrelations(anti-persistence).The non-Gaussian distributions with the power-lawtail allows the generalized Hurst exponent to exceed therange of (0 , , see [19] and h ( q ) is a nontrivial functionof the variable q . It generally exhibits combination ofthe effects caused by the probability distribution and theautocorrelation function. In order to extract informationon both we execute next two additional steps. h. Shuffling the dataset We shuffle the originaltimeserie { x ( i ) } Ni =1 to form (cid:8) x shuf ( i ) (cid:9) Ni =1 using theFisher–Yates algorithm, see [48] for the description andthe historical note. By calculating the fluctuation func-tion of the shuffled timeserie F shufq ( s ) performing stepsthrough to we obtain the shuffled fluctuation function F shufq ( s ) = F q, { x shuf ( i ) } ( s ) , (11)where the averaging is executed for the different real-izations of shuffling. Shuffling of the original timeseriedestroys (if present) autocorrelations and thus the shuf-fled fluctuation function carries information about thedistribution function. By analogy of Eq. 10, the shuffledfluctuation function behaves F shufq ( s ) ∼ s h shuf ( q )+1 , (12)where h shuf ( q ) called shuffled Hurst exponent describesthe scaling behavior of the fluctuation function of theshuffled timeserie. i. Correlation Hurst exponent We define the auto-correlation Hurst exponent h cor ( q ) as follows: h cor ( q ) = h ( q ) − h shuf ( q ) . (13)The autocorrelation Hurst exponent contains infor-mation on the autocorrelation function of the timeserie { x ( i ) } Ni =1 . For H cor > the timeserie is long rangeautocorrelated while for H cor < there is long rangeanti-autocorrelated behavior. In case H cor = 0 the time-serie is either non-autocorrelated or short range autocor-related. j. Surrogate dataset We perform additional test us-ing the surrogate dataset. We intend to verify robustnessof our conclusions with respect to the non-linear effectsand the non-Gaussian features. The surrogate dataset conserves the autocorrelation function. Calculation ofthe surrogate dataset { x sur ( i ) } Ni =1 is obtained using theinverse Fourier transform of b X sur ( ω ) = b X fluc ( ω ) · exp ( i · Q ( ω )) , (14)where Q ( ω ) is a IID random variable with uniform distri-bution within the range [ − π, π ] . After then the MFDFAis used on { x sur ( i ) } Ni =1 using the above steps from 4 to8. We perform samples of the surrogate dataset andthe surrogate fluctuation function is defined by formula F surq ( s ) = F q, { x suri } ( s ) . (15)Analogically, the surrogate fluctuation function follows F surq ( s ) ∼ s h sur ( q )+1 , (16)where h sur ( q ) is the surrogate Hurst exponent. k. Distribution Hurst exponent The distributionHurst exponent is defined by h dist ( q ) = h ( q ) − h sur ( q ) . (17)The distribution Hurst exponent contains informationabout the distribution function and its scaling properties. l. Multifractal spectrum The fundamental proper-ties of the autocorrelation function and the distributionfunction of the time series can be studied using multifrac-tal spectrum f ( α ) that is Legendre transform of scalingfunction τ ( q ) = q · h ( q ) − f ( π ) ≡ q ( π ) · π − τ ( q ( π )) (18) π ≡ dτdq = q · h ′ ( q ) + h ( q ) . (19)Since we deal with three generalized Hurst exponents( h ( q ) , h cor ( q ) , h shuf ( q ) ) we calculate f ( π ) , f cor ( π ) and f shuf ( π ) . In case that the autocorrelation func-tion has monofractal property or the distribution func-tion has single exponent then π = const. For multifractalcase there occurs a distribution of α values for both theautocorrelations as well as distributions, see [49]. Thewidth of the spectrum △ π = π max − π min describes thestrength of multifractality, where π max = max q π ( q ) and π min = min q π ( q ) . For broad spectrum △ π is an indica-tor of strong multifractality while for narrow spectrum itis the indicator of weak multifractality of the time series.Analogically the autocorrelation and distribution mul-tifractal spectrums are defined by equations 18 and 19with h cor ( q ) and h dist ( q ) respectively. They produce f cor ( π ) and f dist ( π ) and the widths of multifractal spec-trums are ∆ π cor = π cormax − π cormin and ∆ π dist = π distmax − π distmin .4 s [h] F q ( s ) Fluctuation function for France q =−10q =−5q =−2q =0 q =2q =5q =10s
Figure 4. The fluctuation function for the MFDFA-4 forFrance for q ∈ {− , − , − , , , } . Fluctuation functionsare shifted by power of between each consecutive q . Thetypical power law behavior is presented. III. RESULTS OF ANALYSIS
We perform the methodology which was described inthe previous section on a dataset of the electric powerload for 35 European countries or the independent elec-tric power systems which were obtained from ENTSO-E[50] . This organization roofs local Transmission Sys-tem Operators (TSO) that control the swift and effectivefunctionality of the backbone electric power grid in theEuropean countries. The dataset spans from January 12008 to December 31st 2012 with the one-hour resolutionfor the longest dataset.Each measured datapoint collects produced electric-ity within a region as it counts for exports, imports andadding balance of pumped-storage hydroplants [51]. Itamounts electric power that is provided to users – citi-zens, industry and services – through local electric dis-tribution companies. We also note that actual demandthat is intended to be satisfied must be lower in longerperiods of time than actual supply, otherwise blackoutcan happen on small or on the larger scale. To preventsignificant economic and physical waste of electric powercompanies producing and transmitting electric power aretrying to correlate both variables to the extent possi-ble. We may also perceive the dataset as a measure ofan activity (economic, social) of a society. Autocorrela-tions of the datasets can be understood as an indicatorof autocorrelations in European countries that share datathrough ENTSO-E. h ( q ) Hurst exponents for Sweden −0.020.000.020.040.060.08 h d i s t ( q ) MFDFA−2MFDFA−3MFDFA−4 MFDFA−5MFDFA−6MFDFA−7 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 q h c o r ( q ) Figure 5. The generalized Hurst exponents h ( q ) , h shuf ( q ) and h cor ( q ) are shown from top to bottom respectively. Theorders of the MFDFA between to are plotted with variouscolors for United Kingdom (excluding Northern Ireland). A. Oscillations in power grid data
Firstly, we show examples of datasets of typical Euro-pean countries at the top of the Fig. 1. Each datasetexhibits oscillations of period a year, a week and a daybut the countries differ in proportions of the oscillationson the time scale. Countries of the south Europe usu-ally exhibit limited one-year oscillations in contrast tothe countries of the northern Europe. Since most of thecountries are of Christian origin, the religious holidaysare usually public holidays, e.g., Christmas and Easter.There is a significant decrease of the electric power loadduring that period. Oscillations can generally spoil theresults of MF-DFA and biasing statistical analysis, see[13, 14]. Performing Step 1 of the method we obtain adataset (signal) that is modulated on the carrying signal.We note that performing cuts during Step 1 have beentested for different exponents β with minimal qualitative5 .0 1.1 1.2 1.3 1.4 1.50.00.20.40.60.81.0 f ( π ) Multifractal spectrum for Cyprus f s h u f ( π ) π f c o r ( π ) MFDFA−2MFDFA−3MFDFA−4 MFDFA−5MFDFA−6MFDFA−7
Figure 6. The multifractal spectrum f ( π ) , the shuffled mul-tifractal spectrum f shuf ( π ) and the multifractal autocorre-lation spectrum f cor ( π ) for Cyprus of the MFDFA ordersbetween 2 and 7 is shown. The typical numerical issues canbe observed at the bottom as the wings. impact on results of following analysis, see [38].We show examples of the modulated signal for Ger-many, Italy and Norway at the bottom of the Fig. 1.The modulated signal (cid:8) X fluc ( i ) (cid:9) Ni =1 is then analyzedand tested for the presence of autocorrelations or typesof distribution. We performed Komogorov-Smirnov testof the signal and rejected null hypothesis that there is theGaussian distribution, e.g. for Norway the p-value standsat extremely small value ∼ − and for Germany at ∼ − so that it fails to reject the null hypothesis. Todemonstrate deviations from the Gaussian distribution,see Fig. 3. The power loads of the countries deviatefrom diagonal line and the deviations show the presenceof the asymmetric probability distribution. Due to thefact that Austria does not fit the pattern of the rest ofthe states we suppose that it is due to diverge portfolioof the power plants. A u s t r i a B e l g i u m B o s n i a − H e r ce go v i n a B u l ga r i a C r oa t i a C y p r u s C zec h R e pub li c D e n m a r k E s t o n i a F i n l a nd F r a n ce G e r m a n y G r eece H un ga r y I ce l a nd I r e l a nd I t a l y L a t v i a L i t hu a n i a L u x e m bu r g M a ce d o n i a M o n t e n e g r o N e d e r l a nd s N o r t h I r e l a nd N o r w a y P o l a nd P o r t u ga l R o m a n i a S e r b i a S l o v a k i a S l o v e n i a Sp a i nS w e d e nS w i t ze r l a nd U n i t e d K i n g d o m State or national authority0.40.50.60.70.8 π c o r m a x Comparison of maximus of autocorrelation multifractal spectrum
MFDFA−2MFDFA−3MFDFA−4MFDFA−5MFDFA−6MFDFA−7
Figure 7. Plot of maxims of the autocorrelation multifractalspectrum for states or national authorities in Europe shownfor various orders of the MFDFA from 2 to 7.
Secondly we perform an autocorrelation test of adataset, calculating Q AC ( m ) and comparing with χ m distribution at level of significance. In the Fig. 2,the test rejects hypothesis of no-autocorrelation withina sample dataset for each European country. It moti-vates taking the following steps of the methodology thataims at the calculation of the MF-DFA to reveal proper-ties of the autocorrelation and distribution function. Viaperforming steps from 4 to 8 and 10 we obtain F q ( s ) , F shufq ( s ) and F surq ( s ) that follows the power laws 10, 12and 16, which is for France demonstrated in the Fig. 4. B. MFDFA
The calculation of the appropriate Hurst exponents h ( q ) , h shuf ( q ) and h sur ( q ) is then a key point for calcu-lation h cor ( q ) and h dist ( q ) (steps 9 and 11) that containinformation on the scaling properties of the autocorrela-tion and distribution function. For illustration, we show h ( q ) , h dist ( q ) and h cor ( q ) in the Fig. 5 for Sweden. Wenote that the generalized Hurst exponent h ( q ) is esti-mated above standard range for the Gaussian distribu-tion [0 , [52]. Next h dist ( q ) is close to indicating thatthe non-linear effects does not spoil Hurst exponent es-timation. Finally, the autocorrelation Hurst exponent h cor ( q ) indicates a presence of multifractality for auto-correlation function.Next, we focus on multifractal spectrums defined byEqs. 18 and 19, the multifractal spectrum f ( π ) , theshuffled multifractal spectrum f shuf ( π ) and the auto-correlation multifractal spectrum f cor ( π ) that are shownin Fig. 6 for Cyprus. All spectrums exhibit a peak-like structure with the most frequent value π max , π shufmax and π cormax respectively, forming the top of the appropri-ate multifractal spectrum. In table I we show π shufmax for6 igure 8. The width of the multifractal shuffled (on the top)and correlation (at the bottom) spectrum calculated by theMFDFA-4 and it is visualized on the map of Europe. the states of Europe and it is mostly close to , indi-cating the most frequent presence of the Gaussian dis-tribution. However, maximums of the autocorrelationmultifractal spectrum π cormax varyes significantly from onecountry to another, see the Fig. 7. π cormax is also out ofstandard range (cid:2) − , (cid:3) for the Brownian processes ex-cept for North Ireland, Luxembourg and Latvia [53] thatis on the edge. Most frequently the process is governed by1 and the stochastic process is persistent. However, theother states are governed by the most frequent processdescribed by equation d ˙ X ( t ) = µ (cid:16) t, ˙ X ( t ) (cid:17) · d t + d W (cid:16) t, ˙ X ( t ) (cid:17) . (20)The maximum of the autocorrelation multifractal spec-trum is then deceased by 1 (Similarly when we skip toexecute step 2 but execute steps 7,8 and 10 as they are).The process is then anti-persistent, similar analysis inRef. [33] also revealed anti-persistent processes within shuf ∆ π c o r Comparison of widths of multifractal spectrums
AustriaBelgiumBosnia−HercegovinaBulgariaCroatiaCyprusCzech RepublicDenmarkEstoniaFinlandFranceGermanyGreeceHungaryIcelandIrelandItalyLatvia LithuaniaLuxemburgMacedoniaMontenegroNederlandsNorth IrelandNorwayPolandPortugalRomaniaSerbiaSlovakiaSloveniaSpainSwedenSwitzerlandUnited Kingdom
Figure 9. Comparison of the widths of the correlated andshuffled multifractal spectrums. The order of the MFDFA isindicated by different symbol. Square, circle and triable downare the orders 4,5 and 6 respectively. the power grids. While the nd order differential equa-tions are more stable with respect to the observation ofdisplacement of particles at the start and at the end thanthe st order equations and the weaker anti-persistencealso decreases deviations. We argue that higher π cormax ismore beneficial for stability of the electric power grids. C. Multifractal spectra
We also report that the wide multifractal spectrumsare present within the dataset, see Tab. I. Since thewidth of the autocorrelation multifractal spectrum △ π cor is wider, we conclude that multifractality of autocor-relation function is also stronger than the distributionfunction. Since the width of the distribution multifractalspectrum △ π dist is far smaller than △ π cor and △ π and π distmax is close to zero, we conclude that the non-lineareffects are limited. Nevertheless they are present andmultifractality is not a side effect but it has a systematicnature in the operation of the electric power system. Weadmit that the width of the distribution and correlationspectrums is usually influenced by the numerical issues,see the bottom of the Fig. 6.
1. Operational properties of power grids in Europe
Focusing on the operational properties of the nationalpower grids, we demonstrate the differences in the Fig.8 by metric of the width of the shuffled and autocor-relation multifractal spectrums. North-South or East-West spatial clustering based on political or climatic isnot present for both spectrums. The widest widths of7 π shuf ∆ π cor ∆ π ∆ π sur ∆ π dist ItalyNetherlandsLithuaniaDenmarkPortugalHungaryMontenegroCyprusCzech RepublicUnited KingdomNorth IrelandGermanySlovakiaFranceIrelandNorwayBelgiumAustriaFinlandRomaniaGreeceBulgariaSwedenSloveniaBosnia − HercegovinaCroatiaSwitzerlandSerbiaMacedoniaLatviaLuxemburgPolandEstoniaSpainIceland − . − . − . . . . . . . Figure 10. The renormalized values of various widths of themultifractal spectra are shown and compared for order 5 ofthe method. On left hand side we present clusterization ofthe states in euclidean metric. the shuffled multifractal spectrums is present for NorthIreland, Spain and a cluster of states in the central Eu-rope, namely Germany, Poland, Czech Republic and Slo-vakia. Interestingly, the cluster is investigated in Ref.[1] from the point of view of the trans-border flows re-veals enormous production of the electric energy from thewind power plants in north Germany and consumptionin south Germany with the enormous cross-border elec-tric power traffic spanning Poland, Czech Republic andSlovakia.The autocorrelation multifractal spectrum is usuallywider than the shuffled one. We observe central block ofcountries with relatively narrow spectrum where coun-tries with wider spectrum are placed randomly. Thewidest spectrum is present for Iceland, Poland, Estoniaand Spain while the narrowest spectrum is present forFrance, Latvia and Bulgaria. We refer reader to the Ta-ble I where results of MFDFA for order 4 are presented.We show the space of the widths of the shuffled andautocorrelation multifractal spectrums in the Fig. 9.Most states operate within an elliptic cluster with centerat (0 . , . a semi-axes (0 . , . with few deviationslike Iceland, Latvia, Spain and Poland which operatewith much wider spectrums. Additionally North Ireland,Latvia, Luxembourg and to some degree additional statesdevelop, due to their location of maximums of the auto-correlation spectrum π cormax close to and considerablewidths of spectrums, Formulas 1 and 20. D. Clusterization of the states
We employed Machine learning tools to search for clus-ters to put out previous observations of clusterization on solid ground, see [54–57]. Since the widths of multifrac-tal spectra are order dependent the clusters presence varyand more than 2 clusters were observed for orders 2, 3and 4 for different methods of clusterization, metric andsubset of the widths of spectra.Most countries of Europe operate power grid systemsin the region with a limited variety of Hurst exponents,see red cluster in Fig. 10. The deviations are present andthat are namely Iceland, Spain, Estonia and Poland(onlyfor the order 5). However, we observed additional pres-ence of Germany, Czech Republic, Slovakia, Portugal andNorth Ireland for different setup of clusterization algo-rithm. For Iceland, portfolio of sources of electricity iswidely different than in the rest of Europe. Presence ofGermany, Czech Republic, Slovakia a Poland in deviat-ing group may be caused by crossborder transfer in cetralEurope, see [1].
E. Note on implementation and robustness ofresults
We note that we performed the analysis for various or-ders of the method with minimal impact on the conclu-sions. In the table I we report the complete multifractalanalysis of European countries including the widths ofall multifractal spectrums as well as the locations of themaxims of all multifractal spectrums for order 4 of theMFDFA.The method have been implemented in
Zarja library[58] with appropriate Python interface. Pythonis then used for post-processing of the datasets. Dur-ing the analysis we employed Matplotlib [59], Basemap,SciPy [60], NumPy and Scikit tools [55, 56]. The analysisusing the MFDFA gives the best results for orders 4, 5and 6. The lower orders lacks effective detrending butelectic power grids are relatively stable in this respect.On the other hand, we also admit that the analysiscan be inaccurate for both large | q | and the high ordersof the MFDFA. It is caused by use of floating point num-bers (IEEE 754) that are inaccurate, see [48], due to thefixed length of the data-type. To solve the problem werefer to the use the fractions or to the variable lengthnumbers, however the use is at the expense of the speed.Implementing MFDFA can solve the problem. IV. CONCLUSIONS
We performed the modification of the MFDFA methodon the dataset obtained from the electric power grid ofthe European countries which consists of two level of de-trending. Dataset analyzed are the electric power loadswith one hour frequency. In the study we focus on theproperties of the autocorrelation function and the prob-ability distribution of the stochastic process governingevolution of the system.8o address whether the properties of the distributionand autocorrelation function are unique or more parame-ters are involved we employ multifractal spectrum of theautocorrelation and distribution function. All countriesadmit multifractality of both the autocorrelation and dis-tribution functions. To be more specific, all datasets ad-mit existence of the non-Gaussian probability distribu-tion due to the width of multifractal spectra as well asthe long-range autocorrelation function.We detected presence of the countries that posses theparameters that Langevin equation for both the coordi-nates and the velocities have to be employed. Generally,we verified the robustness of our analysis performing thesurrogated datasets as well as various orders of detrend-ing and no non-linear effects spoils the analysis. The loworders of detrending seems to be inefficient and too highorders are influenced by numeric instability but the mainmessage is robust.The analysis of the spatial properties of the noisewithin the electric grids among the states does not dif-fer so much. The properties of the distribution func-tion among the states are practically the same but thedifferences are present for the autocorrelation function.Few states possess the most frequent description by thederivatives of the power load and the weak anti-persistentprocess. Few states (Luxembourg, Latvia and North Ire-land) follow the strong persistent process of the electricpower loads and most states follow description by thederivatives and strong anti-persistent process. The widewidth of the multifractal spectrums is mostly scatteredamong Europe except for the cluster of central Europeancountries where the width of the shuffled multifractalspectrum can be understood as an effect of the imbalancedue to installations of the wind power plants in northGermany. The problem should be also addressed moreclosely using different datasets. InterestiThe wide width of the autocorrelation multifractalspectrum is scattered on the map of Europe. Since it ispresent for Iceland we believe that it shall be related tothe hidden connection with the structure of the powerproduction. We emphasize that more analysis of thewider set of datasets is needed.To put the paper in the context of the analysis of thetime series and particularly in the analysis electric elec-tric grids, we exploited widely used methodology andimproved it with the additional tests and detrendingmethod. We performed deep analysis of stochastic prop- erties of the time series for European countries that goesbeyond trends and we address the properties of the noiseand its effective description by Langevin equation.Practically, in the long horizon, the study shall be use-ful for unbiased description of the electric power gridsusing the stochastic processes and the effective handlingof the risk management. Based on the analysis we focuson use of the Fractionally integrated autoregressive con-ditional heteroscedastic processes (FIARCH), see [5, 61]or Autoregressive fractionally integrated moving averageprocess (ARFIMA) [62–64] that are able to generate thetime series with the power-law correlations for predic-tion of the electric power loads. Recently, generalizationof both processes is proposed in [65] for symmetric dis-tribution functions. However, our previous study [38] re-vealed that the distribution function is asymmetric evenmore it may not fit to any element from the set of theLévy stable distributions, see [62, 65, 66]. It supportsa hypothesis that the poles of the Mellin transform, see[39], of the distributions are not stationary.
AUTHOR CONTRIBUTIONS
J.K. obtained and prepared the dataset. J.K and H.L.developed the main ideas of the paper. H.L. preparedthe tool for analysis for UNIX-like system. Both authorsalso performed literature search, the analysis of the timeseries and the visualization of the results. J.K and H.Lboth contributed to the writing of the manuscript andthe work was performed under H.L. leadership. The workdescribed in this paper will be used in J.K.’s Ph.D. thesis.
ACKNOWLEDGMENTS
The analysis exhibits personal view of authors of elec-troenergetics. No authority or a company did influencethe analysis. We acknowledge fruitful discussions withP. Jizba, J. Lavička, E. Lutz, T. Kiss, G. Alber, E. Gil,R. Weron, M. Ausloos, L.Matsuoka and H.E. Stanley.It was also supported by Czech Ministry of EducationRVO68407700. This thesis benefited from the EuropeanUnion’s Horizon 2020 Research and Innovation Staff Ex-change programme under the Marie Sklodowska-Curiegrant agreement No 681228. [1] Z. Boldiš, “Czech electricity grid challenged by Germanwind,”
Europhys. News , vol. 44, no. 4, pp. 16–18, 2013.[2] R. Weron,
Modeling and Forecasting Electricity Loadsand Prices: A Statistical Approach . John Wiley and SonsLtd, 2006.[3] C. Harris,
Electricity markets: Pricing, Structures andEconomics . John Wiley and Sons Ltd, 2006. [4] G. E. Uhlenbeck and L. S. Ornstein, “On the theory ofbrownian motion,”
Phys.Rev. , vol. 36, pp. 823–841, 1930.[5] G. Rangarajan and M. Ding, eds.,
Processes with long-range correlations , vol. 621. Springer-Verlag Berlin Hei-delberg, 2003.[6] R. Mantegna and H. Stanley,
An introduction to econo-physics: correlations and complexity in finance . NewYork, NY, USA: Cambridge University Press, 2000. idths of multifractal spectrum Position of maximaCountry ∆ π ∆ π shuf △ π sur △ π cor △ π dist π max π shufmax π surmax π cormax π distmax Austria 0.370 0.175 0.362 0.216 0.102 1.266 0.494 1.232 0.772 -0.018Belgium 0.383 0.215 0.373 0.199 0.099 1.123 0.499 1.110 0.624 -0.029Bosnia-Herzegovina 0.363 0.190 0.310 0.184 0.095 1.137 0.493 1.116 0.644 -0.008Bulgaria 0.378 0.199 0.335 0.192 0.064 1.121 0.500 1.120 0.621 -0.010Croatia 0.374 0.199 0.304 0.176 0.071 1.169 0.497 1.154 0.672 0.010Cyprus 0.461 0.199 0.412 0.290 0.089 1.255 0.495 1.279 0.760 -0.017Czech Republic 0.426 0.227 0.380 0.209 0.090 1.133 0.501 1.122 0.632 -0.021Denmark a a b c a d a a Excludes oversea dominions. b Excluding North Ireland, Guersey, Jersey, Isle of Man and oversea dominions, includes England, Wales and Scotland. c FYROM Macedonia d Part of Great Britain with separated power grid system.
Table I. Collection of results of the MFDFA analysis of the order .[7] J. Kantelhardt, S. Zschiegner, E. Koscielny-Bunde,S. Havlin, A. Bunde, and H. Stanley, “Multifractal de-trended fluctuation analysis of nonstationary time se-ries,” Physica A , vol. 316, p. 87, 2002.[8] J. Kantelhardt, E. Koscielny-Bunde, H. Rego, S. Havlin,and A. Bunde, “Detecting long-range correlations withdetrended fluctuation analysis,”
Physica A , vol. 295,p. 441, 2001. [9] C.-K. Peng, S. Havlin, H. E. Stanley, andA. L.Goldberger, “Quantification of scaling exponentsand crossover phenomena in nonstationary heartbeattime series,”
Chaos , vol. 5, pp. 82–87, 1995.[10] A. L. Goldberger, L. A. N. Amaral, L. Glass, J. M. Haus-dorff, P. C. Ivanov, R. G. Mark, J. E. Mietus, G. B.Moody, C.-K. Peng, and H. E. Stanley, “PhysioToolkit,and PhysioNet: Components of a new research resourcefor complex physiologic signals,”
Circulation , vol. 101, p. e215–e220, 2000.[11] C.-K. Peng, J. Mietus, J. M. Hausdorff, S. Havlin, H. E.Stanley, and A. L. Goldberger, “Long-range anticorrela-tions and non-gaussian behavior of the heartbeat,” Phys.Rev. Lett , vol. 70 (9), pp. 1343–1346, 1993.[12] H. E. Hurst, “Long term storage capacity of reservoirs,”
Trans. Am. Soc. Civ. Eng. , vol. 116, no. 770, 1951.[13] K. Hu, P. C. Ivanov, Z. Chen, P. Carpena, and H. E.Stanley, “Effect of trends on detrended fluctuation anal-ysis,”
Phys. Rev. E , vol. 64, p. 011114, 2001.[14] Z. Wu, N. E. Huang, S. R. Long, and C. K. Pen,“On the trend, detrending, and variability of nonlinearand nonstationary time series,”
PNAS , vol. 104, no. 38,pp. 14889–14894, 2007.[15] D. Horvatic, H. E. Stanley, and B. Podobnik, “Detrendedcross-correlation analysis for non-stationary time serieswith periodic trends,”
Europhys. Lett. , vol. 94, no. 18007,2011.[16] C. V. Chianca, A. Ticona, and T. J. P. Penna, “Fourier-detrended fluctuation analysis,”
Physica A , vol. 357,pp. 447–454, 2005.[17] X.-Y. Qian, G.-F. Gu, and W.-X. Zhou, “Modified de-trended fluctuation analysis based on empirical modedecomposition for the characterization of anti-persistentprocess,”
Physica A , vol. 390, pp. 4388–4395, 2011.[18] B. Podobnik and H. E. Stanley, “Detrended cross-correlation analysis: A new method for analyzing twononstationary time serie,”
Phys. Rev. Lett. , vol. 100,no. 084102, 2008.[19] G.-F. Gu and W.-X. Zhou, “Detrended fluctuation anal-ysis for fractals and multifractals in higher dimensions,”
Phys. Rev. E , vol. 74, no. 061104, 2007.[20] C.-K. Peng, S. V. Buldyrev, S. Havlin, M. Simons, H. E.Stanley, and A. L. Goldberger, “Mosaic organization ofDNA nucleotides,”
Phys. Rev. E , vol. 49, pp. 1685–1689,1994.[21] C.-K. Peng, S. V. Buldyrev, A. L. Goldberger, S. Havlin,F. Sciortino, M. Simons, and H. E. Stanley, “Long-rangecorrelations in nucleotide sequence,”
Nature , vol. 356,pp. 168–170, 1992.[22] T. Preis, P. Virnau, W. Paul, and J. J. Schneider, “Accel-erated fluctuation analysis by graphic cards and complexpattern formation in financial markets,”
New J. Phys. ,vol. 11, p. 093024, 2009.[23] A. Bunde, S. Havlin, J. W. Kanterhardt, T. Penzel, J. H.Peter, and K. Voigt, “Correlated and uncorrelated regionsin heart-rate fluctuations during sleep,”
Phys. Rev. Lett. ,vol. 85, pp. 3736–3739, 2000.[24] P. Talkner and R. O. Weber, “Power spectrum and de-trended fluctaution analysis: Application to daily tem-peratures,”
Phys. Rev. E , vol. 64, pp. 150–160, 2000.[25] T. Preis, W. Paul, and J. J. Schneider, “Fluctuationpatterns in high-frequency financial asset returns,”
Eu-rophys. Lett. , vol. 82, p. 68005, 2008.[26] J. W. Kantelhardt, Y. Askenazy, P. C. Ivanov, A. Bunde,S. Havlin, T. Penzel, J.-H. Peter, and H. E. Stanley,“Characterization of sleep stages by correlations in themagnitude and sign of heartbeat increments,”
Phys. Rev.E , vol. 051908, 2002.[27] K. Koçak, “Examination of persistence properties of windspeed records using detrended fluctuation analysis,”
En-ergy , vol. 34, pp. 1980–1985, 2009.[28] S. Hosseinabadi and M. Rajabi, “Stochastic and fractalproperties of silicon and porous silicon rough surfaces,”
J. of Phys.:Conf. ser. , vol. 454, no. 012035, 2013.[29] D. Stauer and D. Sornette
Physica A , vol. 252, no. 271,1998.[30] K. Ivanova and M. Ausloos, “Application of the de-trended fluctuation analysis (dfa) method for describingcloud breaking,”
Physica A , vol. 274, no. 1-2, p. 349,1999.[31] R. Weron and A. Przybylowicz, “Hurst analysis of elec-tricity price dynamics,”
Physica A , vol. 283, 2000.[32] R. Weron, “Energy price risk management,”
Physica A ,vol. 285, p. 127, 2000.[33] I. Simonsen, “Measuring anti-correlations in the nordicelectricity spot market by wavelets,”
Physica A , vol. 322,pp. 597–606, 2003.[34] P. Norouzzadeh, W. Dullaert, and B. Rahmani, “Anti-correlation and multifractal features of spain electricityspot market,”
Physica A , vol. 380, pp. 333–342, 2007.[35] S. Chan, K. Tsui, H. Wu, Y. Hou, Y.-C. Wu, and F. Wu,“Load/price forecasting and managing demand responsefor smart grids: Methodologies and challenges,”
IEEESignal Processing Magazine , pp. 68 – 85, 2012.[36] T. Hong, P. Pinson, and S. Fan, “Global energy forecast-ing competition 2012,”
Int. J. Forecasting , vol. 30, no. 2,pp. 357–363, 2014.[37] R. Weron, “Electricity price forecasting: A review of thestate-of-the-art with a look into the future,”
Int. J. Fore-casting , vol. 30, pp. 1030–1081, 2014.[38] J. Kracík and H. Lavička, “Fluctuation analysis of highfrequency electric power load in the czech republic,”
Physica A: Statistical Mechanics and its Applications ,vol. 462, pp. 951 – 961, 2016.[39] A. A. Kilbas and M. Saigo,
H-Transforms: Theory andApplication , vol. Analytical methods and special func-tions of
An International Series of Monographs in Math-ematics . Chapman and Hall, 2004.[40] A. Mathai and H. Hausbolt,
Special Functions for Ap-plied Scientists . Springer, 2008.[41] H. Glaeske, A. Prudnikov, and K.A.Skórnik,
OperationalCalculus and Related Topics , vol. Analytical methods andspecial functions of
An International Series of Mono-graphs in Mathematics . Chapman and Hall, 2006.[42] Y. A. Brychkov,
Handbook of Special Functions, Deriva-tives, Integrals, Series and Other Formulas . CRC Press,2008.[43] R. Beals and R. Wong,
Special Functions . CambridgeUniversity Press, 2010.[44] B. Davies,
Integral Transforms and Their Applications .Springer, 2002.[45] W. Magnus, F. Oberhettinger, and F.G.Tricomi,
HigherTranscendental Functions , vol. III. McGraw-Hill BookCompany, 1953.[46] N. Lebedev,
Special Functions nad Their Applications .Translation from Russian Prentice-Hall, 1965.[47] χ m ( α ) = γ (cid:0) m , α (cid:1) , where γ is incomplete Gamma func-tion.[48] D. Knuth, The Art of Computer Programming-Seminumerical Algorithms (Third Edition) , vol. 2.Addison-Wesley, 1997.[49] D. Harte,
Multifractals: Theory and applications . Chap-man and Hall, 2001.[50] Web page .[51] The datasets does not include electric power productionfor local needs – powerplants for a factory or a house –and also power production within islands disconnected rom the mainland.[52] Then the autocorrelation Hurst exponent defined by Eq.13 is in range (cid:2) − , (cid:3) .[53] Well, Latvia is just at the border.[54] A. Boschetti and L. Massaron, Python Data Science Es-sentials . Packt Publishing, 2015.[55] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel,B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer,R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cour-napeau, M. Brucher, M. Perrot, and E. Duchesnay,“Scikit-learn: Machine Learning in Python,”
Journal ofMachine Learning Research , vol. 12, pp. 2825–2830, Oct.2011.[56] R. Garreta and G. Moncecchi,
Learning scikit-learn: Ma-chine Learning in Python . Packt Publishing, 2013.[57] A. Müller and S. Guido,
Introduction to Machine Learn-ing with Python . O’Reilly Media, 2017.[58] http://sourceforge.net/projects/zarja/.[59] J. D. Hunter, “Matplotlib: A 2d graphics environment,”
Computing In Science & Engineering , vol. 9, no. 3, pp. 90–95, 2007.[60] E. Jones, T. Oliphant, P. Peterson, et al. , “SciPy: Opensource scientific tools for Python,” 2001–. [Online; ac-cessed 2017-01-30].[61] C. Granger and Z.Ding, “Varieties of long-range memorymodels,”
J. Econometr. , vol. 73, pp. 61–77, 1996.[62] G. Samorodnitsky and M. S. Taqqu,
Stable Non-Gaussian Random Processes: Stochastic Models with In-finite Variance.
Chapman and Hall, New York, 1994.[63] J. R. M. Hosking, “Fractional differencing,”
Biometrika ,vol. 68, no. 1, pp. 165–176, 1981.[64] C. Granger and R. Joyeux, “An introduction to long-memory time series models and fractional differencing,”
J. Time Ser. Anal. , vol. 1, no. 15, pp. 15–29, 1980.[65] B. Podobnik, P. C. Ivanov, K. Biljakovic, D.Horvatic,H.E.Stanley, and I.Grosse, “Fractional integrated pro-cess with power-law correlations in variables and mag-nitudes,”
Phys. Rev. E , vol. 026121, 2005.[66] K.-I. Sato,
Lévy Processes and Infinitely Divisible Distri-butions . Cambridge University Press, 1999.. Cambridge University Press, 1999.