Novel Multimodal Data for Enhanced Electricity
Spot Price Forecasting Using A CNN-LSTM
Ensemble Learning Model for the Japan ElectricPower eXchange (JEPX) Spot Market
Ziyang Wang, Masahiro Mae, Ryuji Matsuhashi
Dept. of Electrical Engineering and Information Systems
The University of Tokyo
Tokyo, Japanwang-ziyang@ieee.org, mmae@ieee.org, matu@enesys.t.u-tokyo.ac.jp

Abstract—Electricity price forecasting (EPF) is critical inenergy markets, particularly with the rise of renewable energy
sources (RES) in Japan, which can cause day-ahead spot prices
to drop to nearly zero JPY/kWh, impacting retailer profitability.
This study demonstrates that the Convolutional Neural NetworkLong Short-Term Memory (CNN-LSTM) model outperforms the
LSTM model in both prediction accuracy and computational
efficiency in Japan’s electricity spot market. A novel ensemble
learning strategy enhances both the robustness and accuracyof the EPF model is proposed. Novel multimodal explanatory
variables, including electricity spot price, system price, discretionary cost, sell bid amount, buy bid amount, total contracted
volume, actual power generation, actual solar power generation,
meteorological forecasts and calendar forecasts, alongside therolling features of spot prices, are utilized and verified. Furthermore, a ”policy-versus-policy” approach addresses the zeroinflated regression issue of the zero price prediction is proposed.
Our model, with a comprehensive feature integration, achieves
an RMSE of 5.66 JPY/kWh and an R2 of 0.729 during the test
range from 2022.01.01 to 2022.12.31. The paper also introduces anovel method for estimating confidence intervals using ensemble
learning.
Index Terms—electricity price forecasting (EPF), renewableenergy source (RES), CNN-LSTM, ensemble learning, zero electricity spot price, confidence interval

I. I NTRODUCTION
In recent years, the increased penetration of RES, particularly from wind and solar sources, has introduced an unprecedented phenomenon in the US and EU electricity markets:
negative electricity spot prices [1], [2]. However, the situationin Japan presents a distinct narrative. With the significant
introduction of solar and wind energy into the Japanese power
grid, the electricity wholesale spot markets have witnessed
prices approaching nearly zero, specifically at 0.01 JPY/kWh(referred to as zero prices in the following sections of thispaper). The electricity spot prices in Kyushu region, Japan, is
depicted in Fig. 1(a). A closer look at the intermittent zero
prices is presented in Fig. 1(b). As evident from Fig. 1(a), theemergence of considerable zero prices became prominent fromthe year 2020 onwards, largely attributable to the rapid integration of RES. In this evolving energy landscape, day-ahead

electricity price forecasting (EPF) in electricity spot markets
has emerged as a paramount concern [3]. Statistical models
such as the Autoregressive Moving Average (ARMA) [4],[5] and Autoregressive Integrated Moving Average (ARIMA)[6]–[8] have been widely employed in EPF studies. Whilethese models provide foundational approaches, their inherent
linearity can pose challenges. The increased integration of
RES, along with factors like demand fluctuations, introduces
non-linear trends and sudden price anomalies that traditional
statistical models may struggle to capture accurately.
Moreover, as is shown in Fig. 1(b), after data normalization,
the considerable durations of zeros resulted by the zero prices
in the target variable create a zero-inflated regression problem
in machine learning. However, it is almost impossible formost machine learning models, including Random Forest(RF), Support Vector Regression (SVR), and neural networks,
to continuously output zeros. Traditional solutions for zeroinflated regression problem usually requires two models: oneclassification model to identify zero values and another regression model for non-zero values, which doubles the training
time and cost.
To address the challenges mentioned above, this paper’s
objective is to utilize multimodal data as novel features toenhance electricity spot price forecasting (EPF) in Kyushu region, Japan, using LSTM and CNN-LSTM prediction models.
In addition, a novel ensemble learning approach is leveraged
to further enhance the prediction accuracy of LSTM and CNNLSTM due to the their inherent uncertainties. The performanceof the LSTM and CNN-LSTM in EPF is compared in termsof both prediction accuracy and computation time. Moreover,to shoot the zero-inflated problem in the JEPX spot market, anovel ”policy versus policy” strategy is employed to forecast
the zero prices to half the computation time which traditional
two-stage method requires. Furthermore, a natural logarithm
transformation is utilized to improve the spot price’s Skewness
and Kurtosis to enhance prediction accuracy. Moreover, anovel method for extracting meteorological forecast data by
utilizing Google Maps is introduced. Finally, an implication

Electricity spot price (yen/kWh)Electricity spot price (yen/kWh)

250
200
150
100
50
0

models underwent multiple training iterations. Subsequently,
all individual predictions were aggregated using a simple
averaging method to construct the final ensemble prediction,
as depicted in (1), where N represents the total numberof predictions, and k denotes the index of each individual
prediction.

Kyushu electricity spot price

2016

2017

2018

2019

2020
Time [Yearly]D

2021

2022

2023

Kyushu electricity spot price (Zoomed)N40ŷensemble =

30
20

1 X
ŷk
N

(1)

k=1

10
0

2022-04-01

2022-04-02

2022-04-03

2022-04-04
Time [Daily]E

2022-04-05

2022-04-06

For clarity, the pseudo-code for the ensemble learning
procedure is outlined in Algorithm 1.

2022-04-07

Figure 1. Kyushu electricity spot price [JPY/kWh] (a) andzooming-in zero-inflated prices (b).

Algorithm 1 Ensemble Learning Procedure
1: Perform log transformation of the electricity spot priceusing (9).
of a novel estimation method for the confidence interval of
EPF is demonstrated.
II. M ETHODOLOGY
A. LSTM and CNN-LSTM forecasting models

LSTM Layer(M5 × 1)

LSTM Layer(M5 × 1)

LS
TM

LS
TM

A LSTM model and CNN-LSTM model were designed
and employed for EPF and for comparison using the Python
Tensorflow keras library. The architectures of the LSTM and
CNN-LSTM models are delineated in Fig. 2, the hyperparameters were selected empirically based on optimal performance.
Input Layer(N × F )FC Layer(M6 × 1)

M6

48

LS
TM

F features

Fo
re
ca
st
e

LS
TM

LS
TM

d

LS
TM

M5

pr
ic
e

Past N -step data

M5

(a)
LSTM Layer(M5 × 1)

Max Pooling Layer(M1 /2 × M2 /2)

Convolution Layer(M3 × M4 )M3

FC Layer(M6 × 1)

Output Layer(48 × 1)

M1 /2
Past N -step data

N

48

LS

M2 /2
M4

F features

ec
Fo
r

LS

TM

as

te

d

M2

M6

TM

M5

e

M1

pr
ic

Input Layer(N × F )LS
TM

Convolution Layer(M1 × M2 )data normalization).
the exponential transformation to the predictedvalues using (10) (reverse of the log transformation).
8: Calculate the ensemble prediction ŷensemble using (1).
9: Calculate the zero price for ŷi and ŷensemble using (2).
7: Apply

C. ”Policy-versus-policy” zero prices forecasting strategy

Output Layer(48 × 1)N

2: Normalize the training and test data.
3: for i = 1 to 30 do
4:
Train the model to generate prediction ŷi .
5: end for
6: Restore the predicted values to their original scale (reverse(b)

Figure 2. Schematic of the architectures of the (a) LSTM and(b) CNN-LSTM models.
B. Ensemble learning strategy
Given the inherent variability of neural network models dueto their sensitivity to initial conditions and the stochastic natureof their training, training the same neural network multiple
times and averaging the predictions can mitigate individualmodel errors, leading to enhanced prediction performance, asdifferent models will not make identical errors on the test set[9], [10]. Based on this understanding, an ensemble learning
approach was implemented. The CNN-LSTM and LSTM

In this study, we introduce a novel method to address thezero-inflated regression problem of the EPF in JEPX spot
market. The solution begins with understanding the broader
trends in global electricity spot markets. As highlighted bySeel et al. [1], an abundance of RES can lead to negativeelectricity spot prices in the US and EU. Drawing from this,
we infer that negative pricing is a natural consequence of
RES abundance. In Japan, policy dictates that electricity spot
prices cannot drop below 0.01 JPY/kWh, thereby preventing
them from turning negative. Assuming that the circumstancesleading to negative prices in the US and EU are similar to thosein Japan, it is reasonable to infer that the explanatory variables
in both scenarios would exhibit similar patterns. Feedingthese Japanese explanatory variables into a machine learning
regression model would naturally produce negative prices, as
the model is not constrained by Japan’s minimum pricing
policy. Hence, by leveraging this premise, zero prices can beforecast by translating any negative outputs from the model
to zeros, as indicated in Algorithm 1. This approach effectively functions as a policy-versus-policy forecasting strategy,reflecting real-world conditions. Equation (2) illustrates thezero price calculation procedure.
ŷi = max(0, ŷi )(2)

D. Performance evaluation
Performance metrics, including root mean squared error(RMSE) and the coefficient of determination (R2 ), were utilized for evaluation. The computational formulae for RMSE
and R2 are specified in (3)-(4), respectively.
r Pn
RMSE(y, ŷ) =2
t=1 (yt − yˆt )

t
1 X
yi
w i=t−w+1

(7)

t
X
1
(yi − Meant )2
w − 1 i=t−w+1(8)

Meant =vuu
Stdt = t

(3)

n

Pn
(yt − yˆt )2
R (y, ŷ) = 1 − Pt=1n2
t=1 (yt − ȳ)
27-day-long historical data (05:00 - 05:00]

III. DATA PREPARATIONInterval(19 h)

Actual power generation (kyushu), actual solar power generation(kyushu)This study utilizes multimodal data for enhanced EPF,
including the actual total power generation, actual solar power
generation, spot price rolling features (minimum, maximum,
mean, and standard deviation), system spot price, discretionary
cost, sell bid amount, buy bid amount, total contracted volume,
meteorological forecast data, and calendar forecast data. The
overall data architecture and the corresponding time frame areillustrated in Fig. 3. The input data are segmented by theirtemporal deley into three green blocks. To maintain consistency with the Japan Electric Power Exchange (JEPX) spot
price data, all input data were linearly interpolated to a time
resolution of 30 min. A 7-day moving window was applied
to the input data before being fed into the CNN-LSTM and
LSTM models. In the JEPX spot market, all transactions mustbe finalized by the bidding deadline of 10:00 JST. To ensure an
accurate and comprehensive prediction process, ample time isallocated for the execution and potential refinements of neuralnetwork calculations. Given the computational requirements
and the complexities of the forecasting process, a 5-h buffer
before the deadline has been established. The forecasting time
point is set to 05:00 JST, covering the entire following dayfrom 00:00 JST to 23:30 JST, which includes 48 time frames
in total.
In this study, since photovoltaic (PV) power is the primary
RES in the Kyushu region due to its substantially greater installed capacity compared to wind power, features pertaining towind power are not incorporated into the current investigation.
A. Feature engineering
Feature engineering was conducted using a method thatinvolves calculating rolling statistics (minimum, maximum,
mean, and standard deviation) of the electricity spot pricewith a rolling window. This method is crucial for capturing
temporal patterns and trends, which are essential for timeseries feature extraction. The length of the rolling window
was chosen to be 24 h, i.e., w = 48. The equations for thesecalculations are provided in (5)-(8).
Mint = min(yt−w+1 , yt−w+2 , ..., yt )(5)

Maxt = max(yt−w+1 , yt−w+2 , ..., yt )(6)

JEPX spot
market bidding
deadline

Forecasting point(5 h ahead)(4)

Prediction period[00:00 - 00:00)

7-day-long forecast data ① [00:00 - 00:00)
Spot price (kyushu), spot price rolling features (kyushu), system spot
price, discretionary cost (kyushu), sell bid amount, buy bid amount, total
contracted volume

7-day-long forecast data ② [00:00 - 00:00)
Meteorological forecast data, calendar data

05:00

00:00

00:00

・・・・・・

05:00 10:00

00:00

00:00

Figure 3. Illustration of the data architecture with a 30-min
time interval, highlighting the time delays among different
data.
B. Electricity data
The electricity spot prices (JPY/kWh) at a 30-min resolution, including the Kyushu regional spot price, system spot
price, discretionary cost, sell bid amount, buy bid amount, andtotal contracted volume, from 00:00 JST to 23:30 JST (over a
7-day period), were downloaded from the Japan Electric Power
Exchange (JEPX) [11].
The distribution of the Kyushu region electricity spot prices,
as illustrated in Fig. 4(a), exhibits pronounced Skewness
and Kurtosis, indicating a non-normal distribution. Neuralnetwork models typically assume that input data are normally
distributed or, at least, exhibit symmetry in their distribution since this facilitates the model’s learning process byproviding a standardized scale of input features. Deviations
from normality, such as Skewness and Kurtosis, can introduce
biases in the model’s predictions and affect the efficiency ofthe learning algorithm. To mitigate these effects, a natural
logarithm transformation was applied to the electricity spot
prices, as detailed in Fig. 9. This transformation is a common
technique in statistical normalization that reduces the impact
of Skewness by compressing the scale of the distribution,
thereby enhancing symmetry and reducing the influence ofoutliers. The effectiveness of the logarithmic transformationis quantitatively evidenced by the reduction in Skewness andKurtosis of the price data, as outlined in Table I. Following themodel’s prediction output, an exponential back-transformation,
defined in (10), is applied to convert the forecasted values backto their original scale.
y = loge (y + 1)

(9)

y = ey − 1(10)Frequency

0.3

color range are identified as land. Subsequently, Fig. 5(b)showcases the delineated land areas from Fig. 5(a), with red
dots marking these regions. After the identification of land,
the meteorological data of all the land area were averaged as
input features.

0.2
0.1
0.0

0

50
100
150
200
Electricity spot price [yen/kWh]D

250

Frequency

0.06

D

E

0.04

Figure 5. Kyushu area map (a) and the identified land areas(b). Map data: ©2024 Google, TMap Mobility.

0.02
0.00

0
1
2
3
4
5
Natural logarithm-transformed electricity spot priceE

Figure 4. Distribution of Kyushu region’s electricity spot
prices (a) and the corresponding distribution after a natural
logarithm transformation (b).
TABLE I. Skewness and Kurtosis of the original electricity
spot price and the natural logarithm-transformed electricity
spot price.

Skewness
Kurtosis

Original

Log-transformed

8.77
123.00-0.68
3.81

The actual electricity generation data (MW) at a 30-min
resolution in the Kyushu region, including the actual power
generation and actual solar power generation over a 7-day
period from 05:00 JST to 05:30 JST, were downloaded fromthe Organization for Cross-regional Coordination of Transmission Operators, Japan (OCCTO) [12]. These data were used
as input features and are depicted in the leftmost green blockin Fig. 3.

D. Calendar forecast data
The calendar forecast data, by integrating the cyclic data andculturally significant Japanese holidays, aimed to enhance the
EPF accuracy by accounting for the unique electricity consumption behaviors associated with these special occasions,covering a 7-day period from 00:00 JST to 23:30 JST, wereused as features, as illustrated in the rightmost green block in
Fig. 3.
1) Cyclic data: In this study, we employ sine and cosine transformations with periods corresponding to common
cyclical patterns: 1 day, 1 week, 3 months, and 1 year. Suchtransformations have been acknowledged in foundational timeseries literature as a robust method to encapsulate periodic
patterns in data [16]–[18]. The input features consist of thesine and cosine waves representing 1 day, 1 week, 3 months(seasonal effects), and 1 year. The sine and cosine waves witha period of 1 year is shown in Fig. 6.
sin_yearcos_year

