The Recurrence Interval Difference of Power Load in Heavy/Light Industries of China

: The signiﬁcant ﬂuctuation of industrial electricity consumption has a high impact on power load, which makes the research on recurrence intervals between extreme events of theoretical and practical signiﬁcance. The study uses a high-frequency data of heavy and light industries and employs recurrence interval analysis in different thresholds. We ﬁnd that the reoccurrence interval of volatility can ﬁt with the stretched exponential function and the probability density functions of recurrence intervals in various thresholds shows a scaling behavior. Then, the conditional probability density function and the multifractal detrended ﬂuctuation analysis demonstrate the existence of short-range correlation, long-range correlation, and multifractal properties, respectively. We further construct a hazard function, introduce recurrence intervals into VaR calculation and establish a functional relationship between average recurrence interval and threshold. Following this result, we also shed light on policy discussion for multi-industrial electricity supply management.


Introduction
Over the past few decades, China has snowballed and surpassed Japan in 2010 to become the world's second-largest economy [1]. As for energy consumption, China's electricity consumption has also exceeded the United States in 2011 and ranked number one in the world [2]. In 2016, China's total electricity consumption was 5.92 trillion kWh with an annual growth rate of 5.0%. Figure 1 shows the proportion of electricity consumption of each industry in China. Compared to 2015, electricity consumption of the tertiary industry increased by 11.2% and continued to maintain a high growth rate, indicating that service consumption led the growth of China's economy. The electricity consumption of urban and rural household industry, secondary industry, and manufacturing industry increased by 10.8%, 2.9%, and 2.5%, respectively. For the four major energy-consuming sectors (Steel, Nonferrous Metals, Building Materials, Chemistry), the electricity consumption barely increased. Such non-increase is because it is apparently for equipment manufacturing, emerging technologies, and mass consumer goods industries, reflecting that the readjustment and upgrading of the industrial structure have a positive effect and the power consumption structure is continuously optimized. By the end of 2016, China's total power generation installed capacity was 1.65 billion kW with an annual increase of 8.2%, which further aggravated overcapacity in some areas. The non-fossil energy generation continued to increase, and the utilization of thermal power equipment further reduced to 4165 hours, which is the lowest since 1964. In summary, the overall electricity supply and demand in China was loose with a relative surplus in some areas [3]. Regarding price, the electricity price for industries is negotiated reduced to 4165 hours, which is the lowest since 1964. In summary, the overall electricity supply and demand in China was loose with a relative surplus in some areas [3]. Regarding price, the electricity price for industries is negotiated in most times which varies in different regions. For example, we can see from Figure 2 that the price in developed areas such as Shanghai is higher than that in developing areas like Xinjiang (See http://www.askci.com/news/chanye/20161020/13442871189.shtml).  At the same time, the complexity of power consumption structure was also rising all over the country, which made regional electricity management a huge challenge both for public administrations and for power supply corporations. Therefore, since the year 2015, China has started a new round of electricity reform, which mainly included projects like power transmission, distribution price reform, energy market construction, electricity trading institution establishment and electricity supply reform. All of these projects require China's electricity supply companies to rebuild their organization structure and upgrade management plans to deal with more intense market competition. Accurately predicting the trend and change of regional/sector electricity consumption by energy consumption modeling can help electricity companies optimize their management strategies, reduce operation cost and prevent potential power supply risks, and can also help the nation save resources, avoid waste and reduce infrastructure over-investment.
As shown in Figure 1, China's energy consumption concentrated in the industrial sector. For enterprises belonging to different industry sectors (i.e., heavy industry or light industry), the power loads have quite a different magnitude. As is well known, people are more interested in large electricity fluctuations than small ones. For example, in a particular area, if the power load of light Energies 2018, 11, 106 2 of 20 reduced to 4165 hours, which is the lowest since 1964. In summary, the overall electricity supply and demand in China was loose with a relative surplus in some areas [3]. Regarding price, the electricity price for industries is negotiated in most times which varies in different regions. For example, we can see from Figure 2 that the price in developed areas such as Shanghai is higher than that in developing areas like Xinjiang (See http://www.askci.com/news/chanye/20161020/13442871189.shtml).  At the same time, the complexity of power consumption structure was also rising all over the country, which made regional electricity management a huge challenge both for public administrations and for power supply corporations. Therefore, since the year 2015, China has started a new round of electricity reform, which mainly included projects like power transmission, distribution price reform, energy market construction, electricity trading institution establishment and electricity supply reform. All of these projects require China's electricity supply companies to rebuild their organization structure and upgrade management plans to deal with more intense market competition. Accurately predicting the trend and change of regional/sector electricity consumption by energy consumption modeling can help electricity companies optimize their management strategies, reduce operation cost and prevent potential power supply risks, and can also help the nation save resources, avoid waste and reduce infrastructure over-investment.
As shown in Figure 1, China's energy consumption concentrated in the industrial sector. For enterprises belonging to different industry sectors (i.e., heavy industry or light industry), the power loads have quite a different magnitude. As is well known, people are more interested in large electricity fluctuations than small ones. For example, in a particular area, if the power load of light  At the same time, the complexity of power consumption structure was also rising all over the country, which made regional electricity management a huge challenge both for public administrations and for power supply corporations. Therefore, since the year 2015, China has started a new round of electricity reform, which mainly included projects like power transmission, distribution price reform, energy market construction, electricity trading institution establishment and electricity supply reform. All of these projects require China's electricity supply companies to rebuild their organization structure and upgrade management plans to deal with more intense market competition. Accurately predicting the trend and change of regional/sector electricity consumption by energy consumption modeling can help electricity companies optimize their management strategies, reduce operation cost and prevent potential power supply risks, and can also help the nation save resources, avoid waste and reduce infrastructure over-investment.
As shown in Figure 1, China's energy consumption concentrated in the industrial sector. For enterprises belonging to different industry sectors (i.e., heavy industry or light industry), the power loads have quite a different magnitude. As is well known, people are more interested in large electricity fluctuations than small ones. For example, in a particular area, if the power load of light industry Energies 2018, 11, 106 3 of 20 enterprise reaches its peak, the power supply network in the area will not feel the pressure. However, if a heavy industry enterprise exceeds its peak power, this may pose an enormous risk of grid paralysis or even serious accidents, which would usually spread and cause an immeasurable damage and suffering to the society. Such differences require power supply companies to differentiate their load management strategies for different industries. Hence, analysis and comparison of heavy/light industrial electricity consumption patterns, especially the extreme electricity events will provide a theoretical foundation and practical guidance, which can help power supply companies avoid grave accidents, ensure the stability of power supply operation and keep power transfer equipment safe.
Due to the availability of real-time power load data, forecasting extreme events in electricity market becomes possible. In this paper, we use 15-min high-frequency electricity load data from two companies (one belongs to heavy industry, another belongs to light industry) for the city of Nanjing in China to study the recurrence pattern between fluctuations. The questions we tried to answer are as follows: How to describe the volatility behavior? If a massive volatility occurred, will it happen again or not? When will the next significant fluctuation occur? How to estimate the occurrence probability and the time interval between two extreme events?
Our study, therefore, contributes to the literature in the following perspectives. First, among plenty studies that have done in the electric field, research on the electricity consumption pattern between heavy and light industry is rare. So far as we know, this is the first paper to focus on the industrial difference of power load from a policy perspective. Second, as a major city with a fast-growing economy in China, Nanjing, the capital of Jiangsu province, has an urgent demand for energy especially electricity. Research on Nanjing can be a useful guide for other emerging first-tier cities in China. Most previous studies on Jiangsu area only use daily data [4][5][6][7] while here we use high-frequency data which can help to capture more specific properties. Finally, different from many studies using power law to model rare events in many fields [8][9][10][11][12][13] including electricity markets [14,15], we use the stretched exponential function to analyze the extreme events [16][17][18]. Unlike the power-law which displays a linear relationship in log-log plots, stretched exponential function has a significant curvature in log-log plots which can better describe the natural phenomena in nature and economy [19].
The rest of this paper is arranged as follows: Section 2 reviews current studies. Section 3 describes the method and presents basic statistics of the data set. Section 4 initials an empirical research, including distribution function, scaling properties, memory effect and risk estimation. Section 5 delivers management policy implications discussions for China's power supply company. Section 6 concludes.

