Modeling of water usage by means of ARFIMA–GARCH processes

Accepted Manuscript Modelling of water usage by means of ARFIMA-GARCH processes Janusz Gajda, Grzegorz Bartnicki, Krzysztof Burnecki

PII: DOI: Reference:

S0378-4371(18)31077-X https://doi.org/10.1016/j.physa.2018.08.134 PHYSA 20020

To appear in:

Physica A

Received date : 30 July 2017 Revised date : 9 April 2018 Please cite this article as:, Modelling of water usage by means of ARFIMA-GARCH processes, Physica A (2018), https://doi.org/10.1016/j.physa.2018.08.134 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Highlights (for review)

Highlights    

A problem of modelling and prediction of hot water usage is addressed. Hot water usage data depicts antipersistent behavior and time-dependent variance. A simple identification and validation algorithm for ARFIMA-GARCH is introduced. The presented scheme can be also applied to other phenomena.

*Manuscript Click here to view linked References

Modelling of water usage by means of ARFIMA-GARCH processes Janusz Gajdaa , Grzegorz Bartnickib and Krzysztof Burneckia1 a

Faculty of Pure and Applied Mathematics, Hugo Steinhaus Center, Wroclaw University of Science and Technology, Wyspianskiego 27, 50-370 Wroclaw, Poland b Faculty of Environmental Engineering, Wroclaw University of Science and Technology, Wyspianskiego 27, 50-370 Wroclaw, Poland

Abstract This paper addresses an important problem of modeling and prediction of phenomena with antipersistent behavior and variance changing in time. As a proper stochastic model we propose an autoregressive fractionally integrated moving average (ARFIMA) process with generalized autoregressive conditional heteroskedasticity (GARCH) noise. First, we introduce a simple identification and validation algorithm for such model. Second, we apply the algorithm to weekday data of hot water usage at urban residential blocks. We extract the deterministic sinusoidal component from the data and fit successfully the ARFIMA-GARCH model to the stochastic part. The goodness of fit is checked by examining model errors and prediction performance. All analyses are performed by the rigorous statistical procedure. The proposed model allows for real-time accurate predictions and when implemented at a hot water supply level will lead to a better optimization of the control system and energy efficiency use. Keywords: ARFIMA model, GARCH model, Long-range dependence, Water usage 1. Introduction Derivatives and integrals of fractional orders can be used to describe random phenomena that can be characterized by long (power-like) memory or 1

Corresponding author: [email protected]

Preprint submitted to Physica A

April 9, 2018

self-similarity [1]. Long memory (or long-range dependence) is a property of certain stationary stochastic processes describing phenomena where the events that are distant still influence each other exceptionally strong. It has been associated historically with slow decay of correlations and a certain type of scaling that is connected to self-similar processes [2, 3]. A link between the power-law scaling of a long-range correlated series and the Shannon entropy was studied in Ref. [4]. The fractional Brownian motion (FBM) and autoregressive fractionally integrated moving average (ARFIMA) process serve as basic stochastic models for fractional dynamics [5]. The discrete-time ARFIMA with noise belonging to the domain of attraction of Gaussian law is stationary and when aggregated, in the limit, converges to fractional Brownian motion (FBM). In contrast to FBM, it allows for different light- and heavy-tailed distributions, and both long (power-like decay) and short (exponential decay) dependencies [5]. Moreover, as a stationary process, it provides prediction tools [6]. A systematic inference of the parameters of ARFIMA models has been studied in Ref. [7]. Power-like decay of the autocorrelation function leads to the notion of either persistence (positive dependence related to the classical definition of long range dependence) present, e.g., in finance [8, 9, 10], temperature records [11], river flows [12], DNA [13], heart-beat intervals [14], earthquakes [15] or telecommunication [16], or antipersistence (negative dependence) observed, e.g., in the biological anomalous diffusion data [17, 5]. The long memory property represented by ARFIMA processes can be combined with dependence in volatility (conditional variance). Bollerslev [18] generalized the autoregressive conditional heteroskedastic model, named ARCH, proposed in the seminal Nobel Prize winning paper by Engle [19]. The generalized autoregressive conditional heteroskedasticity (GARCH) process takes into account heteroskedastic effects of the time series process typically observed in form of clustering of volatilities [3]. Baillie et al. [20] introduced ARFIMA driven by GARCH noise, namely ARFIMA-GARCH processes. This concept was further rigorously studied by Ling and Li [21]. Such process serves as a basic model for data showing both long memory and short-range volatility effects [21, 6]. We also note that in Ref. [22] the authors modelled coupled long-range power-law correlated multivariate signals by a two-component ARFIMA stochastic process and a two-component fractionally integrated ARCH process. Water is considered as a most significant natural resource and a key factor of location and growth of communities. Rapidly growing urbanization makes 2

especially cities and other residential areas to be the main consumers of water resources. Thus planning of water distribution and forecasting of water demand plays a key role in water resource management especially now when many countries in the world is expected to struggle with water shortages [23, 24]. Probabilistic modeling and forecasting of water demand is now a well developing field of study. And due to the importance of the problem a proper usage and development of new techniques of water demand modeling plays a crucial role. Probably the most common methods adapt linear regressions, simple time series analysis [25, 26] and Box and Jenkins models [27]. Bakker et al. [28] compared three models: adaptive heuristic model, transfer model and multiple linear regression model with and without weather input. Recently the technique of artificial neural networks (ANN) is often used in the area of water demand modeling, for instance, Bougadis et al. [29] forecasted peak water demand with three different ANNs together with regression and seven time series models. Later Adamowski [30] considered similar problem and employed multiple linear regressions, autoregressive integrated moving average (ARIMA) processes and ANNs. A comprehensive review of literature in the context of ANNS one can find in [31] where authors consider a problem of forecasting of daily water demand for Al-Khobar city with coupled general regression neural networks and time series techniques. An interesting approach related to prediction of water demand was presented by Dotzauer [32], where the heat demand model was considered in relation to the outdoor temperature and the social behavior of the consumers. Based on communications theory, time series analysis seeks to identify regularities in the behavior of a system over time, trying to quantify cyclic or chaotic (strange) attractors [33]. An explanatory stochastic model for water and energy consumption based explicitly on time series ARIMA approach was considered in [34]. Another time series approach based on regression model was considered in [35] in the context of forecasting electricity load demand. We note that ARIMA or more general autoregressive fractionally integrated moving average (ARFIMA) models are frequently used in modeling of many real life phenomena which reflect the property of long memory (power-like decay of the autocorrelation function), see, e.g. [5] where telecommunication [16], astrophysical [36] and biological data [37] were studied, [38] for applications in climate science, and the monograph by Beran [6] and references therein. In the biophysical context ARFIMA processes with a negative memory 3

parameter are closely related to the notion of subdiffusion [5, 37]. ARFIMAGARCH models with variance changing in time can be useful in description of varying diffusion coefficient which results in so-called transient anomalous diffusion [39]. In Ref. [39] the authors showed that the receptor undergoes changes of diffusivity. This behavior was attributed to the occurrence of very long trapping events. The distribution of diffusion coefficient was considered to have a power-law behavior and the transit time (i.e., the time a particle moves with a given diffusion coefficient) was also considered to be random. Changing diffusivity is also related to ergodicity breaking. In Ref. [40] the ergodicity breaking was found to take place due to transient confinement effects since the molecules alternate between free diffusion and confined motion. ARFIMA combined with GARCH can describe both power-law decay of the autocorrelation function with finite-lag effects (ARFIMA part) and varying diffusion exponent (GARCH part). In this paper we introduce a simple identification and validation scheme for the ARFIMA-GARCH model. We note that the identification and prediction formulas can be easily implemented in any software. The scheme is illustrated on the hot water usage data from 2014 from two buildings with local gas boiler installation. We show that the data, after deseasonalization, depicts long memory (power-like) antipersistent behavior and high variability changing in time. In order to capture well these characteristics, the ARFIMA-GARCH process is identified and validated. The crucial advantage of this stationary model are ready-to-use prediction tools which allow for automated control of the heating equipment. It is also worth to mention that most articles in the literature deal with prediction of annual and monthly water demand for large areas like the cities and only a few address the problem of modeling and predictions of the daily water usage of single or several buildings. Therefore, this work provides new insights in terms of efficient modelling and prediction of the hot water consumption in individual buildings. 2. ARFIMA-GARCH identification, validation and prediction 2.1. ARFIMA model The ARFIMA time series models are fractionally differenced autoregressive moving average (ARMA) processes [41, 42, 43]. The model generalizes three broad classes of time series processes, namely the autoregressive (AR), the integrated (I), and the moving average (MA) [43]. In the finite second 4

moment case ARMA models posses the property of exponentially decaying correlation function. The ARFIMA generalization allows for a power-law correlation structure which leads to the notion of long memory dependence. Such an approach is widely used in modeling the dynamical behavior of various complex systems [6, 5]. The basic building blocks of the ARFIMA model are: AR(p) process X(n) = a1 X(n − 1) + a2 X(n − 2) + . . . + ap X(n − p) + ǫt and MA(q) process X(n) = ǫn − b1 ǫn−1 − b2 ǫn−2 − bq ǫn−q , where the sequence (ǫj ) is called the white noise process (sequence) [43]. We assume that ǫj ’s are independent and identically distributed (i.i.d.). In the literature the white noise term is often related to uncorrelated random variables but we don’t distinguish here between uncorrelation and independence. The noise sequence may have a Gaussian or non-Gaussian distribution. We only assume that the random variables belong to the domain of attraction of a Gaussian law. The process AR(p) is a causal or future-independent function of noise and stands for the regression. The explanatory variables are the observations prior to our current observation. For the MA(q) process X(n) it appears that X(s) and X(t) are independent whenever |t − s| > q. The dependence (correlation) is only q lag long. The time series process X(t) is an ARMA(p, q) process if it is stationary and satisfies (for every t) linear difference equation with constant coefficients: Φp (∇)∆d X(n) = Θq (∇)ǫn ,

,

(1)

where n = 0, ±1, . . . and ∇ stands for the backward shift operator, i.e. ∇(X(n)) = X(n − 1), and ∆ is the difference operator, i.e. ∆X(n) = X(n) − X(n − 1). In the above description the two polynomials Φ and Θ have the form Φp (∇) = 1 − a1 ∇ − a2 ∇2 − . . . − ap ∇p , Θq (∇) = 1 − b1 ∇ − b2 ∇2 − . . . − bq ∇q .

(2)

The former corresponds to the autoregressive (AR) part while the latter to the moving average (MA) part. ARMA models introduce short memory of the process since their autocorrelation function ρ at lag h converges exponentially to zero as h → ∞. They provide a general framework for studying stationary processes and simple prediction algorithms to describe future behavior of the dynamics. A generalization of this class, which incorporates a wide range of nonstationary 5

series, is provided by ARIMA processes, i.e. processes that reduce to ARMA processes when differenced finitely many times. Finally, stationary processes with much more slowly decreasing autocorrelation function, i.e. power-like are called ARFIMA processes. They were originally introduced by Granger and Joyeux [41] and Hosking [42]. The ARFIMA(p, d, q) time series process is defined as a stationary solution of equations Φp (∇)∆d X(n) = Θq (∇)ǫn , (3) where n = 0, ±1, . . .. The name fractional corresponds to the memory parameter d < 1/2 (or fractional differencing exponent). The fractional difference operator ∆d has the binomial expansion ∆d = (1 − B)d =

∞ X

πj (d)B j ,

(4)

j=0

where the coefficients πj (d)’s follows from the Taylor expansion of the function f (z) = (1 − z)d and can be represented by the gamma function Γ: πj (d) =

Γ(j − d) , Γ(−d)Γ(j + 1)

j = 0, 1, . . .

(5)

Moreover, if all roots of the polynomial Φp lie outside the unit circle, the ARFIMA(p, d, q) process defined by (3) is stationary and has a causal moving average (linear) form ∞ X X(n) = cj (d)ǫn−j , (6) j=1

where cj (d)’s are defined by the equation ∞

Θq (z)(1 − z)−d X = cj (d)z j , Φp (z) j=0

|z| < 1,

(7)

for details see [44]. Such processes are asymptotically self-similar with the Hurst parameter H = d + 1/2. The process X(n) is well-defined for −∞ < d < 1−1/2. The rate of decay of the autocovariance function ρn = EX(n)X(0)− EX(n)EX(0) for the ARFIMA model is n2d−1 . Moreover, for d > 0 we have 6

P∞

|ρn | = ∞. This serves as a classical definition of long memory [6]. However, in this paper we call all ARFIMA processes with d 6= 0 long-term dependent to distinguish between ARMA models with exponentially decaying autocorrelation function and ARFIMA models with the power-law decay. If d < 0 we arrive at negative (but still power-like) correlations which gives rise to the notion of antipersistence. We can estimate parameters of the ARFIMA process by using a Whittle estimator [44, 45]. We define the (p + q + 1)-dimensional vector n=0

β 0 = (a1 , a2 , . . . , ap , b1 , b2 , . . . , bq , d), where a1 , a2 , . . . , ap and b1 , b2 , . . . , bq are the coefficients of the polynomials Φp and Θq respectively. For a sample {X1 , X2 , . . . , XN } we denote the normalized periodogram by P 2 N −iλN t=1 Xt e I˜N (λ) = , −π ≤ λ ≤ π, (8) PN 2 t=1 Xt and let

γ(0) =

∞ X

c2j (d).

(9)

j=0

We also introduce the power transfer function 2 Θq (e−iλ , β) g(λ, β) = . −iλ −iλ d(β) Φp (e , β)(1 − e )

Finally, the estimator β N of the true parameter vector β 0 is defined as β N = arg minβ∈E

Z

π −π

I˜N (λ) dλ. g(λ, β)

(10)

The main advantage of this method is that we estimate all parameters together by simple minimization of one quantity. For the considered class of ARFIMA models the linear predictor based on the finite past Xn , . . . , X0 is defined as [44] n X ˆ n+k = X fj Xn−j , (11) j=0

7

where fj = − ′

and ξj (d) s are given by

k−1 X

cj (d)ξj+k−t(d),

t=0

∞

Φp (z)(1 − z)d X = ξj (d)z j , Θq (z) j=0

|z| < 1.

If the residuals of the fitted ARFIMA process seem to be uncorrelated with common variance, this can indicate correctness of the ARFIMA model. If positive autocorrelations in the squared and absolute values are observed, then classical formulations (such as ARFIMA models) are inappropriate. The fact that large absolute returns tend to be followed by large absolute returns (whatever the sign of the variations) is hardly compatible with the assumption of constant conditional variance. This phenomenon is called conditional heteroscedasticity. 2.2. GARCH model Time series models with a time-varying conditional variance were first proposed by Engle [19]. A generalized autoregressive conditional heteroscedasticity (GARCH(k, l)) model is defined by [18]: p (12) ǫt = ht ηt , ht = α0 +

k X

αi ǫ2t−i +

i=1

l X

βj ht−j ,

(13)

j=1

where ηt is a sequence of i.i.d. random variables. The typical examples of the noise are normal and student’s t distributions. In most applications it is enough to consider the order of the model as (1, 1) [46], so k = 1 and l = 1. For this process, in order to estimate its parameters, one can use the maximum likelihood method or follow again the Whittle approach applying the fact that a GARCH(k,l) process (Eqs. (12), (13)) can be rewritten in the ARMA(k,g) form, where g = max(k, l), in the following way ǫ2t −

g X i=1

ζi ǫ2t−i = α0 + νt − 8

l X j=1

βi νt−i ,

(14)

where ζi = αi + βi . If i ∈ (l, k] set βi = 0, and if i ∈ (k, l], set αi = 0 and
νt = ǫ2t − ht = ht ηt2 − 1 .