1.0
0.5
0.0

C. Meteorological forecast data

0.5

The meteorological forecast data of Kyushu region, including the wind speed, air temperature, relative humidity, cloudiness, precipitation, solar radiation, were downloaded fromthe Japan Meteorological Business Support Center (JMBSC)[13], as illustrated in the rightmost green block in Fig. 3.
The Kyushu region map was a rectangular screenshot from
Google Maps Styling Wizard [14]. The identification of theland area from the Kyushu region map was conducted using
the Python OpenCV library [15] applied to the screenshot
image of Kyushu region from Google Maps. Fig. 5(a) displays
the original screenshot of the map for the Kyushu region.
In this step, pixels that are not within the defined blue

1.0
2016 2017 2018 2019 2020 2021 2022 2023
Time [year]

Figure 6. Sine and cosine waves with a period of 1 year.
2) Holiday data: The holidays and non-holidays wereencoded into numerical labels 1 and 0, respectively.
E. Prediction approach
A One-time prediction schematic, shown in Fig. 7(b), isused to validate the proposed features as it requires muchless training time compared with the day-by-day predictionschematic shown in Fig. 7(a). In this scenario, the training
data spaned from April 1, 2016 — the date marking the full
liberalization of the electricity retail market — to December
31, 2021. The testing data encompass the period from January
1, 2022 to December 31, 2022. All the features proposed
in Section III has been validated for prediction accuracy
enhancement.
After validating the proposed features using the test data
which spans the entire 2022 year, the day-by-day prediction
approach was conducted using all proposed features to explore
the actual prediction accuracy in real-application scenarios, as
illustrated in Fig. 7(a), which shares the same test range inthe one-time prediction schematic. The training of the models
was conducted using two NVIDIA Quadro RTX 8000 GPUs
using the Python keras package in Windows OS.

TABLE III. Ensemble learning metrics for the day-by-day
prediction approach using all proposed features.
Metrics

Data set
Training data

The proposed ensemble learning approach effectively reduces
the variability inherent in individual model predictions. Although some individual predictions may diverge markedly
from the actual electricity spot prices, the ensemble method,
described in (1), averages these forecasts to produce a morereliable final prediction. The model is capable of predicting
zero prices, indicating the validity of the ”policy-versuspolicy” zero prices forecasting strategy. The ensemble learning
metrics for the day-by-day prediction approach using all the
proposed features is illustrated in Table III. By using the dayby-day prediction approach, the prediction accuracy is much
improved compared with the one-time prediction approach.

RMSE
5.66

Test data

Training data

R2
0.729

Test data

Training data

Test data
・・・・・・・・・・・・

Day i

Day i+1

Day i+2

Day i+3

B. Implications for EPF’s confidence interval

・・・・・・(a)Data set
Training data

Test data

・・・・・・

Day i

Day i+1

・・・・・・

Day i+365(b)

Figure 7. Day-by-day prediction schematic (a) and one-time
prediction schematic (b).
IV. R ESULTS AND DISCUSSION
A. Ensemble learning results
An ensemble learning technique was applied to each feature
set to generate ensemble predictions. Table II presents theprediction accuracy and the computation time of the proposed
LSTM and CNN-LSTM models. According to Table II, the
CNN-LSTM only has half the computation time of the LSTM
for same training times of ensemble learning, while stillboasts higher prediction accuracy over the LSTM using lesscomputation time.
TABLE II. Prediction accuracy and computation time comparison.
Model

Computation time

Ensemble times

R2

RMSE

LSTM
LSTM
CNN-LSTM
CNN-LSTM

59 min
120 min60 min40 min

15
30
30
20

0.57820.5821
0.58250.5824

7.062
7.028
7.027
7.031

The prediction results utilizing all the features using theday-by-day prediction approach is depicted in Fig. 8. It isnoteworthy from Fig. 8 that each training iteration results in aunique individual prediction on the test set, underscoring theinherent uncertainty in the neural network training process.

In Fig. 8, it is noteworthy that when the ensemble learning
prediction accuracy is high, the 30 individual predictions
tend to highly overlap each other and also demonstrate highprediction accuracy, as shown for the period 2022.05.092022.05.14 in Fig. 8(b). Conversely, when the ensemble learning fails to predict the actual value accurately, the 30 individual
predictions tend to diverge from each other without overlaping.
The standard deviation of the 30 individual predictions at thesame time point is calculated and used to compute the Person
correlation efficient with the ensemble learning prediction
error at the same time point. The correlation efficient value
is 0.462, indicating a moderate correlation. By utilizing theabove-mentioned phenomenon, a novel confidence interval
estimation method for EPF can be developed in our future
work.
V. C ONCLUSIONS AND FUTURE WORKThis study has proposed an innovative EPF framework thatutilizes multimodal data augmented by an ensemble learningtechnique. The CNN-LSTM model is superior over the LSTM
model in both prediction accuracy and computation time. The
logarithm transform of the pre-processing for the electricity
spot price data has been shown crucial for EPF in the Japan
electricity spot market. A ”policy-versus-policy” strategy has
been proposed to solve the zero-inflated regression problem,
which halfs the computation time compared with traditional
two-stage method. By using the day-by-day prediction approach, the ensemble learning method achieved achieved a
RMSE and (R2 ) of 5.66 JPY/kWh and 0.728 over the test
range of the full year 2022. An implication for a novel
estimation method of the EPF’s confidence interval using theensemble learning approach has been proposed.