Literature Review
Currently, plenty of models have been applied to forecast electricity load, and most of them can be classified into three categories: regression, data series analysis, and neural network.
Regression analysis considered power load as a dependent variable of other factors, such as weather or holidays [20], and tried to construct a function relationship between electricity load and other influence factors [21]. Such results were confirmed by the reality in some residential areas [22]. With technology evolving, a decomposition model was applied, and many vital factors were pinned in the evolution progress [23]. Also, GDP [24], income [24], seasonal variations [25], etc. were all confirmed to be relevant to electricity load forecast by regression analyses. In the 21st century, the emergence of smart technology could help researchers grab high-frequency data from a personal level and add human behavior into a regression model to improve the forecast accuracy [26,27]. Briefly, regression analysis aimed to describe the quantitative relationship between the observed variables in statistics, however, it could not capture the spatiotemporal variation and is sometimes restricted by the data volume.
Time series analysis considered that the electricity load pattern was a time series signal, which could be used to predict the future load with historical data [28]. With technology development, various time series models have been proposed [29]. Carmine et al. [30]  detected a periodic pattern, which could be used to improve electricity management. Dong et al. [31] proposed a hybrid model to predict residential load hour and day ahead with five different algorithms and the results improved prediction accuracy. Liang and Liang [32] proposed a mathematical hybrid method to analyze electricity demand in China and predicted the demand evolution in the next five years from 2016 to 2020. Though the predictable time spectrum of time series analyses varies from one hour to several years, the claim to the accuracy of historical data was very high. The algorithm might be complex and unstable in some cases. Additionally, when applying in short-term power load forecast, it was not sensitive to weather factors and not able to solve the inaccuracy problem caused by meteorological factors.
Artificial neural networks (ANNs) could estimate future loads with previous data without the presumption of the functional relationship between electricity load and other relative variables [29]. Compared with traditional methods, ANN was better at dealing with nonlinear relationships [33]. At present, the neural network was well applied in various fields including electric load forecasting [34][35][36][37]. Moreover, the integration of ANN and other methods has become a research hotspot. ANN was integrated with fuzzy logic [38], genetic algorithm [39], wavelet analysis [40,41], chaos theory [42], grey system [43], etc. However, ANN has its limitations: it is hard to avoid learning deficiency or over-fit phenomena, and the convergence speed is slow and easy to fall into local minima.
In this paper, we will use a different method: recurrence interval analysis (RIA), which belongs to time series analysis, to investigate the fluctuation characteristics of heavy/light industries. RIA has been widely used to analyze extreme events in many fields, for example, climate [44], earthquake activities [45,46], heartbeat monitoring [47,48] and financial volatility [49,50]. These events occur with a high magnitude and low probability. For a long time in the past, extreme events were assumed to be spontaneous, namely mutually uncorrelated. Nevertheless, researchers show that extreme events do not occur independently; contrarily, they gather and happen in a relatively short period. RIA focuses on the electricity fluctuation rather than the power load, which is good at depicting the volatility pattern and estimating the risk. Assuming occurrence probability of future events is constant and dependent on past events, RIA can estimate the probability of the events to happen again in the future. Here, Recurrence interval means the time interval between two consecutive events beyond (under) one certain threshold, which can be positive or negative. It is also powerful to characterize the fluctuation behaviors in different magnitudes [51] and sometimes even to construct a relationship among them.

