- Email: [email protected]
Data preprocessing and quantile regression for probabilistic load forecasting in the GEFCom2017 final match
International Journal of Forecasting 35 (2019) 1460–1468
Contents lists available at ScienceDirect
International Journal of Forecasting journal homepage: www.elsevier.com/locate/ijforecast
Data preprocessing and quantile regression for probabilistic load forecasting in the GEFCom2017 final match ∗
Isao Kanda , J.M. Quintana Veguillas Japan Meteorological Corporation, Crystal Tower 17F, 1-2-27 Shiromi, Chuo, Osaka, 540-6017, Japan
article
info
Keywords: Electricity load GEFCom2017 Generalized additive method Probabilistic forecast Quantile regression
a b s t r a c t Team QUINKAN competed in the GEFCom2017 final match of hierarchical probabilistic load forecasting by adopting the quantile regression method using the R package quantreg. The weather stations were clustered into 11 groups, from which an optimal one was chosen for each load meter using the boosting method. The load meter records were cleaned and/or supplemented by various methods in order to secure robust quantile predictions. The variation in the regression formulas was kept as small as possible by introducing measures for suppressing prediction instability, although special formulas were employed for loading meters that were of an industrial nature. Several procedures were applied to help improve the accuracy, such as the smoothing of season transitions, coarse graining of the relative humidity, the use of load-oriented day-type definition, the averaging of weather data, and outlier removal. © 2019 The Authors. Published by Elsevier B.V. on behalf of International Institute of Forecasters. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
1. Introduction
2. Preprocessing
The task of the final match of the GEFCom2017 (Hong, Xie, & Black, 2019) was to forecast the hourly electricity load values in 2012 at 161 meters in an unknown area of the United States. The forecasts were required to be probabilistic ones given at the 10th, 25th, 50th, 75th, and 90th percentiles, and the 161 meters belonged to a threelevel hierarchical network, shown in Fig. 1. Past hourly records of the electricity load were given from 1 Jan. 2005 to 31 Dec. 2011. There were many missing periods, with an extreme case being for meter number 0499, which had data only from around the end of Nov. 2011. Hourly temperature and relative humidity data were given at 27 meteorological stations from 1 Jan. 2005 to 31 Dec. 2011 without any missing hours. No information was provided on the relationship between the load meters and the meteorological stations.
2.1. Meteorological data
∗ Corresponding author. E-mail address: [email protected] (I. Kanda).
2.1.1. Exclusion based on a visual inspection The temporal trends of the temperature (T ) and the relative humidity (q) were plotted and inspected by eye. Nine stations (indices 1, 2, 5, 8, 9, 15, 16, 19, and 25) were excluded from further analysis because of frequent occurrences of unnaturally long periods with constant values, singular spikes, or odd behaviors of q (such as being proportional to T for almost a year). 2.1.2. Station clustering The records of some stations were found to be similar. For simplicity and robustness, meteorological stations were clustered using the hierarchical clustering function hclust in R. The following two measures were employed as the distance for clustering: d0,ij = ⏐T i − T j ⏐ ,
⏐
⏐
(1)
https://doi.org/10.1016/j.ijforecast.2019.02.005 0169-2070/© 2019 The Authors. Published by Elsevier B.V. on behalf of International Institute of Forecasters. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468
1461
Fig. 1. Hierarchical network of load meters. The numbers in parentheses indicate the numbers of meters under the various nodes. Table 1 MSEs for original and averaged meteorological data. Meter
Day-type
Season
Original
Averaged
0168 0045 0045 0419
Weekday Weekday Weekday Weekday
Cold Cold Hot Hot
7 649 799 132 444 170 241 743 201
7 215 479 116 215 155 146 667 453
d1,ij
⏐) (⏐ = var ⏐Ti,I − Tj,I ⏐ ,
(2)
maintained properly. Several stations’ records exhibited such symptoms, often in the form of an upper limit on q of around 90% even on presumably rainy days. We also note that the precision of q is typically around 1%, which is relatively larger than that of T (∼0.1 ◦ C). Accounting for these properties, the q values were transformed using the following formula:
q←
⎧ ⎨5
(0 ≤ q < 10),
85
(85 < q),
(3)
5 [q/5] (10 ≤ q < 85),
⎩ where T i is the mean temperature of station i and Ti,I is the temperature of station i at sequential hourly index I. The distance measure d0,ij distinguishes climatic zones due for example to latitude differences, and d1,ij distinguishes phase differences with respect to synoptic or mesoscale weather systems due for example to longitude differences. Meteorological stations were clustered initially based on d0,ij , then each resulting cluster was divided into sub-clusters based on d1,ij . Fig. 2 shows the clustering results. 2.1.3. Averaging It was inferred that the load data at 1 h were the accumulated values from 0 h to 1 h, whereas the weather data at 1 h were the instantaneous values at 1 h. Thus, the weather conditions affecting the load at 1 h are the average of the values at 0 h and 1 h. We tested this hypothesis by conducting crossvalidation of the deterministic forecast using the generalized additive model (GAM; Hastie, 1992) based on the fit formula employed in the main probabilistic forecast (Section 5.2) for a few load meters (matching meteorological station clusters were chosen by the method described later in Section 3). The cross-validation was done by year (e.g., training = (2005, . . . , 2009, 2011), test = (2010)). The meteorological data corresponding to the load data DI were either the original (TI , qI ) or ( ) the averaged 21 (TI + TI −1 ), 21 (qI + qI −1 ) . Table 1 lists the resulting mean squared errors (MSEs). We found a clear reduction in MSEs when the meteorological data were averaged. Hence, 21 (TI + TI −1 ) and 21 (qI + qI −1 ) were used as the meteorological conditions at the sequential hourly index I. 2.1.4. Transformation of relative humidity Relative humidity sensors exhibit large uncertainties at low and high humidities, as is reflected by the less strict accuracy requirements below 10% and above 90% in the US Climate Reference Network (USCRN, 2007). The uncertainty increases if the sensors are not installed or
where [ ] gives the quotient of integer division. This transformation sets upper and lower limits, and coarsens the resolution of q from 1 to 5. 2.2. Season classification The year was divided into cold, mild, and hot seasons based on the diurnal pattern of the electricity load. In summer, a load peak occurred at around 17 h (in Daylight Saving Time), when the air conditioning demand was high due to the high air temperature. In winter, two peaks occurred, at around 6 h and 20 h (in local standard time), due to heating demand during periods of high human activity and low air temperature. In summer, the load rose as the temperature increased, whereas in winter, the load rose as the temperature decreased. As a result of the distinctly different characteristics of the load behaviors, we constructed separate forecasting models for different seasons. We defined electric-load seasons by examining the diurnal patterns of the load in detail, and considered only weekday records in order to keep the focus on seasonal variation. Writing the electricity load for the day of the year (J) and hour of the day (H) as DJ (H), the normalized diurnal pattern was defined as
ˆ J (H) = D
DJ (H) − minH {DJ (H)} maxH {DJ (H)} − minH {DJ (H)}
.
(4)
ˆ J (H), the Based on a visual inspection of the plots of D summer and winter months were selected, as a preliminary season division, as: Summer: June–September, Winter = November–April .
ˆ J (H) in the summer and winter The averages of D months were defined as ∑ ˆ A (H) = 1 ˆ J (H), D D (5) NA J ∈ΩA
1462
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468
Fig. 2. Clustering result for weather stations. The numbers in parentheses indicate the indices of the original weather stations.
where ΩA is a set of J in A (summer or winter), and NA is the number of elements in ΩA . Then, the normalized mean patterns in summer and winter were defined as
{
ˆ A (H) − minH Dˆ A (H) D FA (H) =
{
}
{
}
ˆ A (H) − minH Dˆ A (H) maxH D
}.
(6)
ˆ J (H) Fig. 3 shows FA (H) for 2015. Finally, the distances of D from FA (H) were defined as dJ ,A =
23 ∑
|Dˆ J (H) − FA (H)|.
(7)
h=0
If dJ ,summer < dJ ,winter , the diurnal pattern is more similar to the summer average than to the winter average, and vice versa. These quantities were calculated for each year. The electric load seasons were defined according to the daily mean temperature as follows. Fig. 4 shows dJ ,A against the daily mean temperature T J for meter 0054 in 2005. We observe that dJ ,summer (red) decreases with T J , whereas dJ ,winter (cyan) increases with T J . The points where dJ ,summer and dJ ,winter cross can be regarded as changes in seasons. An inspection of the plots for different meters and different years allowed us to define (rather subjectively) three seasons based on T J : Cold : T J < 50, Mild : 50 ≤ T J < 70, Hot : T J ≥ 70. For future days (prediction days), T J was predicted using GAM, available in R, by means of the following formula: T J ∼ I + (m = 1, 2, 3, 4),
(8)
where the notation (m = 1, 2, 3, 4) represents the sinusoidal functions described in Section 5.2.
2.3. Load data 2.3.1. Gross classification The sample included various types of load meters. A visual inspection of the temporal trends enabled the meters to be classified initially into three types: residential, industrial, and disused meters (Fig. 5). The residential meters exhibited diurnal and seasonal patterns due to human activities (see Figs. 3 or 8 for meter 0054). Some of the residential meters showed behaviors indicative of large universities or research institutes, because the load tended to decrease in summer. The industrial meters had distinctly different patterns from the residential ones: the diurnal or seasonal variation was often weak, and large differences existed among the industrial meters (see Fig. 8 for meter 0469). The disused meters were all discontinued at different points through the training period, and were expected to remain so in the forecast year 2012. 2.3.2. Cleaning Singular values from the load data were removed in multiple steps: (1) Manual removal of periods with clearly different levels, and automated removal of exceptionally low or high values, (2) Automated removal of spikes using the median filter, (3) Further manual removal of odd data such as the exactly duplicate periods of 0246 in July 2005, and any spikes that escaped the filter in step (2). 2.3.3. Gap filling Some meters did not have sufficient data to enable the production of robust quantile forecasts. For such meters, data from similar meters were inserted in the missing periods by the following procedure, which corresponds to the reconstruction of merged or split meters.
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468
1463
Fig. 3. Normalized mean diurnal patterns FA (H) of the electricity load of meter 0054 on weekdays in 2005. Red: summer; cyan: winter. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 4. Functional distance dJ ,A from the summer pattern (red) and winter pattern (cyan) of meter 0054 on weekdays in 2005, plotted against the daily mean temperature. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(1) The meter data were scaled as
ˆi = D
Di − Di
σ (Di )
,
(9)
where σ (Di ) is the standard deviation and Di is the mean for the period for which the meter of which the missing periods are to be supplemented has a valid continuous record. (2) Similar meters were sought by the hierarchical clustering function hclust in R, with the distance defined as
(
)
ˆ i,I − Dˆ j,I . dij = var D
(10)
A matching meter was chosen from the neighborhood of the resulting cluster dendrogram, taking into account that the match must be from the same midlevel node of the load network hierarchy (Fig. 1),
must have the same associated meteorological station cluster (described in Section 3), must have a similar D level, and must have a long and steady record. The search for a matching meter was also conducted on disused meters. Exceptions were allowed in some conditions depending on the degree of satisfaction of the remaining conditions. (3) The data in the missing periods of meter i were generated from the matching meter j as Di,I =
σ (Di,I ) ˆ Dj,I + Di . σ (Dj,I )
(11)
An example of the filtering and insertion procedures is shown in Fig. 6. For meter 0214, the data toward the end of 2009 were excluded because of the large drop in the level. A similar meter 0215 was identified in the same hierarchical group E022. The data in the period in which the meter 0214 became empty were inserted from meter 0215 using the rescaling formula in Eq. (11).
1464
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468
Fig. 5. Examples of (a) a residential meter, (b) an industrial meter, and (c) a disused meter.
Fig. 6. Example of the filtering and prepension of load data.
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468 Table 2 Cross-validation mean MSE values for meter 0054 with the different weather-station clusters. Meter
Weather-station cluster
MSE
0054 0054 0054 0054 0054 0054 0054 0054 0054 0054 0054
ws_1_1 ws_2_1 ws_2_2 ws_2_3 ws_3_1 ws_3_2 ws_3_3 ws_4_1 ws_4_2 ws_5_1 ws_5_2
975 971 517 068 1 029 371 589 437 1 030 616 746 558 1 171 720 340 271 709 251 709 676 911 032
5. Probabilistic forecast 5.1. Outline For the sake of simplicity, we denote time information (year, month, day, and hour) as t, and treat D, T , and q as functions of t instead of using the subscripts I and J as in the previous sections. It is assumed that t is from a particular season and a particular type of day. Our goal is to evaluate D(τ |t) at quantile τ conditional on t. As random variables, we assume that T depends on t, q depends on T and t, and D depends on T , q, and t. Then, the probability density function p(D | t) can be written as
∫ ∫
p(D | q, T , t) p(q | T , t) p(T | t) dq dT
p(D|t) = T
3. Matching between meter and weather-station clusters
q
∫ ∫ = T
For each residential meter, deterministic forecasts were produced using the generalized boosting method (package gbm of R; Ridgeway, 2017) with T and q from different meteorological station clusters, and cross-validation MSEs were evaluated on a yearly basis (e.g., train = (2005, 2006, . . . , 2009, 2011), test = (2010)). The boosting method was chosen because the resulting MSE was smallest among all of the forecast methods tested so far, although the stability of the quantile predictions, which is the goal of the final match’s task, was extremely poor. For simplicity, we evaluated weekdays in the hot season. The matching meteorological station cluster was selected as the one with the minimum MSE (cross-validation mean). For most meters, the best matching meteorological station cluster had a distinctly smaller MSE than the secondbest. Table 2 shows an example where weather station cluster 4_1 was chosen for meter 0054. We observe that the MSE of 340271 for 4_1 is substantially smaller than the second-best MSE of 517068 for 2_1.
For residential meters, based on human activities, the day types were defined as follows: Day-type Monday Friday Saturday Sunday Weekday
Definition Workday after a non-workday Workday before a non-workday, and not monday Non-workday that is not sunday Non-workday before a workday Workday that is not monday or friday
For industrial meters, there were periods when the load became distinctly small, e.g., early July and late December. These periods were labeled ‘‘mtn’’, standing for maintenance. Days of type ‘‘mtn’’ were identified manually in the past load data from 2005 to 2011, and the ‘‘mtn’’ days in July and December 2012 could be inferred readily based on the week numbers.
q
∂τ D ⏐⏐ ∂τ q ⏐⏐ ∂τ T dq dT , ⏐ ⏐ ∂ D T ,q ∂ q T ∂ T
(12)
where τ D , τ q , and τ T are the cumulative distribution functions for D, q, and T , respectively. By definition, we obtain
τ (D | t) =
D
∫
p(D′ | t) dD′ −∞ 1
∫
1
∫
τ D (D | T , q, t) dτ q dτ T .
= 0
(13)
0
Approximating the integral by the first-order discretization, we have
τ (D | t) =
∑
τ D (D | Tk , qjk , t)(τjq+1 − τjq )(τkT+1 − τkT ), (14)
j,k
where Tk and qjk are quantile values at the middle of inq q tervals [τkT , τkT+1 ] and [τj , τj+1 ], respectively. Specifically,
∫
Tk
p(T | t) dT = −∞
∫
qjk
−∞
4. Day-type definition
1465
1 2
(τkT + τkT+1 ),
p(q | Tk , t) dT =
1 2
q
(15)
q
(τj + τj+1 ).
(16) q
In the actual calculation, τkT and τj were taken in 0, 0.1, 0.2, · · ·, 1.0. Hence, Tk and qjk were evaluated at ten quantiles 0.05, 0.15, · · ·, 0.95. The R package quantreg was used to calculate the quantile values in the following steps:
• for a given t, Tk was calculated; • for each Tk , qjk was calculated; • for each (Tk , qjk ), quantile values Dijk at τiD in 0.05, 0.15, · · ·, 0.95 were calculated. Furthermore, the calculated Dijk (τiD ) was inverted to obtain τ D (D | Tk , qjk , t) on the right-hand-side of Eq. (14). Then, τ (D | t) was obtained by taking the sum of Dijk (τiD ) over j and k. Finally, D(τ | t) was obtained by inverting τ (D | t). The numerical techniques used in the summation and the function inversion are described in Section 5.3. 5.2. Regression formulas We obtained quantile forecasts using quantreg by employing different fit formulas for different quantities
1466
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468
or types of meter. We note that the quantreg forecasts were conducted separately for different seasons and types of day. In order to achieve periodicity, sinusoidal functions were used in the fit formulas (Harben & Giasemidis, 2016). For the day of the year J and hour of the day H, the following terms were used:
[sin, cos](2π mJ /365),
(17)
[sin, cos](2π nH /24),
(18)
where m and n are natural numbers in most cases, but n = 1/2 was also used to account for non-periodic patterns (see Eq. (24)). In addition, the following cross terms of J and H were introduced:
[sin, cos](2π J /365) × [sin, cos](2π H /24).
(19)
When the sinusoidal functions were used, sin and cos were used together for each m or n. For simplicity, the terms for J and H are denoted by the values of m and n, and the cross terms are denoted by cs. T and q were modeled by T ∼ I + (m = 1, 2, 3) + (n = 1, 2, 3) + cs,
(20)
q ∼ I + (m = 1, 2, 3) + (n = 1, 2, 3) + cs + T + T + T . (21) 2
3
For most residential meters, the load was modeled by D ∼ I + (m = 1, 2, 3) + (n = 1, 2, 3, 4, 5, 6, 7, 8)
+ cs + T + T 2 + q.
(22)
When the quantile forecast was unstable, simpler versions were employed. The most frequently-used simpler formula was D ∼ I + (m = 1, 2) + (n = 1, 2, 3, 4, 5, 6) + T + q.
(23)
For industrial meters, even simpler formulas were employed. For the non-maintenance periods of meters 0208 and 0469, D ∼ (m = 1, 2) + (n = 1/2, 1, 2, 3, 4) + T + q,
(24)
while for the maintenance days of meters 0208 and 0469, D ∼ (m = 1) + (n = 1, 2, 3) + T + q.
(25)
5.3. Numerical techniques
where 50 ± 34.1% are the percentiles of the standard normal distribution N(0, 1) at ± one standard deviation, respectively. Then, Ds was set by Ds − median(D)
−1 = τN(0 (27) ,1) ([0.05, 0.15, . . . , 0.95]), σ where τN(0,1) is the cumulative distribution function for N(0, 1).
5.3.2. Function inversion Given D(τ ), the values of τ at D = Ds were obtained by the following method. Because D(τ ) is not linear and tends toward limiting values as τ approaches 0 or 1, linear interpolation does not work well near τ = 0 or 1 unless Ds is distributed densely near τ = 0 or 1. By gross approximation, we assume
τ∼
) 1( tanh D′ + 1 , 2
(28)
where D′ is some linear function of D. Then, if we write x = tanh−1 (2τ − 1),
(29)
′
D and x, and hence D and x, have approximately linear relationships. We can then regress D ∼ β0 + β1 x or any other simple formula such as natural splines to obtain τ at D = Ds . This inversion technique was adopted for calculating τjk (Ds ) (= τ D (D | Tk , qjk , t) on the right-handside of Eq. (14) from Dijk (τiD ), and for calculating D(τ | t) from τ (D | t) on the left-hand side of Eq. (14). Note that, unlike in the preceding steps, the final target D(τ | t) was evaluated at τ = 0.1, 0.25, 0.5, 0.75, and 0.9. 5.3.3. Season transition The season transition was made smooth by introducing buffer regions between the seasons defined in Section 2.2. As Fig. 7 shows, the model coefficients in the regression formulas in Section 5.2 were interpolated linearly between the seasons in the temperature ranges ±5 ◦ F around the thresholds 50 ◦ F and 70 ◦ F. 5.3.4. Final filtering Even after applying the measures described so far, some singular quantile predictions still occurred. Because they did not occur very often, they were detected by a visual inspection of the trend plots and removed by the median filter with an appropriate set of parameters. 5.4. Prediction results
5.3.1. Summation of Dijk (τ Summing over j and k in Eq. (14) required Dijk for τiD to be inverted to τ D as a function of D. To do this, Dijk (τiD ) was interpolated at a set of values Ds , and the resulting τjk (Ds ) was summed over j and k. We assumed that Dijk obeys a normal distribution, and selected Ds from the quantiles around the median. We used the median rather than the mean because there were often outlier values. The standard deviation σ was determined by D i )
⏐} {⏐ σ = min ⏐percentile(Dijk , 50 ± 34.1%) − median(Dijk )⏐ , (26)
The hourly electricity loads in 2012 at the 161 meters were predicted using the methods described above. For the calculation, we used a workstation with dual CPU (Intel Xeon E5-2683 v4 2.1 GHz), each with 16 cores. Using the parallel computation functions of R, the annual forecast of 161 meters took about 6.7 h. Fig. 8 shows the prediction results of a residential meter (0054) and an industrial meter (0469) in January and August 2012. We observe that the difference between the 10th and 90th percentiles for meter 0054 remains almost constant through the day in January, while that in August varies diurnally: small at night and large during the day. It should
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468
1467
Fig. 7. Season transition.
Fig. 8. Prediction results for a residential meter (0054) and an industrial meter (0469) in January and August 2012. The borders of the color bands indicate the 10th, 25th, 50th, 75th, and 90th percentiles, from the bottom to the top. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
be noted that the diurnal variation in the 10th–90th percentile difference in August is similar to that obtained in the GEFCom2017 qualifying match for the prediction from January to April in ISO-NE of the USA. However, the observed constancy of the 10th–90th percentile difference in January was confirmed to be consistent with the constancy of the load variance at different hours of the day. The industrial meter 0469 consumes large amounts of electricity throughout the day on days of types ‘weekday’, ‘monday’, and ‘friday’ (as defined in Section 4), and occasionally ‘saturday’. This irregularity on ‘saturday’ is the cause of the larger variance on ‘saturday’ than on ‘sunday’. The moderately large variance in ‘friday’ is due to the occasional low load on ‘friday’ (though usually not as low as on ‘saturday’ and ‘sunday’). 6. Discussion Although we were fortunate enough to win the final match, there are several aspects of our methods that we find unsatisfactory. (1) Use of boosting We did use the boosting method to find the best matches between the load meters and the meteorological stations, because boosting resulted in a mean
square error that was almost half the size of that obtained from GAM using the same fit formula as in the quantile regression (see Section 3). However, a naive application of boosting to the quantile regression task resulted in an irregular ordering of quantile values at most of the data hours. We suppose that this behavior is due to the random-number generation in boosting; however, since we did not have either the knowledge or the skills to overcome this problem at the time of the competition, we simply dismissed this method. A substantial improvement in accuracy could be achieved if quantile regression based on boosting could be conducted in a stable manner. (2) Effective variables We conducted so-called feature engineering by redefining meteorological variables (Sections 2.1.3 and 2.1.4). The ideas behind these operations are derived from meteorological knowledge, but we expect that more such effective features could be inferred from external wisdom such as biometeorology. One such feature that came to our notice recently is the contribution of the wind speed (Xie & Hong, 2017). (3) Regression formulas The regression formulas used for quantreg were determined in an ad-hoc manner, e.g. the tuning
1468
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468
of the m and n values in Eqs. (17) and (18). We most regret the use of cyclic functions for the diurnal pattern in Eq. (18). Because the models were constructed for different types of days, it is unlikely that the demand at the beginning of a day will always be equal to that at the end of the day. Although the idea of using separate models for each hour (e.g. Dordonnat, Pichavant, & Pierrot, 2016) was dismissed for fear of losing the stability of quantile forecasts due to the relatively small amount of training data, it might be worth trying. (4) Automation Our method involves a fair number of manual operations: data cleaning (Section 2.3.2), gap filling (Section 2.3.3), matching between meters and meteorological stations (Section 3), and selecting optimal regression formulas (Section 5.2). The manual operations were considered acceptable because the number of exceptional cases was manageable. However, we find that it would be feasible to automate some of the processes because some patterns of exceptions have been clarified (e.g. typical ways in which meters are split or combined).
Acknowledgments We are grateful to JMC president M. Suzuki for his long-standing support of our challenges. References Dordonnat, V., Pichavant, A., & Pierrot, A. (2016). GEFCom2014 probabilistic electric load forecasting using time series and semiparametric regression models. International Journal of Forecasting, 32, 1005–1011. Harben, S., & Giasemidis, G. (2016). A hybrid model of kernel density estimation and quantile regression for GEFCom2014 probabilistic load forecasting. International Journal of Forecasting, 32, 1017–1022. Hastie, T. J. (1992). Generalized additive models. In J. M. Chambers, & T. J. Hastie (Eds.), Statistical models in S. Wadsworth & Brooks/Cole. Hong, T., Xie, J., & Black, J. (2019). Global energy forecasting competition 2017: Hierarchical probabilistic load forecasting. International Journal of Forecasting. Ridgeway, G. (2017). gbm: Generalized Boosted Regression Models. R package version 2.1.3. USCRN (2007). Functional requirements document, revision 1., United States Climate Reference Network, National Environmental Satellite, Data, and Information Service. Xie, J., & Hong, T. (2017). Wind speed for load forecasting models. Sustainability, 9, 795.
Contents lists available at ScienceDirect
International Journal of Forecasting journal homepage: www.elsevier.com/locate/ijforecast
Data preprocessing and quantile regression for probabilistic load forecasting in the GEFCom2017 final match ∗
Isao Kanda , J.M. Quintana Veguillas Japan Meteorological Corporation, Crystal Tower 17F, 1-2-27 Shiromi, Chuo, Osaka, 540-6017, Japan
article
info
Keywords: Electricity load GEFCom2017 Generalized additive method Probabilistic forecast Quantile regression
a b s t r a c t Team QUINKAN competed in the GEFCom2017 final match of hierarchical probabilistic load forecasting by adopting the quantile regression method using the R package quantreg. The weather stations were clustered into 11 groups, from which an optimal one was chosen for each load meter using the boosting method. The load meter records were cleaned and/or supplemented by various methods in order to secure robust quantile predictions. The variation in the regression formulas was kept as small as possible by introducing measures for suppressing prediction instability, although special formulas were employed for loading meters that were of an industrial nature. Several procedures were applied to help improve the accuracy, such as the smoothing of season transitions, coarse graining of the relative humidity, the use of load-oriented day-type definition, the averaging of weather data, and outlier removal. © 2019 The Authors. Published by Elsevier B.V. on behalf of International Institute of Forecasters. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
1. Introduction
2. Preprocessing
The task of the final match of the GEFCom2017 (Hong, Xie, & Black, 2019) was to forecast the hourly electricity load values in 2012 at 161 meters in an unknown area of the United States. The forecasts were required to be probabilistic ones given at the 10th, 25th, 50th, 75th, and 90th percentiles, and the 161 meters belonged to a threelevel hierarchical network, shown in Fig. 1. Past hourly records of the electricity load were given from 1 Jan. 2005 to 31 Dec. 2011. There were many missing periods, with an extreme case being for meter number 0499, which had data only from around the end of Nov. 2011. Hourly temperature and relative humidity data were given at 27 meteorological stations from 1 Jan. 2005 to 31 Dec. 2011 without any missing hours. No information was provided on the relationship between the load meters and the meteorological stations.
2.1. Meteorological data
∗ Corresponding author. E-mail address: [email protected] (I. Kanda).
2.1.1. Exclusion based on a visual inspection The temporal trends of the temperature (T ) and the relative humidity (q) were plotted and inspected by eye. Nine stations (indices 1, 2, 5, 8, 9, 15, 16, 19, and 25) were excluded from further analysis because of frequent occurrences of unnaturally long periods with constant values, singular spikes, or odd behaviors of q (such as being proportional to T for almost a year). 2.1.2. Station clustering The records of some stations were found to be similar. For simplicity and robustness, meteorological stations were clustered using the hierarchical clustering function hclust in R. The following two measures were employed as the distance for clustering: d0,ij = ⏐T i − T j ⏐ ,
⏐
⏐
(1)
https://doi.org/10.1016/j.ijforecast.2019.02.005 0169-2070/© 2019 The Authors. Published by Elsevier B.V. on behalf of International Institute of Forecasters. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468
1461
Fig. 1. Hierarchical network of load meters. The numbers in parentheses indicate the numbers of meters under the various nodes. Table 1 MSEs for original and averaged meteorological data. Meter
Day-type
Season
Original
Averaged
0168 0045 0045 0419
Weekday Weekday Weekday Weekday
Cold Cold Hot Hot
7 649 799 132 444 170 241 743 201
7 215 479 116 215 155 146 667 453
d1,ij
⏐) (⏐ = var ⏐Ti,I − Tj,I ⏐ ,
(2)
maintained properly. Several stations’ records exhibited such symptoms, often in the form of an upper limit on q of around 90% even on presumably rainy days. We also note that the precision of q is typically around 1%, which is relatively larger than that of T (∼0.1 ◦ C). Accounting for these properties, the q values were transformed using the following formula:
q←
⎧ ⎨5
(0 ≤ q < 10),
85
(85 < q),
(3)
5 [q/5] (10 ≤ q < 85),
⎩ where T i is the mean temperature of station i and Ti,I is the temperature of station i at sequential hourly index I. The distance measure d0,ij distinguishes climatic zones due for example to latitude differences, and d1,ij distinguishes phase differences with respect to synoptic or mesoscale weather systems due for example to longitude differences. Meteorological stations were clustered initially based on d0,ij , then each resulting cluster was divided into sub-clusters based on d1,ij . Fig. 2 shows the clustering results. 2.1.3. Averaging It was inferred that the load data at 1 h were the accumulated values from 0 h to 1 h, whereas the weather data at 1 h were the instantaneous values at 1 h. Thus, the weather conditions affecting the load at 1 h are the average of the values at 0 h and 1 h. We tested this hypothesis by conducting crossvalidation of the deterministic forecast using the generalized additive model (GAM; Hastie, 1992) based on the fit formula employed in the main probabilistic forecast (Section 5.2) for a few load meters (matching meteorological station clusters were chosen by the method described later in Section 3). The cross-validation was done by year (e.g., training = (2005, . . . , 2009, 2011), test = (2010)). The meteorological data corresponding to the load data DI were either the original (TI , qI ) or ( ) the averaged 21 (TI + TI −1 ), 21 (qI + qI −1 ) . Table 1 lists the resulting mean squared errors (MSEs). We found a clear reduction in MSEs when the meteorological data were averaged. Hence, 21 (TI + TI −1 ) and 21 (qI + qI −1 ) were used as the meteorological conditions at the sequential hourly index I. 2.1.4. Transformation of relative humidity Relative humidity sensors exhibit large uncertainties at low and high humidities, as is reflected by the less strict accuracy requirements below 10% and above 90% in the US Climate Reference Network (USCRN, 2007). The uncertainty increases if the sensors are not installed or
where [ ] gives the quotient of integer division. This transformation sets upper and lower limits, and coarsens the resolution of q from 1 to 5. 2.2. Season classification The year was divided into cold, mild, and hot seasons based on the diurnal pattern of the electricity load. In summer, a load peak occurred at around 17 h (in Daylight Saving Time), when the air conditioning demand was high due to the high air temperature. In winter, two peaks occurred, at around 6 h and 20 h (in local standard time), due to heating demand during periods of high human activity and low air temperature. In summer, the load rose as the temperature increased, whereas in winter, the load rose as the temperature decreased. As a result of the distinctly different characteristics of the load behaviors, we constructed separate forecasting models for different seasons. We defined electric-load seasons by examining the diurnal patterns of the load in detail, and considered only weekday records in order to keep the focus on seasonal variation. Writing the electricity load for the day of the year (J) and hour of the day (H) as DJ (H), the normalized diurnal pattern was defined as
ˆ J (H) = D
DJ (H) − minH {DJ (H)} maxH {DJ (H)} − minH {DJ (H)}
.
(4)
ˆ J (H), the Based on a visual inspection of the plots of D summer and winter months were selected, as a preliminary season division, as: Summer: June–September, Winter = November–April .
ˆ J (H) in the summer and winter The averages of D months were defined as ∑ ˆ A (H) = 1 ˆ J (H), D D (5) NA J ∈ΩA
1462
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468
Fig. 2. Clustering result for weather stations. The numbers in parentheses indicate the indices of the original weather stations.
where ΩA is a set of J in A (summer or winter), and NA is the number of elements in ΩA . Then, the normalized mean patterns in summer and winter were defined as
{
ˆ A (H) − minH Dˆ A (H) D FA (H) =
{
}
{
}
ˆ A (H) − minH Dˆ A (H) maxH D
}.
(6)
ˆ J (H) Fig. 3 shows FA (H) for 2015. Finally, the distances of D from FA (H) were defined as dJ ,A =
23 ∑
|Dˆ J (H) − FA (H)|.
(7)
h=0
If dJ ,summer < dJ ,winter , the diurnal pattern is more similar to the summer average than to the winter average, and vice versa. These quantities were calculated for each year. The electric load seasons were defined according to the daily mean temperature as follows. Fig. 4 shows dJ ,A against the daily mean temperature T J for meter 0054 in 2005. We observe that dJ ,summer (red) decreases with T J , whereas dJ ,winter (cyan) increases with T J . The points where dJ ,summer and dJ ,winter cross can be regarded as changes in seasons. An inspection of the plots for different meters and different years allowed us to define (rather subjectively) three seasons based on T J : Cold : T J < 50, Mild : 50 ≤ T J < 70, Hot : T J ≥ 70. For future days (prediction days), T J was predicted using GAM, available in R, by means of the following formula: T J ∼ I + (m = 1, 2, 3, 4),
(8)
where the notation (m = 1, 2, 3, 4) represents the sinusoidal functions described in Section 5.2.
2.3. Load data 2.3.1. Gross classification The sample included various types of load meters. A visual inspection of the temporal trends enabled the meters to be classified initially into three types: residential, industrial, and disused meters (Fig. 5). The residential meters exhibited diurnal and seasonal patterns due to human activities (see Figs. 3 or 8 for meter 0054). Some of the residential meters showed behaviors indicative of large universities or research institutes, because the load tended to decrease in summer. The industrial meters had distinctly different patterns from the residential ones: the diurnal or seasonal variation was often weak, and large differences existed among the industrial meters (see Fig. 8 for meter 0469). The disused meters were all discontinued at different points through the training period, and were expected to remain so in the forecast year 2012. 2.3.2. Cleaning Singular values from the load data were removed in multiple steps: (1) Manual removal of periods with clearly different levels, and automated removal of exceptionally low or high values, (2) Automated removal of spikes using the median filter, (3) Further manual removal of odd data such as the exactly duplicate periods of 0246 in July 2005, and any spikes that escaped the filter in step (2). 2.3.3. Gap filling Some meters did not have sufficient data to enable the production of robust quantile forecasts. For such meters, data from similar meters were inserted in the missing periods by the following procedure, which corresponds to the reconstruction of merged or split meters.
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468
1463
Fig. 3. Normalized mean diurnal patterns FA (H) of the electricity load of meter 0054 on weekdays in 2005. Red: summer; cyan: winter. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 4. Functional distance dJ ,A from the summer pattern (red) and winter pattern (cyan) of meter 0054 on weekdays in 2005, plotted against the daily mean temperature. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(1) The meter data were scaled as
ˆi = D
Di − Di
σ (Di )
,
(9)
where σ (Di ) is the standard deviation and Di is the mean for the period for which the meter of which the missing periods are to be supplemented has a valid continuous record. (2) Similar meters were sought by the hierarchical clustering function hclust in R, with the distance defined as
(
)
ˆ i,I − Dˆ j,I . dij = var D
(10)
A matching meter was chosen from the neighborhood of the resulting cluster dendrogram, taking into account that the match must be from the same midlevel node of the load network hierarchy (Fig. 1),
must have the same associated meteorological station cluster (described in Section 3), must have a similar D level, and must have a long and steady record. The search for a matching meter was also conducted on disused meters. Exceptions were allowed in some conditions depending on the degree of satisfaction of the remaining conditions. (3) The data in the missing periods of meter i were generated from the matching meter j as Di,I =
σ (Di,I ) ˆ Dj,I + Di . σ (Dj,I )
(11)
An example of the filtering and insertion procedures is shown in Fig. 6. For meter 0214, the data toward the end of 2009 were excluded because of the large drop in the level. A similar meter 0215 was identified in the same hierarchical group E022. The data in the period in which the meter 0214 became empty were inserted from meter 0215 using the rescaling formula in Eq. (11).
1464
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468
Fig. 5. Examples of (a) a residential meter, (b) an industrial meter, and (c) a disused meter.
Fig. 6. Example of the filtering and prepension of load data.
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468 Table 2 Cross-validation mean MSE values for meter 0054 with the different weather-station clusters. Meter
Weather-station cluster
MSE
0054 0054 0054 0054 0054 0054 0054 0054 0054 0054 0054
ws_1_1 ws_2_1 ws_2_2 ws_2_3 ws_3_1 ws_3_2 ws_3_3 ws_4_1 ws_4_2 ws_5_1 ws_5_2
975 971 517 068 1 029 371 589 437 1 030 616 746 558 1 171 720 340 271 709 251 709 676 911 032
5. Probabilistic forecast 5.1. Outline For the sake of simplicity, we denote time information (year, month, day, and hour) as t, and treat D, T , and q as functions of t instead of using the subscripts I and J as in the previous sections. It is assumed that t is from a particular season and a particular type of day. Our goal is to evaluate D(τ |t) at quantile τ conditional on t. As random variables, we assume that T depends on t, q depends on T and t, and D depends on T , q, and t. Then, the probability density function p(D | t) can be written as
∫ ∫
p(D | q, T , t) p(q | T , t) p(T | t) dq dT
p(D|t) = T
3. Matching between meter and weather-station clusters
q
∫ ∫ = T
For each residential meter, deterministic forecasts were produced using the generalized boosting method (package gbm of R; Ridgeway, 2017) with T and q from different meteorological station clusters, and cross-validation MSEs were evaluated on a yearly basis (e.g., train = (2005, 2006, . . . , 2009, 2011), test = (2010)). The boosting method was chosen because the resulting MSE was smallest among all of the forecast methods tested so far, although the stability of the quantile predictions, which is the goal of the final match’s task, was extremely poor. For simplicity, we evaluated weekdays in the hot season. The matching meteorological station cluster was selected as the one with the minimum MSE (cross-validation mean). For most meters, the best matching meteorological station cluster had a distinctly smaller MSE than the secondbest. Table 2 shows an example where weather station cluster 4_1 was chosen for meter 0054. We observe that the MSE of 340271 for 4_1 is substantially smaller than the second-best MSE of 517068 for 2_1.
For residential meters, based on human activities, the day types were defined as follows: Day-type Monday Friday Saturday Sunday Weekday
Definition Workday after a non-workday Workday before a non-workday, and not monday Non-workday that is not sunday Non-workday before a workday Workday that is not monday or friday
For industrial meters, there were periods when the load became distinctly small, e.g., early July and late December. These periods were labeled ‘‘mtn’’, standing for maintenance. Days of type ‘‘mtn’’ were identified manually in the past load data from 2005 to 2011, and the ‘‘mtn’’ days in July and December 2012 could be inferred readily based on the week numbers.
q
∂τ D ⏐⏐ ∂τ q ⏐⏐ ∂τ T dq dT , ⏐ ⏐ ∂ D T ,q ∂ q T ∂ T
(12)
where τ D , τ q , and τ T are the cumulative distribution functions for D, q, and T , respectively. By definition, we obtain
τ (D | t) =
D
∫
p(D′ | t) dD′ −∞ 1
∫
1
∫
τ D (D | T , q, t) dτ q dτ T .
= 0
(13)
0
Approximating the integral by the first-order discretization, we have
τ (D | t) =
∑
τ D (D | Tk , qjk , t)(τjq+1 − τjq )(τkT+1 − τkT ), (14)
j,k
where Tk and qjk are quantile values at the middle of inq q tervals [τkT , τkT+1 ] and [τj , τj+1 ], respectively. Specifically,
∫
Tk
p(T | t) dT = −∞
∫
qjk
−∞
4. Day-type definition
1465
1 2
(τkT + τkT+1 ),
p(q | Tk , t) dT =
1 2
q
(15)
q
(τj + τj+1 ).
(16) q
In the actual calculation, τkT and τj were taken in 0, 0.1, 0.2, · · ·, 1.0. Hence, Tk and qjk were evaluated at ten quantiles 0.05, 0.15, · · ·, 0.95. The R package quantreg was used to calculate the quantile values in the following steps:
• for a given t, Tk was calculated; • for each Tk , qjk was calculated; • for each (Tk , qjk ), quantile values Dijk at τiD in 0.05, 0.15, · · ·, 0.95 were calculated. Furthermore, the calculated Dijk (τiD ) was inverted to obtain τ D (D | Tk , qjk , t) on the right-hand-side of Eq. (14). Then, τ (D | t) was obtained by taking the sum of Dijk (τiD ) over j and k. Finally, D(τ | t) was obtained by inverting τ (D | t). The numerical techniques used in the summation and the function inversion are described in Section 5.3. 5.2. Regression formulas We obtained quantile forecasts using quantreg by employing different fit formulas for different quantities
1466
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468
or types of meter. We note that the quantreg forecasts were conducted separately for different seasons and types of day. In order to achieve periodicity, sinusoidal functions were used in the fit formulas (Harben & Giasemidis, 2016). For the day of the year J and hour of the day H, the following terms were used:
[sin, cos](2π mJ /365),
(17)
[sin, cos](2π nH /24),
(18)
where m and n are natural numbers in most cases, but n = 1/2 was also used to account for non-periodic patterns (see Eq. (24)). In addition, the following cross terms of J and H were introduced:
[sin, cos](2π J /365) × [sin, cos](2π H /24).
(19)
When the sinusoidal functions were used, sin and cos were used together for each m or n. For simplicity, the terms for J and H are denoted by the values of m and n, and the cross terms are denoted by cs. T and q were modeled by T ∼ I + (m = 1, 2, 3) + (n = 1, 2, 3) + cs,
(20)
q ∼ I + (m = 1, 2, 3) + (n = 1, 2, 3) + cs + T + T + T . (21) 2
3
For most residential meters, the load was modeled by D ∼ I + (m = 1, 2, 3) + (n = 1, 2, 3, 4, 5, 6, 7, 8)
+ cs + T + T 2 + q.
(22)
When the quantile forecast was unstable, simpler versions were employed. The most frequently-used simpler formula was D ∼ I + (m = 1, 2) + (n = 1, 2, 3, 4, 5, 6) + T + q.
(23)
For industrial meters, even simpler formulas were employed. For the non-maintenance periods of meters 0208 and 0469, D ∼ (m = 1, 2) + (n = 1/2, 1, 2, 3, 4) + T + q,
(24)
while for the maintenance days of meters 0208 and 0469, D ∼ (m = 1) + (n = 1, 2, 3) + T + q.
(25)
5.3. Numerical techniques
where 50 ± 34.1% are the percentiles of the standard normal distribution N(0, 1) at ± one standard deviation, respectively. Then, Ds was set by Ds − median(D)
−1 = τN(0 (27) ,1) ([0.05, 0.15, . . . , 0.95]), σ where τN(0,1) is the cumulative distribution function for N(0, 1).
5.3.2. Function inversion Given D(τ ), the values of τ at D = Ds were obtained by the following method. Because D(τ ) is not linear and tends toward limiting values as τ approaches 0 or 1, linear interpolation does not work well near τ = 0 or 1 unless Ds is distributed densely near τ = 0 or 1. By gross approximation, we assume
τ∼
) 1( tanh D′ + 1 , 2
(28)
where D′ is some linear function of D. Then, if we write x = tanh−1 (2τ − 1),
(29)
′
D and x, and hence D and x, have approximately linear relationships. We can then regress D ∼ β0 + β1 x or any other simple formula such as natural splines to obtain τ at D = Ds . This inversion technique was adopted for calculating τjk (Ds ) (= τ D (D | Tk , qjk , t) on the right-handside of Eq. (14) from Dijk (τiD ), and for calculating D(τ | t) from τ (D | t) on the left-hand side of Eq. (14). Note that, unlike in the preceding steps, the final target D(τ | t) was evaluated at τ = 0.1, 0.25, 0.5, 0.75, and 0.9. 5.3.3. Season transition The season transition was made smooth by introducing buffer regions between the seasons defined in Section 2.2. As Fig. 7 shows, the model coefficients in the regression formulas in Section 5.2 were interpolated linearly between the seasons in the temperature ranges ±5 ◦ F around the thresholds 50 ◦ F and 70 ◦ F. 5.3.4. Final filtering Even after applying the measures described so far, some singular quantile predictions still occurred. Because they did not occur very often, they were detected by a visual inspection of the trend plots and removed by the median filter with an appropriate set of parameters. 5.4. Prediction results
5.3.1. Summation of Dijk (τ Summing over j and k in Eq. (14) required Dijk for τiD to be inverted to τ D as a function of D. To do this, Dijk (τiD ) was interpolated at a set of values Ds , and the resulting τjk (Ds ) was summed over j and k. We assumed that Dijk obeys a normal distribution, and selected Ds from the quantiles around the median. We used the median rather than the mean because there were often outlier values. The standard deviation σ was determined by D i )
⏐} {⏐ σ = min ⏐percentile(Dijk , 50 ± 34.1%) − median(Dijk )⏐ , (26)
The hourly electricity loads in 2012 at the 161 meters were predicted using the methods described above. For the calculation, we used a workstation with dual CPU (Intel Xeon E5-2683 v4 2.1 GHz), each with 16 cores. Using the parallel computation functions of R, the annual forecast of 161 meters took about 6.7 h. Fig. 8 shows the prediction results of a residential meter (0054) and an industrial meter (0469) in January and August 2012. We observe that the difference between the 10th and 90th percentiles for meter 0054 remains almost constant through the day in January, while that in August varies diurnally: small at night and large during the day. It should
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468
1467
Fig. 7. Season transition.
Fig. 8. Prediction results for a residential meter (0054) and an industrial meter (0469) in January and August 2012. The borders of the color bands indicate the 10th, 25th, 50th, 75th, and 90th percentiles, from the bottom to the top. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
be noted that the diurnal variation in the 10th–90th percentile difference in August is similar to that obtained in the GEFCom2017 qualifying match for the prediction from January to April in ISO-NE of the USA. However, the observed constancy of the 10th–90th percentile difference in January was confirmed to be consistent with the constancy of the load variance at different hours of the day. The industrial meter 0469 consumes large amounts of electricity throughout the day on days of types ‘weekday’, ‘monday’, and ‘friday’ (as defined in Section 4), and occasionally ‘saturday’. This irregularity on ‘saturday’ is the cause of the larger variance on ‘saturday’ than on ‘sunday’. The moderately large variance in ‘friday’ is due to the occasional low load on ‘friday’ (though usually not as low as on ‘saturday’ and ‘sunday’). 6. Discussion Although we were fortunate enough to win the final match, there are several aspects of our methods that we find unsatisfactory. (1) Use of boosting We did use the boosting method to find the best matches between the load meters and the meteorological stations, because boosting resulted in a mean
square error that was almost half the size of that obtained from GAM using the same fit formula as in the quantile regression (see Section 3). However, a naive application of boosting to the quantile regression task resulted in an irregular ordering of quantile values at most of the data hours. We suppose that this behavior is due to the random-number generation in boosting; however, since we did not have either the knowledge or the skills to overcome this problem at the time of the competition, we simply dismissed this method. A substantial improvement in accuracy could be achieved if quantile regression based on boosting could be conducted in a stable manner. (2) Effective variables We conducted so-called feature engineering by redefining meteorological variables (Sections 2.1.3 and 2.1.4). The ideas behind these operations are derived from meteorological knowledge, but we expect that more such effective features could be inferred from external wisdom such as biometeorology. One such feature that came to our notice recently is the contribution of the wind speed (Xie & Hong, 2017). (3) Regression formulas The regression formulas used for quantreg were determined in an ad-hoc manner, e.g. the tuning
1468
I. Kanda and J.M.Q. Veguillas / International Journal of Forecasting 35 (2019) 1460–1468
of the m and n values in Eqs. (17) and (18). We most regret the use of cyclic functions for the diurnal pattern in Eq. (18). Because the models were constructed for different types of days, it is unlikely that the demand at the beginning of a day will always be equal to that at the end of the day. Although the idea of using separate models for each hour (e.g. Dordonnat, Pichavant, & Pierrot, 2016) was dismissed for fear of losing the stability of quantile forecasts due to the relatively small amount of training data, it might be worth trying. (4) Automation Our method involves a fair number of manual operations: data cleaning (Section 2.3.2), gap filling (Section 2.3.3), matching between meters and meteorological stations (Section 3), and selecting optimal regression formulas (Section 5.2). The manual operations were considered acceptable because the number of exceptional cases was manageable. However, we find that it would be feasible to automate some of the processes because some patterns of exceptions have been clarified (e.g. typical ways in which meters are split or combined).
Acknowledgments We are grateful to JMC president M. Suzuki for his long-standing support of our challenges. References Dordonnat, V., Pichavant, A., & Pierrot, A. (2016). GEFCom2014 probabilistic electric load forecasting using time series and semiparametric regression models. International Journal of Forecasting, 32, 1005–1011. Harben, S., & Giasemidis, G. (2016). A hybrid model of kernel density estimation and quantile regression for GEFCom2014 probabilistic load forecasting. International Journal of Forecasting, 32, 1017–1022. Hastie, T. J. (1992). Generalized additive models. In J. M. Chambers, & T. J. Hastie (Eds.), Statistical models in S. Wadsworth & Brooks/Cole. Hong, T., Xie, J., & Black, J. (2019). Global energy forecasting competition 2017: Hierarchical probabilistic load forecasting. International Journal of Forecasting. Ridgeway, G. (2017). gbm: Generalized Boosted Regression Models. R package version 2.1.3. USCRN (2007). Functional requirements document, revision 1., United States Climate Reference Network, National Environmental Satellite, Data, and Information Service. Xie, J., & Hong, T. (2017). Wind speed for load forecasting models. Sustainability, 9, 795.