Spot price (yen/kWh)Spot price (yen/kWh)Individual Predictions

Ensemble Prediction

Actual

40
20
0

75
50
25
0

2022-05-01

2022-05-05

2022-05-09

2022-05-13

2022-05-17
Time [Daily]D

2022-09-01

2022-09-05

2022-09-09

2022-09-13

2022-09-17
Time [Daily]E

2022-05-21

2022-09-21

2022-05-25

2022-09-25

2022-05-29

2022-06-01

2022-09-292022-10-01

Figure 8. Individual and ensemble predictions for the test range of May and September in 2022, using all the proposedfeatures and day-by-day prediction approach.

R EFERENCES[1] J. Seel, D. Millstein, A. Mills, M. Bolinger, and R. Wiser, “Plentiful
electricity turns wholesale prices negative,” Advances in Applied Energy,
vol. 4, p. 100073, Nov. 2021.
[2] A. Shiri, M. Afshar, A. Rahimi-Kian, and B. Maham, “Electricity priceforecasting using support vector machines by considering oil and natural
gas price impacts,” in 2015 IEEE International Conference on Smart
Energy Grid Engineering (SEGE). ieeexplore.ieee.org, Aug. 2015, pp.
1–5.
[3] K. Özen and D. Yıldırım, “Application of bagging in day-ahead electricity price forecasting and factor augmentation,” Energy Econ., vol. 103,
p. 105573, Nov. 2021.
[4] K. He, Y. Xu, Y. Zou, and L. Tang, “Electricity price forecasts using
a curvelet denoising based approach,” Physica A: Statistical Mechanics
and its Applications, vol. 425, pp. 1–9, May 2015.
[5] Z. Yang, L. Ce, and L. Lian, “Electricity price forecasting by a hybrid
model, combining wavelet transform, ARMA and kernel-based extreme
learning machine methods,” Appl. Energy, vol. 190, pp. 291–305, Mar.
2017.
[6] N. Chaâbane, “A hybrid ARFIMA and neural network model forelectricity price prediction,” Int. J. Electr. Power Energy Syst., vol. 55,
pp. 187–194, Feb. 2014.
[7] A. J. Conejo, M. A. Plazas, R. Espinola, and A. B. Molina, “Day-ahead
electricity price forecasting using the wavelet transform and ARIMA
models,” IEEE Trans. Power Syst., vol. 20, no. 2, pp. 1035–1042, May
2005.
[8] G. P. Girish, “Spot electricity price forecasting in indian electricity
market using autoregressive-GARCH models,” Energy Strategy Reviews,
vol. 11-12, pp. 52–57, Jun. 2016.
[9] I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning. MIT Press,
Nov. 2016.
[10] J. Chen, G.-Q. Zeng, W. Zhou, W. Du, and K.-D. Lu, “Wind speed
forecasting using nonlinear-learning ensemble of deep learning timeseries prediction and extremal optimization,” Energy Convers. Manage.,
vol. 165, pp. 681–695, Jun. 2018.
[11] J. E. P. Exchange, “Day ahead market,” https://www.jepx.jp/en/
electricpower/market-data/spot/, 2023, accessed: 2023-8-19.
[12] J. Organization for Cross-regional Coordination of Transmission Operators, “Menu,” https://occtonet3.occto.or.jp/public/dfw/RP11/OCCTO/
SD/LOGIN login, 2023, accessed: 2023-8-1.
[13] J. M. B. S. Center, “Numerical weather prediction model gpvmsm,” http://www.jmbsc.or.jp/jp/online/file/f-online10200.html, 2023,
accessed: 2023-7-15.
[14] “Google maps,” https://www.google.com/maps/@36.2932467,137.
3408308,6z?entry=ttu, accessed: 2024-5-3.
[15] O. team, “Opencv library,” https://opencv.org/, 2023, accessed: 2023-819.

[16] G. E. P. Box, G. M. Jenkins, G. C. Reinsel, and G. M. Ljung, Time
Series Analysis: Forecasting and Control. John Wiley & Sons, May
2015.
[17] S. J. Taylor and B. Letham, “Forecasting at scale,” Am. Stat., vol. 72,
no. 1, pp. 37–45, Jan. 2018.
[18] R. J. Hyndman and G. Athanasopoulos, Forecasting: principles andpractice. OTexts, May 2018.