Data Description
The data in this paper is extracted from the transformer of each corresponding area in 15-min frequency, and we chose two enterprises, which represented heavy and light industries, respectively. In the national economy, heavy and light industries refer to the sectors which provide capital and consumption goods, respectively. Heavy industry includes steel industry, metallurgical industry, machinery, energy (electricity, oil, coal, natural gas, etc.), chemistry and materials science. These industries provide technical equipment, power and raw materials for all branches of the national economy. Light industry mainly includes foodstuff, textile, leather, papermaking, daily chemical industry, culture, education and sporting goods industry. According to the data from National Bureau of Statistics of China, the finished goods of heavy and light industry were 2.608 and 1.341 trillion yuan in 2015, respectively (See http://data.stats.gov.cn/easyquery.htm?cn=C01&zb=A0E010D&sj=2015). In this paper, we select a steel/textile enterprise to represent the heavy/light industry. The sample period was from 1 January 2016 through 31 December 2016 and finally we ended up with 35,316 electricity load observations. The returns of time series are calculated by taking differences of the load as: where l(t) is the electricity load (unit: kW) of the tth time and ∆t = 15 because the data is 15-min frequent. Figure 3 shows the returns of heavy/light industries and Figure 4 is the frequency distribution of fluctuation. The statistics are summarized in Table 1.  Table 1 show that the returns are not normally distributed. One has a sharp peak, and another has platy kurtosis. Though their average and standard deviation values have magnitude difference, the two industries also have common properties, like skewness. From Figure 3, it is clear that there are zero values in Figure 4b, which means in a specified period, the light industry enterprise was not in operation. When combining with Table 1, we can see that even a large company could not operate continually, which is the reason why the minimum values are both zero and why we used different values to describe fluctuation. where ( ) is the electricity load (unit: kW) of the th time and ∆ = 15 because the data is 15-min frequent. Figure 3 shows the returns of heavy/light industries and Figure 4 is the frequency distribution of fluctuation. The statistics are summarized in Table 1. Table 1 show that the returns are not normally distributed. One has a sharp peak, and another has platy kurtosis. Though their average and standard deviation values have magnitude difference, the two industries also have common properties, like skewness. From Figure  3, it is clear that there are zero values in Figure 4b, which means in a specified period, the light industry enterprise was not in operation. When combining with Table 1, we can see that even a large company could not operate continually, which is the reason why the minimum values are both zero and why we used different values to describe fluctuation.   The x-ray of Figure 5 is the value of threshold , and y-ray is the number of extreme events of ;

Figures 3 and 4 and
measures the volatility of the normalized ( ). Due to magnitude discrepancy of electricity consumption between heavy and light industries, we normalize the fluctuation for comparison purposes and the results are shown in Figure 5. It can be clearly seen that the number decreases with an increasing , indicating that a larger fluctuation has less occurring probability, which is in accordance with reality. The slope is getting smaller when q increases for each curve, and there is an intersection point in the figure, indicating that small fluctuations more likely happen in heavy industry and larger volatilities more likely happen in light industry. where ( ) is the electricity load (unit: kW) of the th time and ∆ = 15 because the data is 15-min frequent. Figure 3 shows the returns of heavy/light industries and Figure 4 is the frequency distribution of fluctuation. The statistics are summarized in Table 1. Table 1 show that the returns are not normally distributed. One has a sharp peak, and another has platy kurtosis. Though their average and standard deviation values have magnitude difference, the two industries also have common properties, like skewness. From Figure  3, it is clear that there are zero values in Figure 4b, which means in a specified period, the light industry enterprise was not in operation. When combining with Table 1, we can see that even a large company could not operate continually, which is the reason why the minimum values are both zero and why we used different values to describe fluctuation.   The x-ray of Figure 5 is the value of threshold , and y-ray is the number of extreme events of ;

Figures 3 and 4 and
measures the volatility of the normalized ( ). Due to magnitude discrepancy of electricity consumption between heavy and light industries, we normalize the fluctuation for comparison purposes and the results are shown in Figure 5. It can be clearly seen that the number decreases with an increasing , indicating that a larger fluctuation has less occurring probability, which is in accordance with reality. The slope is getting smaller when q increases for each curve, and there is an intersection point in the figure, indicating that small fluctuations more likely happen in heavy industry and larger volatilities more likely happen in light industry.  The x-ray of Figure 5 is the value of threshold q, and y-ray is the number of extreme events of q; q measures the volatility of the normalized r(t). Due to magnitude discrepancy of electricity consumption between heavy and light industries, we normalize the fluctuation for comparison purposes and the results are shown in Figure 5. It can be clearly seen that the number decreases with an increasing q, indicating that a larger fluctuation has less occurring probability, which is in accordance with reality. The slope is getting smaller when q increases for each curve, and there is an intersection point in the figure, indicating that small fluctuations more likely happen in heavy industry and larger volatilities more likely happen in light industry.

Probability Density Function
Consider as the recurrence interval when threshold is , the overwhelming consensus is that the recurrence intervals of volatility are distributed as a stretched exponential.
In Equation (2), ̅ is the mean recurrence interval that depends on the threshold , where the probability distribution function of recurrence interval is ( ), and , , are the parameters [16][17][18], where α and are the coefficients related to the threshold , while γ is the correlation exponent reflecting the long term memory of the recurrence intervals.
We normalize the returns of each series by dividing its standard deviation: ( ) = ( ) [ ( ) 2 − 2 ( )] 1/2 , where [ ( ) 2 − 2 ( )] 1/2 is the standard deviation of ( ). For each threshold , we can obtain a data set of recurrence interval , and then get the probability density function of . The here is positive ( > 0) because our interest is the overload operation of the electrical system. The empirical parameters of the probability distribution function ( ), calculated by maximum likelihood estimation are shown in Table 2. Note that the coefficients of the exponent γ of large enterprises are in the range between 0.448 and 0.684, while the corresponding exponents of the light enterprises are much smaller within the interval (0.163, 0.171). It is obvious that the long-term memory effect in the recurrence intervals of electricity usage in large enterprises dies out faster than that in light companies.

Probability Density Function
Consider τ as the recurrence interval when threshold is q, the overwhelming consensus is that the recurrence intervals of volatility are distributed as a stretched exponential.
In Equation (2), τ is the mean recurrence interval that depends on the threshold q, where the probability distribution function of recurrence interval τ is P q (τ), and α, β, γ are the parameters [16][17][18], where α and β are the coefficients related to the threshold q, while γ is the correlation exponent reflecting the long term memory of the recurrence intervals.
We normalize the returns of each series by dividing its standard deviation: is the standard deviation of r(t). For each threshold q, we can obtain a data set of recurrence interval τ, and then get the probability density function of τ. The q here is positive (q > 0) because our interest is the overload operation of the electrical system. The empirical parameters of the probability distribution function P q (τ), calculated by maximum likelihood estimation are shown in Table 2. Note that the coefficients of the exponent γ of large enterprises are in the range between 0.448 and 0.684, while the corresponding exponents of the light enterprises are much smaller within the interval (0.163, 0.171). It is obvious that the long-term memory effect in the recurrence intervals of electricity usage in large enterprises dies out faster than that in light companies.  Figure 6 shows that first, for large fluctuations, P q (τ) decays slowly, indicating that with q rises, recurrence intervals will be longer, in accordance with the fact that large fluctuations have more long intervals and less short intervals than small fluctuations. Second, when threshold q is fixed, P q (τ) decreases with increasing recurrence interval τ. This implies that when the previous fluctuation occurred t units of time ago, if ∆t 1 < ∆t 2 , the probability P q (τ) of the next fluctuation to occur after ∆t 1 is larger than ∆t 2 . The result suggests that if the enterprise suffers a huge electricity load rise, the enterprise should be aware that the probability that if another large load rise occurs within a shorter time interval it is likely to be larger with a longer time interval. Moreover, when recurrence interval is fixed, P q (τ) increases with q. If an extreme event occurs, the probability that another extreme event occurs will be higher for a larger threshold q. Last, comparing Figure 6a,b, we can find for a same threshold q, P q (τ) decreased faster along with τ for heavy enterprises than low enterprises, indicating that the probability of large fluctuations of heavy enterprise is less than light enterprise, which in accordance with the fact that heavy industry is less risk-taking to large volatility than light industry.
Energies 2018, 11,106 7 of 20 Figure 6 shows that first, for large fluctuations, ( ) decays slowly, indicating that with rises, recurrence intervals will be longer, in accordance with the fact that large fluctuations have more long intervals and less short intervals than small fluctuations. Second, when threshold is fixed, ( ) decreases with increasing recurrence interval . This implies that when the previous fluctuation occurred t units of time ago, if ∆ 1 < ∆ 2 , the probability ( ) of the next fluctuation to occur after ∆ 1 is larger than ∆ 2 . The result suggests that if the enterprise suffers a huge electricity load rise, the enterprise should be aware that the probability that if another large load rise occurs within a shorter time interval it is likely to be larger with a longer time interval. Moreover, when recurrence interval is fixed, ( ) increases with . If an extreme event occurs, the probability that another extreme event occurs will be higher for a larger threshold . Last, comparing Figure 6a,b, we can find for a same threshold , ( ) decreased faster along with for heavy enterprises than low enterprises, indicating that the probability of large fluctuations of heavy enterprise is less than light enterprise, which in accordance with the fact that heavy industry is less risk-taking to large volatility than light industry. From Figure 6, it can be seen that with rising, recurrence intervals will be longer, in agreement with the fact that large fluctuations have more long intervals and less short intervals than small fluctuations, which means for large fluctuations, there is a higher probability that the time intervals between two consecutive events will increase.
Besides, by the observations in Table 2 and Figure 6, we find all the probability density function curves have similar shapes, which make us wonder whether there is a scaling behavior between these probability density functions? To examine that, we use Yamasaki et al.'s [52] method. Not only discovering the scaling behavior of events of different thresholds, but it can also efficiently solve the problem of insufficient data.
where / ̅ is scaled recurrence interval and ( ) ̅ is scaled PDFs. When threshold changes, ̅ changes and there is ( ̅ )/( ) > 0 , indicating that with volatility increases, the mean time of recurrence interval increases, which is consistent with the fact. If ( / ̅ ) is independent to , then there will be a unique function f(x), for different threshold q.
Namely, the scaled probability distribution ( / ̅ ) will converge to the single curve ( / ̅ ) and recurrence intervals have a scaling behavior. To test that, we picture the scatter diagram of ( / ̅ ) as the function of / ̅ in Figure 7. It can be seen that for different thresholds , ( ) ̅ converge to one curve regardless of which enterprise it belongs, indicating that there exists a scaling behavior and we can deduce the behavior characteristics of large fluctuations by those of small fluctuations. From Figure 6, it can be seen that with q rising, recurrence intervals will be longer, in agreement with the fact that large fluctuations have more long intervals and less short intervals than small fluctuations, which means for large fluctuations, there is a higher probability that the time intervals between two consecutive events will increase.
Besides, by the observations in Table 2 and Figure 6, we find all the probability density function curves have similar shapes, which make us wonder whether there is a scaling behavior between these probability density functions? To examine that, we use Yamasaki et al.'s [52] method. Not only discovering the scaling behavior of events of different thresholds, but it can also efficiently solve the problem of insufficient data.
where τ/τ is scaled recurrence interval and P q (τ)τ is scaled PDFs. When threshold q changes, τ changes and there is (dτ)/(dq) > 0, indicating that with volatility increases, the mean time of recurrence interval increases, which is consistent with the fact. If f q (τ/τ) is independent to q, then there will be a unique function f (x), for different threshold q.
Namely, the scaled probability distribution f q (τ/τ) will converge to the single curve f (τ/τ) and recurrence intervals have a scaling behavior. To test that, we picture the scatter diagram of f q (τ/τ) as the function of τ/τ in Figure 7. It can be seen that for different thresholds q, P q (τ)τ converge to one curve regardless of which enterprise it belongs, indicating that there exists a scaling behavior and we can deduce the behavior characteristics of large fluctuations by those of small fluctuations.

Short-Term Memory
Short-term memory means that small tends to follow small 0 , and big tends to follow big 0 . The conditional probability density functions ( | 0 ) [52] can be used to reflect the occurrence probability of a recurrence interval immediately after the recurrence interval 0 . If there is no short-term memory, ( | 0 ) is independent of 0 . However, due to the insufficiency of interval sample, it is impossible to calculate ( | 0 ) for a single value 0 . In order to obtain more data, we adopt a method based on the idea of coarse-graining [53][54][55] and calculate ( | 0 ) for values of 0 in a certain interval instead of a single value 0 . Specifically, for a given threshold q, all the recurrence intervals in the set T are sorted in an increasing order, and then we divide the set T into four nonoverlapping subsets with the same size, where 1 contains the smallest quarter while 4 contains the largest. Then we estimate the conditional probability density functions: If there is no short-term memory, we should find Figure 8 shows the results of ( | 0 ) ̅ as a function of / ̅ for 0 in the smallest subset 1 (filled symbols) and the largest subset 4 (open symbols). Obviously, no matter for heavy or light enterprises, ( | ) ≠ ( | ). We also find that ( | 0 ∈ 1 ) is bigger than ( | 0 ∈ 4 ) for small / ̅ , while for big / ̅ , ( | 0 ∈ 1 ) is smaller than ( | 0 ∈ 4 ) in Figure 8a,b. There exists very short-term memory in the recurrence intervals.

Short-Term Memory
Short-term memory means that small τ tends to follow small τ 0 , and big τ tends to follow big τ 0 . The conditional probability density functions P q (τ|τ 0 ) [52] can be used to reflect the occurrence probability of a recurrence interval τ immediately after the recurrence interval τ 0 . If there is no short-term memory, P q (τ|τ 0 ) is independent of τ 0 . However, due to the insufficiency of interval sample, it is impossible to calculate P q (τ|τ 0 ) for a single value τ 0 . In order to obtain more data, we adopt a method based on the idea of coarse-graining [53][54][55] and calculate P q (τ|τ 0 ) for values of τ 0 in a certain interval instead of a single value τ 0 . Specifically, for a given threshold q, all the recurrence intervals in the set T are sorted in an increasing order, and then we divide the set T into four non-overlapping subsets with the same size, where T 1 contains the smallest quarter while T 4 contains the largest. Then we estimate the conditional probability density functions: If there is no short-term memory, we should find P q (τ|T i ) = P q (τ|T i ), i = j (7) Figure 8 shows the results of P q (τ|τ 0 )τ as a function of τ/τ for τ 0 in the smallest subset T 1 (filled symbols) and the largest subset T 4 (open symbols). Obviously, no matter for heavy or light enterprises, P q (τ|T i ) = P q (τ|T i ). We also find that P q (τ|τ 0 ∈ T 1 ) is bigger than P q (τ|τ 0 ∈ T 4 ) for small τ/τ, while for big τ/τ, P q (τ|τ 0 ∈ T 1 ) is smaller than P q (τ|τ 0 ∈ T 4 ) in Figure 8a,b. There exists very short-term memory in the recurrence intervals.

Short-Term Memory
Short-term memory means that small tends to follow small 0 , and big tends to follow big 0 . The conditional probability density functions ( | 0 ) [52] can be used to reflect the occurrence probability of a recurrence interval immediately after the recurrence interval 0 . If there is no short-term memory, ( | 0 ) is independent of 0 . However, due to the insufficiency of interval sample, it is impossible to calculate ( | 0 ) for a single value 0 . In order to obtain more data, we adopt a method based on the idea of coarse-graining [53][54][55] and calculate ( | 0 ) for values of 0 in a certain interval instead of a single value 0 . Specifically, for a given threshold q, all the recurrence intervals in the set T are sorted in an increasing order, and then we divide the set T into four nonoverlapping subsets with the same size, where 1 contains the smallest quarter while 4 contains the largest. Then we estimate the conditional probability density functions: If there is no short-term memory, we should find Figure 8 shows the results of ( | 0 ) ̅ as a function of / ̅ for 0 in the smallest subset 1 (filled symbols) and the largest subset 4 (open symbols). Obviously, no matter for heavy or light enterprises, ( | ) ≠ ( | ). We also find that ( | 0 ∈ 1 ) is bigger than ( | 0 ∈ 4 ) for small / ̅ , while for big / ̅ , ( | 0 ∈ 1 ) is smaller than ( | 0 ∈ 4 ) in Figure 8a,b. There exists very short-term memory in the recurrence intervals.

Long-Term Memory
Experience shows that scaling behavior and short-term memory is usually a foreshadowing of the long-term memory [56,57]. To figure out whether there exists long-term memory in electricity use of the high/light industries, we adopt the MF-DFA (multifractal detrended fluctuation analysis) method [56,[58][59][60].
The conventional DFA method invented by Peng et al. (1994) [61], which is used to investigate the statistical self-affinity in time series analysis, but limited to scale the second order statistical moment of one-dimensional fractal time series by computing Hurst exponent.
In mathematics, a fractal is an abstract object used to describe and simulate naturally occurring objects. Artificially created fractals can exhibit similar patterns at increasingly small scales. So far several models applied in electricity market are proposed based on fractal theory [62][63][64][65]. With the study objects becoming increasingly diverse and complicated, fractal with a single exponent (the fractal dimension) is not capable of describing the dynamics in reality, like coastlines length, stock market time series, heartbeat dynamics, real-world scenes, etc. A continuous spectrum of exponents (singularity spectrum) is needed and then emerges the multifractal [66][67][68][69]. A multifractal system is a generalization of a fractal system and can be discovered in nature, which we can apply in a variety of practical situations such as electricity demand [70,71].
In 2002, Kantelhardt et al. [72] combined multifractal with DFA and put forward the MF-DFA, which can describe the multifractal characteristics of time sequence and computes the H(p) for all p-order statistical moments. As one practical method to test whether a non-stationary time series has multifractal characteristics, MF-DFA is widely applied in the financial markets [68], molecular biology [73], disaster prevention and control [69], power pricing analysis [74,75] and electricity market [71,76]. Assume the sample X = {x k , k = 1, 2, . . . , N}, the specific steps of this algorithm are as follows.
(1) Calculate the accumulated deviation N times of the original data and construct new time series: where x is the average of X.
(4) Calculate p-order volatility function of sequence Energies 2018, 11, 106 10 of 20 when p = 2, it is standard DFA. (5) Based on Equation (10), we know F p (s) is positively correlated with s, with log-log coordinate, we have F p (s) ∼ s h p (12) where h p is Hurst exponent, for non-stationary time series, only when 0.5 < h p < 1, this series has a long-term correlation, indicating that the system has the fluctuation pattern in long-term evolution. When h p is the function of p, X(t) has multiple fractal characteristics.
In Figure 9, there are ten subfigures, and each subfigure has four sub-subfigures, which show the results of MF-DFA, respectively. It can be seen that the p-order Hurst exponent of each line is higher than 0.5 in a certain area, suggesting that long-term correlations and multifractal characteristics exist in the recurrence intervals. When h p < 0.5; it means the volatility is of anti-continuity.
Energies 2018, 11, 106 10 of 20 (5) Based on Equation (10), we know ( ) is positively correlated with , with log-log coordinate, we have (s)~ℎ (12) where ℎ is Hurst exponent, for non-stationary time series, only when 0.5 < ℎ < 1, this series has a long-term correlation, indicating that the system has the fluctuation pattern in long-term evolution. When ℎ is the function of , X(t) has multiple fractal characteristics.
In Figure 9, there are ten subfigures, and each subfigure has four sub-subfigures, which show the results of MF-DFA, respectively. It can be seen that the p-order Hurst exponent of each line is higher than 0.5 in a certain area, suggesting that long-term correlations and multifractal characteristics exist in the recurrence intervals. When ℎ < 0.5; it means the volatility is of anticontinuity.

Risk Estimation
The hazard probability function (∆ | ) is an important method to estimate risk in recurrence interval analysis. Considering the fact that units of time have passed since last large volatility is greater than occurred, what fascinates us is the probability that the next large volatility greater than will occur within ∆ units of time. The hazard probability function is as follows: we can calculate the theoretical value of ( ) with the parameters given in Table 2 if the distribution ( ) fits with a stretched exponential index, and calculate the empirical value by rewriting (∆ | ) into:

Risk Estimation
The hazard probability function W q (∆t|t) is an important method to estimate risk in recurrence interval analysis. Considering the fact that t units of time have passed since last large volatility is greater than q occurred, what fascinates us is the probability that the next large volatility greater than q will occur within ∆t units of time. The hazard probability function is as follows: (13) we can calculate the theoretical value of P q (τ) with the parameters given in Table 2 if the distribution P q (τ) fits with a stretched exponential index, and calculate the empirical value by rewriting W q (∆t|t) into: where "count τ q > t " is the number of recurrence intervals greater than t and "count t < τ q ≤ t + ∆t " is the number of recurrence intervals greater than t and less than t + ∆t for a given q. Figure 10 is the hazard function; the symbols are empirical values and the curves are the theoretical values of W q (∆t = 15|t). We can observe that they are fitted very nicely and the discrepancy between empirical and theoretical curves decreases when t increases. Furthermore, W q (∆t = 15|t) decreases with increasing t, confirming that recurrence intervals exhibit clustering behaviors and long-term memory between volatilities and theoretical values will overestimate the risk in a short time period. For a given threshold q, we can calculate the recurrence probability of extreme events.
For risk estimation, value at risk (VaR) is widely applied. Here we use loss probability density function to define the risk at rising q: q −∞ P(R)dR = P * , where P(R) is the probability density function of normalized series R(t), and P * is the rise probability. The mean recurrence interval τ q is: Energies 2018, 11, 106 14 of 20 where " ( > )" is the number of recurrence intervals greater than and " ( < ≤ + ∆ )" is the number of recurrence intervals greater than and less than + ∆ for a given . Figure 10 is the hazard function; the symbols are empirical values and the curves are the theoretical values of (∆ = 15| ) . We can observe that they are fitted very nicely and the discrepancy between empirical and theoretical curves decreases when increases. Furthermore, (∆ = 15| ) decreases with increasing , confirming that recurrence intervals exhibit clustering behaviors and long-term memory between volatilities and theoretical values will overestimate the risk in a short time period. For a given threshold , we can calculate the recurrence probability of extreme events. For risk estimation, value at risk (VaR) is widely applied. Here we use loss probability density function to define the risk at rising : ∫ ( ) −∞ = * , where ( ) is the probability density function of normalized series ( ), and * is the rise probability. The mean recurrence interval ̅ is: is the number of recurrence intervals below the threshold , ∑ , =1 is the total number of recurrence returns. Then the relation between mean recurrence interval and VaR can be expressed as: number of ( ) below total number of ( ) Equation (16) means that one/ ̅ defines the rise probability of which is depicted in Figure  11. If one wants to know the risk level at 1% probability of rise they can find the q for 1/ ̅ = 1%, which is the VaR one is looking for. N q is the number of recurrence intervals below the threshold q, ∑ τ q i=1 τ q,i is the total number of recurrence returns. Then the relation between mean recurrence interval and VaR can be expressed as: Equation (16) means that one/τ q defines the rise probability of q which is depicted in Figure 11. If one wants to know the risk level at 1% probability of rise they can find the q for 1/τ q = 1%, which is the VaR one is looking for. For risk estimation, value at risk (VaR) is widely applied. Here we use loss probability density function to define the risk at rising : ∫ ( ) −∞ = * , where ( ) is the probability density function of normalized series ( ), and * is the rise probability. The mean recurrence interval ̅ is: is the number of recurrence intervals below the threshold , ∑ , =1 is the total number of recurrence returns. Then the relation between mean recurrence interval and VaR can be expressed as: number of ( ) below total number of ( ) Equation (16) means that one/ ̅ defines the rise probability of which is depicted in Figure  11. If one wants to know the risk level at 1% probability of rise they can find the q for 1/ ̅ = 1%, which is the VaR one is looking for.

Discussion
In the context of China's new round of power reform, according to the empirical analysis above, we propose several suggestions for power supply companies in various regions of China to promote future management.
Firstly, based on the analysis of probability density function and risk estimation, we can see that the power load of light industry is more likely to fluctuate than for heavy industry in China. Such fluctuation difference means that the power companies need more intensive and meticulous maintenance of the electricity load in the light industry, which is supported by all analyses including the risk analysis. However, since the power load consumed by the heavy industry is much higher than that of the light industry at the same time interval, the management of the heavy industry should be more biased towards safety checks of the pre-process and the normative constraint of the overall operation flow instead of increasing the times of maintenance of the electricity load.
Secondly, power enterprises can make load and customer management based on the analysis of load data. Power companies can increase the communication with their clients, plan the line redundancy, enlarge load redundancy of the existing heavy industry and eliminate the risk of equipment damage caused by capacity loads. Also, the electric power enterprise should incorporate the improvement of electricity utilization rate into management assessment, through the introduction of clean energy and distributed energy, and improve energy supply structure, to cope with the effects of the large power enterprises' sudden peak with least environment negative externalities.
Thirdly, existing power companies should pay more attention to the fluctuation characteristics of different industries to realize optimal management expenditure. As we demonstrated in Section 4, the enterprises of two distinct industries in this study have some similarities in fluctuation pattern, which makes it possible to use similar management strategies for maintaining related equipment. Moreover, the difference between the two enterprises' load peaks allows the power companies to design different maintenance teams to eliminate frequency to optimize the cost of management and maintenance.
Finally, an electric power company can make use of the correlation between small and large fluctuations, especially small fluctuations, and provide them as a value-added services to the existing large power customer to increase corporate earnings in this competitive market. For example, using the current fluctuation trends, the power enterprises can fully expect the changes of large fluctuations according to the characteristics of small fluctuations in the relevant industries. Based on pre-judgment, power companies can design different energy storage or line maintenance and support strategies, pack these strategies into different packages, and offer them to energy consumers in the form of value-added services. It can improve the efficiency of regional energy utilization, optimize local network and equipment management of electric circuitry, and also can bring extra income for the electric power company and improve the survivability of the electric power enterprises under the market reform.

Conclusions
The paper uses recurrence interval analysis to investigate the property of recurrence intervals of electricity fluctuations of heavy/low industries in China for different thresholds using 15-min high-frequency data, attempting to understand the behaviors of large volatilities in different industries.
First, we observe the distribution functions of recurrence and find the results that for different thresholds, the probability density function fits the stretched exponential function. We found that there is scaling behavior in the distributions for different thresholds after the recurrence intervals are scaled with the mean recurrence interval. Then we use conditional probability density function and MF-DFA to investigate the short-term, long-term memory and multi-fractural characteristics separately, indicating that the clusters of recurrence intervals of volatilities are not only caused by present conditions, but also long-term correlations. Later, we apply recurrence interval analysis to evaluate the risk of heavy/light industries, which provide relatively accurate estimates of hazard and forge a link between loss possibility and VaR. Finally, we propose management suggestions to Chinese energy supply companies based on all previous analyses.
Of course, this paper still has deficiencies. For example, other forecasting methods could be used or more enterprises from heavy/light industry could be chosen to test the robustness of research results. These deficiencies could be perfected by subsequent researchers to help energy company administrators evaluate risk.