Thus, we can estimate GARCH parameters by the same technique as for ARMA (a special case of ARFIMA) processes. Under the assumption that ht is a strictly stationary process and V ar(ǫ21 ) < ∞ (which is the considered here case) the sequence νt constitutes a white noise process. Observe that this procedure is simple and robust as before, since we just maximize periodogram function with respect to the transfer function of the ARMA model. We finally note that the procedure fails for a class of stable GARCH processes [47]. 2.3. ARFIMA-GARCH model If we combine ARFIMA with GARCH processes by assuming the noise of the former is modelled by the latter we obtain a powerful composite ARFIMA(p, d, q)-GARCH(k, l) model that takes into account short and long memory, and the variance changing randomly in time [20, 21]. This process will serve as a basic model for the hot water usage data. We estimate parameters of this composite process sequentially, namely first we estimate the parameters of the ARFIMA model, then we extract its residuals for which we estimate GARCH parameters. This two-step procedure is similar to the one presented in Ref. [43] (see Example 10.3.3 therein) and implemented in the package ITSM [43]. Another approach is to estimate ARFIMA-GARCH parameters simultaneously [20, 21] but this is much more computationally expensive, which we try to avoid in the context of automation software. We also believe that the introduced methodology is more robust than the simultaneous estimation of the parameters, which is implemented, e.g., in the R package ’rugarch’ (https://cran.r-project.org/package=rugarch), where convergence problems can arise. Prediction of ARFIMA-GARCH processes is the same as prediction of ARFIMA processes, since formula (11) does not depend on the noise sequence. To predict one observation ahead for the GARCH(1,1), we simply write ˆ t+1 = α h ˆ0 + α ˆ 1 ǫ2t + βˆ1 ht , where we use estimated parameters α ˆ0, α ˆ 1 , βˆ1 and the last values in the series 2 ǫt and conditional variance ht .

9

2.4. Identification and validation algorithm In order to fit the ARFIMA-GARCH model to real empirical data, we propose the following procedure, which is illustrated in Fig. 1.

Figure 1: Identification and validation algorithm of ARFIMA-GARCH processes.

ARFIMA-GARCH identification and validation algorithm S1 If the data depict non-stationary behavior by a visible trend or seasonality, remove them from the data, see, e.g., [43]. S2 Having stationary data, check for the property of long memory. It can be done by means of the mean-squared displacement (MSD). Its linear form (at + b) in the double logarithmic scale with slope a of the line not equal to 1 suggests the long range (power-like) correlations. The MSD also provides a rough estimate of the long memory parameter d = (a − 1)/2 [48, 5] (see also Section 3). 10

S3 Estimate the order of the ARFIMA model, namely p and q parameters. To this end fractionally difference the data to obtain a standard ARMA(p, q) process [16, 49]. Next, for the ARMA process apply common model order criteria, like the bias-corrected Akaike [43]. To estimate the order of the model, one can also use visual techniques based on autocorrelation and partial autocorrelation functions. Finally, we note that since in many applications both p and q are not grater than 2, one can just fix p = 2 and q = 2 without a thorough investigation. S4 Given the order, fit the ARFIMA parameters by means of the Whittle estimator given by formula (10). S5 Check if residuals (one-step prediction errors) form a white noise, i.e. a sequence of i.i.d. random variables. This can done by a visual inspection of autocorrelation functions (ACF)’s of both residuals and squared residuals or/and by rigorous statistical tests like the LjungBox or McLeod-Li [48] (see also Section 3). If the residuals resemble i.i.d. sequence, proceed to Step S9. S6 If the residuals fail the visual or statistical tests and they show correlations in squared values, fit the GARCH process to residuals by applying the Whittle estimator to the squared values, cf. formula (14). S7 Calculate the residuals of the GARCH model as one-step prediction errors. S8 Check if residuals form a white noise. This can done analogously as in Step S5. If the residuals resemble an i.i.d. sequence, proceed to the next step. S9 Fit the distribution to the residuals (Gaussian, student’s t, normal inverse Gaussian, other) and check the goodness of fit by visual comparison of the tails of the empirical and fitted analytical distribution functions or/and by rigorous statistical tests [48] (see also Section 3). If the above steps fail for the analyzed data, the ARFIMA-GARCH process should be excluded as a potential model. For the fitted ARFIMAGARCH model on can calculate forecasts and prediction intervals. In order to illustrate these quantities on the original data one should remember to invert the transformations from Step S1. 3. Investigated system The analyzed system consists of 2 buildings built and settled in the years 2006-2008. In total, 162 premises are located in the buildings. The heat is 11

produced in a local gas boiler installation and is supplied for heating and domestic hot water purposes. All apartments are equipped with a heat consumption measuring system (heat meter with wing flow meter) and wing water meter for the hot tap water consumption measurement. The measuring system is complemented by the hot water consumption measurement located in the gas boiler room. Water meter with pulse transmitter is mounted on the supply side of the hot water cylinder. Indications of the water meter are recorded at 15 seconds intervals, and pulse transmitter has a frequency of 1dm3 /impulse. In Figs. 2 and 3 we present illustrations of the local gas boiler installation.

Figure 2: Local gas boiler installation – real photo.

We encountered several issues with the recorded data. First, the data contain missing values which have to carefully handled. Second, they contain holidays which should be treated separately. Since our investigation of the data revealed different daily and monthly deterministic patterns for the purpose of analyzes in the paper we take into consideration water usage for one day of a month, namely Wednesdays in June, 2014. We selected one day of month which is free of missing data and holidays and treat it is a representative sample. As a preliminary step we aggregated the recorded data on 15-minute intervals in order to smooth possible short-time variations. The data are depicted in Fig. 4. Notice that high water usage values correspond to mornings and evenings. We will now follow the steps of the identification and validation algorithm presented in Section 2.4. 12

Figure 3: Local gas boiler installation – schematic picture. 4.5 Data Deterministic component

Water usage in June 2014

4 3.5 3 2.5 2 1.5 1 0.5 0 0AM

12:0

4) (06-0

1) 9PM (06-1 11:5 0AM 12:0

8) 9PM (06-1 0AM 12:0

11:5

5) 9PM (06-2 0AM 12:0

11:5

PM

9 11:5

Figure 4: Water usage empirical data with fitted deterministic part given by sum of two sinuses with 24 and 12 hours periods.

Step S1. We can notice a visible sinusoidal behavior of the series. Since two periods were determined (24 and 12 hours), we subtracted the following sinusoidal function from the data (see the red solid line in Fig. 4): p1 + p2 sin(2π/d1 (t + p3 )) + p4 sin(2π/d2 (t + p5 )), where p1 = 1.1, p2 = −0.93, p3 = 1.75, p4 = −0.41, p5 = −13.42 and di , i = 13

2.5

Water usage in June 2014 after trend removal

2 1.5 1 0.5 0 -0.5 -1 -1.5 0AM

12:0

4) (06-0

1) 9PM (06-1 11:5 0AM 12:0

8) 9PM (06-1 0AM 12:0

11:5

5) 9PM (06-2 0AM 12:0

11:5

PM

9 11:5

Figure 5: Stochastic part of the considered data, i.e. observations with removed deterministic part.

1, 2 are period lengths so d1 = 12 and d2 = 24 hours. As a result we obtain the series presented in Fig. 5. For the deseasonalized data we plot their autocorrelation function (ACF) in Fig. 6. 1 0.8

Autocorrelation

0.6 0.4 0.2 0 -0.2 -0.4 0

10

20

30

40

50

Lag

Figure 6: Autocorrelation function of the stochastic component of the data with confidence bounds.

14

The ACF at lag k is defined as

rk =

NP −k t=1

(Xt − X)(Xt+k − X) N P

t=1

,

for k = 0, 1, . . .

(15)

(Xt − X)

and serves as a measure of the linear dependence between observations with time lag k [43]. For example, if we consider i.i.d. data, then rk = 0 for all k 6= 0. When we investigate moving average processes of order q (MA(q)), the autocorrelation function is zero for lags beyond q, thus it can be useful for estimating the q parameter. As we can notice in Fig. 6, there is a significant dependence in the data which suggests to model them with fractionally differenced processes. Step S2. In order to confirm appropriateness of a long memory model, we estimate the memory parameter d by using the sample (time-averaged) meansquared displacement (MSD) statistics [50]. Consider a discrete time series {Yn , n = 1, . . . , N} of length N. Then, the sample (time-averaged) MSD can be calculated as N −τ 1 X MN (τ ) = (Yk+τ − Yk )2 . N − τ k=1

(16)

If N becomes large and τ small, the sample comes from an aggregate ARFIMA process with noise belonging to the domain of attraction of Gaussian law, then MN (τ ) ∼ τ 2d+1 [50] (distribution of the sample MSD of a general discrete Gaussian process was studied in Ref. [51]). We present the results in Fig. 7, and since d is essentially different from 0 we confirm that our series possesses long memory. Steps S3-S4. We set the order of the ARFIMA model to (2, 2) and evaluate ARFIMA parameters with the use of Whittle’s estimator (see Section 2). The orders of p and q were selected by sequentially analyzing information criteria for the models of orders p = 0, 1, 2 and q = 0, 1, 2. The final model is chosen as the one which minimizes the Akaike criterion [43]. We obtain the following ARFIMA model ∆−0.3 (X(t) − 1.43X(t − 1) + 0.54X(t − 2)) = ǫt − 0.87ǫt−1 + 0.31ǫt−2 . (17) 15

0.4 0.2

ln(MSD(t))

0 −0.2 −0.4 −0.6 −0.8 −1 −1.2 0

0.5

1

1.5

2

2.5

ln(t)

Figure 7: MSD statistic (blue circles) in the log-log scale calculated for the stochastic component of the data. The estimated slope of the fitted black thick line corresponds to d = −0.3, which suggests long memory (power-like) antipersistent behavior. The red dashed lines correspond to the confidence interval for the MSD statistic calculated by means of Monte Carlo simulations from the fitted ARFIMA-GARCH model. The empirical MSD fits well into the confidence region.

In order to check the goodness of fit of the model we calculate residuals (ǫj ) as one-step ahead prediction errors. Steps S5-S6. For the residuals we calculate the ACF for both raw and squared values. We observe in Fig. 8 that for squared values there is a substantial dependence at lags 1, 2 and 3. In addition, we performed the Ljung-Box and McLeod-Li tests which take into account magnitude of correlation as a group. The Ljung-Box statistic for a given maximum lag K is given by [52]:  2  ρ1 ρ22 ρ2K Q = N(N + 2) , (18) + + ...+ N −1 N −2 N −K where ρi is a squared autocorrelation at lag i and K is a number of autocorrelation lags included in the statistic. The McLeod-Li statistic is defined as the Ljung-Box but calculated on the squared values of the process. The Q statistic is χ2 distributed as N → ∞. This allows to test the null hypothesis that there is no serial correlation in the data. It appears that the Ljung-Box test does not reject the randomness hypothesis (for K = 22 the p-value equals 0.26), whereas the McLeod-Li strongly rejects it (for K = 22 the p-value equals 0). 16

We decide to apply the GARCH model for ǫj ’s to capture this phenomenon. Thus, we fit the model to the noise by means of Whittle’s estimator (see Section 2). The fitted GARCH process is given by the following formula p ǫt = ht ηt , (19)

where

ht = 0.23 + 0.12ǫ2t−1 + 0.74ht−1 .

(20)

1 0.5 0 Lag

-0.5

Autocorrelation squares

Autocorrelation

Steps S7-S8. We extract residuals for this model and calculate, as previously, ACF for both raw and squared data. The results are presented in Fig. 9 and they confirm the lack of dependence between values of residuals. Moreover, both the Ljung-Box and McLeod-Li tests do not reject the randomness hypothesis (for K = 22 the p-value of the former test equals 0.16 and for the latter is 0.78). From these results we infer that there is no serial correlation in the data. In sum, we successfully fitted ARFIMA-GARCH model, which is given by Eqs. (17) and (20).

1

0

10

20

30

40

50

0

10

20

30

40

50

0.5 0 Lag

-0.5

Figure 8: Autocorrelation function of residuals of the fitted ARFIMA model. We can clearly observe that residuals exhibit significant dependence especially when looking at the autocorrelation of squared values at the first three lags. This fact suggests fitting GARCH model to residuals.

Step S9. We now investigate the distribution of the residuals. Following [53] we first checked if the underlying distribution belongs to the domain of attraction of the Gaussian or non-Gaussian stable distribution by examining its rate of convergence. Since the results suggested the Gaussian domain of 17

0.5 0 Lag

-0.5

Autocorrelation squares

Autocorrelation

1

1

0

10

20

30

40

50

0

10

20

30

40

50

0.5 0 Lag

-0.5

Figure 9: Autocorrelation function of residuals of the fitted ARFIMA-GARCH model. We notice that residuals of the model exhibit small serial correlations and they can be considered as a white noise.

attraction we consider three typical light-tailed probability distributions for the residuals of the ARFIMA-GARCH model, namely Gaussian, t locationscale and normal inverse Gaussian (NIG). The t location-scale distribution generalizes ordinary student’s t distribution. The probability density function (PDF) of this distribution is given as ft (x) =

− 12

n
σB n2 , 21

2

1+

((x − µ)/σ) n

!− n+1 2

,

where B (·, ·) is the beta function and n denotes the degree of freedom parameter. The distribution is useful for modeling data with heavier tails than the normal. It approaches the Gaussian distribution as n approaches infinity, and smaller values of n yield heavier tails. A random variable X is said to have a NIG distribution if it has PDF p αδ δ√α2 −β 2 +β(x−µ) K1 (α δ 2 + (x − µ)2 ) p , −∞ < x < ∞. (21) e fN IG (x) = π δ 2 + (x − µ)2 The NIG distribution, introduced in Ref. [54], is described by four parameters (α, β, δ, µ), where α stands for tail heaviness, β for asymmetry, δ is the scale parameter, and µ is the location. The normalizing constant Kλ (t) in (21) is the modified Bessel function of the third kind with index λ, also known 18

Table 1: Values of the estimated parameters of the Gaussian, t location-scale and NIG laws.

Gaussian µ=0

σ = 0.47 t location-scale µ=0 σ = 0.39 n = 6.35 NIG α = 3.65 β = 1.67 δ = 0.59 µ = −0.31 as the MacDonald function. The NIG distribution is more flexible than the Gaussian since it allows for fat-tails and skewness. The normal distribution arises as a special case by setting β = 0, δ = σ 2 α, and letting α → ∞. The fitted parameters of the distributions are presented in Table 1. In order to check goodness of fit for the three distributions considered here we follow the testing procedure from [48]. We use two classes of goodness of fit statistics, namely the Kolmogorov-Smirnov and Cramer-von Mises [55, 56]. The Kolmogorov-Smirnov statistics D is defined as D = sup |FN (x) − F (x)|.

(22)

x

Moreover, it can be written equivalently as a maximum of two non-negative supremum statistics D + and D − : D + = sup{FN (x) − F (x)} and D − = sup{F (x) − FN (x)}. x

(23)

x

Within the same class we consider also the Kuiper statistic given by V = D+ + D− .

(24)

For the Cramer-von Mises family of statistics CM we have the following formula Z ∞ CM = N (FN (x) − F (x))2 ψ(x)dF (x). (25) −∞

Here ψ(x) is a weight function. In particular for ψ(x) = 1 we obtain the W 2 Cram´er-von Mises statistic and for ψ(x) = [F (x)(1 − F (x))]−1 we arrive at the A2 Anderson and Darling. 19

Table 2: Values of the test statistics and corresponding p-values (in parentheses) under the assumption of the Gaussian, t location-scale and NIG laws.

V W2 A2 Gaussian 0.08 0.12 0.69 3.67 (0) (0) (0) (0) t location-scale 0.05 0.1 0.27 0 (0) (0) (0) (0) NIG 0.03 0.06 0.05 0.35 (0.24) (0.11) (0.19) (0.12) D

The results of the testing procedure are depicted in Table 2. The clear winner is the NIG distribution with all p-values above 0.1. We also calculated the confidence interval for the MSD of the fitted ARFIMA-GARCH process and checked that the empirical MSD falls into the confidence region, see Fig. 7. This justifies the proposed model. Finally, to confirm the goodness of fit of the identified model we perform the procedure of backtesting. It consists in assessing the accuracy of fitted models predictors using existing historical data. The backtesting process starts by selecting a threshold date within a time span covered by the historical data. Then, for the threshold, the historical data are truncated at the threshold, the model is fitted to the truncated data (training dataset) and the original untruncated data (testing dataset) are compared with the forecasts calculated from the model and a prediction interval. Prediction intervals can be constructed from the noise term of the model or nonparametrically from the prediction errors. The results for both detrended and original data are presented in Fig. 10 where the testing dataset was chosen as the first Wednesday in June 2015. All forecasts are one-step-ahead. The prediction interval was calculated from the empirical distribution of one-step-ahead prediction errors by computing the appropriate sample quantiles [57]. We can see that roughly 80% of the true observations lie within the prediction bounds and the predicted values are close to them. We now concentrate on the original data (with added trend). To exam20

ine the point forecasting performance, we calculated the root mean square error (RMSE) for the ARFIMA-GARCH model. We found that RMSE is 0.48, which is acceptable. In order to check the interval forecasting quality we performed the Christoffersen tests of coverage and independence [58]. For the 80% prediction interval, p-values for the unconditional coverage, independence and conditional coverage likelihood ratio tests were (0.22, 0.07, 0.12), so the null hypotheses were clearly retained at the 5% significance level. We have also checked the out-of-sample performance for special cases of the ARFIMA-GARCH model, namely ARMA-GARCH (without the long memory part) and AR (no long memory and GARCH parts) processes. The obtained RMSEs were 0.57 and 0.53 respectively, which are higher than for the ARFIMA-GARCH model. The unconditional coverage, independence and conditional coverage likelihood ratio tests for ARMA-GARCH and AR processes provided the following p-values: (0.04, 0.07, 0.08) and (0.03, 0.11, 0.04), respectively. This confirms the choice of the fitted ARFIMA-GARCH model.

2 1 0 -1 3

Water usage, data without trend [m ] prediction interval

Empirical Data

ARFIMA-GARCH forecast

4 3 2 1 0

3

Water usage, data with trend [m ] 12:00AM

6:00AM

12:00PM

6:00PM

11:59PM

Figure 10: Water usage data (blue line) with one-step prediction (red line) and 80% prediction intervals (shaded area) of the fitted ARFIMA-GARCH model for the data with deterministic trend removed (top panel) and for original data (bottom panel). We can see that roughly 80% true observations lie between the prediction bounds which confirms the goodness of fit of the model.

To sum up, we conclude that the hot water usage data for Wednesdays in June can be well described by the ARFIMA-GARCH model. Simpler 21

structures like AR or ARMA-GARCH provide worst fits both in terms of in-sample fits and out-of-sample performance. We have also performed preliminary analyzes for different days for different months and we found that the deterministic pattern of the dynamic changes but the underlying stochastic dynamic is of the same type, namely ARFIMA-GARCH. This also agrees with the intuition that the stochastic pattern should be similar for different days only with the model parameters changing. Therefore, mathematically, the overall process can be treated as a piecewise ARFIMA-GARCH (each day of a week with different parameters). 4. Conclusions In this paper we considered an important issue of modeling and prediction phenomena which posses the property of long-range (power-like) dependence and randomly changing variance. As an appropriate model, we proposed the ARFIMA-GARCH process. For the general ARFIMA-GARCH model described in Section 2 we introduced an easy-to-use identification and validation algorithm. We notice that the identification part is self-adaptive [33], namely having fixed the order of the model, ARFIMA-GARCH parameters can be fitted periodically from the historical data with a fixed time window length. Estimation and prediction formulas are given in an elementary form and can be easily implemented in any software. This two-step estimation procedure is similar to the one used in the package ITSM [43]. We are aware of the fact that the introduced identification algorithm is not as accurate as the simultaneous estimation of all ARFIMA-GARCH parameters presented in Refs. [20, 21] and implemented, e.g., in the R package ’rugarch’ (https://cran.r-project.org/package=rugarch) but we believe it is simpler and more robust without possible convergence issues. As a result it allows for reliable real-time predictions. We applied the introduced methodology to hot water usage data collected at urban residential blocks. We extracted the deterministic sinusoidal component from the weekday data of a Wroclaw residential block from June 2014 and showed that the data, after deseasonalization, depicts strong antipersistence and high variability. As a consequence, we proposed an ARFIMAGARCH process as a proper stochastic time series model which captures well the discovered characteristics and as a stationary process enables an efficient prediction of future values. We fitted successfully the ARFIMA(2, 2)GARCH(1, 1) model to the stochastic part of the data. The goodness of 22

fit was checked by examining model errors and prediction performance. All analyses were performed with the use of rigorous statistical tools. We also checked that simpler structures like ARMA-GARCH (no long memory part) and AR (no long memory and GARCH parts) provide worst results. In the paper, we treat Wednesdays as a representative sample. Our preliminary analyzes have shown that different days (within different months) differ with respect to the deterministic pattern but the underlying stochastic dynamic is of the same type, namely ARFIMA-GARCH. Therefore, mathematically, the overall process of hot water usage data can be treated as a piecewise ARFIMA-GARCH. We note that an effective forecasting of the hot water consumption in short time intervals would enable to develop control and adjust algorithms, which would optimize fuel consumption. The implementation of new control algorithms based on the prediction of hot water demand would enable elimination of the unnecessary hot water preparation system launching, reduction of heat losses from hot water tanks and reduction of fuel consumption. Currently used control algorithms are designed to keep the set parameter (in this case, the hot water temperature) within the limits specified by hysteresis of temperature set point. While there is no analysis if the demand for hot water will occur in the nearest future; as a result the hot water preparation system keeps being turned on unnecessarily. This creates heat losses reducing energy efficiency of the heat supply system. Finally, we note that the model studied here is closely related to the notion of transient diffusion dynamics observed in biological experiments which manifest through randomly varying diffusivity [39]. This behavior is attributed to the occurrence of very long trapping events. In Ref. [40] the ergodicity breaking was found to take place due to transient confinement effects since the molecules alternate between free diffusion and confined motion. Hence, we also believe that the introduced methodology can be applied to a large class of anomalous diffusion data both power-law decay of the autocorrelation function with finite-lag effects (ARFIMA part) and varying diffusion exponent (GARCH part). Acknowledgments K.B. and J.G. gratefully acknowledge the financial support from the NCN Maestro Grant No. 2012/06/A/ST1/00258.

23

References [1] Klafter J, Lim S and Metzler R, 2012 Fractional Dynamics: Recent Advances. Singapore: World Scientific; 2012. [2] Beran J. Statistics for Long-memory Processes. New York: Chapman & Hall; 1994. [3] Mantegna RN, Stanley HE. An Introduction to Econophysics: Correlations and Complexity in Finance. Cambridge, UK: Cambridge University Press; 2000. [4] Carbone A, Stanley HE. Scaling properties and entropy of longrange correlated time series. Physica A. 2007; 384: 21-24. doi.org/10.1016/j.physa.2007.04.105. [5] Burnecki K, Weron A. Algorithms for testing of fractional dynamics: a practical guide to ARFIMA modelling. J. Stat. Mech. 2014; P10036. doi:10.1088/1742-5468/2014/10/P10036. [6] Beran J, Feng Y, Ghosh S and Kulik R. Long-Memory Processes. Probabilistic Properties and Statistical Methods. Berlin: Springer-Verlag; 2013. [7] Graves T, Franzke CLE, Watkins NW, Gramacy RB, Tindale E, Systematic inference of the long-range dependence and heavy-tail distribution parameters of ARFIMA models. Physica A. 2017; 473: 60-71. doi.org/10.1016/j.physa.2017.01.028. [8] Liu Y, Cizeau P, Meyer M, Peng CK, Stanley HE. 1997 Correlations in economic time series. Physica A. 1997; 245: 437–440. doi.org/10.1016/S0378-4371(97)00368-3 [9] Weron R. Estimating long-range dependence: finite sample properties and confidence intervals Physica A. 2002; 312: 285-299. doi.org/10.1016/S0378-4371(02)00961-5. [10] Muchnik L, Bunde A, Havlin S. Long term memory in extreme returns of financial time series. Physica A. 2009; 388: 4145-4150. doi.org/10.1016/j.physa.2009.05.046.

24

[11] Koscielny-Bunde Eva, Bunde A, Havlin S, Roman HE, Goldreich Y, Schellnhuber HJ Indication of a universal persistence law governing atmospheric variability. Phys. Rev. Lett. 1998; 81: 729–732. doi.org/10.1103/PhysRevLett.81.729. [12] Eichner JF, Kantelhardt JW, Bunde A, Havlin S. Extreme value statistics in records with long-term persistence. Phys. Rev. E. 2006; 73: 016130. doi:10.1103/PhysRevE.73.016130. [13] Peng CK, Buldyrev SV, Goldberger AL, Havlin S, Sciortino F, Simons M, Stanley HE. Long-range correlations in nucleotide sequences. Nature. 1992; 356: 168–170. doi:10.1038/356168a0. [14] Peng CK, Mietus J, Hausdorff JM, Havlin S, Stanley HE, Goldberger AL. Long-range anticorrelations and non-Gaussian behavior of the heartbeat. Phys. Rev. Lett. 1993; 70 1343–1346. doi:10.1103/PhysRevLett.70.1343. [15] Livina VN, Havlin S, Bunde A Memory in the occurrence of earthquakes. Phys. Rev. Lett. 2005; 95: 208501. doi: 10.1103/PhysRevLett.95.208501. [16] Burnecki K, Sikora G. Identification and validation of stable ARFIMA processes with application to UMTS data. Chaos, Solitons & Fractals. 2017; 102: 456-466. doi.org/10.1016/j.chaos.2017.03.059. [17] Metzler R, Jeon J-H, Cherstvy AG, Barkai E. Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking. Phys. Chem. Chem. Phys. 2014; 16: 24128. doi: 10.1039/C4CP03465A. [18] Bollerslev T. Generalized autoregressive conditional heteroskedasticity. J. Econometrics. 1986; 31: 307-327. doi.org/10.1016/03044076(86)90063-1. [19] Engle RF. Autoregressive Conditional Heteroscedasticity with Estimates of the Variance of United Kingdom Inflation. Econometrica. 1982; 50: 987-1007. doi: 10.2307/1912773.

25

[20] Baillie R, Chung, C-F, Tieslau, MA. Analysing inflation by the fractionally integrated ARFIMA-GARCH model. J. Appl. Econometrics 1996; 11:23-40. [21] Ling S, Li WK. On fractionally integrated autoregressive moving-average time series models with conditional heteroscedasticity. J. Amer. Statist. Assoc. 1997; 92: 1184-1194. doi:10.2307/2965585. [22] Podobnik B, Horvatic D, Lam Ng A, Stanley HE, Ivanov PCh. Modeling long-range cross-correlations in two-component ARFIMA and FIGARCH processes. Physica A. 2008; 387: 3954-3959. doi.org/10.1016/j.physa.2008.01.062. [23] Fan L, Liu G, Wang F, Ritsema CJ, Gissen V. Domestic water consumption under intermittent and continuous modes of water supply. Water Resour. Manag. 2014; 28: 853-865. doi: 10.1007/s11269-014-0520-7. [24] Gain AK and Wada Y. Assessment of future water scarcity at different spatial and temporal scales of the Brahmaputra River Basin. Water Resour. Manag. 2014; 28: 999-1012. doi 10.1007/s11269-014-0530-5 [25] Jain A, Varshney K, Joshi U. Short-term water demand forecast modeling at IIT Kanpur using artificial neural networks. Water Resour. Manag. 2001; 15: 299–321. doi: 10.1023/A:1014415503476. [26] Gato S, Jayasuriya N, Roberts P. Forecasting residential water demand: Case study. J. Water Resour. Plan. Manag. 2007; 133: 309-319. dx.doi.org/10.1061/(ASCE)07339496(2007)133:4(309)sthash.eOCmbDnT.dpuf. [27] Maidment D, Miaou S, Crawford M. Transfer function models of daily urban water use. Water Resour. Res. 1985; 21: 425-432. doi: 10.1029/WR021i004p00425. [28] Bakker M, van Duist H, van Schagen K, Vreeburg J, Rietveld L. Improving the performance of water demand forecasting models by using weather input. Procedia Eng. 2014; 70: 93-102. doi.org/10.1016/j.proeng.2014.02.012.

26

[29] Bougadis J, Adamowski K, Diduch R. Short-term municipal water demand forecasting. Hydrol. Process. 2005; 19: 137-148. doi: 10.1002/hyp.5763. [30] Adamowski J. Peak daily water demand forecast modeling using artificial neural networks. J. Water Resour. Plan. Manag. 2008; 134: 119-128. dx.doi.org/10.1061/(ASCE)0733-9496(2008)134:2(119). [31] Al-Zagarini MA, Abo-Monsar A. Urban residential water demand prediction based on artificial neural networks and time series models. Water Resour. Manag. 2015; 29: 3651-3662. doi: 10.1007/s11269-015-1021-z. [32] Dotzauer E. Simple model for prediction of loads in district-heating systems. Appl. Energy. 2002; 73: 277-284. doi.org/10.1016/S03062619(02)00078-8. [33] Mac´ıas-Escriv´a FD, Haber R, del Toro R, Hernandez V. Selfadaptive systems: A survey of current approaches, research challenges and applications. Expert Syst. Appl. 2013; 40: 7267-7279. doi.org/10.1016/j.eswa.2013.07.033. [34] Lowry G, Bianeyin FU, Shah N. Seasonal autoregressive modelling of water and fuel consumptions in buildings. Appl. Energy. 2007; 84: 542552. doi.org/10.1016/j.apenergy.2006.02.003. [35] Vaghefi A, Jafari MA, Bisse E, Lu Y, Brouwer J. Modeling and forecasting of cooling and electricity load demand. Appl. Energy. 2014; 136: 186-196. doi.org/10.1016/j.apenergy.2014.09.004. [36] Stanislavsky A, Burnecki K, Magdziarz M, Weron A. FARIMA Modeling of solar flare activity from empirical time series of soft X-ray solar emission. Astrophys. J. 2009; 693: 1877-1882. doi:10.1088/0004637X/693/2/1877. [37] Sikora G, Kepten, R, Weron A, Balcerek M, Burnecki K. An efficient algorithm for extracting the magnitude of the measurement error for fractional dynamics. Physical Chemistry Chemical Physics. 2017; 19: 26566-26581.

27

[38] Graves T, Gramacy RB, Franzke CL, Watkins NW. Efficient Bayesian inference for natural time series using ARFIMA processes. Nonlin. Processes Geophys. 2015; 22: 679-700. doi:10.5194/npg-22-679-2015. [39] Manzo C, Torreno-Pina JA, Massignan P, Lapeyre, Jr., GJ, Lewenstein M, Garcia Parajo, MF. Weak Ergodicity breaking of receptor motion in living cells stemming from random diffusivity. Phys. Rev. X. 2015; 5: 011021. doi: 10.1103/PhysRevX.5.011021. [40] Weron A, Burnecki K, Akin EJ, Sole L, Balcerek M, Tamkun MM, Krapf D. Ergodicity breaking on the neuronal surface emerges from random switching between diffusive states. Scientific Reports 2017; 7: 5404. [41] Granger CWJ, Joyeux R. An introduction to long memory time series models and fractional differencing. J. Time Ser. Anal. 1980; 1: 15-29. doi: 10.1111/j.1467-9892.1980.tb00297.x. [42] Hosking JRM. Fractional differencing. Biometrika. 1981; 68: 165-176. doi.org/10.1093/biomet/68.1.165 [43] Brockwell PJ, Davis RA. Introduction to Time Series and Forecasting. New York: Springer; 2002. [44] Kokoszka PS. Prediction of infinite variance fractional ARIMA. Probab. Math. Statist. 1996; 16: 65-83. [45] Burnecki K, Sikora G. Estimation of FARIMA parameters in the case of negative memory and stable noise. IEEE Trans. Signal Process. 2013; 11: 2825-2835. [46] Lv X, Shan X. Modeling natural gas market volatility using GARCH with different distributions. Physica A. 2013; 392: 5685-5699. doi.org/10.1016/j.physa.2013.07.038. [47] Mikosh T, Straumann D. Whittle estimation in a heavy-tailed GARCH(1,1) model. Stoch. Proc. Appl. 2002; 100: 187-222. doi.org/10.1016/S0304-4149(02)00097-2. [48] Burnecki K, Gajda J, Sikora G. Stability and lack of memory of the returns of the Hang Seng index. Physica A. 2011; 390: 3136-3146. doi.org/10.1016/j.physa.2011.04.025 28

[49] Xiu J, Jin Yao. Empirical study of ARFIMA model based on fractional differencing. Physica A. 2007; 377: 138-154. doi.org/10.1016/j.physa.2006.11.030. [50] Burnecki K, Weron A. Fractional L´evy stable motion can model subdiffusive dynamics. Phys. Rev. E. 2010; 82: 021130. doi:10.1103/PhysRevE.82.021130. [51] Grebenkov DS, Probability distribution of the time-averaged meansquare displacement of a Gaussian process. Phys. Rev. E 2011; 84: 031124. doi.org/10.1103/PhysRevE.84.031124. [52] Ljung GM, Box GEP. On a measure of lack of fit in time series models. Biometrika. 1978; 65: 297-303. doi.org/10.1093/biomet/65.2.297. [53] Burnecki K, Wyloma´ nska A, Chechkin A. Discriminating between lightand heavy-tailed distributions with limit theorem. PLoS ONE. 2015; 10: e0145604 doi.org/10.1371/journal.pone.0145604. [54] Barndorff-Nielsen OE. Normal Inverse Gaussian Processes and the Modelling of Stock Returns. Research Report 300. Department of Theoretical Statistics, University of Aarhus; 1995. [55] D’Agostino RB, Stephens MA. Goodness-of-Fit Techniques. New York: Marcel Dekker; 1986. [56] Burnecki K, Janczura J, Weron R. Building loss models. In Cizek P, H¨ardle W, Weron R, editors. Statistical Tools for Finance and Insurance. Berlin: Springer; 2011 pp. 293-328. [57] Weron R, Misiorek A. Forecasting spot electricity prices: A comparison of parametric and semiparametric time series models. International Journal of Forecasting. 2008; 24: 744-763. [58] Christoffersen P. Evaluating Interval Forecasts International Economic Review 1998; 39: 841-862.

29