=Paper= {{Paper |id=Vol-3039/paper26 |storemode=property |title=Water Consumption Periodic Autoregressive Time Series Iterative Forecasting |pdfUrl=https://ceur-ws.org/Vol-3039/paper26.pdf |volume=Vol-3039 |authors=Taras Mykhailovych,Mykhailo Fryz,Iaroslav Lytvynenko |dblpUrl=https://dblp.org/rec/conf/ittap/MykhailovychFL21 }} ==Water Consumption Periodic Autoregressive Time Series Iterative Forecasting== https://ceur-ws.org/Vol-3039/paper26.pdf
Water Consumption Periodic Autoregressive Time Series
Iterative Forecasting
Taras Mykhailovycha, Mykhailo Fryza, Iaroslav Lytvynenkoa
a
    Ternopil ‘Ivan Puluj’ National Technical University, Rus’ka st., 56, Ternopil, 46001, Ukraine


                 Abstract
                 Problems of urban water provision enterprises, such as: geographical zoning, pressure
                 normalization, emergency detection, were stated. Water consumption of such an urban water
                 provision enterprise is being studied by the authors. Water consumption mathematical model
                 was presented and justified as cyclostationary. Hourly water consumption time series were
                 justified in terms of linear forecasting. Water consumption operative interval forecasting
                 method and its existing authors’ information technology were presented and revised.
                 Problems of existing information technology were stated. Improvements to information
                 technology were proposed. The advantages of the iterative forecasting method were
                 explained. Streaming algorithm of water consumption forecasting was proposed and
                 presented. Algorithm implementation and deployment preferences were proposed. Next steps
                 of the research were defined.

                 Keywords 1
                 Water consumption, cyclostationary process, operative linear forecasting, time series,
                 periodic autoregressive, prediction interval, streaming algorithm, iterative method

1. Introduction
   The practical value of the urban water consumption forecasting technology is providing efficient
management of clean water supply at the pressure, needed by customers [1, 2]. As per [3], an
information technology, which implements hourly water consumption operative interval forecasting
method, may solve the significant problems of urban water provision enterprises, such as:
geographical zoning, pressure normalization, emergency detection. It may be mostly useful in scaling
ahead the water pipeline systems and their emergency states detection.
   Subjects of study are as following.
    Water consumption — a stochastic process: an oscillation of water flow (measured in liters per
   hour) in time
    Hourly water consumption — a stochastic time series: hourly water volume (measured in liters)
    Unit water consumption — one short-term consumption session (i.e. temporary valve opening)
   by any impartible user (i.e. a person). Example of unit water consumption is “Mr. John Smith
   washes his hands after returning home at 9pm, Sep 18, 2021”
   First two items above describe the consumption from urban water provision enterprise in a broad
sense: by some user: a person, an apartment, a building, a zone (geographic or logistic group of
buildings) etc., and of any scope (any combination of consumption sessions by that user).
   Object of study is operative (one step ahead) interval forecasting. Resulting prediction intervals are
considered as highly likely expected water consumption bounds and are “red lines” of theoretically
normal operation of the water provision system. Falling out from the prediction interval is a display of
probable emergency, which may be false positive due to the human factor, yet another reason to
ensure the sanity of the related system node.
1
ITTAP’2021: 1nd International Workshop on Information Technologies: Theoretical and Applied Problems, November 16–18, 2021,
Ternopil, Ukraine
EMAIL: tarquas@gmail.com (T. Mykhailovych); mykh.fryz@gmail.com (M. Fryz); d_e_l@i.ua (I. Lytvynenko)
ORCID: 0000-0002-8138-3642 (T. Mykhailovych); 0000-0002-8720-6479 (M. Fryz); 0000-0001-7311-4103 (I. Lytvynenko)
              ©️ 2021 Copyright for this paper by its authors.
              Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
              CEUR Workshop Proceedings (CEUR-WS.org)
    In [3] the authors described an information technology: it was developed for specific needs of an
urban water provision enterprise. It allows the collection and processing of data from water
consumption monitoring hardware.
    The most computation-heavy part of information technology is PAR (periodic autoregressive
model) parameters estimation method, which requires inversions and multiplications of very big
matrices on a regular basis (every hour). This computation improves the PAR parameters estimates
after acquiring each new hourly water consumption value. The updated parameter estimates are then
used to compute prediction intervals for the next hours.
    Recurring processing of all available historical data for PAR parameters estimation have
introduced a problem: the computation became slower for each next cycle of estimation, and kept
requiring more and more runtime and computing resources. This is so, because the calculation was
done using all available historical data from scratch, without reusing the results from previous
computations.
    This paper proposes developing a method to update the PAR parameters estimates, without need to
recompute them from a whole range of historical data. This will be achieved by developing an one-
pass streaming algorithm for real-time PAR parameter estimation upon the iterative method of least
squares updating, proposed in [4].
    Another problem is current deployment of mentioned information technology on a custom
dedicated server, which is not reliable, because it has following limitations.
     Requires specific computational hardware (a video card)
     Requires extreme electric energy consumption (for video card)
     Requires 24/7 unchoked operation
     Is not horizontally scalable
     Is not backed by a replica server
    A new water consumption operative interval forecasting method, upon updating streaming
algorithm implementation (which is low-memory and low-CPU consuming), allows to take advantage
of cloud environments and move the deployment to a serverless computing platform, in order to do
the following.
     Conserve the energy resources
     Conserve the costs on hardware wearout and development operations
     Scale the computations horizontally by users
     Have reliable 24/7 service and backing from the platform
    An approach of deployment transition to such a platform is also proposed in this paper.

2. Related Works
     Recently known water consumption models, short-term forecasting methods, and current state of
research in the field were described in a digest [5] of considered publications. They were analyzed by
authors in [3], in respect to stated in [3] problems of urban water provisioning on an enterprise, and
are considered as incomplete for the following reasons.
      Not being stochastic
      Not being cyclostationary
      Having insufficient theoretical background
     In [3], the authors defined the water consumption as an additive stochastic mathematical model of
unit water consumptions. Each unit water consumption ξ̆ k (t ), k ∈ℤ is approximated by scaled
boxcar [6] function ξ̆ k (t )=α k Π τ , τ +~l (t ), defined by three parameters: start time τ k , τ k ∈ℝ ,
                                        k   k   k
                                                                     ~ ~
 τ m ≠ τ n ∀ m≠n, average flow α k , α k ∈ ( 0 ,∞ ) , and duration l k , l k ∈ ( 0 ,∞ ) . According to this,
water consumption process ξ (t ), t ∈ℝ, is a sum of unit water consumptions:
                              ∞            ∞                     ∞
                                                                               ~                      (1)
                      ξ (t )= ∑ ξ̆ (t )= ∑ α Π
                                   k                k    ~ (t )= ∑ α ϕ (t − τ ; l ),
                                                                                 k   k   k
                                                        τ k , τk + l k
                          k =−∞             k=−∞                         k =−∞
where
                                  ϕ (s ; l)= 1 , if s∈ [ 0 , l ) , l ∈ ( 0 , ∞ ) .
                                                                                                      (2)
                                                {
                                             0 , otherwise
   Process (1) was justified in [3] as conditional linear random process (CLRP) in terms of its
properties [7], and was represented∞ in its form as a Stieltjes integral, compatible with (1):
                                               ~                                               (3)
                            ξ (t )= ϕ (t− τ ; l )d π ( τ ) ,t ∈(−∞ ,∞),
                                             ∫                 π (t )             α
                                             −∞
                 ~
where ϕ (t− τ ; l π (t) ) — random function, a kernel [7]; π (t ) — nonhomogeneous simple Poisson
process with unit jumps in time moments τ π (t) , which values are used as indexes of unit water
consumptions, π (t )∈ ℤ; π α (t ) — generating process: compound Poisson process, having jumps of
size α π( t) at the same time moments as π (t ), probability P( π (0 )=0 )=1 , variance
 Var ( πα (t ))<∞ .
   According to [3], finite-dimension characteristic functions of process (3):
                                                                           n                          n        ∞
                                                                       i ∑ u ξ (t )
                                                                                  k   k
                                                                                                          ~
                                                                                        i ∑ u ∫ ϕ (t − τ ; l
                                                                                                           k )d π ( τ )
                                                                                                                    k
                                                                                                                          (4)
                                                                                                                        π (t k )   α

                   f ξ (u 1 ,u 2 ,... ,u n ; t 1 ,t 2 ,... , t n )=E e    k =1
                                                                                    =Ee              k=1       −∞
                                                                                                                        ,
                                                                            n
                                                       ∞ ∞ ∞            ix ∑ u k ϕ ( t k − τ ; y )
         ln f ξ (u 1 , u 2 , ... ,u n ; t 1 ,t 2 ,... ,t n )=∫ ∫∫ (e       k =1
                                                                                                 −1)d F α ( x)d F l ( y) λ ( τ )d τ ,
                                                       −∞ 0 0

are T-periodic by their time arguments:
            f ξ (u 1 ,u 2 ,... ,u n ; t 1 ,t 2 ,... , t n )= f ξ (u 1 ,u 2 ,... ,u n ; t 1 +T ,t 2 +T ,... ,t n +T ),n ∈ℕ, (5)
since intensity       λ ( τ )    of unit water consumers connection to a water provision system, is
(theoretically) daily periodic: λ ( τ )= λ ( τ +T ), T =24 h. Thus, process (3) is cyclostationary [7].
   In [3], hourly water consumption ξ t , t ∈ℤ was represented as periodic autoregressive (PAR)
time series (with period T =24 h), adequate to (3) in terms of linear forecasting approach [8]:
                                 t                                                          p
                                                                                                                           (6)
                          ξ t = ∫ ξ ( s ) ds ,t ∈ℤ , ξ t =mt + ξt , ξ t =∑ a k ,t ξ t− k +bt ηt ,
                                t−1                                                           k =1
where mt — periodic mean of hourly water consumption: periodic (with period T ) time series,
expectation E ξ t+kT =mt , t ∈[1 ,T ], k ∈ℤ and mt =m t+T ∀t ; ξ t — centered PAR time series with
period T , E ξ t =0 ,∀ t ; a k , t — a-parameters of PAR: weights of dependency from sample
 {ξ k∣k ∈[t− p , t−1]}, containing previous p number of values before ξ t , a k , t =a k ,t +nT ∀ t ,
 k ∈[1 , p], n∈ ℤ; p — PAR model order: number of a-parameters, p ∈ℕ ; b t — b-parameter of
PAR: periodic standard deviation of generating time series (a white noise b t ηt , which is a stochastic
component of PAR model), b t =b t +T ∀ t , E b t =0 , Var (b t )<∞; ηt — basic white noise, E ηt =0 ,
 Var ( ηt )=1.
    A method of hourly water consumption operative forecasting was proposed in [3]. According to
this method, an operative forecast (one hour ahead prediction value) is a point estimate ξ t of
stochastic value ξ t , and is determined from: recent p number of centered values
 {ξ k∣k ∈[t− p , t−1]} of hourly water consumption historical data, mean estimates m^ t , and a-
parameters estimates a^ k ,t :
                                                                p
                                                                                                  (7)
                                             ^ t + ξ^ t , ξ^ t =∑ a^ k ,t ξ t− k .
                                      ξ^ t = m
                                                                                  k =1
     For each period index t , mean estimates m          ^ t are determined from historical data as:
                                          q−1
                                        1                                                             (8)
                                  m^ t = ∑ ξ t+k T ,t =[ 1 ,T ] , m  ^ t =m^ t+nT , n∈ℤ ,
                                        q k =0
                                                           n −t −1                                    (9)
                                                      q= ξ      [T
                                                                         ,                ]
where q — size of historical data periodic sample for period t ; n ξ — size of all historical data;
 [ x ] — floor integer of x .
     The a-parameters estimates a^ k ,t are determined from centered (by periodic mean estimates)
historical data ξ t = ξt − m
                           ^ t ∀ t , by solving the matrix equation:
                                        A^ t =( W Tt Wt )−1 W Tt X t , t=[1 ,T ] ,                   (10)
                                           T
                                         ^ t =( a^ 1, p+t a^ 2 , p+t ... a^ p , p+t ) ,
                                        A
                                            X Tt =( ξ p+t ξ p+t +T                ... ξ p+ t+ ( q−1 ) T ) ,
                                         ξt + p−1               ξ t+ p−2              ...      ξt         (11)




    Mean estimates m
                            W t=
                                  (     ξ t+T + p− 1
                                            ...
                                                               ξ t+T + p− 2
                                                                   ...
                                      ξt +( q− 1) T + p−1 ξ t+ ( q−1 ) T + p−2
                                                                                      ...
                                                                                      ⋱
                                                                                            ξ t+T .
                                                                                              ...
                                                                                                    )
                                                                                      ... ξt +( q− 1) T
                       ^ t and a-parameters estimates a^ k ,t are improved every hour for corresponding
period index t .
    The b-parameter estimates b^ t are used once to determine the best-fit PAR order p by applying
the method of minimal description length (MDL) [9], and was calculated from centered (by mean
estimate) historical data and a-parameters estimates a^ k ,t :
                            q−1            n                        2
                                                                                                        (12)
               ^b p+t = 1 ∑ ξ p+ x T +t −∑ a^ k , p+t ξ p +x T +t− k ,t=[ 1 ,T ] , b^ t =b^ t+nT , n∈ℤ.
                       q−1 x=0(          k =1
                                                                                  )
    A method of hourly water consumption operative interval forecasting was proposed in [3].
According to it, each prediction interval is determined from prediction errors εt ,t ∈ℤ. Each
prediction error ε is a difference of hourly water consumption ξ from its most recent forecast ξ^
                   t                                                                          t                  t
(which is the only one in operative forecasting):
                                                εt =ξ t − ξ^ t .                                  (13)
      Analysis of (13) using known "portmanteau" test [10] has shown that its empirical distribution is
not Gaussian. Therefore, the interval forecasting is based on methods of non-parametric statistics.
Thus, for some preset confidence probability γ , prediction interval [ L t , H t ] of water consumption
 ξ t , is:
                                            L t = ξ^ t +h 1−γ ;                                   (14)
                                                                    ε,
                                                                          2

                                                       H t = ξ^ t + h 1+ γ ,
                                                                     ε,
                                                                              2
where Lt — lower bound of prediction interval; H t — upper bound of prediction interval; h ε , x —
x-quantile estimate of time series εt distribution (known to be symmetric):
                                                    ε ([ n x+1 ] )− ε( [ n ( 1−x )+1 ])             (15)
                                           h ε , x=        ε              ε
                                                                                        ,
                                                                   2
where ε( k ) — k-th order statistic of time series εt (array of all its values sorted ascending); n ε —
size of ε( k ) ; [ x ] — floor integer of x.
    The above hourly water consumption operative interval forecasting method was developed into
information technology and tested on an urban water provision enterprise. Two versions of software,
implementing the PAR parameter high-performance estimation parallel algorithm, had been
developed and split-tested. Those versions differ by the following parallel computation technologies.
     OpenCL — allows more precise computations (double-precision floating point), but requires
    specific drivers (i.e. CUDA for nVidia) and Node.js OpenCL bindings (“node-opencl” [11] was
    used)
     OpenGL — utilizes vertex shaders and GLSL and is widely available without specific
    computational drivers. It may use WebGL interface of modern web-browsers (this is handy to
    offload the computations onto the web clients), and headless Node.js module [12]; yet this method
    is limited in precision (usually only single-precision floating point computations are available)
    Despite the simplicity of OpenGL-based software integration, the authors considered it’s precision
as insufficient for big ranges of historical data, used for PAR parameters estimation, so this version of
software has been deprecated in favor of its OpenCL-based version.

3. Proposed methodology
   The improvement, described in this paper, concerns updating (8) and (10) without recalculation
from scratch. This is achieved by using iterative methods of real-time estimates updating,
programmed into a one-pass streaming algorithm, which is based on iterative data processing, without
need to store them all (only a small recent fixed-length subset of them is stored, which is a sliding
window) [13]. The OpenCL-based algorithm will be used once for initial a-parameter estimation (10);
all subsequent updating will be done using the proposed iterative method.

3.1. Mean estimates updating
   Mean estimation method, used in water consumption forecasting information technology, can be
depicted as follows: for period iteration q (which is algorithm iteration t +(q−1)T ), equation (8) is
represented as:
                                                 S                                             (16)
                                          ^ t = t , t ∈[1 ,T ],
                                          m
                                                  q
where S t — cumulative sum of data samples by period index t , defined as:
                                                    q−1
                                                                                               (17)
                                              S t =∑ ξt +k T .
                                                                     k =0
    Persistent variables t , q , and {S t∣t ∈[1 , T ]} are defined and preserved between algorithm
iterations in a fail-safe storage (i.e. “Redis” key-value storage), and in runtime cache. Such estimation
algorithm consumes only 2 integer, and T double precision floating point values of memory and
storage.

3.2. The a-parameters estimates updating
    According to [4], a-parameters estimates a^ k , t , k ∈[1 , p], t ∈[1 ,T ], can be improved, using least
squares updating method, which can be depicted in matrix form, for period iteration q (not less than
 p ), as:
                      ^ t ,q+1 =A
                      A             ^ t ,q +Pt ,q Q t ,q (Q Tt , q Pt , q Q t , q +1)−1 ( ξ p+t +qT −Q Tt ,q A
                                                                                                             ^ t , q ), (18)
                                                                         T                         T
                                                                                             −1
                               Pt , q+ 1=Pt ,q −Pt , q Q t , q (Q t , q Pt , q Q t , q +1) Q t , q Pt ,q ,
                                         Q Tt ,q =( ξ t +qT + p−1 ξ t+qT + p−2 ... ξ t+qT ) ,
                                                           q≥ p , t ∈[1 ,T ],
        ^
where A t ,q — column matrix of a-parameters estimates (10) before an update (for iteration q of
period index t ); A  ^ t ,q+1 — updated column matrix of a-parameters estimates (for iteration q +1 of
period index t ); Q t ,q — column matrix of size p×1 ; Pt ,q — symmetric square matrix of size
 p× p .
             ^ t ,q and P are defined as recurrence relations. According to [4], their initial values
    In (18), A               t ,q
can be calculated from scratch by using the historical data, using matrix equations, compatible with
(10):
                                                         ^ t , q =P t , q W Tt ,q Xt ,q ,
                                                         A                                                              (19)
                                                                                      −1
                                                         Pt ,q =( WTt , q W t ,q ) ,
                                                  ξt + p−1              ξt + p− 2          ...      ξt
                               W t ,q =
                                          (      ξt +T + p−1
                                                     ...
                                                                       ξt +T + p− 2
                                                                           ...
                                              ξt + ( q−1 ) T + p−1 ξt +( q− 1 ) T + p− 2
                                       X Tt , q =( ξ p+t ξ p+t +T ... ξ p+t +( q−1 ) T ) .
                                                                                           ...
                                                                                           ⋱
                                                                                                  ξt +T ,
                                                                                                   ...
                                                                                                          )
                                                                                           ... ξt + ( q−1 ) T


   According to (18), A     ^ t ,q can be additively updated. So, its current state can be preserved in
                         ^
persistent variables {A t∣t ∈[1 ,T ]}. In each algorithm iteration, deltas
                           ΔA^ t , q=Pt , q Q t ,q (QTt , q Pt , q Q t ,q +1)−1 ( ξ p+ t+ qT −Q Tt ,q A
                                                                                                      ^ t ,q) (20)
                                                                                          ^
may be calculated for period index t into a temporary matrix Δ A , and used to update matrix A t in-          ^
place. Variables {A ^ t∣t ∈[1 ,T ]} consume p T double precision floating point values of memory and
storage.
    The same, additive updating is also applied to Pt ,q , current state of which can be preserved in
persistent variables {Pt∣t ∈[1 ,T ]}. In each algorithm iteration, deltas
                             Δ P t ,q =−Pt , q Q t ,q (QTt ,q P t , q Qt ,q +1)−1 Q Tt ,q Pt , q                (21)
may be calculated for period index t into a temporary matrix Δ P , and used to update matrix Pt in-
place. Variables {Pt∣t ∈[1 ,T ]} consume p 2 T double precision floating point values of memory and
storage. Also, they can be compressed by deduplicating the symmetric entries of matrices, to consume
 T p (1+ p)
              double precision floating point values.
      2
    As it concluded from [4], the current state of matrices Q t ,q can be represented by a single sliding
window [13]: on each algorithm iteration, one element is appended to its rear (this value is the current
hourly water consumption value ξ ), and one element is removed from its front (the oldest value in it).
An iterable queue of fixed length p is enough to preserve the state of this sliding window between
algorithm iterations. This queue can be defined in persistent variable Q (which is considered in
calculations as column matrix, where the bottom entry of the matrix is the front element of the queue).
It consumes p double precision floating point values of memory and storage.
    Persistent variables t , Q , {Pt∣t ∈[1 ,T ]}, and {A                    ^ t∣t ∈[1 ,T ]} are defined and preserved
between algorithm iterations in a fail-safe storage, and in runtime cache. Such estimation algorithm
                                                       p (3+ p )
consumes 1 integer and T ( p 2 + p+1) (or T (                         + 1) with deduplication of symmetric entries)
                                                             2
double precision floating point values of memory and storage.

3.3. Streaming algorithm of operative forecasting method
      The a-parameters estimation algorithm uses mean estimates, and the operative forecast from
previous iteration. So, streaming algorithms for mean and a-parameters estimation are combined
along with an operative forecasting algorithm. The generalized algorithm is optimized by computation
complexity to minimize the number of floating-point divisions and multiplications, thus conserving
the cost of CPU usage in serverless computing environments, for which the algorithm is targeted. It’s
a one-pass streaming algorithm, which means that it inputs each hourly water consumption value only
once. This value is latched to a sliding window Q , and persists there only during next p iterations.
Thus, algorithm iterations of the improved forecasting method don’t require the whole array of
historical data values anymore.
      Streaming algorithm of water consumption operative interval forecasting has two phases [14]:
setup and iteration. Setup phase is performed at start; then, after its successful completion, algorithm
iterations (recurring runs of iteration phase) are performed. The algorithm doesn’t need a wrap-up
phase, as persistent variables are stored after each iteration to a fail-safe storage. Algorithm can be
applied to any water consumption user, for which the scope of available contiguous historical data
 ξ t ,T ∈[1 , n ξ ] allows inversion in matrix Pt , q calculation (19), i.e. matrix (WTt ,q Wt , q ) is not
singular.
      Subsequent iteration phase runs for the same user must not overlap or fail: each next run must
happen only after successful completion of its predecessor. So these runs must not be triggered
directly by hourly water consumption value arrival, but backed by fail-safe execution queue (i.e.
“Google PubSub”); so that arrival event pushes the input value to this queue and then the queued
values get processed strictly sequentially by the platform. If the run fails, it should be retried until the
success. Such architecture may lead to the chokes, especially when backing services experiencing
long-term connectivity issues. To minimize such risks, all backing services of information technology
(i.e. fail-safe storage) must persist on the same platform, region, and network; and forecast demanding
client should either be backed by outgoing queue, or do not require contiguous forecast delivery
(treating delivery failures as success, allowing algorithm to continue).
      This algorithm may be horizontally scaled by water consumption users, as their water
consumptions (and thus, parameters and forecasts) are independent.
3.3.1.      Setup phase
    This phase must be performed prior to recurring runs of the iteration phase for each user. It
organizes the operating state of computing environment, and prepares the persistent variables by the
following actions.
     Reserving space for them in fail-safe storage
     Reserving space for them in runtime cache (if setup and iterative phase runtimes are sharing
    operative memory within one execution space)
     Setting their initial values
     Saving them to runtime cache and fail-safe storage
     Using them for initial estimation and forecasting (optional)
    Setup phase steps of water consumption streaming algorithm are as follows.
    1. Input initial historical data array ξ t ,T ∈[1 , n ξ ]
    2. Define persistent variables t , q , {S t∣t ∈[1 , T ]}, {m   ^ t∣t ∈[1 ,T ]}, Q , {Pt∣t ∈[1 ,T ]},
       ^                     ^
     {A t∣t ∈[1 ,T ]}, and ξ ; and reserve them in storage
    3. Calculate period index t :=(n ξ −1) mod T +1
                                                  n −t−1
    4. Calculate period iteration index q := ξ      [T          ]
    5. Calculate cumulative sums:
              q−1

      {
   a. S u =∑ ξ t +k T∣u∈[1 , t ] ,
              k=0
              q− 2
                                     }
      { ∑
   b. S u =
              k=0
                     ξ t +k T∣u∈[t +1 ,T ]}
   6. Calculate mean estimates:
               S
      {
   a. m  ^ k = k ∣k ∈[1 ,t ] ,
                q              }
                 S
      {
   b. m  ^ k = k ∣k ∈[t +1 ,T ]
               q−1                    }
   7. Define matrices {W t∣t ∈[1 ,T ] } and column matrices {X t∣t ∈[1 ,T ] } as views of historical data
   array ξ t ,t ∈ ℤ, using functions:
    W t ( y , x )=ξ t+ p− x+ ( y −1) T = ξt + p−x +( y−1) T − m ^ t + p−x ,
    X t ( y )=ξ p+t+( y−1) T =ξ p+t+( y −1) T − m   ^ t+ p− x ,
   where W t ( y , x ) — row y column x entry of matrix W t ; X t ( y ) — row y entry of X t
   8. Define matrices {W Tt ∣t ∈[1 ,T ] } as views of matrices {W t∣t ∈[1 ,T ] }, using function
    V t ( y , x)= Wt ( x , y),
   where V t ( y , x) — row y column x entry of W Tt
   9. Define variables {Sk∣k ∈[1 ,T ]}, and calculate them, using a parallel computing algorithm for
   matrix multiplication [12], as:
   a. {Sk =WTk , q W k ,q∣k ∈[1 ,t ]},
   b. {Sk =WTk , q−1 W k ,q−1∣k ∈[ t+ 1 ,T ]},
   where W k ,n = Wk ( y , x )∣x∈[1, p] , y ∈[1 ,n] ; W Tk ,n =Vk ( y , x )∣x ∈[1 , p], y∈[1, n]
   10.        Calculate matrices {P k =S−1            k ∣k ∈[1 ,T ]}, using a parallel computing algorithm for matrix
   inversion [12].
   11.        Create queue Q , and fill it with elements {ξ t+qT +k =ξ t+qT +k −m                ^ t +k∣k ∈[0 , p−1]}, inserted
   in order of growing k
   12.        Calculate a-parameters estimates:
   a. {A  ^ k =P k W Tk ,q X k ,q∣k ∈[ 1 ,t ] },
   b. {A  ^ k =P k W Tk ,q−1 X k , q− 1∣k ∈[t +1 ,T ]},
    where X k ,m = Xk ( y )∣y ∈[1 ,m]
    13.      Calculate operative forecast ξ^ =m  ^ t +Q T A
                                                          ^ t . Output it to the client
    14.      Receive feedback from client, that output was either successfully consumed or failed. This
    step is optional
    15.      Store modified persistent variables t , q , {S t∣t ∈[1 , T ]}, {m          ^ t∣t ∈[1 ,T ]}, Q ,
                         ^                    ^
     {Pt∣t ∈[1 ,T ]}, {A t∣t ∈[1 ,T ]}, and ξ into runtime cache and fail-safe storage
    16.      Allow iterator phase runs
    Variables S u , m ^ k , Sk , and A ^ k are calculated in two passes, because they can have different q
within period index intervals [1 , t ] and [t +1 ,T ], for historical data not aligned by period T . If
historical data are aligned by period T , steps 5.b, 6.b, 9.b, and 12.b are not performed.
    Step 14 is optional, and allows the algorithm to continue after a failure. This is so, because the
failure to deliver this initial forecast won’t violate contiguity of subsequent forecasts (which may be
required by the client).
    There is no concern about computation complexity of setup phase, as this phase is performed once
per algorithm run (thus once for each water consumption user), and affords more budget for
computation.

3.3.2.     Iteration phase
    After acquiring each next hourly water consumption value ξ from the sensor, the iteration phase
is performed. Persistent variables get loaded by either of the following ways.
     From fail-safe storage — for first run on new execution context (and then are preserved in
    runtime cache)
     From runtime cache — for recurring runs on the same execution context
    As the iteration phase is recurring in real-time, it’s memory consumption and computing
complexity are minimized as much as possible.
    Iteration phase steps of water consumption streaming algorithm are as follows:
    1. Input hourly water consumption value ξ
    2. Load persistent variables t , q , {S t∣t ∈[1 , T ]}, {m      ^ t∣t ∈[1 ,T ]}, Q , {Pt∣t ∈[1 ,T ]},
     {A^ ∣t ∈[1 ,T ]}, and ξ^
       t
   3. If t=T : reset period index t :=1, increase period iteration index q :=q+1; else: increase
   period index t :=t +1
   4. Update cumulative sum for period index t : S t :=S t + ξ
                                                                   S
   5. Update mean estimate for period index t : m        ^ t := t . This needs one floating-point division
                                                                   q
   6. Define and calculate variable R=P t Q . Multiplication of these matrices may be done using the
   classic algorithm. It needs p 2 multiplications of floating-point values
   7. Define and calculate variable Y=R (Q T R+1)−1 . Calculation of this expression needs 2 p
   multiplications and one division (to calculate reciprocate number) of floating-point values
   8. Define and calculate variable ε= ξ −ξ^ k . Its value is operative forecasting error
   9. Define and calculate variable Δ A=Y^     ε . This needs p floating-point multiplications
   10.      Update a-parameters estimate A   ^ t := A ^ t +Δ A   ^
   11.      Define and calculate variable Δ P=−Y R T . This needs p 2 floating-point multiplications
   12.      Update matrix Pt :=Pt +Δ P
   13.      Push ξ to rear of queue Q , and remove one element from its front
   14.      Calculate operative forecast ξ^ =m   ^ t +Q T A  ^ t . This needs p floating-point multiplications.
   Output it to the client
   15.      Receive acknowledgment (feedback of success) from client, that output was successfully
   consumed. If it has failed, terminate the algorithm run
   16.      Store modified persistent variables t , q , {S t∣t ∈[1 , T ]}, {m              ^ t∣t ∈[1 ,T ]}, Q ,
                       ^                   ^
    {Pt∣t ∈[1 ,T ]}, {A t∣t ∈[1 ,T ]}, and ξ to runtime cache and fail-safe storage
    Computation of iteration phase needs 2 p 2 +4 p multiplications and two divisions of floating-
point values.
    Step 15 preserves the client’s concern about contiguity of subsequent forecasts. If the client
doesn’t have such a concern, this step does not terminate the algorithm run on forecast delivery
failure.
    Operative forecasting error ε (calculated on step 8) will be useful in further development of this
algorithm to support interval forecasting.

3.4. Algorithm implementation and deployment
   In ECMAScript (JavaScript), which is used as a high-level programming language in water
consumption forecasting information technology, the streaming algorithms can be effectively
implemented using asynchronous iterators [15]. The information technology is proposed to be
deployed using Cloud Functions on Google Cloud Platform [16], selected as best-fit serverless
computing environment. It also has a rich variety of backing services, which may be used as in-
platform fail-safe storage (“Datastore”, “Firestore” etc.).

4. Conclusion
    Since memory consumption of iteration phase is constantly low, computing operations are simple
(not so many multiplications and couple of divisions), and computing complexity is low, the
algorithm iteration phase of information technology is suitable to be deployed on wide range of
computing environments: dedicated servers, serverless clouds, and even embedded devices. For
instance, in the studied water provision enterprise [3], this phase can be integrated directly into Fargo
Maestro 100 GSM modems, which are used to transfer the hourly water consumption data.
    Setup phase (even with any number of subsequent iterations) may be done externally, and then
initial values of persistent variables may be copied into fail-safe memory of embedded device before
running iteration phase on it. This may be used for off-line automated emergency detection, and for
sealing the potentially damaged nodes of the water provision system.
    Authors plan to improve this algorithm to support water consumption operative interval
forecasting. This will require implementing a non-parametric quantile estimation method into both
phases.

5. References
[1] A. Candelieri, D. Soldi and F. Archetti, Short-term forecasting of hourly water consumption by
    using automatic metering readers data, Procedia Engineering 119 (2015) 844–853.
    doi:10.1016/j.proeng.2015.08.948.
[2] E. Pacchin, S. Alvisi and M. Franchini, A Short-Term Water Demand Forecasting Model Using a
    Moving Window on Previously Observed Data, Water 9 (2017) 1–15. doi:10.3390/w9030172.
[3] T. Mykhailovych, M. Fryz, Model and Information Technology for Hourly Water Consumption
    Interval Forecasting, in: Proceedings of the 15th International Conference on Advanced Trends
    in Radioelectronics, Telecommunications and Computer Engineering, TCSET 2020, IEEE, Lviv-
    Slavske, Ukraine, 2020, pp. 341–346. doi:10.1109/TCSET49122.2020.235452.
[4] A. Kapustinskas, A. Nemura, Identification of Linear Random Processes, Vilnius „Mokslas"
    Publishers, 1983.
[5] R. Benítez, C. Ortiz-Caraballo, J. C. Preciado, J. M. Conejero, F. S. Figueroa, Alvaro Rubio-
    Largo, A short-term data based water consumption prediction approach, Energies 2019, Vol. 12,
    2359 (2019), doi:10.3390/en12122359.
[6] E. W. Weisstein, Boxcar Function. From MathWorld — A Wolfram Web Resource, 2013. URL:
    https://mathworld.wolfram.com/BoxcarFunction.html .
[7] M. Fryz, B. Mlynko, Properties of Stationarity and Cyclostationarity of Conditional Linear
    Random Processes, in: Proceedings of the 15th International Conference on Advanced Trends in
     Radioelectronics, Telecommunications and Computer Engineering, TCSET 2020, IEEE, Lviv-
     Slavske, Ukraine, 2020, pp. 166–171. doi:10.1109/TCSET49122.2020.235415.
[8] W. A. Gardner, A. Napolitano, and L. Paura, Cyclostationarity: Half a century of research, ISSN
     0165-1684: Signal Processing, Elsevier, vol.86, no. 4, 639–697, 2006.
[9] S. L. Marple, Digital Spectral Analysis with Applications, New Jersey: Prentice Hall, 1987.
[10] A. I. McLeod, Diagnostic checking periodic autoregression models with application, ISSN 0143-
     9782: Journal of Time Series Analysis, vol. 15, no. 2, 221–233, 1995.
[11] Mikeseven, Node OpenCL, 2020. URL: https://github.com/mikeseven/node-opencl .
[12] T. Mykhailovych, GPU-EZ: Easy GLSL parallel computing on GPU.JS, 2019. URL:
     https://github.com/tarquas/gpu-ez .
[13] N. Alon, Y. Matias, M. Szegedy, The space complexity of approximating the frequency
     moments, Journal of Computer and System Sciences, 58 (1999) 137–147.
     doi:10.1006/jcss.1997.1545.
[14] M. Wolf, High-Performance Embedded Computing, ISSN 978-0-12-410511-9: Morgan
     Kaufmann, Elsevier, Second Edition, 139–200, 2014.
[15] I. Kantor, Async iteration and generators, 2021. URL: https://javascript.info/async-iterators-
     generators .
[16] Google Inc., Cloud Functions, 2021. URL: https://cloud.google.com/functions .