Weather forecast information is typically observed at intervals of one day and not 30 min (e.g., the maximum temperature or average humidity). For this reason, we assume that the weather forecast information does not depend on j. Let a vector of weather information be s _i . We express b _ij as a function of s _i . Here, we assume two structures as follows:

• It is well known that the relationship between weather variables and consumption is nonlinear. For example, the relationship between maximum temperature and consumption is approximated by a quadratic function (e.g., 30) because air conditioners are used on both hot and cold days. For this reason, it is assumed that b _ij is expressed as some nonlinear function of s _i .

Under the above considerations, we propose expressing b _ij as follows: b i j = ∑ m = 1 M β m ( j ) g m ( s i ) , (3) where g _m ( s _i ) (m = 1, … , M) are basis functions given beforehand, β _m (j) are functions of regression coefficients, and M is the number of basis functions.

We also use the basis expansion for β _m (j): β m ( j ) = ∑ q = 1 Q γ q m h q ( j ) , (4) where h _q (j) (q = 1, … , Q) are basis functions given beforehand and γ _qm are the elements of the coefficient matrix Γ = (γ _qm ). Substituting (4) into (3) results in the following: b i j = ∑ m = 1 M ∑ q = 1 Q γ q m h q ( j ) g m ( s i ) . (5)

Because h _q (j) and g _m ( s _i ) are known functions, the parameters concerning b _ij are γ _qm .

Since the effect of weather is assumed to be smooth according to both j and s _i , we use basis functions h _q (j) and g _m ( s _i ), which produce a smooth function, such as B-spline and the radial basis function (RBF). Our numerical study adopts the B-spline and cyclic B-spline functions as basis functions of g _m ( s _i ) and h _q (j), respectively. For detail, please refer to Section 4.1. The RBF could be used but it includes a tuning parameter that highly affects the forecast error, resulting in low forecast accuracy in our experience.

Since μ _ij is the effect of past consumption, one can assume that μ _ij is expressed as a linear combination of past consumption y ( i − t − L α ) j , i.e., μ i j = ∑ t = 1 T α j t y ( i − t − L α ) j + ∑ u = 1 U β j u y i ( j − u − L β ) , (6) where T and U are positive integers, which denote how far we trace back through the data and α _jt (t = 1, … , T) and β _ju (u = 1, … , U) are positive values given beforehand. The α _jt correspond to the effects of past consumptions for the same time interval on previous days, while β _jt are the coefficients for different time intervals on the same day. For the day-ahead market, we assume that β _ju ≡ 0. Here, L _α and L _β are nonnegative integers that describe the lags; these values change according to the closing time of transactions. For example, transactions of the day-ahead market in the JEPX close at 17:00 on the previous day of the delivery of electricity. Thus, for instance, when the trade of the electricity demand is in the 17:30–18:00 interval on Sep 20th, we cannot use the electricity demand in the same time interval on Sep 19th due to the trading hours of the market, resulting in L _α = 1.

The design of T and U depends on the aim of the analysis. For example, if the analyst should use only the latest demand, T and U are selected as (T, U) = (1, 1). The usage of only one latest demand value could often yield unstable forecast value; in this case, larger T and U would provide better results. However, too large T and U can result in overfitting. Thus, we do not generally know the best T and U to achieve the best forecast accuracy. In our numerical experiments in Sections 4 and 5, all tuning parameters, including T and U, are selected by cross-validation. In our experience, the cross-validation generally provides good forecast accuracy.

In practice, however, it is assumed that past consumption also depends on past weather, such as daily temperature. For example, suppose that it was exceptionally hot yesterday and it is cooler today. In this case, it is not desirable to directly use past consumption as the predictor; it is better to remove the effect of past temperature from past consumption. In other words, we can use y ( i − t − L α ) j − b ( i − t − L α ) j and y i ( j − u − L β ) − b i ( j − u − L β ) instead of y ( i − t − L α ) j , and y i ( j − u − L β ) , respectively. As a result, μ _ij is expressed as follows: μ i j = ∑ t = 1 T α j t ( y ( i − t − L α ) j − b ( i − t − L α ) j ) + ∑ u = 1 U β j u ( y i ( j − u − L β ) − b i ( j − u − L β ) ) . (7)

Substituting (5) into (7) results in the following: μ i j = ∑ t = 1 T α j t y ( i − t − L α ) j − ∑ t = 1 T ∑ m = 1 M ∑ q = 1 Q α j t γ q m h q ( j ) g m ( s i − t − L α ) + ∑ u = 1 U β j u y i ( j − u − L β ) − ∑ u = 1 U ∑ m = 1 M ∑ q = 1 Q β j u γ q m h q ( j − u − L β ) g m ( s i ) . (8)

The appropriate values of α _jt and β _ju are chosen by several approaches. A simple method is α _jt = 1/T and β _ju = 1/U, which implies μ _ij is the sample mean of the past consumption. Another method is based on the AR(1) structure, i.e., α j t = ρ α t and β j u = ρ β u , where ρ _α and ρ _β satisfy ∑ t = 1 T ρ α t = 1 and ∑ u = 1 U ρ β u = 1 , respectively. Note that ∑ t = 1 T ρ t = ρ ( 1 − ρ T ) / ( 1 − ρ ) , so ∑ t = 1 T ρ t = 1 is equivalent to ρ ^T+1 − 2ρ + 1 = 0, whose numerical solution is easily obtained.

By combining the expressions of b _ij in Eq. 5 and μ _ij in Eq. 8, the model Eq. 2 is expressed as follows: y i j = 1 T ∑ t = 1 T y ( i − t − L α ) j − 1 T ∑ t = 1 T ∑ m = 1 M ∑ q = 1 Q γ q m h q ( j ) g m ( s i − t − L α ) + 1 U ∑ u = 1 U y i ( j − u − L β ) − 1 U ∑ u = 1 U ∑ m = 1 M ∑ q = 1 Q γ q m h q ( j − u − L β ) g m ( s i ) + ∑ m = 1 M ∑ q = 1 Q γ q m h q ( j ) g m ( s i ) + ε i j . (9)

The model Eq. 9 is equivalent to the linear regression model y ̃ = X γ + ε , (10) where γ = vec(Γ) and ɛ = vec( E ) with E = (ɛ _ij ). Here, X and y ̃ are considered as the design matrix and the response vector, respectively. The definitions of y ̃ and X are given in Section S1 in the Supplementary Material.

To estimate the regression coefficient vector γ , one can use the least squares estimation (LSE) min γ ‖ y ̃ − X γ ‖ 2 2 .

In our experience, however, the elements of least squares estimate γ ̂ often become negative. In such cases, the estimate of b _ij is negative because the basis functions h _q (j) and g _m ( s _i ) generally take positive values. When b _ij < 0, one can interpret the weather effect is negative. Nevertheless, the one-day curve of b ̂ i j turns out to be counterintuitive; the effect of weather negatively increases as the electricity demand increases. In other words, b _ij is negatively large at working hours and small at midnight. As a result, μ _ij becomes extremely large in the working time. We observe this phenomenon on the analysis of three datasets presented in this study; one of them is the well-known Global Energy Forecasting Competition 2014 (GEFCom 2014) data. Therefore, this phenomenon could occur in other datasets.

A clear interpretation is realized when the weather effect b _ij is nonnegative. To achieve this, we employ the nonnegative least squares (NNLS) estimation, in which we minimize the loss function under a constraint the regression coefficients are nonnegative: min γ ‖ y ̃ − X γ ‖ 2 2 s u b j e c t t o γ ≥ 0 . (11)

The optimization problem Eq. 11 is a special case of quadratic programming with nonnegativity constraints [e.g., (Franc et al., 2005)]. As a result, the NNLS problem becomes a convex optimization problem. Several efficient algorithms to obtain the solution in Eq. 11 have been proposed in the literature [e.g., (Lawson and Hanson, 1995; Bro and De Jong, 1997; Timotheou, 2016)].

We add the ridge penalty (Hoerl and Kennard, 1970) to the loss function of the NNLS estimation: min γ ‖ y ̃ − X γ ‖ 2 2 + λ ‖ γ ‖ 2 2 s u b j e c t t o γ ≥ 0 , (12) where λ > 0 is a regularization parameter. In our experience, the ridge penalization improves the forecast accuracy, and also produces a smoother function of b ̂ i j than the unpenalized NNLS, which makes the interpretation of the weather effect easier.

For the day-ahead forecast, we forecast the demand on the next day, y ̂ ( i + 1 ) j , for the given NNLS estimate γ ̂ and weather information s _i+1. The forecast value y ̂ ( i + 1 ) j is expressed as follows: y ̂ ( i + 1 ) j = μ ̂ ( i + 1 ) j + b ̂ ( i + 1 ) j .

Here, b ̂ ( i + 1 ) j = ∑ m = 1 M ∑ q = 1 Q γ ̂ q m h q ( j ) g m ( s i + 1 ) and μ ̂ ( i + 1 ) j = ∑ t = 1 T α j t ( y ( i + 1 − t − L α ) j − b ̂ ( i + 1 − t − L α ) j ) . On the intraday forecast, we may use information about the demand on that day so that μ ̂ ( i + 1 ) j is expressed as μ ̂ ( i + 1 ) j = ∑ t = 1 T α j t ( y ( i + 1 − t − L α ) j − b ̂ ( i + 1 − t − L α ) j ) + ∑ u = 1 U β j u ( y ( i + 1 ) ( j − u − L β ) − b ̂ ( i + 1 ) ( j − u − L β ) ) .

Construction of a forecast interval based on Eq. 11 or Eq. 12 is not easy due to the constraints of the parameter. To derive the forecast interval, we employ a two-stage procedure; first, we estimate the parameter via NNLS to extract variables that correspond to nonzero coefficients. Then, we employ the least squares estimation based on the variables selected in the first step. With this procedure, we should derive the forecast interval after model selection. To achieve this result, the post-selection inference (Lee and Taylor, 2014; Lee et al., 2016) is employed. The post-selection inference for the NNLS estimation is detailed in Section S2 in the Supplementary Material.

We use maximum temperature in Tokyo as the weather variable s _i , available at Japan Meteorological Agency (https://www.jma.go.jp/jma/indexe.html).

Figure 1 shows the relationship between the maximum temperature and demand. We depict this relationship in different time intervals for Sunday, Monday, and Wednesday. The curves are depicted by fitting the ordinary least squares estimation with the cubic B-spline function (e.g., Hastie et al., 2008) with equally-spaced knots. The number of basis functions is M = 5. The B-spline function of k-degree is constructed by linear combination of M basis functions with M + 2k knots, say t ₁, … , t _M+2k . For given n data points in increasing order, say x ₁, … , x _n , we set the knots as t ₁ ≤ t ₂ ≤ ⋯ ≤ t _k = x ₁ ≤ ⋯ ≤ x _n = t _M+k+1 ≤ ⋯ ≤ t _M+2k . The B-spline function is the density function of a convolution of uniform random variables (Hastie et al., 2008), and it is constructed by the following recurrence relation (De Boor, 1972): B j ( x ; 0 ) = 1 [ t j , t j + 1 ) ( x ) , B j ( x ; m ) = x − t j t j + m − t j B j , m − 1 ( x ) + t j + m + 1 − x t j + m + 1 − t j + 1 B j + 1 , m − 1 ( x ) , where 1 A ( x ) is the indicator function, that is, 1 A ( x ) = 1 when x ∈ A and 1 A ( x ) = 0 otherwise. We adopt the cubic B-spline (k = 3). The left panel of Figure 2 depicts the cubic B-spline function with M = 10. For all settings, the curve-fitting works well, which implies the usage of the B-spline function as a basis function g( s _i ) would be reasonable.

The curve shape for 9:00–10:00 is different from that for 15:00–16:00, which suggests that the regression coefficients must be different among time intervals. Meanwhile, the curve shapes for 15:00–16:00 and 18:00–19:00 are similar. In this case, it is reasonable to assume that the regression coefficients for 15:00–16:00 are similar to those for 18:00–19:00. The regression coefficients are assumed to change depending on the time interval yet to be smooth according to the time interval. To achieve a smooth function of β _m (j), we also use the cubic B-spline as a basis function of h _q (j). As the ordinary B-spline function cannot produce a smooth curve around the boundary (i.e., 23:00–0:00 and 0:00–1:00), we employ the cyclic B-spline function, where the basis functions wrap at the first and last knot locations. The cyclic B-spline function is implemented in the cSplineDes function in the mgcv package in R. The right panel of Figure 2 shows the cyclic B-spline function of degree 3 with M = 10.

We also observe that the curve shapes are different among each day of the week. Therefore, we construct the statistical models by day of the week separately; seven statistical models are constructed. To forecast the demand, we select a model that matches the day of the week. All national holidays are regarded as Sunday; therefore, the number of observations on Sunday is larger than on other weekdays.

We employ our proposed method based on two estimation procedures: ridge estimation min γ ‖ y ̃ − X γ ‖ 2 2 + λ ‖ γ ‖ 2 2 , and NNLS estimation with the ridge penalty in Eq. 12. We label these estimation procedures as “LSE” and “NNLS,” respectively. For both LSE and NNLS, we prepare a wide variety of statistical models by changing the tuning parameters. Here, we review the role of each tuning parameter and present their candidates as follows:

• Q: the number of basis functions in the varying coefficient model. As Q increases, the electricity fluctuations caused by weather becomes large in time interval j (j = 1, … , J). As we will show in Section 4.3, the value of Q may not highly affect the forecast accuracy. Therefore, it would be sufficient to set the candidates of Q as Q = 5, 10.

In addition to the above candidates of models, we consider two types of α _jt : sample mean of the past demand (i.e., α _jt = 1/T) and AR(1) structure. Details of the AR(1) structure are presented at the end of Section 2.2. As a result, the total number of candidates of the models is 336 (= 2 × 2 × 2 × (20 + 1) × 2). These candidates are summarized in Table 1.

The forecast accuracy of our proposed method depends on the tuning parameters presented above. We investigate the impact of tuning parameters on the mean absolute percentage error (MAPE) defined as M A P E = 1 J n ∑ i = 1 n ∑ j = 1 J | y i j − y ̂ i j | y i j . (13)

Figure 3 shows the relationship between regularization parameter λ and MAPE. The relationships are investigated for all combinations of (Q, M, T).

The result shows that the AR(1) structure on α _jt produces smaller MAPE than the mean structure (α _jt = 1/T) for all candidates of ridge parameter λ. For AR(1) model, the performance for T = 4 is better than that for T = 2. We observe that a large value of λ may result in poor forecast accuracy; meanwhile, a small amount of λ generally results in good accuracy. The result for λ = 0 is not displayed here because the MAPE can be extremely large due to the non-convergence of parameters; the ridge penalty helps avoid such a non-convergence. In summary, the AR(1) structure on α _jt , T = 4, and a small value of λ will result in a small MAPE for this dataset.

We compare the performance of our proposed method with the following popular machine learning techniques: random forest (RF), support vector machine (SVM), Least absolute shrinkage and selection operator (Lasso), LightGBM (LGBM; (Ke et al., 2017)), and Long Short Term Memory [LSTM; (Hochreiter and Schmidhuber, 1997; Van Houdt et al., 2020)]. The RF, SVM, LGBM, and LSTM involve nonlinear functions, allowing us to capture complex structures. We use R packages randomForest, ksvm, glmnet, and lgbm to implement RF, SVM, Lasso, and LGBM, respectively. The LSTM is implemented by tensorflow.keras library in python. The LGBM is based on the gradient boosting decision tree and achieves good forecast accuracy in various fields of research, including energy (e.g., (Ju et al., 2019)). The LSTM is the state-of-the-art technique for analyzing time series data using artificial neural networks (ANN). For these machine learning techniques, the forecast is made by the electricity demand of past T days and maximum temperature, that is, y ̂ i j = f ( s i , y ( i − 1 ) j , … , y ( i − T + 1 ) j ) .

The characteristics of our proposed method and the machine learning techniques are compared in Table 2.

A detailed description of Table 2 is as follows:

• With small sample sizes, the machine learning methods that capture nonlinear structure provide good accuracy except for the LSTM; in general, the ANN can often yield poor accuracy with small sample sizes (e.g., Tange et al., 2017). The accuracy of the lasso is not sufficiently but relatively high because the lasso cannot capture the nonlinear structure. When the number of observations is large, the LGBM and LSTM will perform perfectly well because these methods capture complex nonlinear structures. However, the sample sizes of the dataset used in this Section (and also other Sections) may not be large.

In practice, we need to select a set of tuning parameters among candidates. The candidates of the tuning parameters for the proposed method are presented in Section 4.3. For machine learning methods, the candidates are detailed in Section S3 in the Supplementary Material. For all methods, we must select an appropriate set of tuning parameters. These tuning parameters are generally selected with cross-validation. In this study, we choose these tuning parameters to minimize the MAPE in the past year. The values of tuning parameters are changed every month.

Table 3 shows the monthly MAPE for our proposed method and existing methods from April 2019 to March 2020. The result shows that the proposed method performs better than existing machine learning techniques in total. The SVM achieves the best performance in specific months but sometimes provides much larger MSE than other methods, resulting in worse performance than our proposed method in total. The RF and LGBM yield similar performance, and these methods produce larger MAPE than our proposed method. The lasso cannot capture the nonlinear relationship between temperature and demand as in Figure 1, resulting in poor performance in total. The LSTM provides the worst performance, probably due to the small number of observations: it is from n = 132 to n = 302. In general, the ANN does not perform well with such small sample sizes; (Tange et al., 2017) showed that the ANN performed worse than SVM for small sample sizes, and several similar numerical results are found in the literature (Bećirović and Ćosović, 2016; Dehalwar et al., 2016; Maldonado et al., 2019). We observe that the NNLS and LSE result in similar values of MAPE. These similar results are obtained probably due to the identifiability of the parameter. Our proposed procedure decomposes y _t as y _t ≈ μ _t + b _t , and this decomposition is not unique. Therefore, the NNLS and LSE produce different values of μ ̂ t and b ̂ t but similar values of μ ̂ t + b ̂ t .

We also compute the MAPE for well-known Global Energy Forecasting Competition 2014 (GEFCom 2014) data (Hong and Fan, 2016), and obtain a similar result as that shown in Table 3. For detail, please refer to Section S4 in the Supplementary Material.

With our proposed method, the estimated model can be interpreted by decomposing the forecast value by the effects of temperature and past demand: y ̂ i j = b ̂ i j + μ ̂ i j . The values, b ̂ i j and μ ̂ i j , for NNLS and LSE from October 1st to 14th, 2019, are depicted in Figure 4.

Although the forecast accuracy of the LSE estimation is similar to that of the NNLS estimation, as shown in Table 3, the results of the decomposition y ̂ i j = b ̂ i j + μ ̂ i j for these two methods are unalike in terms of values. The LSE often results in negative values of b ̂ i j , and then μ ̂ i j becomes substantially larger than the actual demand. In particular, the value of b ̂ i j at working hours is negatively larger than that at midnight; on the other hand, the electricity is used primarily during working hours. As a result, the behavior of b ̂ i j is counterintuitive; b ̂ i j becomes negatively larger as y ̂ i j increases. Thus, interpreting the effect of weather turns out to be difficult with LSE. This issue occurs because there are no restrictions on the sign of b ̂ i j . We observe a similar behavior of b ̂ i j on other datasets, including GEFCom 2014.

In contrast, the effect of weather is appropriately captured by the NNLS estimation. Indeed, in many cases, the value of b ̂ i j at working hours is positively larger than that at midnight. The constraint on nonnegativeness of γ significantly improves the interpretation of the weather effect while still maintaining excellent forecast accuracy.

To further investigate the effect of temperature, we depict b ̂ i j on Tuesday when maximum temperatures are 5°C (cold day), 20°C (cool day), and 35°C (hot day), which is shown in Figure 5. For existing methods, it is difficult to depict b ̂ i j because these methods do not assume the existence of b ̂ i j . For existing methods, instead of b ̂ i j , the weather effect is calculated depending on the input: we forecast the demand with specific temperature (5°C, 20°C, and 35°C) and average annual temperature (21.7°C in this case), and subtract the forecast value with average annual temperature from that with a specific temperature. We depict the weather effect on October 8th, 2019, and October 15th, 2019, for existing methods; both of these days are Tuesday.

The results of our proposed procedure show that b ̂ i j is highly dependent on the temperature: the weather effect may be substantial for cold and hot days due to air conditioner use. For NNLS, the weather effect is always positive, which allows clear interpretation compared with LSE. Both NNLS and LSE result in smooth curves due to the smooth basis function of the regression coefficients. We observe that the weather effect is stable unless M is large and λ is excessively small.

The existing machine learning methods, RF, SVM, LGBM, and LSTM, are unstable and highly depend on the input. With the Lasso, the impact of temperature is almost zero, which implies that the weather effect cannot be captured. Our proposed method results in smoother curves than existing methods thanks to the B-spline basis expansion of regression coefficients in Eq. 4. As a result, our proposed method is more suitable for interpreting the weather effect than existing methods.

In Japan, the most frequently-used daily information about COVID-19 is the daily number of infections, available at https://www3.nhk.or.jp/news/special/coronavirus/data-all/. However, this variable turns out to be relatively unstable. In order to use more stable information about COVID-19, we may use the moving average of the number of infections; that is, the average number of infections in the past several days.

The research facility used in this study promotes the staff to work from home to prevent the spread of COVID-19. We expect that the electricity usage of the research facility decreases as the number of infected people of the COVID-19 increases; therefore, the COVID-19 information may affect the forecast results. Figure 6 shows the relationship between the number of daily infections and demand, and between the average number of infections in the past 14 days and demand. We depict these plots on working days. The curve fitting is done by polynomial regression with cubic function. The B-spline may not be suitable in this case due to the limited number of observations.

The number of daily infections includes outliers; thus resulting in unstable curve fitting. For example, when the number of infections is relatively large, the demand increases as the number of infections increases. This phenomenon is counterintuitive because this facility encourages working at home as much as possible to prevent infections. On the other hand, the average number of infections in the past 14 days is negatively correlated with the demand, and the curve fitting of the cubic function works well. For this reason, we use the average number of infections in the past 14 days as daily information about COVID-19. As a result, the daily variable s _i becomes a two-dimensional vector that consists of the maximum temperature and the average number of infections in the past 14 days. We use the cubic B-spline function as a basis function of maximum temperature, and the cubic function as a basis function of the average number of infections in the past 14 days.

In the previous real data example, we construct the statistical models by day of the week separately. However, in this case, the number of observations affected by COVID-19 is excessively small. Thus, only two statistical models are constructed based on working days (from Monday to Friday) and holidays (Saturday, Sunday, and national holidays). We consider the problem of the day-ahead forecast from March 1st, 2020 to May 30th, 2020 (data in 2017–2019 are only used for training). The training data consist of all demand data up to the previous day of the forecast day.

We investigate how the input of COVID-19 information improves the forecast accuracy. The NNLS does not produce negative regression coefficients but we observe that the electricity usage decreases as the average number of infections in the past 14 days increases, as shown in Figure 6. Thus, the electricity fluctuation caused by COVID-19, say b i j covid 19 , should be negative. As the basis function, we use the cubic function of negative number of infections in the past 14 days; that is, b i j covid 19 = β 1 ( j ) ( − s i covid 19 ) 3 + β 2 ( j ) ( − s i covid 19 ) 2 + β 3 ( j ) ( − s i covid 19 ) , where s i covid 19 is the number of infections in the past 14 days.

Figure 7 shows the forecast values with our proposed method using NNLS and the actual demand. We depict two forecast values: one uses the information about COVID-19 and the other does not. The result shows that the usage of the COVID-19 data slightly improves the accuracy; for example, on the 10th and 13th.

Figure 8 shows b ̂ i j of both temperature and COVID-19. The results show that after April 10th, the negative effect of COVID-19 is observed. The government declared a state of emergency on April 7th, and following this, the usage pattern of the electricity changed. The change in electricity pattern on the 8th and 9th may not be captured due to excessively small sample sizes related to the number of infections. However, after the 9th, the effect of COVID-19 is captured. (insert Figure 8 here).

Table 4 shows MAPE for our proposed method and existing methods on March, April, and May 2020. We also compare the performance of two estimation procedures: one uses the information about COVID-19 and the other does not. In March, the result shows that the proposed method performs the best and Lasso yields the worst performance. This is because the effect of temperature can be appropriately captured as shown in the previous example. The effect of COVID-19 is not crucial in March because the performance becomes worse when the COVID-19 information is included.

In April and May, the information about COVID-19 significantly improves the performance for almost all methods. In particular, our proposed method and Lasso both produce small MAPE values, and interestingly, Lasso performs slightly better than our proposed method. There are two reasons to explain the superiority of Lasso. First, the relationship between the number of infections and the demand is approximated by a linear function, as shown in Figure 6. Second, the temperature effect is not crucial because the temperature does not change much in April and May. The nonlinear machine learning techniques, that is, SVM, RF, LGBM, and LSTM yield large MAPE in April, probably due to the small number of observations related to COVID-19; these techniques result in overfitting. In particular, LSTM produces much worse MAPE than other methods. Such poor results may be due to the instability of the electricity demand; the neural networks can be highly sensitive to noise [e.g., (Goodfellow et al., 2014)]. Indeed, we observe that the electricity demand of this dataset turns out to be unstable compared with the datasets from Tokyo electric power company holdings and Global Energy Forecasting Competition 2014. To improve the accuracy of the LSTM, we need to change the input/output by incorporating much more information of the data (e.g., not daily but hourly weather and hourly COVID-19 data).