Hidden Markov Mixture of Gaussian Process Functional Regression: Utilizing Multi-Scale Structure for Time-Series Forecasting
Authors: Tao Li
Jinwen Ma
Venue: Submitted to ICLR 2023 (2023)

Abstract: The mixture of Gaussian process functional regressions (GPFRs) assumes that there are a batch of time-series or sample curves which are  generated by independent random processes with different temporal structures. However, in the real situations, these structures are actually transferred in a random manner from a long time scale. Therefore, the assumption of independent curves is not true in practice. In order to get rid of this limitation, we propose the hidden Markov based GPFR mixture model (HM-GPFR) by describing these curves with both fine and coarse level temporal structures. Specifically, the temporal structure is described by the Gaussian process model at the fine level and hidden Markov process at the coarse level. The whole model can be regarded as a random process with state switching dynamics. To further enhance the robustness of the model, we also give a priori to the model parameters and develop Bayesian hidden Markov based GPFR mixture model (BHM-GPFR). Experimental results demonstrate that the proposed methods have both high prediction accuracy and good interpretability.


Under review as a conference paper at ICLR 2023
HIDDEN MARKOV MIXTURE OF GAUSSIAN PROCESS FUNCTIONAL REGRESSION: UTILIZING MULTI-SCALE STRUCTURE FOR TIME-SERIES FORECASTING
Anonymous authors Paper under double-blind review
ABSTRACT
The mixture of Gaussian process functional regressions (GPFRs) assumes that there are a batch of time-seriss or sample curves which are generated by independent random processes with different temporal structures. However, in the real situations, these structures are actually transferred a random manner from a long time scale. Therefore, the assumption of independent curves is not true practice. In order get rid of this limitation, we propose the hidden Markov based GPFR mixture model (HM-GPFR) by describing these curves with both fine and coarse level temporal structures. Specifically, the temporal structure is described by the Gaussian process model at the fine level and hidden Markov process at the coarse level. The whole model can be regarded further enhance the robustness of the model, we also give a priori the model parameters and develop Bayesian hidden Markov based GPFR mixture model (BHM-GPFR) Experimental results demonstrate that the proposed methods have both high prediction accuracy and good interpretability.
random process with state switching dynamics. To
INTRODUCTION
The time-series considered this paper has the multi-scale structure: the coarse level and the fine level. We have observations (y)...., . Vr) where each y: (.. -Ut, itself is a time-series of length L. The whole time-series is arranged as
21, 11.2............................ LI 22 12.2
22 ................................ 17,1 21. MT.L.
(I)
The subscripts of {3r )T=1 called coarse level indices, while the subscripts of (ye.i are called fine level indices. Throughout this paper, we take the electricity load dataset as a concrete example. The electricity load dataset consists of T 365 consecutive daily records, and in each day there L = 96 samples recorded every quarter-hour In this example, the coarse level indices denote "day", while the fine level indices correspond to the time resolution of 15 minutes. The aim is to forecast both short-term and long-term electricity loads based on historical records. There may be partial observations yr 11,1 1T+1.M with M L, so the entire observed time-series has the form
11,1: V1,2. ... V1.L 12,1 V2.2. V2.L : FT,1, VT. VT 2T.1.1.................1 UT-1,M
(2)
The task is predict future response ye.,,, where t, 2 +,11 < < are positive integers.
The coarse level and fine level provide different structural information about the data generation process. In the coarse level, each 3/. can be regarded
time-series. and there is certain cluster structure (Shi & Wang 2008 Wu & Ma 2018) underlying these time-series {y) ) we can divide {y.)²= into groups such that time-series within each group share a similar evolving trend. Back to the electricity load dataset, such groups correspond to different electricity consumption patterns. We use 2 to denote the cluster label of y:. the fine level, observations { can be regarded as realization of stochastic process, and the properties of the stochastic process determined by the cluster label žz-
The mixture of Gaussian processes functional regression
model Shi & Wang 2008 Shi & Choil 2011) is powerful for analyzing functional data or batch data, and is applicable to the multi-scale time-series forecasting task. Mix-GPER assumes there are K Gaussian processes

---

Under review as a conference paper at ICLR 2023
Flectricity load (million kW
10
8
Up
200
300
400
500
600
1000
Time
Figure 1: An illustration of multi-scale time-series.
functional regression (GPFR) (hi et al. 2007 components, and associated with each Yr there is latent variable 2, indicating y: generated by which GPFR component. Since GPFR is good at capturing temporal dependency, this model successfully utilizes the structure information in the fine level. However, the temporal information in the coarse level is totally ignored since mix-GPFR assumes {z+]== are i..d.
In this work, we propose model the temporal dependency in the coarse level by the hidden Markov model, which characterizes the switching dynamics of 2i...., by the transition probability matrix. We refer the proposed model HM-GPFR. Mix-GPFR is able to effectively predict Br+1.M+1...... UT+LL when M > O. To predict the responses yr++,, we must determine the cluster label 2T+1 based on observations 1T+ 1.+ otherwise we do not know 4/++1 is governed by which evolving pattern. If there is no observation at day T + 1 (i.e., M 0), then mix-GPFR fails to identify the stochastic process that generates 1/+1- For the same reason, mix-GPFR is not suitable for long-term forecasting (t, > T 1). On the other hand, HM-GPFR is able to infer 2, for any t, based on the transition probabilities of the hidden Markov model even M = 0. Therefore, HM-GPFR makes of coarse level temporal information and solves the cold start problem in mix-GPFR, Besides, when a new day's records yr+11 have been fully observed, one needs to re- trairan mix- model to utilize yT+ while HM-GPFF can adjust the parameters incrementally without retraining the model.
RELATED WORKS
Gaussian process Rasmussen Williams 2006) is a powerful non non-parametric Bayesian model. In (Girard et al. 2002 Brahin-Belhouari & Bermak 2004 Girard & Murray-Smite 2005), GP has been applied for time-series forecasting. Shi et al. proposed the GPFR model to process batch data (Shi et al.) 2007). To effectively model multi-modal data, the mixture structure is further introduced
to GPFR and the mix-GPFR model was proposed (Shi & Wang 2008 Shi & Choi 2011). In (Wu & Ma 2018 Li al. 2019 Cao al., 202T), GP related methods for electricity load prediction have been evaluated thoroughly. However, in these works daily records are treated as i.i.d. samples, and the temporal information the coarse level ignored.
Multi-scale time-series was proposed in (Ferreira et al. 2006 Ferreira & Lee 2007b a), and further developments in this direction have been achieved in recent years. The time-series considered in this work is different from the multi-scale time-series since at the coarse level there is no aggregated observation from the samples at the fine level. In this paper, we mainly emphasize the multi-scale structure of the time-series.
PRELIMINARIES
3.1 HIDDEN MARKOV MODEL
For a sequence of observations {ye {y the hidden Markov model (HMM) RRaiiner Juang 1986 Elliott et al., 2008) assumes there is a hidden state variable 2; associated with yt- The sequence of hidden states {x,]T_ forms a homogeneous Markov process. Usually, {z, (=)) are categorical variables taking values in (1,..., K.. and the transition dynamics is governed by P(z,, 2|2-11
2

---

Under review as a conference paper at ICLR 2023
= Pka. There are K groups of parameters k and 2, k indicates that the observation yo is generated by P(y: 0z). The goal of learning is to identify the parameters and infer the posterior distribution of hidden states {2, ) Usually. the Baum-Welch algorithm (Baum & Petrie 1966 Baum et al. 1970) is utilized to learn the HMM, which can be regarded as specifically designed EM algorithm based the wwwwdddbckwward algorithm. Once the model has been trained, we able to simulate future behavior of the system.
3.2 GAUSSIAN PROCESS FUNCTIONAL REGRESSIONS
Gaussian process stochastic process that any finiteddimensional distribution of samples is a multivariate Gaussian distribution. The property of a Gaussian process is determined by the mean function and the covariance function. We write the mean function as y((-) and the covariance function as c(.,)). Suppose that we have dataset D = { y:) The relationship between input and output is connected by a function V, i.e., (x,) = Di. Let x [211 22,,22...... ................................
,r y
(yn 12, .......... yrl", then we assume y|æ N(u, C) where f = [u(ær), p((r)........ (r ) and C2 = c(2,, 2,). In machine learning, the mean function and the covariance function are usually parameterized. Here. we use the squared exponential covariance function (Rasmussen & Williams
2006; Shi Choi 2011 Wu & Ma 2018 c(z,, Ty: 0) 0f exp
(-02 (-)--))) + fzgby. where 2
dy is the Kronecker delta function and 0 00 B2. 0z].
linear The mean function is modeled
combination of B-spline basis functions Shi et al., 2007 Shi & Choi 2011). Suppose that we have
D B B-spline basis functions {x) Let (xx) = d= badd(2)) and & be an L. x D matrix with Lid = (z(z) b = [bi. by bD] then y x ~ N($b, C). From the function perspective, this model can be denoted as 3/ (x) GPFR( O).
We can use the Gaussian process to model the multi-scale time-series considered in this paper, and the key-point is transform the multi-scale time-series to a batch dataset. For each coarse level index t, we can construct a dataset Dj = {( ) where 2t,; is the sampling time of i-th sample in time-series yt- Let 3 be the function underlying dataset D, i.c., 3 (Iz,z) = yt.i, then these (Di)T=
GPFR(z; b, 0). Without loss of generality, we may assume Ft,5 i, and thus a. = da(i), Ct2 = j(i.j: do not depend on the coarse level index f. Therefore, is equivalent to assume y 17" are independently and identically distributed as ((@b, C). To learn the parameters b and 0 we apply the Type-Il maximum likelihood estimation technique (Rasmussen & Williams| 2006 Shi & Choi 2011).
can be regarded as independent realizations of a GPFR, which assumes 3/ (x=) i.d.
As for prediction, given a new record {zz,, H and we want to predict the corresponding output y/t,, .. at .t. where M < i. S L., from the definition of Gaussian process we immediately know that Ut... also obeys a Gaussian distribution (Rasmussen & Williams 2006). Let
x. -lsta,,,............. st.,M] y. [y...,,...
Mt..M]"
μ. = [u((t..1) ..., u(xt M)]" IT, [C.lis = (.
then the mean of yt..i, is p(Ie,..)) + c(F....,, C, (y. - -..). and the variance of yr..s. is c(zi, d.: (.: c( xx) C. c(E,,Tt,,,, Note that if M = O, the prediction is simply given by (H(2z,, c(Te, ,t, )), which equals to the prior distribution of ye.,,, and fails to utilize the temporal dependency with recent observations. In the electricity load prediction example, this means we can only effectively predict new day's 's electricity loads when already have the first few observations of this day. In practice, however, it is very common to predict a new day's 's electricity loads from scratch.
3.3 THE MIXTURE OF GAUSSIAN PROCESS FUNCTIONAL REGRESSIONS
GPFR implicitly assumes that all {3r)7 are generated by the same stochastic process, which is not the case in practice. In real applications, they may be generated from different signal sources, thus single GPFR is not flexible enough to model all the time series, especially when there are variety of evolving trends. Take the electricity load dataset for example, the records corresponding to winter and summer are very likely to have significantly different trends and shapes. To solve this problem. Shi et al.() suggested to introduce the mixture structure to GPFR, and proposed the mixture of Gaussian process functional regressions (mix-GPFR). In mix-GPFR, there are K GPFR
3

---

Under review as a conference paper at ICLR 2023
components with different parameters {bk. 0x] and the mixing proportion of the k-th GPFR component is Ty- Intuitively, there are K different signal sources or evolving patterns in mix-GPFR to describe temporal data with different temporal properties. For each yt, there is associated latent indicator variable 2, E (1,2..... K), and 21 = k indicates ye is generated by the k-th GPFR component. The generation process of mix-GPFR is
i.id Categorical(Ti, . TK),
(S)
(x) = ~ GPFR(z; by. 0x)
Let C₁ L L be the covariance matrix calculated by 0g. i.e., [Cals j(i,j;0x) then the above equation is equivalent to y/ ~ (bb,, C₁).
Due the existence of latent variables, the parameter learning of mix-GPFR involves the EM algorithm (Dempster et al. 1977, Shi & Wang/ 2008). As for prediction, K GPFR components of miX-GPFR first make predictions individually, then we weight these predictions based on the posterior probability P(=,, klyr. b Oz). Note that if M = O then P(z,, klys._ bb, 0₂. wa. which equals to the mixing proportions and also fails to utilize recent observations. Therefore, mix-GPFR also suffers from the cold start problem.
4 PROPOSED METHODS
4.1 HIDDEN MARKOV BASED GAUSSIAN PROCESS FUNCTIONAL REGRESSION MIXTURE MODEL
Similar to mix-GPFR, the hidden Markov based Gaussian process functional regression mixture model also assumes the time-series is generated by K signal sources. The key difference is that the signal source may switch between consecutive observations in the time resolution of the coarse level. The temporal structure in the coarse level is characterized by the transition dynamics of {2+] and the temporal dependency in the fine level is captured by Gaussian processes. Precisely,
21 Categorical 11; n2, ........... TK), P( = 1 - Pai 1 2.3....,
(6)
(xx))2 ~ b, 0x) 1.2........
Here, π [TT, ....... TK] is the initial state distribution, and P - [Pu]]x is the transition probability matrix. We refer to this model as HM-GPFR. In GPFR and miX-GPFR, the observations {yr)== are modeled as independent and exchangeable realizations of stochastic processes, thus the temporal structure in the coarse level is destroyed. However, in HM-GPFR, consecutive yt-1, Be are connected by the transition dynamics their corresponding latent variables 2,-1, Zt, which is more suitable for time-series data. For example, if today's 's electricity loads are very high, then it is unlikely that tomorrow's 's electricity loads are extremely low.
The learning algorithm for HM-GPFR is based on the EM algorithm, and derive the algo- rithm in the appendix. After the parameters have been learned, we assign the latent variable & arg max N' (k) and regard { as deterministic. For prediction, we consider two cases: t. = T and f, > + 1. When to = T + 1, the latent variable 27+1 is determined by both the conditional transition probability 2r+1 12r and partial observations yr- 1. More precisely, suppose द = 1, then
(7) where the square brackets denote slicing operation. If M O. then Wi- Pte. is determined by the last hidden state and transition dynamics, which more accurate than mix-GPFR. Suppose the prediction of the k-th component is y² then the final prediction is given by EKnwey
W/ = P = k kTT, = θ) x PuNN 2[1 M.:bb. C[1 M, M)), :M).
We next consider the case f, > T the main difference is the posterior distribution of 21, this case, we need to use the transition probability matrix recursively. First, we calculate the distribution of 2 according to Equation 77. Then by the Markov property, we know
P(z.. = K|T, yT 2т 4:
P(Zr+1 =
= O) T- -T
(8)
The final prediction is also given by ΣK Whye (k) =Ek,
wa wi..
4

---

Under review as a conference paper at ICLR 2023
Pa
00
A 1,2
1.... K
21
P2
NN
f1
1/2
UN
0z: 7k, by
0c. 7kk
bi
1726. Eny
A 1,2 ....K K
k= 1,2...
(a) HM-GPFR
(b) BHM-GPFR
Figure 2: Probabilistic graphical models of HM-GPFR and BHM-GPFR.
4.2 BAYESIAN HIDDEN MARKOV BASED GAUSSIAN PROCESS FUNCTIONAL REGRESSION MIXTURE MODEI.
One drawback of HM-GPFR is that there are too many parameters and thus has the risk of overfitting. In this section, we further develop a fully Bayesian treatment of HM-GPFR. We place a Gaussian prior N ((m,, Σx) the coefficients of B-spline functions (b&) k= For the transition probabilities, let P [Ps1 Pk...... P |T be the probabilities from state ic to other states, then assume ри: obeys Dirichlet prior Dir(au..... Ou). The generation process of Bayesian HM-GPFR
by ~ N((m, Σx), k ==
K
Pk Dir(ag: ................................ a0), 1 1,2 K
21 ~ Categorical (arriical((11111)) (2, TK),
(9)
P( 1|2:-1111 k) PaI, 1 2,3,...,
{[()|< 2 = k ~~ GPFR(( by 02), t 1,2....,
The detailed learning algorithm is presented in the appendix. After learning, we set the latent variables to their maximum posteriori (MAP) estimates 52. Specifically, by m Pst Σm= 3km
&g.
arg max =1,..1 1, Ye (k:). The rest of prediction the same as HM-GPFR.
EXPERIMENTAL RESULTS
5.1 EXPERIMENT SETTINGS
In this section, we use the electricity load dataset issued by the State Grid of China for city in northwest China. The dataset records electricity loads every 15 minutes, thus there are 96 records per day. Using the electricity load records of 2010 for training, we pre- dict the subsequent S-step electricity loads in a time-series prediction fashion, where 1, 2. 3, 4, 10. 20, 30, 50, 80, 100, 200, 500, 1000. This setting allows both short-term and long-term predictions to be evaluated. For more comprehensive and accurate assessment of the performance, we roll the time series by 100 rounds. Based on the electricity loads of 2010, the r-th round also puts the first (r 1) records of 2011 I into the training set. In each round, we predict the subsequent S-step electricity loads. In r-th round, suppose the ground-truths are y1, ....... ..... ys and the predictions are 91. V2..... Us. we use the Mean Absolute Percentage Errors (MAPEs) to evaluate the prediction results. Specifically, MAPE,
DS 1v-- For overall evaluation, report the average [8ol
of 100 MAPES to obtain MAPE TDU randomness, we repeat the algorithms for 10 runs and report the average results.
2r MAPE,. Since the algorithms influenced by
100
We compare HM-GPER and BHM-GPFR with other time-series forecasting methods. Specifically,
5

---

Under review as a conference paper at ICLR 2023
Table 1: MAPE of various methods on the electricity loads dataset under different step lengths and parameter settings.
Mehsod
Iregh
10
20
30
50
90
100
200
3043
500
10007
-3055
5.135
2.55%
1.31%
6.9%%
8.85%
13-7555
13:955
15.05
15 13AMKX
S7
13:54%
10: 4.-85%
1.01%
1355
1.21%
2.13%
2.57%
4.00%
6.90%
8.63%
13.15%
13.23%
10.3%%
17 17.56
17.6%
15.91%
16-39%
5.3055
5.74%
641355
9.3335
11.85%
14.85%
15.2%
15.5855
16.35%
15. 4848
17.35%
16.05%
MA
2.585%
3
4.44%
0.813
6.838
11.1.3
14.2156
15.27%
18.05%
18
13-27%
18.02%
1,01%
2.345
L,355
3.125
2.02%
1.45%
6.85%
5.0%
13.25%
13.55%
14.01%
17.5%%
16.77%
15.35%
ARMA
3.01%
2.445
17255
2003
4.3%%
6.57%
8.5255
13.15%
13.05%
141.%%
17.2%
37 7744%
15-71%
16.34%
12
3,00%
.305
1.705
21475
3.07%
4.30%
6.5%
5.4%
19.99%
13.85%
13.975
17 17
17.45
15-455
0.985
311
1.21%
2.188
2574
1.58%
7.3%%
9085
11,99%
16.45%
15339
16
16675
17:055
17:45%
ARIMA
1.01%5
-355
1.7555
2.1
25553
1.585
1.355
1:25
13:155
14335%
15
58 385033
56
15.5%
16.37%
1.01%
1-K5
1,335
2.55%
2.85
4.88%
3.61%
11
15-11%
23.05%
21.335
28 100
24 383
21.7%
25.33%
1.66%
2.085
4.06%
6.119
7.67%
SARMA
0.555
..88%
L.125
1.005
1,6744
4,945
9,1255
0.305
1.545
5.54%
0315
7.5854
1.12 12
0.525
-477%
1.358
1.823
1.725
2.625
4.065
6.165
0.355
7.1754
5.115
10 103
10005
1.905
7.505
19 21018
13.44%
132275
15
16
17
10
10
17
17.555
LSTM
12
234
12.3975
11,455
11.49 13.335
31.45%
31.40 1.40%
12.00 8508
11.5%
11.8%
12.055 0.055
11.0%
13 1973
1.533
11 1355
35%
15.2555
17.30%
16 6595%
77.2%
18 00033
15.8%
17 53.5%%
11.415
1. 08
10.1%
10.1133
10.117
30.10%
10.2%
10.49%
11 17275
11 11467
13945
13.00%
13.00%
16 1054
17 17.02%
17.98%
18.37%
6.31%
11.55
FNN
1- 08
1. 24
12
0.557
015855
..35%
1.15
1963
12295
1.37%
11155
1402%
1-333
1.5835
1185
1.2075
1.852
5.366
5.155
0.2355
5.345
$:500
411355
5
13368
T.1255
1613
11:455
11.55%
59954
10 30
30 55
18 NE3
36
3:20
155
15:055
19:015
1722
17.87%
19.8%%
5.54%
1.71%
2.035
5.55%
10.254
14.55
15.77
18 18047
10
SVR
4-11 12
1000
1.005
1355%
755
1.555
1905
1.0036
2.1755
1.8555
3.725%
1.03
0.55%%
4.835%
0.0155
6385%
5.50%
13.05
17.85%
13-435
10.24
24 05.3
25.359
25.355
15-3655
23.65%
16-13%
1.25%
485
1.415
1.48%
1.96%
7.88%
0.1%
5.31%
7.9%
11-13%
15.05%
12.39%
9.87%
6 6.5%
5.0%%
4.K
0.97%
205
8.387%
4.00%
6.4255
7.5475
10.25%
13.45%
16.005
1617%%
17
17:10%
15.81%
mm 4.K
0.50%
.25%
L845
190
2.5225
4,565
0.35%
1.6155
19.185
13.45%
10
18.00 13.885
17.50%
17.16%
15.83%
0.57%
-20%%
LA5S
1.003
2.855
4.005
0.45%
7.57%
10:255
15.054
15.077
18 18-105
17.555
17.66%
10.42%
0.58%
3,60%
5.82%
13.4455
14.12%
15.35%
15-78%
18.41%
10
EGPM
K
0.90%
19%
1.51%
1.11%
1,8%
1.85%
5.5%
1.85%
3.67%
5.00%
7.95%
8.700
16.15%
13.31%
15 15.0%%
15 15.35%
19 1000
19
19 16.
19
33: 13-1755
15. 15.155
15.13%
21.515
0.97%
1422
185%
7.97%
5.81%
0.51%
12/02%
11.35%
17 115
15:455
21.45
1.025
5.42%
7.0355
9.055
11.45%
12.0%%
53
10 12.00%
48 K
1.02%
1.25%
1,225
1.425
1.555
3.85%
5.63%
7.4455
0.605
12.35%
14.42%
54 12%
12 05.85%
35.55%
45.86%
L- 65.8
10
3.025%
1.535
1.703
1.55%
3.945
0.7%
7.435
0.38%
12.20%
14.35%
53.01
32.75
35.96%
43.09%
MIEGIM
K-33
0.785
0.91%
1.11%
1.962
1.4%
7.10%
3.17%
134%
5.85%
6.0%%
9.27%
10
7.3%
1.55%
DIM-GIEE
0.606
0491%
0377
10336
1.00%
1.436
3.89%
9.6
330%
1.184
1.8%
10 15%
19.30%
12.36%
13.0%
RM GIFR
3.15
1.078
8.51%
4.07%
10
4.805
Classical times-series forecasting methods: auto-regressive (AR), moving average (MA), auto-regressive moving average (ARMA), auto-regressive integrated moving average (ARIMA), seasonal auto-regressive moving average (SARMA).
Machine learning methods: long short-term memory (LSTM), feedforward neural net- work (FNN), support vector regression (SVR), enhanced Gaussian process mixture model (EGPM).
GPFR related methods: the mixture of Gaussian process functional regressions (mix- GPFR), the mixture of Gaussian processes with nonparametric mean functions (miX-GPNM), Dirichlet process based mixture of Gaussian process functional regressions (DPM-GPFR).
Detailed parameter settings of comparison methods are shown in the appendix. The main parameters in HM-GPFR and BHM-GPFR are the number of components K and the number of B-spline basis functions D, and we set = 5, D = 30.
5.2 PERFORMANCE EVALUATION AND MODEL EXPLANATION
The prediction results of various methods on the electricity load dataset are shown in table the table, we can see that the prediction accuracy of classical time-series forecast methods decreases significantly as we increase the prediction step. Among them, SARMA outperforms AR, MA, ARMA, and ARIMA, because SARMA takes the periodicity of data into consideration and can fit data more effectively. The results of machine learning methods LSTM, NN, SVR, and EGPM also have similar phenomena, that is, when S is small, the prediction accuracy is high, and when S is large, the prediction accuracy is low. This observation indicates that these methods are not suitable for long-term prediction. In addition, machine learning methods are also sensitive to the settings of parameters. For example, the results of FNN and SVR are better when L = 4, which is close to SARMA, while the long-term prediction accuracy of EGPM decreases significantly when L is relatively large. It is challenging appropriately set hyper-parameetrss in practice. When making
T From
6

---

Under review as a conference paper at ICLR 2023
03417
08337
0.333
5000C
DYTTSS
0.390
COON 7
02309
0LE000
D0833
2000
(a) HM-GPFR
(b) BHM-GPFR
Figure 3: Evolving law of electricity loads and transition dynamies learned by HM-GPFR and BHM-GPFR.
long-term prediction, classical time-series prediction methods and machine learning methods need to recursively predict the subsequent values based on estimated values, which will cause the accumulation and amplification of errors. On the other hand, Gelaerelated methods first make predictions according to the mean function, then finely correct these predictions based on observed data. The mean function part can better describe the evolution law of data, which enables us to historical information and structural information in data more effectively. MiX-GPER, mix-GPNM, and 3PM-GPFR obtain similar results in long-terg prediction compared with SARMA, and can even achieve the best results in short-term prediction. This observation demonstrates the effectiveness of GPFR-related methods. However, these methods cannot deal with long-term prediction tasks well due to the "cold start" problem. Overall, the performances of the proposed HM-GPER and BHM-GPFR are more comprehensive. For medium-tern and short-term prediction, the results of HM-GPFR and BHM-GPFR are slightly worse than those of mix-GPFR. mix-GPNM, and DPM-GPFR, but they still enjoy significant advantages compared with other comparison methods. In terms of long-term forecasting, HMM-PPRR and BHM-GPP outperform miX-GPFR, which shows that considering the multi-scale temporal structure between daily electricity load time- series can effectively improve the accuracy of long-term forecasting. In addition, BHM-GPFR is generally better than HM-GPFR, which shows that giving prior distributions to the parameters and learning in a fully Bayesian way can further increase the robustness of the model and improve the prediction accuracy.
and PPMR,
HM-GPFR and 9HM-GPFR have strong interpretability. Specifically, the estimated values of hidden variables obtained after training { ) divide the daily electricity load records into K categories according to the evolution law. Each evolution pattern can be represented by the mean function of GPFR component, and these evolution patterns transfer each other with certain probabilities. The transfer law is characterized by the transfer probability matrix in the model. In fig. 3 we visualize the evolution patterns and transfer laws learned by HM-GPFR and BHM-GPFR. We call the evolution law corresponding to the mean function represented by the orange curve (at the top of the figure) mode 1, and call the five evolution modes as mode mode 5 respectively in clockwise order. Combined with the practical application background, some meaningful laws can be found according to the results of learned models. Examples are as follows:
The electricity load of mode is the lowest. Besides, mode is relatively stable: when the system is in this evolution pattern, then it will stay in this state in the next step with a probability of about 0.5. In case of state transition, the probability of transferring to the mode with second lowest load (mode 2 in Figure 3al and mode 3 in Figure (b) is high, while the probability of transferring the mode with highest load (mode 5 in Figure 3a and mode
2, mode 5 in Figure 3b) is relatively low:
The evolution laws of mode 2 and mode 5 in fig. 3b are very similar, but the probabilities of transferring to other modes both of them can be regarded as the mode with the highest load. When the system is in the mode with the highest load (mode 5 in Figure 3a and mode 2, mode 5 in Figure 35), the
different. From the perspective of electricity load alone,
7

---

Under review as a conference paper at ICLR 2023
10
Electrieity (millinn kW
hon
500
1000
1500
2000
2500
Time
(a)
Electricity load (million kW-b)
M m / (1)
mm
2000
Time
(b) GEE-
Figure 4: One-step-ahead rolling prediction results of HM-GPFR and BHM-GPFR.
probability of remaining in this state in the next step is the same as that of transferring the mode with the lowest (mode 1);
When the system is in the mode with the second- -iihhstt load (mode 3 in fig. 3a and mode 4 in fig. 3b), the probability of remaining in this state in the next step is low, while the probabilities of transferring to the modes with the lowest load and the highest load are high.
These laws are helpful for to understand the algorithm, have certain guiding significance for production practice, and can also be further analyzed in combination with expert knowledge.
The case S lin table I is the most common in practical application, that is, one-step-ahest ahead rolling forecast. As discussed in section[4.] when making rolling prediction, HM-GPFR and BHM- GPFR can dynamically adjust the model incrementally after collecting new data without retraining the model. The results of the one-step-ahead rolling prediction of HM-GPFR and BHM-GPFR on the electricity load dataset are shown in fig. can be seen that the predicted values of HM-GPFF and BHM-GPFR are very close to the ground-truths, indicating that they are effective for rolling prediction. the figure, the color of each point
weighted average of the colors corresponding to each mode in fig. 3 according to the weight wk. Note that there color changes in some electricity load curves in fig. 4a and fig. 4b) Taking the time-series in fig. 4al in the range of about I 1100- 1200 as example, when there
few observation data that day, HM-GPFR believes that the electricity load evolution pattern of that day is more likely to belong to mode 3. With the gradual increase of observation data, the model tends think that the electricity load evolution pattern of that day belongs to mode 5, and then tends to mode 3 again. This shows that HM-GPFR and BHM-GPER can timely adjust the value of 2,, according the latest information during rolling prediction.
5.3 ABLATION STUDY
In this section, we mainly compare HM-GPFR, BHM-GPFR with mix-GPFR, mix-GPNM, and DPM-GPFR to explore the impact of introducing coarse-grained temporal structure on the prediction performance. The MAPEs reported in Table 2 are averaged with respect to 1..... ... 100, while in this section we pay special attention to the case of r 1. In this case, the observed data is the electricity load records in 2010, and there is partial observations M = 0 in eq. (2). Therefore, mix-G GiXR mix miX-G and DPM-GPFR will encounter the cold-start problem. Table reports the MAPE of these methods at different prediction steps when 1. It can
January 1, 2011 (i.c.,
8

---

Under review as a conference paper at ICLR 2023
Table 2: MAPE of GP related methods under the cold-start setting 1).
Methon
Step longts
2
19
20
30
50
50
100
200
300
500
1000
mit-GRI GHT
6.725
6.75%
6.98%
T.0%%
1.215
7.125
T.18%
1.00%
5.06%
5.00%
8.49%
11,005
10,85%
1.50%
7.4%
mix- BIX-
6.7256
6,785
6.96%
1.05%
1.29
7,1195
1.08%
1.00%
8.0655
5.00%%
8.4855
11.065
10.8596
7.00%
7.4%
DPM GEER
11.79%
11.81%
12,035
12.14%
12.30%
12,42%
12.00%
12.53%
10,7356
12.97%
12.88%
15.33%
11,955
8.4%
6.4175
HIM-OPER BHMM-CIPR
6.41%
6.635
6.61%
6.60%
6.78%
6.6%
6.21%
0.385
8.415
5.16%
1.5055
11,505
10.45%
7.44%
6.71%
0.585
4.6155
4.84%
5.8%%
5.0756
4.9895
5. 1.18
5.5758
7:33
704456
0.7096
10.67%
10.3596
7.58%
7.24%
Electricity load (million kW
10-
training Ground-trunh
8
Predictions of M'mix--PER
Predictions mix-GPNM
Predictions of DPM-
Predictions
Predictions of BHM-GPER
1000
Time
Figure 5: Multi-step prediction results of mix-GPPR, miX-GPNM, DPM-GPFR, HM-GPFR and BHM-GPFR.
be seen from the table that the prediction accuracy HM-GPFR and BHM-GPFR higher than that of mix-GPFR, mix-GPNM, and DPM-GPER at almost every step, which shows that coarse-grained temporal information is helpful improve the prediction performance, and the use of Markov chain to model the transfer law of electricity load evolution patterns make effective use of coarse-grained temporal information.
Figure further shows the results of multi-step prediction of these methods the electricity load dataset. Here also the case of cold start" (r = 1), and we predict the electricity loads in the next 10 days (960 time points in total). It can be seen from the figure that these methods effectively utilize the periodic structure in the time-series, and the prediction results show periodicity, but the prediction results of HM-GPFR and BHM-GPER are slightly different from other methods. Due the problem of "cold start", the predictions of miX-GPFR, miX-GPNM, and DPM-GPER for each +1. while HM-GPFR and BHM-GPFR will use
day are the same, i.e., NN++1 DN coarse-grained temporal information when making predictions, and then adjust the predicted values of each day. Based on the predicted values of other methods, it can be seen from the figure that the predicted values of HM-GPFR and BHM-GPP on the first day are higher, and with the increase in step size, the predicted values will tend to the weighted average value of the mean function of each GPFR component.
+10
CONCLUSION
In this paper, we have proposed the concept of multi-scale time series. Multi-scale time series have two granularity temporal structures. We established the HM-GPFR model for multi-scale time series forecasting and designed an effective learning algorithm. In addition, we also gave priori to the parameters in the model, and obtain more robust BHM-GPFR model. Compared with conventional GPFR-related related methods (mix-GPFR, miX-GPNM. DPM-GPFR), the proposed method can effectively use the temporal information of both fine level and coarse level, alleviate the "cold start" problem, and has good performance in short-tern prediction and long-tern prediction. HM-GPFR and CHMM-GPFR not only achieve high prediction accuracy but also have good interpretability. Combined with the actual problem background and domain knowledge, by the model.
explain the state transition law learned
REFERENCES
Leonard E Baum and Ted Petric. Statistical inference for probabilistic functions of finite state markov chains. The Annals of Mathematical Statistics, 1554-1554-1563 1966.

---

Under review as a conference paper at ICLR 2023
Leonard E Baum, Ted Petrie, George Soules, and Norman Weiss. A maximization technique occurring in the statistical analysis of probabilistic functions of markov chains. The Annals of Mathematical Statistics, 41(1):164-171, 1970.
Sofiane Brahim-Belhouari and Amine Bermak. Gaussian process for nonstationary time series prediction. Computational Statisties Data Analysis, 47(4):705-712 2004.
Di Cao. Junbo Zhao, Weihao Hu, Yingchen Zhang, Qishu Liao, Zhe Chen, and Frede Blaabjerg. Robust deep gaussian process-basee probabilistic electrical load forecasting against anomalous events. IEEE Transactions on Industrial Informatics, 18(2): (5(2):1142-1) 142-112211511 1251, 2021.
Arthur P Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39(1): 1-22, 1977.
Robert J Elliott, Lakhdar Aggoun, and John Moore. Hidden Markov models: estimation and control, volume 29. Springer Science & Business Media, 2008
Marco Antonio Rosa Ferreira and Herbert KH Lec. Multiscale modeling: a Bayesian perspective, volume 2. Springer, 2007a.
Marco AR Ferreira and Herbert KH Lee. Multiscale time series. Multiscale Modeling: Bayesian Perspective, pp. 113- 143, 2007b.
Marco AR Ferreira, David M Higdon, Herbert KH Lee, and Mike West. Multi-scale and hidden resolution time series models. Bayesian Analysis, 1(4):947-967, 2006.
Agathe Girard and Roderick Murray-Smith Gaussian processes: prediction at a noisy input and application to iterative multiple-step ahead forecasting of time-series. In Switching and Learning in Feedback Systems, pp. 158-184. Springer, 2005.
Agathe Girard, Carl Rasmussen, Joaquin Q Candela, and Roderick Murray-Smite Gaussian process priors with uncertain inputs: application to multiple-step ahead time series forecasting. In Advances in Neural Information Processing Systems (NeurlPs), pp. 552,, 2002.
Ling-Ling Ling Li, Jin Sun, Ching-Hsin Wang, Ya-Tong Zhou, and Kuo-Ping Lin. Enhanced gaussian process mixture model for short-term electric load forecasting. Information Sciences, 477:386-398, 2019.
Lawrence Rabiner and Biinghwang Juang. An introduction hidden markov models. IEEE ASSP Magazine, 3(1):4-16, 1986.
Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian processes for machine learning. The MIT Press, 2006.
Jian Qing Shi and Taeryon Choi. Gaussian process regression analysis for functional data. Chapman and Hall/CRC, 2011.
JQ Shi and B Wang. Curve prediction and clustering with mixtures of gaussian process functional regression models. Statistics and Computing, 18(3):267-283, 2008
JQ Shi, B Wang, Roderick Murray-Smith, and DM Titterington. Gaussian process functional regression modeling for batch data. Biometrics, 63(3):114-723, 2007.
Di Wu and Jinwen Ma. A two-layer mixture model of gaussian process functional regressions and its meme em algorithm. IEEE Transactions on Neural Networks and Learning Systems (T-NNLS), (99):1-11, 2018.
10

---

Under review as a conference paper at ICLR 2023
A LEARNING ALGORITHMS OF THE PROPOSED METHODS
A.1 HM-GPFR
Due to the existence of latent variables {z apply the EM algorithm to learn the HM-GPFR model. We write = {D:) to denote observations, Q =Ppx) ( by, 0s ) to denote
all parameters, and Ω = { f to denote all latent variables. First, the complete data log-likelihood
K
T-1 K
[(0,T,92)
Σ
[( =
log T
ΣΣΣ
[ = 1, = =)) log Pa
km1
(10)
I(a k)) log by, 0x) Σ Σ
(y:
km
In the E-step of the EM algorithm, we need calculate the expectation of Equation TTPp with respect to the posterior distribution of latent variables Ω to obtain the Q-function. However, it is not necessary explicitly calculate P θ), which categorical distribution with KWW possible values, and
suffices to obtain P( = 1, 2 T[:; () and (= :|T; We first introduce some notations as follows:
(k()) (31 1 12
yr 2t O)
(k) = P 11+2. 3/1 21 = k; O),
(11)
Ye = P(2 (), (k,,1) P(= = k =1+1 = IT: 9).
The key-point is calculate (*k)) and (k) Note that
Th (k) = P(2 = k|T )) x P( = k, T; = (k()))((k),
(12)
(kk () = a, (k) (y Tb,, (1)
Therefore, the problem boils down calculate or (k) and (k)). We can derive them recursively based on the forward-backward algorithm. According to the definition of oe we have
K
a1 (A) = (yy Dbg; C₁)
(k)) = Σ
(1))
N(y Dbx. Ct).
(13)
Similarly, according the definition of (k).,, we have
(r())
(k) Σ PRWNyeer Tb. C1)
(14)
To summary, in the E-step first Equations (3) and 14 to calculate (x(k), (kk) recursively
based on current parameters, then calculate "Ve E 1) according to Equation (12) The Q-function given by
K
2
Σ 21 (k) log T& k=1
ΣΣΣ (k, log Par+
ΣΣ
r (k) Nog,N(;; bbb. Ca)- (15)
In the M-step, we need to maximize Q with respect to parameters. The parameters { K and (P can be optimized in closed form,
(k) ΣK (1)
EE Ee(k,
Th
PR
(16)
Σ Σm= (((km
The parameters (ba. 0K) k= cannot be solved in closed form, and we apply the gradient ascent
algorithm to optimize Q ( with gradients
(k) ((C +bb)(yy
bbb) "C C=1)
8Ck
00, 2
00,
(17)
a
T mkkk T 4bb).
ab:
11

---

Under review as a conference paper at ICLR 2023
The complete algorithm is summarized in Algorithm When the partial observations yr+11...... 9T+1,M become complete as we collect more data, we can adjust the parameters mnrr- mentally without retraining the model. This is achieved by continuing EM iterations with current parameters until the iteration converges again.
Algorithm 1: The EM algorithm for learning HM-GPFR.
Initialize parameters θ; while not converged do
// E-step a1 (k: = TrN((11 Tbg, Cs);
for 1 2,..., do
for k = 1,,2.... K do
ax(k) = Ei= (1) NYy: Pbk. Ck);
end end
ST = 1:
for T 1,..., do
for k = 1,2, 1....... K do
3.(k) E/= Pu.N(uee Tb, C.) (1);
11
end
12
end
for t 1,..., T do
14
for k = 1,, ...,.... do (k)) ar ( for = 1,2...., K do
16
(k, () as (k) PNN Pb; Ci) (l);
end
Normalize (k, 1)
20
end
Normalize {7((k) end
21
22
// M-step
23
for k = 1, 2........ K do
21 (k)
Th& Σ (1))*
for = 1, 11,2....
...,, K do
T-
(k..)
Psi = T- Σx (i)km)
end
28
Using gradient ascent algorithm to optimize ((@) with respect to 02 and by according to Equation [1]];
29
Update C₁ with new parameters 0%:
30
end
end
12

---

Under review as a conference paper at ICLR 2023
A.2 BHM-GPFR
We still use the EM algorithm to learn the parameters of the BHM-GPFR model. However, this case more complicated since there are more latent variables. The complete data log-likelihood is
K
K
K
[(0;T,9)
log N(bb ma; Σb)
coo 1) log Pkt
k=
1=1
T-1 K
[(21 k)g log nī + ΣΣΣ
1( 1, 12)=k)ooo ko Pst
(18)
T K
+
I(2 = log log (y: bbk. Ck).
k=1
Compared with Equation (1O), the first two terms extra due to the prior distributions. In the E-step of EM algorithm, we need to take expectation of Equation (T8) with respect to the posterior distribution of latent variables. However, the posterior of Ω is intractable since (bx) ] {px) )= and {=+)== correlated. We use the variational inference method and try find optimal approximation of [ ()) with simple form. We adopt the mean-field family approximation, which factorizes the joint distribution of S2 product of several independent distributions,
K
K
(
II Q(bR) k=1
II Q(p:)((z). 1=1
(19)
Similar to the HM-GPFR case, (=)) is categorical distribution with KT possible values, but we do not need to calculate ((z) explicitly and only need to calculate Ye(k) Q(z = k) and (,k,() l) = Q(+41 2 k). According the variational inference theory, we iterate Q(b.).
((pc) and Q(z) alternately until convergence.
For Q(ba).
Q(be)) x exp FTK Q(p ((0:
exx Bg(z)
log log:V(bb; mp, ET) Σ 1l = k) Nog.Noo 2b,, Cx)
1 =(bs. (ymm))) ET,"" (bx ms)
(20)
exp
log EEb
2
Te(k)
(-1 log C+l
ž (yr - pbx) C (3/ 2bb))
By completing the square, obtain the approximate posterior of by is N(m:: Σx) with
T
Σx =
ΣB
ne (k) Φ C l
m Eg
Σ m + Σ
[kk CT C.. y/
(21)
For @(px).
((px) exp Ii ((0))()))(CC)))))))))))))))))))) ((0::)))
T-1 K
=exp EC(z)
(ao - 1) log Pat +
I(a411 = 1, = k)) log Pri
1=1
[=1
1=1
K
(22)
exp Σ
a0
log Pet +
Σ
(k, log Pi
K
p + (k(
13

---

Under review as a conference paper at ICLR 2023
Therefore, the approximate posterior of py is Dir(axa...... ................................ okk) with agd au ET= (k)
For Q2)),
(()) exp R Th% G(b.) ((v))
k
T-1 K
K
= exp 1 TT. ((bb) ((pe)
1[21 = log TTk +
[( = 4, = =k)=
=k)log log Pkt
T( k) log (yr Dbk, C+t)
t=1 k=
K
T- ΣΣΣ =1
exp
T(21 kt) log Th
(2++1=10 1,
Rop flog Pkd]
T
K
11 = k (lgg bby C+:]]
(23)
Note that this equation has exactly the same form Equation (10), thus we can the forward- backward algorithm to obtain (a()) and (k( 1). To see this, let
K
)
Pa = exp papp PRr] exp
(((aw))
asst)
(-) 1
Py: my Ex. 0x) = exp Ee(b.) [lo. ( Pb:, C&)l = NNyi Pm.z exp
ex)
( tr tr DC,
x0))
(24)
then Equation 23 can be rewritten as
T-1 K ΣΣΣ I 1(2-11 1,, log Pst
log (2))
Σ I(2) log TR. +
(25)
T K
1( log P(y m. Σx, 0n)
To obtain (k) and (ekk; l), we run the Belc-Welch algorithm with sufficient statistics Tk. Pki, P(y my Σg, 0).
Taking expectation of Equation 18 with respect to the approximate posterior Q(2), the Q function
K (
K Σ
T K
ms. Σb)] +
71 (k) log Te +
M>>
7e
Doog Pbx C+)
1=1
k=1
(log. N Vmm;i mé, Ex) ] E5 +
71 (k) log n
km
km/
T
(k) ((ye om,, C+)
tr Σx FC
(26)
Maximizing 2(0) with respect to Th. my and Eg, we obtain
(k)
1 Σ Σx + (mg ma) m²
Σn
m
m
Σ (1)
K
(27)
14

---

Under review as a conference paper at ICLR 2023
The parameters {0x) cannot be solved in closed form, and apply the gradient ascent algorithm to optimize Q(@). The gradient of ())) with respect to 0 is
(Q(@)
St, C OCk (
(k)tr "r
St.. (yr
(3 km;)) + 1-""E E₁ Ck.
00
00%
(28)
The complete algorithm is summarized in Algorithm 2
B MORED EXPERIMENTAL RESULTS
B.1 DETAILED EXPERIMENT SETTINGS
For AR, MA, ARMA, ARIMA, and SARMA, we set the model order L. in {4, 8, 12]. For SARMA, the seasonal length is set to be 96 since there 96 records per day, which implicitly assumes that the overall time-series exhibits periodicity in days. LSTM, NN, SVR, and EGPM transform the time-series prediction problem into regression problem, i.e., use the latest L observations predict the output at the next point and then use the regression method to train and predict. In the experiment, set L in {4, 12, 24, 48). The neural network in the FNN has two hidden layers with 10 and neurons, respectively. The kernel function in SVR is the Gaussian kernel whose scale parameters are adaptively selected by cross-validation. The number of components for EGPM is in (3, 5, 10). In addition, we use the recursive method ? for multi-step prediction. For mix- MPFF, mix- and DPM-GPFR, we first convert the time-series data into curve datasets and then use these methods to make predictions. The number of components K in mix-GPFR and mix-GPNN is set to 5 and the number of B-spline basis functions D in mix-GPFR and DPM-GPFR set to 30.
B.2 CLUSTERING STRUCTURE
Estimated values of latent variable & also indicate the evolution mode corresponding to the data of the i-th day. Figure 6 visualizes some training data with different colors indicating different evolution modes, so we can intuitively see the multi-scale structure the electricity load time- series. According the learned transition probability, we can obtain the stationary distribution of Markov chain (21, ......... 2N). which is (0.4825, 0.2026, 0.0513, 0.1124, 0.1513] in HM-GPFR, and [0. 4501, 0.0427, 0.2992, 0.1381, 0.0700 T in BHM-GPFR. The proportion of each mode in fig. is roughly consistent with the stationary distribution.
B.3 MULTI-STEP PREDICTION UNDER COLD-START SETTING
In order to more clearly see the role of Markov chain structure of hidden variables in the cold start setting, in Figure 7 and Figure 8 we show the predicted values of HM-GPFR and BHM- GFPR for electricity load in the next five days .N+ UN45 and distributions of latent variables AN+1..... 45 conditioned on y = It can be seen from the figure that HM-GPFR and BHM- GPFR have different predictions for each day's electricity load, which will be adjusted according to the transition probability of evolution law. For example, in Figure 7 when &n = 1, the power load that day is low, and the predicted value of HM-GPFR on the (N + 1)-th day is also low. When halzx = 5, the electricity load on that day is higher, and the predicted value of HM-GPFR on the (N + 1) day is also higher. Figure has similar phenomenon. In addition, it can be seen that with the increase of i., IP(,.) quickly converges to the stable distribution of the Markov chain, and the predicted value ÿi. also tends to be the weighted average of the mean function each GPFR component. In conclusion, these phenomena demonstrate that HM-GPFR and BHM-GPFR can effectively use the coarse- -raiined temporal structure adjust the prediction of each day.
B.4 SENSITIVITY OF PARAMEER
There are two main parameeas
HM-GPFR and BHM- GPFR: the number of B-spline basis functions D and the number of GPFR components K. Here mainly focus the selection of K. We vary in (3, 4, 5, 6, 7,8, 7,8,9 9, 10, 15, 20, 30, 50). train HM-GPFR and BHM-GPER respectively, and report the results in table 3
For HM-GPFR, its prediction performance tends to deteriorate with
15

---

Under review as a conference paper at ICLR 2023
Algorithm 2: The Variational EM algorithm for learning BHM-GPFR.
Initialize parameters θ;
while not converged do
/ Variational E-step
Initialize variational parameters {mg: Eg, 4,1, a....... a&k ke while not converged do
// Calculate surrogate parametezs.
for k= 1,2,.... ..... K do
PP(;;; m2 Σg. 0k) (y: Imgg C&) exp (-2ttrI,, I));
for 1 = 1.2..... K do
Pri = exp (((i)) 2 (E an)):
5
end end
// Forward-backware a algorithm (1(k) TV P (yr bu, Ox);
for = ...... do
10
for k = -1,2.... do
(k) E Ei=1 1())()) Ply Ex. 0x);
end end
13
14
3T 1:
for 1 T - -1,..., do
16
for kt = 1,2.... do
3e(k) = Σ Pw (1:+1; mi, E,, 0.)3++1 (I);
end end
20
for
1 ,,...,TTe do for k: = 1,2 2....... do (k) x ar (k) for 1,2, ...,, K do
21
(k, az ( m, E,, 01)3++1(1)
end
Normalize (k, )) end
27
Normalize {>x)))]] end
29
// Update posterior Q(b+) and Q(ps).
30
for k = ..,..... K do
Σn Eg Σt
(T) '
31
xT TT yr m, E₁ E m E (
33
33
for 1,2..... ... K do
asl = ap + ΣT- (k,
34
end
36
end
37
end
// M-step
38
E. 1 K
K k=
Eg + (m, mm) mm m) my
1 Σ K k=
m
39
for k = 1,2,
K do
(k)
TT = (1)
Using gradient ascent algorithm to optimize Q(@) with respect to O2 according to
Equation 28
Update C₁ with parameters 0k:
43
end
end
16

---

Under review as a conference paper at ICLR 2023
Ellctricity koad (million KW h)
10
int V d hol bro
2500
3500
4000
Time
HM-GPFR
Illctriciy lond (million kW h)
10
root nottoooo Certor oooonn Mo / 1 / P
4000
Time
(b) BHM-GPFR
Figure 6: Training time-series are divided into different evolving laws based on the learning results of HM-GPFR and BHM-GPER.
Table 3: Sensitivity of HM-GPFR and BHM- with respect to the number of components K.
Method
Step length S
10
20
30
300
100
200
300
300
1400
3
0.84%
1.015%
1.195
1.85%
1.5256
2.35%
8.80%
4.82%
6.00%
S.48%
9.63%
30.71%
9.4356
6.75%
6.735
0.0%%
1.09%
1.22%
1.44%
1.6115
4.0235
5.00%
6.6%
5.57%
9.68%
10.71%
0.4355
6.78%
6.73%
0.98%
1.125
1.355
1.48%
1.6956
2.51%
4.07%
5.18%
6.799%
8.855
9.85%
10.76%
9.4955
6.82%
6.77%
1.00%%
1.32%
1.57%
1881%
2.0455
3.135
4.77%
5.80%
7.3%
9.0455
10.0395
10.95
9.595%
6.88%
6.8%
1.00%
1.27%
1. 1005
1.7%
1.9156
2.87%
4.5%%
5.56%
7-28%
9.175
10.1.58
10.88%
9.5854
6.87%
6.795
HM-GIPER
8
0.97%
1.15
1.185%
1.33%
1 .016
1. 1.79%
1.54%
2.0055
1.70%
2.57%
3.063
4.14%
4.77%
5. 3.225
6.009%
7. 7.40%%
9.35
S:5556
10.3095
9.92%
10.81%
10.09%
9.5756
6.88%
6.88%
6.795
6.5%
10
1.18%
1.41%
605
1.885
2.11%
3.22%
6.085
7.535
0.5990
10.4598
11.17%
9.71%
0.05%
6.83%
1.945
4.17%
3.29%
6.32%
7.09%
6.36%
1.31%
1.545
1.00%
1.74
5.12%
1.00%
5.38%
4.225
3.69%
8.30%
5.07%
6.14%
7. -09%
7.00%
0.55%
1003555
10.7890
9.555
7.01%
7.01%
1905
2.25%
2.5%%
3.7%
6.62
1.20%
103555
11.007 1
12.01%
10.00%
7.12%
6.91%
1.105
1,305%
1.525
2.35%
4.55%
6,11%
30.15
9.47%
6.842
0.78%
1.00%
1.18%
..0%
3.82%
4.8%
8-25%
9.0455
6.08%
6
0.77%
0.025
1.0%%
1.88%
1.30
500%
5.88%
50%
1.00%
8.047
9.855
0,2155
6.045
7.135
0.8%
0.08%
1.1%
1.25%
1,5955
2.00%
3.17 1.77
97%
5.22%
8.522
9.0155
0.395
7.01%
7.185
7
0.70%
11.0%%
1.15
1,225
1.53%
1.94%%
3.019
3.79
5.12%
7.999%
9.70%
9.3458
7.18%
T,89F
RHM-CIPER
8
9
0.78%
0.78%
0.94%5
0.03%
10855
1005%
1.19%
1148
1.5115
1.29%
1.8%%
188%
2.94%
2.865
3.61%
5.71%
6.03%
1.2%%
6.74%
6.895
7.77%
T.29%
9.73%
9.7%
9.4955
0.525
7.46%
7.535%
7.8%
1.7%
10
0.82%
0.08%
1 11155
1365
1.4%
2.11%
3.29%
4.00%
5.379
7.01%
8.045
9.9455
9.895
7.80%
6.125
15
0.70%
11.4%%
1025
1.88%
1.296
1.84%
2.80%
3.64%
6.805
7.7%
9.89%
9.9654
8.25%
8.5%
20
0.79%
0.34%
1025
1.17%
1.288
83%
2.83%
3.6%
1.88%
6.51%
7.5%
9.963
10.325
8.88%
9.5155
30
0.85
0 1.95%
1075
1.88
1999
805
2.82%
5.58%
1 885
8.52%
7.58%
10.405
9.07%
9.52%
90
0.83%
0.0 1.86%
1.1%
1.22%
1339
88%
2.9%
3.68
1. 196%
6.3%
7.46%
10.145
10.71%
9.45%
0.95
the increase of K. In short-term prediction, MAPE increases significantly, while MAPE changes less long-term prediction. With the increase of K, the number of parameters in the model also increases, and the model tends to suffer from over-fitting, For BHM-GPFR, with the increase of K, its long-term prediction performance decreases significantly, while the medium-term and short-term prediction results do not change much. This shows that BHM-GPFR can prevent overfitting to certain extent after introducing prior distributions to parameters. In addition, we also note that when K 10, the difference between the results corresponding to different K is not significant, which is more reasonable choice. From the perspective of application, we set K 5 in the experiment, which take both expression ability and interpretability of the model into consideration.
17

---

Under review as a conference paper at ICLR 2023
2N|
Mean function of N-th day
Predicted values of N L, ................................
Electricity load (million kW
Klean fanction of With
Electricity (million
1.00
0.54
0.49
0.48
0.48
0.48
0.00
0.17
0.220
0.20
0.220
0.20
0.00
0.04
0.05
0.05
0.05
0.05
0.00
0.13
0.11
0.11
0.11
0.00
0.12
0.15
0.15
0.15
0.15
10
15
20
10
15
20
Eleenieity (calliion kW-h))
10-
Xiex fauction of With
Eleariatty load (million kW-h)
10:
Frttersonot Lahir
Predictoss of efl,ntee Prodiatioss afr,ntes
0.00
0.55
0.49
0.48
0.48
0.48
1.00
0.10
0.20
0.20
0.20
0.20
2
7
0.00
0.15
0.04
0..5
0.05
0.05
0.00
0.05
0.11
0.11
0.11
0.11
0.00
0.15
0.15
0.15
0.15
0.15
Time (h)
Time (h)
Electricity (mallion kW-h)
Moen of Medider
Electriitty (million kW-h)
0.00
0.59
0.48
0.48
0.48
0.48
0.00
0.20
0.19
0.20
0.20
0.20
3
.000
0.00
0.06
0.05
0.05
0.05
0.00
0.02
0.13
0.11
0.11
0.11
0.00
0.20
0.15
0.15
0.15
0.15
10
15
20
10
15
20
Eleetricity (milliion KW (W)
10
Nizz tascticn of 24th day
Electriity load (milliin kW-h)
10
Pediitios no
of
Predction ar1,
0.00
0.42
0.48
0.48
0.48
0.48
0.00
0.33
0.19
0.20
0.20
0.20
0.00
0.00
0.07
0.05
0.05
0.05
1.00
0.08
0.11
0.11
0.11
0.11
0.00
0.17
0.15
0.15
0.15
0.15
10
15
20
15
20
Time
Time
10-
Eleatriity (mnillion kW-h) h)
10
Prodution sfieeded afe
Electriiit (eniliion
of
0.00
0.22
0.45
0.48
0.48
0.48
0.00
0.36
0.22
0.20
0.20
0.20
51
0.00
0.00
0.06
0.05
0.05
0.05
0.00
0.21
0.11
0.11
0.11
0.11
1000
0.22
0.16
0.15
0.15
0.15
15 20
Time
Time (h)
Figure 7: Estimated values 1N+1. 4N45 and distributions of WN, where 1, 1, 5.
2N45 of HM-GPER under
18

---

Under review as a conference paper at ICLR 2023
2N|
Mean function of N-th day
Predicted values of N L,
Electricity load (million kW-))
Electricity load (million kWh)).
201
1.00
0.52
0.47
0.45
0.45
0.45
0.00
0.05
0.04
0.04
0.00
0.04
0.00
0.29
0.29
0.30
0.30
0.30
0.00
0.14
0.14
0.14
0.14
0.00
0.03
0.06
0.07
0.07
0.07
Morn Fanction N:B day
10
15
20
10
15
20
Eleenieity (calliion kW-h)
10-
Eleariatty load (million kW-h)
10-
Predittoss of efl Prodiatioss wll niee
Prckction Lonkig
0.00
0.12
0.39
0.44
0.45
0.45
1.00
0.12
0.00
0.04
0.00
0.04
2
0.00
0.53
0.32
0..0
0.30
0.30
0.00
0.12
0.17
0.14
0.14
0.14
0.00
0.12
0.08
0.08
0.07
0.07
N-B
Time (h)
Time (h)
Electricity (mallion h)
Electriitty (million kW-h)
0.00
0.48
0.47
0.45
0.45
0.45
0.00
0.02
0.04
0.04
0.04
0.04
3
.000
0.25
0.28
0.30
0.30
0.30
0.00
0.22
0.12
0.14
0.14
0.14
0.00
0.04
0.08
0.07
0.07
0.07
Mon Functire ofe Nith day
10
15
20
10
15
20
Eleetricity (milliion KW (W)
10
10
Electriity load (million kW
Medittoo no
of
Predction af1,
0.00
0.44
0.39
0.45
0.45
0.45
0.00
0.04
0.05
0.04
0.04
0.04
0.00
0.24
0.36
0.29
0.30
0.30
1.00
0.04
0.15
0.14
0.14
0.14
0.00
0.24
0.06
0.07
0.07
0.07
fanction day
10
15
20
15
20
Time
Time
10
Eleatriity (mnillion kW-h)
10
Prodution sfieeded afa,
Electriiit (eniliion
0.00
0.07
0.42
0.44
0.45
0.45
0.00
0.07
0,003
0.04
0.04
0.04
0.00
0.58
0.,29
0.30
0.30
0.30
0.00
0.20
0.17
0.14
0.14
0.14
1000
0.07
0.09
0.08
0.07
0.07
Man fanction of Nots day
15 20
Time
Time (h)
Figure 8: Estimated values 9N+1..... under where == 3, 4,5.
4N and distributions of 2N NZ of BHM-GPFR
19


Reviews:
[{"review_id": "bcBp1RtTgA", "reviewer": "Reviewer_hJkz", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposes a mixture Gaussian process functional regression model coupled with a hidden Markov chain at the coarse level. A Gaussian process regression is applied to the intra-day time series and then inter-day correlations are modeled through the hidden Markov chain. A Bayes version (BHM-GPFR) of the proposed model (HM-GPFR) is proposed to further improve the robustness of the method. The effectiveness of the proposed is demonstrated through an application to the electricity load data.", "summary_of_the_review": "The paper is well-written and easy to follow. The novelty of the paper is a bit limited and the numerical support of the proposed method is not enough convincing.", "strength_and_weaknesses": "Strength: the paper is well-written and easy to follow.\nWeakness: the proposed method is only a straightforward extension of the existing literature, e.g., [1]. The only difference is to introduce of the correlation between the indicator variables via a Markov chain. The benefits of this modification, however, are not well justified.\n\n1. The performance of the proposed method is only evaluated through one data set and results are not very impressive compared to existing methods (mix-GPFR, DPM-GPFR). In addition, the real data has too many irregularities and it is hard to grasp the idea of what are the advantages and limitations of the proposed method. In my opinion, this can be better achieved through carefully designed synthetic experiments, especially considering that the work is purely empirical.\n\n2. In section 5.1, it is probably better to consider SARIMA model compared to the SARMA model because the process is clearly nonstationary even after taking into account the seasonality.\n\n3. In section 5.3, HM-GPFR and BHM-GPFR perform better when $r=1$. Can you study how the performances change as $r$ gradually increases? This perhaps will give more insights into why the proposed method work.\n\n4. In equation (6), the transition probability only depends on the previous state, which may not be very suitable for grid data. It may be a better idea to include other coarse-level covariates such as predicted temperature, day of the week, etc.\n\n5. The claim on page 2, \"while HM-GPFR can adjust the parameters incrementally without retraining the model \", is not correct. The model parameters still need to be updated when new data are available.\n\nReference:\n\n[1] Shi, Jian Qing, and B. Wang. \"Curve prediction and clustering with mixtures of Gaussian process functional regression models.\" Statistics and Computing 18, no. 3 (2008): 267-283.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "WvcShLIHoIA", "reviewer": "Reviewer_CSZt", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The submission proposed the hidden Markov based GPFR mixture model (HM-GPFR) by describing the time series data from the perspective of both fine and coarse level. The time series data at fine level is characterized by Gaussian process model and and that at coarse level is characterized by a hidden Markov process. To avoid overfitting, the work place priori on the model parameters and develop a Bayesian version. ", "summary_of_the_review": "The major concern of the submission is novelty. All key ideas of the proposed method: gaussian process regression, hidden Markov model, and EM (forward-backward) algorithm, are all existing in previous works. There is nothing new in the methodology in this submission. Another minor conceren is that the submission is not self-contained very well; it is confusing to see some notations without definition, e.g., $\\gamma_t(k)$ above Eq.(7) and below Eq. (9), $m_k$ and $a_{kl}$ below Eq.(9). I guess these may be defined in the appendix, but the main text should be written in a self-contained way. ", "strength_and_weaknesses": "Strength: The paper is easy to follow and the derivation of inference seems correct but I did not check it carefully.\n\nWeakness: The submission is lacking in novelty and not self-contained. ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "M_vkJKgW0h", "reviewer": "Reviewer_tgKE", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper proposes an interesting method for multi-scale time series based on Hidden Markov mixture models and Gaussian processes. Overall, the model is capable of fitting the data with switching dynamics, and experiments demonstrate good performance compared to the baselines. ", "summary_of_the_review": "The paper proposes an interesting model but, overall, lacks reproducibility details which are necessary for a publication. The settings of the baseline methods are not discussed in sufficient detail. Based on this, I cannot recommend acceptance at the current moment. ", "strength_and_weaknesses": "Strengths: \n-- The model combines some existing ideas from the literature into one complex model which is able to deal with multi-scale time series. \n-- A good level of interpretability is next to good performance.\nWeaknesses: \n-- The structure of the model is fairly complex, and modes of failure are not completely clear (i.e. what are the limitations of the proposed approach). \n-- I would expect the inference method to be fairly slow and slower than competitors; there is, unfortunately, no comparison regarding computational costs. \n-- Details regarding the specifications of the competitors are not clear. For example, did ARIMA include a seasonal component, how exactly was LSTM specified and trained, what was the choice of the kernels in GP-based approaches etc.? \n\nMinor: it would be useful if Table 1 and Table 3 highlighted the best-performing methods to get an overall picture better. \n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject