# Geometric interpretation of phyllotaxis transition

Takuya Okabe  Faculty of Engineering, Shizuoka University, 3-5-1 Johoku,
Hamamatsu 432-8561, Japan<EMAIL_ADDRESS>
###### Abstract

The original problem of phyllotaxis was focused on the regular arrangements of
leaves on mature stems represented by common fractions such as 1/2, 1/3, 2/5,
3/8, 5/13, etc. The phyllotaxis fraction is not fixed for each plant but it
may undergo stepwise transitions during ontogeny, despite contrasting
observation that the arrangement of leaf primordia at shoot apical meristems
changes continuously. No explanation has been given so far for the mechanism
of the phyllotaxis transition, excepting suggestion resorting to genetic
programs operating at some specific stages. Here it is pointed out that
varying length of the leaf trace acts as an important factor to control the
transition by analyzing Larson’s diagram of the procambial system of young
cottonwood plants. The transition is interpreted as a necessary consequence of
geometric constraints that the leaf traces cannot be fitted into a fractional
pattern unless their length is shorter than the denominator times the
internode.

###### keywords:

Schimper-Braun’s law, Populus deltoides, Fibonacci numbers, golden ratio,
phyllotaxy

##  1 INTRODUCTION

The spiral arrangement of leaves on a stem, phyllotaxis, is represented by the
fraction of the circumference of the stem traversed by the spiral in passing
from one leaf to the next. Braun and Schimper discovered that leaves are lined
up in ranks parallel to the stem so that the fraction is literally represented
by a common fraction (Braun (1835)). In a 3/8 phyllotaxis, for instance, every
eighth leaf comes over one below it after three turns of the spiral, so that
eight straight ranks are visible along the stem (Fig. 1). To this day, the
following list of the most common fractions has been circulated in books and
websites for non-specialists; 1/2 for elm, lime and linden, 1/3 for beech and
hazel, 2/5 for oak, cherry, apple, holly and plum, 3/8 for poplar, rose and
pear, 5/13 for almond, etc. Such correspondence tables seem to have existed
already in the middle of the nineteenth century (Henfrey (1870)). Sometimes
willow is listed in 5/13 (Coxeter (1961)) and sometimes in 3/8 (Adam (2006)).
As a matter of fact, it had been remarked since early times that even an
individual plant sometimes makes transitions between different fractions
(Braun (1835)). Therefore, the phyllotaxis fraction is not a determined
characteristic of each species. Most notably, Larson (1980) has revealed the
manner in which the vascular system is rearranged through the phyllotaxis
transition by mapping arrangement of the leaf traces, the portion of vascular
bundles of leaves that resides in the stem (Fig. 2). By contrast, it has been
commonly accepted that the arrangement of leaf primordia in the bud or at the
shoot apical meristem does not conform to any fractional number; primordia do
not appear as radial rows. The angular divergence between successive primordia
stays close to a unique “ideal” angle of 137.5∘, which is the golden mean
0.3820 of 360∘ (Church (1920); Richards (1951)). The golden mean is a
mathematical limit of the sequence of the above fractions. Much attention has
been paid to mathematics of these numbers and to mechanisms of leaf primordia
formation by which the ideal angle is achieved and regulated (Adler et al.
(1997)). In contrast, the phyllotaxis transition has been left unexplained
without attracting interest from researchers. The fractional phyllotaxis and
the phyllotaxis transition are two sides of the same problem. A key point
noted in the present study is that phyllotaxis in the bud changes continuously
whereas phyllotaxis on the mature stem changes stepwise. Evidently, the
arrangement of primordia at the apical growing point must be a determining
factor of leaf arrangement on the mature stem. But what causes the transition
from 2/5 to 3/8, for instance? There must be another factor. Larson (1980) has
suggested that the transition is programmed in the plant to occur at specific
stages of ontogeny. In another view, phyllotaxis of primordium formation and
vascularization are controlled by some higher-level system (Romberger et al.
(1993)). Kuhlemeier (2007) remarks that virtually nothing is known about the
molecular mechanisms that underlie the transitions between different spiral
systems (e.g. 3/8 to 5/13 patterns), except that larger meristems seem to have
higher Fibonacci numbers. The present paper aims at pointing out that this
phenomenon is consistently explained by considering geometry of growing leaf
traces. By means of a full quantitative analysis, which surely is not standard
in this field of research, it is shown (1) that the positions of the
phyllotaxis transition are located based on the angular positions at which
leaf traces exit from the vascular cylinder and (2) that the phyllotaxis
transition is caused as a necessary consequence of change in size of leaf
traces relative to internode length.

Figure 1: A young poplar in a 3/8 phyllotaxis with eight vertical ranks of
leaves.

A caveat: this paper deals with the mechanism of the phyllotaxis transition in
spiral systems. Although a common keyword of “phyllotaxis” may suggest, it
should not be confused with mechanisms of phyllotactic primordium formation,
for which considerable advances have been made over the last decade
(Kuhlemeier (2007)). Models of the latter category deal with continuous
changes of apical meristems, while they are not concerned with the fractional
expression of a phyllotactic pattern.

##  2 MATERIAL AND METHODS

Figure 2:  Diagram of the procambial system of a typical cottonwood plant
compiled by Larson (1980). The vascular cylinder is displayed as if unrolled
and laid flat. The ordinate is the Leaf Plastochron Index (LPI) for leaves,
whereas the abscissa corresponds to the angular coordinate in a full turn
about the stem axis. Each leaf is entered by three traces; a central
($\times$), right ($\blacktriangle$), and left ($\triangle$) traces.
Phyllotaxis orders, 1/2, 1/3, 2/5, 3/8 and 5/13, are indicated by the right
vertical axis.  Figure 3:  Vertical length of leaf traces in Fig. 2 is plotted
against the Leaf Plastochron Index. Orders of phyllotaxis, 2/5, 3/8 and 5/13,
are denoted at the same position as indicated in Fig. 2. See Fig. 2 for trace
designations. The phyllotaxis transitions from 2/5 to 3/8 and from 3/8 to 5/13
are caused by the traces indicated by a solid arrow and a dashed arrow,
respectively. Horizontal dashed lines at 5 and 8, separating the different
phyllotaxis regimes, are drawn for reference sake (see Results and
Discussion).

The analysis is based on a diagram of the procambial system of a cottonwood
plant (Populus deltoides) reconstructed by Larson (1980), which is reproduced
in Fig. 2. The vascular cylinder is displayed as if unrolled and laid flat.
The ordinate is the leaf plastochron index (LPI) for numbering leaves. Each
leaf has three traces: central, right and left traces exit the vascular
cylinder at positions denoted by symbols $\times$, $\triangle$ and
$\blacktriangle$, respectively. For more details, see Larson (1980) and
references cited therein.

By virtue of the convention to take LPI as the ordinate, the vertical scale of
the diagram is normalized such that differences in height between two
successive leaves, or internodes, are a unit of length in the vertical
direction. Therefore, the vertical component of a line segment in Fig. 2
represents not its actual length but an effective length relative to the
internode length, namely the length measured in internode units. Accordingly,
vertical lengths of the leaf traces in internode units are directly evaluated
by applying a digitizing ruler to the diagram. The lower end points of right
and left traces are located without ambiguity, for they are connected to other
types of traces. Lengths of central traces are fixed by decomposing the whole
pattern into clusters consisting of three adjacent traces, the right trace of
leaf $n$, the central trace of leaf $n+2$ and the left trace of leaf $n-1$,
where $n$ for LPI is an integer. For illustration purposes, the clusters from
$n=0$ to 7 are shown in Fig. 5. The length of the leaf traces thus obtained is
plotted against LPI in Fig. 3, where reading errors are of the order of an
internode at most. On the other hand, divergence angles are evaluated from the
horizontal, angular coordinates of the exit points of the leaf traces denoted
by the symbols. Evaluated angles show rhythmic, systematic variations around
the “ideal” angle, which are typically observed in a quantitative analysis
(Okabe (2012a)). The systematic variations, which are correlated with the
angular positions, are suppressed apparently by deleting from the diagram a
blank rectangular strip of a narrow width along the left ordinate. The width
of the strip is determined so as to minimize the standard deviation of
divergence angles for ten youngest central traces (LPI less than $-6$).
Results for the divergence angle thus corrected are shown in Fig. 4. The width
ratio of the deleted strip is 0.023, whereby the standard deviation of the
divergence angle is suppressed from 4.2∘ to 0.89∘. The correction is made for
ease of understanding implications of Fig. 4. It does not affect the results
discussed below qualitatively.

##  3 RESULTS

Figure 4:  Divergence angle between the traces at LPI $n$ and $n+1$ is plotted
against LPI $n$. See Fig. 2 for trace designations. The inset shows the angle
between the central traces at LPI $n$ and $n+5$: for ideal patterns of 2/5,
3/8 and 5/13 orders, it should be zero, 45 and 28 degrees, respectively. The
positions at which the phyllotaxis transitions from 2/5 to 3/8 and from 3/8 to
5/13 occur are indicated by solid and dashed arrows, respectively.

Results in Fig. 3 indicate that the three traces grow in length steadily all
alike. This behavior is consistent with other quantities reported by Larson
(1980). According to Fig. 3, three traces at LPI 2 have almost the same length
of about five internodes. This means that the traces at LPI 2 extend down to
the height of LPI 7 (see Fig. 2). As indicated by the right ordinate of Fig.
2, the phyllotaxis order fraction changes from 2/5 through 3/8 to 5/13 as we
climb up the stem, or as LPI decreases. The order fractions are indicated at
the top of Fig. 3 at the same LPI coordinates as in Fig. 2 by Larson (1980).
Horizontal dashed lines at five and eight internodes in Fig. 3 are drawn to
separate regimes of different phyllotaxis orders (see below).

Divergence angle between the leaf traces at LPI $n$ and $n+1$ is plotted
against LPI $n$ in Fig. 4. The inset shows the angle between the central
traces at LPI $n$ and $n+5$, that is, the net angle of inclination of
5-parastichies. The angle should become zero, 45 and 28 degrees for 2/5, 3/8
and 5/13 ideal patterns, respectively; for instance, the ideal angle for 3/8
is $360\times 3/8\times 5=360\times 2-45$, which is congruent to $-45$. The
inset of Fig. 4 clearly indicates stepwise transitions between the three
distinct fraction regimes. Thus, the positions of the phyllotaxis transition
are located based on the exit points, or the bases, of the leaf traces, i.e.,
without inspecting internal changes in the vascular structure. This crucial
property for us is brought to light probably because the subject plants are
grown under controlled uniform conditions (Larson (1980).

The cause of the transition is traced back by close inspection based on the
quantitative results. The transition from 3/8 to 5/13, indicated by a dashed
arrow in the inset of Fig. 4, is caused by an irregular decrease of divergence
angle at LPI 2, which therefore is also indicated by a dashed arrow in the
main figure. Similarly, a solid arrow in the inset indicates the transition
from 2/5 to 3/8, which is ascribed to an increase of divergence angle at LPI
7, a solid arrow in the main figure.

In connection with the latter transition, it is worth a remark that divergence
angle in the 2/5 phyllotaxis regime is not held at an ideal constant value of
144∘ $(360\times 2/5=144)$. This means that five ranks (orthostichies) of a
real 2/5 pattern is not equally spaced. According to Fig. 4, two full turns
($360^{\circ}\times 2$) of the 2/5 pattern is divided roughly into unequal
parts of $140^{\circ}\times 4+160^{\circ}\times 1$, instead of a regular
spacing with $144^{\circ}\times 5$. The figure shows that the irregular shift
at LPI 7 is shared with LPI 12. This is just as expected for the 2/5
arrangement $(7+5=12)$. Similarly, a cycle of $360^{\circ}\times 3$ of a 3/8
pattern is divided into $137^{\circ}\times 7+120^{\circ}\times 1$, as
indicated by the dashed arrow. Thus, the quantitative analysis reveals that
divergence angle on a stem is a secondary property, as it is very unlikely
that the exceptional angles of $160^{\circ}$ and $120^{\circ}$ are intrinsic
to the plant. Accordingly, the pattern of Fig. 2 is to be viewed as a result
of secondary processes.

The final step is to identify the cause of the irregular shift in divergence
angle. For brevity, let the central ($\times$), right ($\blacktriangle$) and
left ($\triangle$) trace of LPI $n$ be denoted as $n$C, $n$R and $n$L,
respectively. At the transition from 2/5 to 3/8, the irregular shift of 7C
(solid arrow) is accompanied by 5R on the right side (Fig. 4). Inspection of
Fig. 2 reveals that this collective shift is caused as a result of 0R
intervening between 3R and 5R. This is illustrated in Fig. 5, a schematic
excerpt from Fig. 2. For this reason, the trace 0R is indicated by a solid
arrow in Fig. 3 as the very cause of the phyllotaxis transition from 2/5 to
3/8. A horizontal line at five is drawn in Fig. 3 to indicate that 0R
interferes with 5R if only the former length exceeds $5-0=5$ internodes (cf.
Fig. 5). Indeed, the filled triangle at the solid arrow in Fig. 3 lies well
above the horizontal line at five internodes. Similarly, the trace causing the
transition from 3/8 to 5/13 is indicated by a dashed arrow in Fig. 3, where
the upper horizontal line at eight is drawn as a threshold length for the
transition.

Figure 5:  A schematic excerpt from Fig. 2 near the transition from 2/5 to
3/8. The right trace at LPI 5 is denoted as 5R, and 7C signifies the central
trace at LPI 7. According to Fig. 4, the traces 5R and 7C are shifted slightly
to the right as compared to the preceding (lower) traces. This figure shows
that the shift (solid arrow) is caused by a longer trace 0R intervening
between 3R and 5R, thereby the transition is initiated. If the length of 0R
were shorter than five internodes, 0R should have been aligned with 5R to keep
the 2/5 order, as the preceding 1R is with 6R.

##  4 DISCUSSION

Figure 6:  Distinct patterns of 2/5 and 3/8 orders on the mature stem
cylinders (bottom) result from similar arrangements of leaf primordia at the
shoot apical meristems (top). Each dot represents a leaf node, through which a
dotted line is drawn to demarcate internodes. Leaf traces (line segments) are
aligned in a 2/5 phyllotaxis pattern if their length is greater than three and
less than five internodes (left), whereas a 3/8 phyllotaxis results if it is
greater than five and less than eight internodes (right).

The trace length represented in internode units is an important geometric
factor as it imposes constraints on possible fractional patterns to be
realized (Okabe (2011, 2012b)). The geometric effect invoked in interpreting
the above results is schematically illustrated in Fig. 6. During growth of the
shoot apex, the precursor of the vascular system associated with a leaf
primordium develops to become the leaf trace (Romberger et al. (1993)). Both
lengths of the leaf trace and internode change continuously during growth. In
this view, the phyllotaxis order to be realized on the mature stem is
determined depending on the trace length represented in internode units. Top
two patterns in Fig. 6 represent the initial arrangements of the most common
case of 137.5∘ angular divergence, whereas their growth rates in the radial
direction are different, i.e., the two patterns are characterized with
different plastochron ratios (Richards (1951)). At this point, the difference
is not a qualitative but a quantitative one. Leaf traces (line segments) in
the left pattern traverse about four internodes, while those in the right
pattern traverse about six internodes. As the figure shows, a qualitative
difference is brought about through the quantitative difference in the trace
length: geometric constraints tend to achieve either a 2/5 or 3/8 phyllotaxis
on the mature stem depending on whether the trace length is shorter or longer
than five internodes. Experimentally, the difference would be judged on
whether leaf trace 0 is identified or not in the cross section at the level of
node 5. In a similar manner, the threshold value of the trace length for the
transition from 3/8 to 5/13 order is eight internodes, the denominator of the
lower order fraction. These are indicated by horizontal dashed lines in Fig.
3, as already remarked. In terms of the threshold length, the phyllotaxis
transition is interpreted consistently without resort to operations of
elaborate mechanisms like genetic programs; the phyllotaxis transitions are
caused because changing length of the leaf trace happens to cross the
threshold values of five and eight internodes. As illustrated in Fig. 6, this
view provides a consistent explanation of empirical observation that large
meristems result in arrangements of higher order fractions. Unequal
distribution of divergence angle noted in Fig. 4 is circumstantial evidence of
a unique intrinsic angle close to 137.5∘ and secondary distortions therefrom.

Fibonacci numbers 5 and 8 enter as the threshold values because leaf 0 appears
close to leaves 5 and 8 (Fig. 6). This in turn is a mathematical consequence
of the golden-mean angle at the apex. Mathematically, any number can be
approximated by a common fraction with any assigned degree of accuracy. The
greater the denominator, the better the approximation. A fast-converging
sequence of approximate fractions is uniquely determined for a given number
(Hardy and Wright (1979)). For the golden mean 0.3820, it is 1/2, 1/3, 2/5,
3/8, 5/13, etc., the main sequence of phyllotaxis. The denominators of these
approximate fractions comprise the index differences of nearby leaves, namely
the threshold values for the trace length. Thus, the geometric interpretation
predicts a correlation between the fraction index of the phyllotaxis order and
the trace length represented in internode units. Indeed, it has been remarked
as a general rule; the higher phyllotaxis order is associated with the longer
leaf trace (Girolami (1953); Esau (1965)). To conclude, the plant in its
maturity, as it were, achieves rational approximations to a divergence angle
at the apex in conformity with the leaf-trace length in internode units.

Mathematically, the golden mean is the worst “approximable” real number in the
sense that the sequence of the approximate fractions converges most badly. In
this connection, it has been commonly mentioned, despite objections, that the
golden-mean divergence is advantageous because leaves are distributed most
evenly to sunlight. When viewed in the context of this study, the golden-mean
divergence distributes the leaf traces most efficiently to coordinate the
vascular system. Although “phyllotaxis” is the arrangement of leaves on a stem
according to dictionaries, the geometric view on the leaf-trace organization,
the arrangement of leaves in a stem, may shed a new light on the long-standing
problem of phyllotaxis in vascular plants.

## References

  * Adam (2006) Adam, J., 2006. Mathematics in Nature: Modeling Patterns in the Natural World. Princeton University Press.
  * Adler et al. (1997) Adler, I., Barabé, D., Jean, R. V., 1997. A history of the study of phyllotaxis. Annals of Botany 80, 231–244.
  * Braun (1835) Braun, A., 1835. Dr. Carl Schimper’s Vorträge über die Möglichkeit eines wissenschaftlichen Verständnisses der Blattstellung, nebst Andeutung der hauptsächlichen Blattstellungsgesetze und Insbesondere der Neuentdeckten Gesetze der Aneinanderreihung von Cyclen Verschiedene. Flora 18, 145–191.
  * Church (1920) Church, A. H., 1920. On the interpretation of phenomena of phyllotaxis. Botanical memoirs. Hafner Pub. Co.
  * Coxeter (1961) Coxeter, H. S. M., 1961. Introduction to Geometry. Wiley, New York and London.
  * Esau (1965) Esau, K., 1965. Vascular differentiation in plants. New York: Holt, Rinehart and Winston.
  * Girolami (1953) Girolami, G., 1953. Relation between phyllotaxis and primary vascular organization in linum. American Journal of Botany 40, 618–625.
  * Hardy and Wright (1979) Hardy, G., Wright, E., 1979. An Introduction to the Theory of Numbers. Oxford Science Publications. Clarendon Press.
  * Henfrey (1870) Henfrey, A., 1870. An elementary course of botany: structural, physiological and systematic. J. Van Voorst.
  * Kuhlemeier (2007) Kuhlemeier, C., 2007. Phyllotaxis. TRENDS in Plant Science 12, 143–150.
  * Larson (1980) Larson, P. R., 1980. Interrelations between phyllotaxis, leaf development and the primary-secondary vascular transition in Populus deltoides. Annals of Botany 46, 757–769.
  * Okabe (2011) Okabe, T., 2011. Physical phenomenology of phyllotaxis. Journal of Theoretical Biology 280, 63–75.
  * Okabe (2012a) Okabe, T., 2012a. Systematic variations in divergence angle. Journal of Theoretical Biology 313, 20–41.
  * Okabe (2012b) Okabe, T., 2012b. Vascular phyllotaxis transition and an evolutionary mechanism of phyllotaxis. http://arxiv.org/abs/1207.2838.
  * Richards (1951) Richards, F. J., 1951. Phyllotaxis: its quantitative expression and relation to growth in the apex. Philosophical Transactions of the Royal Society of London. Series B 225, 509–564.
  * Romberger et al. (1993) Romberger, J., Hejnowicz, Z., Hill, J., 1993. Plant structure: function and development : a treatise on anatomy and vegetative development, with special reference to woody plants. Springer-Verlag.


# SSDNet: State Space Decomposition Neural Network for Time Series Forecasting

Yang Lin  School of Computer Science
The University of Sydney
Sydney, Australia
<EMAIL_ADDRESS>Irena Koprinska  School of Computer Science
The University of Sydney
Sydney, Australia
<EMAIL_ADDRESS>Mashud Rana  Data61
CSIRO
Sydney, Australia
<EMAIL_ADDRESS>
###### Abstract

In this paper, we present SSDNet, a novel deep learning approach for time
series forecasting. SSDNet combines the Transformer architecture with state
space models to provide probabilistic and interpretable forecasts, including
trend and seasonality components and previous time steps important for the
prediction. The Transformer architecture is used to learn the temporal
patterns and estimate the parameters of the state space model directly and
efficiently, without the need for Kalman filters. We comprehensively evaluate
the performance of SSDNet on five data sets, showing that SSDNet is an
effective method in terms of accuracy and speed, outperforming state-of-the-
art deep learning and statistical methods, and able to provide meaningful
trend and seasonality components.

###### Index Terms:

time series forecasting, time series decomposition, state space model, deep
learning

##  I Introduction

Time series forecasting is an important task in many practical applications,
e.g. predicting electricity demand, stock prices, immune response and disease
progression over time.

Statistical methods such as ARIMA and exponential smoothing are well
established for time series forecasting and are a part of the more general
framework of State Space Models (SSMs) [1]. SSMs are considered interpretable
models as components such as trend and seasonality can be extracted and used
for decision making and explanation. However, they are not able to infer
shared patterns from a set of related time series as each time series is
fitted independently [2, 3].

Deep learning methods have been recently investigated as a promising
alternative, due to their ability to learn from raw data with minimum domain
knowledge and to extract complex patterns, including shared patterns across
related time series. Prominent examples include DeepAR [4], a probabilistic
forecasting model based on Long Short Term Memory (LSTM) neural networks,
DeepSSM [5] which combines SSM with LSTM, Deep Factor models [6] which combine
LSTM with a local probabilistic model, and LogSparse Transformer [3] and
Informer [7], which are modifications of the Transformer architecture [8] for
time series forecasting. However, deep learning models are difficult to
interpret; generating interpretable forecasts is crucial for adopting these
systems in practical applications.

Approaches for producing interpretable forecasts using deep learning models
have been recently proposed. For example, Guo et al. [9] and Li et al. [10]
explored the structure of LSTM and DeepSSM attempting to learn the variable
importance for prediction. Oreshkin et al. [2] developed the deep learning
architecture N-BEATS. They showed that its interpretable configuration
(N-BEATS-I) is able to successfully learn and output the trend and seasonality
components of the forecast, which is useful for forecasting practitioners.
However, N-BEATS is designed for univariate time series, it provides point
forecasts and has limited ability to model trends with varying slope because
it uses a deterministic trend model. In this paper, we present a new
forecasting approach to address these issues.

The main contributions of this work are as follows:

  1. 1.

We present a new forecasting approach, called State Space Decomposition Neural
Network (SSDNet), which combines the Transformer deep learning architecture
with a SSM. SSDNet combines the advantages of deep learning (learning from raw
data without intensive feature engineering and ability to infer shared
patterns from related time series) with the interpretability of SSM models. It
employs the Transformer architecture to learn the temporal pattern and
directly estimate the parameters of SSM. To facilitate interpretability, we
used a fixed form SSM to provide trend and seasonality components and the
attention mechanism of the Transformer to identify which parts of the past
history are most important for the prediction.

  2. 2.

We evaluate the performance of SSDNet on five time series forecasting tasks.
The results show that SSDNet achieved higher accuracy than the state-of-the-
art deep learning models DeepAR, DeepSSM, LogSparse Transformer, Informer and
N-BEATS and the statistical models SARIMAX and Prophet. SSDNet was also able
to provide interpretable results by showing nonlinear trend and seasonality
components. To evaluate the effectiveness of the fixed form SSM part, we
conduct an ablation study by replacing the Transformer with LSTM, and this
architecture also showed competitive results.

##  II Problem Formulation

We consider three tasks: 1) solar power forecasting, 2) electricity demand
forecasting and 3) exchange rate forecasting. Accurate forecasting of the
generated solar power and the electricity demand is needed for optimal
scheduling of generators and large-scale integration of solar into the
electricity grid. Over-prediction may lead to wasting energy, under-prediction
may result in blackouts. Exchange rate forecasting is used to derive future
monetary value, earn profits and avoid risks in international business
environments.

###  II-A Data Sets

Figure 1: Partial autocorrelation graphs for (a) Sanyo and Hanergy, (b) Solar,
(c) Electricity and (d) Exchange data sets

Figure 2: Time series additive decomposition for (a) Sanyo, (b) Hanergy, (c) Solar, (d) Electricity data sets and (e) one series from Exchange data set |  start date |  end date |  granularity |  $L_{d}$ |  $N$ |  $n_{T}$ |  $n_{C}$ |  $T_{l}$ |  $T_{h}$
---|---|---|---|---|---|---|---|---|---
Sanyo |  01/01/2011 |  31/12/2016 |  30 minutes |  20 |  1 |  4 |  3 |  20 |  20
Hanergy |  01/01/2011 |  31/12/2017 |  30 minutes |  20 |  1 |  4 |  3 |  20 |  20
Solar |  01/01/2006 |  31/08/2006 |  1 hour |  24 |  137 |  0 |  3 |  24 |  24
Electricity |  01/01/2011 |  07/09/2014 |  1 hour |  24 |  370 |  0 |  4 |  168 |  24
Exchange |  01/01/1990 |  31/12/2016 |  1 day |  1 |  8 |  0 |  3 |  30 |  20
TABLE I: Dataset statistics. $L_{d}$ \- number of steps per day, $N$ \- number
of series, $n_{T}$ \- number of time-based features, $n_{C}$ \- number of
calendar features, $T_{l}$ \- length of input series, $T_{h}$ \- length of
forecasting horizon.

We use five data sets: Sanyo [11], Hanergy [12], Solar [13], Electricity [14]
and Exchange [15]; their statistics are summarised in Table I.

Sanyo and Hanergy contain solar power generation data from two PV plants in
Australia - from 01/01/2011 to 31/12/2016 (6 years) for Hanergy and 01/01/2011
to 31/12/2017 (7 years) for Sanyo. Only the data between 7am and 5pm is
considered and was aggregated at half-hourly intervals. For both datasets,
weather and weather forecast data was also collected (see [16] for more
details) and used as covariates. Solar contains solar power data from 137 PV
plants in Alabama, USA, from 01/01/2006 to 31/08/2006. Electricity contains
electricity consumption data for 370 households from 01/01/2011 to 07/09/2014.
The Solar and Electricity data is aggregated into 1-hour intervals. Exchange
contains daily exchange rate data for 8 countries from 01/01/1990 to
31/12/2016. Following [3, 16, 17], calendar features are added according to
the granularity of the datasets: Sanyo and Hanergy use month, hour-of-the-day
and minute-of-the-hour, Solar uses month, hour-of-the-day and age, Electricity
uses month, day-of-the-week, hour-of-the-day and age and Exchange uses month,
day-of-the-week and age.

Fig. 1 shows the partial autocorrelation plots for each data set. Fig. 1 (a)
shows the results for both Sanyo and Hanergy since they have a similar partial
autocorrelation pattern; both datasets include only 1 time series. Fig. 1 (b),
(c) and (d) show the results for Solar, Electricity and Exchange datasets
respectively. All plots are for 8 time series for consistency since 8 is the
common number of time series in the three datasets, see Table I.

Fig. 2 shows the additive decomposition results for a period of 7 days for
Sanyo, Hanergy, Solar and Electricity. For Exchange we just show one time
series as this dataset doesn’t have a clear seasonality and additive
decomposition methods cannot be applied.

We can see that all time series except Exchange exhibit seasonality and
repetitive patterns with high partial autocorrelation, and sometimes there are
also significant short-term trend variations. For the Exchange set, there is
no clear repetitive pattern and its partial autocorrelation at lag 1 is the
only significant one.

All data was normalized to have zero mean and unit variance.

###  II-B Problem Statement

Given is: 1) a set of $N$ univariate time series (solar, electricity or
exchange series) $\\{\mathbf{Y}_{i,1:T_{l}}\\}^{N}_{i=1}$, where
$\mathbf{Y}_{i,1:T_{l}}\coloneqq[{y}_{i,1},{y}_{i,2},...,{y}_{i,T_{l}}]$,
$T_{l}$ is the input sequence length, and ${y}_{i,t}\in\Re$ is the value of
the $i$th time series (generated PV solar power or consumed electricity) at
time $t$; 2) a set of associated time-based multi-dimensional covariate
vectors $\\{\mathbf{X}_{i,1:T_{l}+T_{h}}\\}_{i=1}^{N}$, where $T_{h}$ denotes
the length of the forecasting horizon. The covariates for the Sanyo and
Hanergy datasets include: weather data
$\\{\mathbf{W1}_{i,1:T_{l}}\\}_{i=1}^{N}$, weather forecasts
$\\{\mathbf{WF}_{i,T_{l}+1:T_{l}+T_{h}}\\}_{i=1}^{N}$ and calendar features
$\\{\mathbf{Z}_{i,1:T_{l}+T_{h}}\\}_{i=1}^{N}$, while the covariates for the
Solar, Electricity and Exchange datasets include only calendar features. Our
goal is to predict the future values of the time series
$\\{\mathbf{Y}_{i,T_{l}+1:T_{l}+T_{h}}\\}^{N}_{i=1}$, i.e. the PV power or
electricity usage for the next $T_{h}$ time steps after $T_{l}$.

Specifically, SSDNet produces the probability distribution of the future
values, given the past history:

$\begin{split}p\left(\mathbf{Y}_{i,T_{l}+1:T_{l}+T_{h}}\mid\mathbf{Y}_{i,1:T_{l}},\mathbf{X}_{i,1:T_{l}+T_{h}};\Phi\right)\\\
=\prod_{t=T_{l}+1}^{T_{l}+T_{h}}p\left({y}_{i,t}\mid\mathbf{Y}_{i,1:t-1},\mathbf{X}_{i,1:t};\Phi\right),\end{split}$
(1)

where $\Phi$ denotes the parameters of SSDNet, and the input of SSDNet at step
$t$ is the concatenation of ${y}_{i,t-1}$ and ${x}_{i,t}$. The models are
applicable to all time series, so the subscript $i$ will be omitted in the
rest of the paper for simplicity.

##  III Related Work

DeepSSM [5] combines SSMs with Recurrent Neural Networks (RNNs). A RNN is used
to generate the parameters of a linear Gaussian SSM for time series
forecasting and Kalman filter is used to derive posterior knowledge. DeepSSM
is considered interpretable because it can recover the SSM parameters and they
can be inspected. The Deep Factor model DF-LDS [6] is a variation of DeepSSM
that can handle non-Gaussian observations and it also uses Kalman filters.

Although Kalman filters have been successfully used in linear dynamic systems
for decades, the need for covariance matrix inversions is computationally
expensive and limits their applicability for processing large data sets [18,
19]. In contrast, Eleftheriadis et al. [20] proposed Gaussian Process State
Space Model (GPSSM) merging Gaussian process with a bi-directional RNN to
approximate the posterior of a non-linear system. However, GPSSM has high
complexity of inference ($O(T)$ per time step) which is a disadvantage for
forecasting multi-horizon large data sets.

N-BEATS [2] is a novel deep neural architecture using backward and forward
residual links and stacks of fully connected layers for univariate point
forecasting. Each N-BEATS layer contains a number of blocks which produce
partial forecasts; these forecasts are aggregated at stack level and then
network level in hierarchical fashion. The final forecast is the sum of the
partial forecasts of all blocks. N-BEATS can produce interpretable forecasting
results without expert knowledge by employing a fixed form of block function
to generate prediction with decomposition components. The interpretable
version (N-BEATS-I) uses the deterministic polynomial trend model and assumes
a stationary trend with constant slope. However, the local trend of time
series over a short period (one or a few days) could vary significantly and
inconsistently - e.g. the solar power generation or electricity consumption
could decrease suddenly due to temporary cloud coverage or irregular
activities of users (see Fig. 2). Hence, a stochastic model which assumes a
nonstationary trend with varying slope would be more suitable for such cases
[21] and this is what we use in SSDNet.

In summary, compared to DeepSSM, the proposed SSDNet: i) employs the
Transformer architecture to learn temporal patterns, and estimates the
parameters of SSM and the probability term directly instead of using Kalman
filters; ii) uses the SSM with fixed and non-trainable transition matrix at
the decoder to generate decomposition results within the forecasting horizon,
while DeepSSM uses a more flexible SSM along all time steps; iii) combines
Mean Absolute Error (MAE) and Negative Log-Likelihood (NLL) together as loss
function to achieve accurate point and probabilistic forecasting results
simultaneously instead of using NLL only; iv) shows better interpretability by
using a stochastic decomposition technique in the form of SSM with an
innovation term to model nonlinear trend and seasonality components and
attention mappings to show the importance of historical time steps for
forecasting. Compared to N-BEATS, SSDNet can model trend with nonlinear slope,
provide probabilistic forecasts and also better utilizes the covariates.

##  IV SSDNet

The aim of SSDNet is to model temporal patterns effectively in order to
provide accurate probabilistic and interpretable forecasts.

###  IV-A Network Architecture

Fig. 3 illustrates the architecture of SSDNet. SSDNet is based on the encoder-
decoder framework and uses two key components: SSM and Transformer. SSM (in
its fixed form) provides interpretable forecasting results, and the
Transformer learns temporal patterns and estimates the SSM parameters from its
decoder.

The feedforward process of SSDNet includes two steps. First, the Transformer
of SSDNet processes the historical time steps and generates the latent
components to estimate the parameters of SSM and the variance of the
forecasted distribution. Second, SSM takes the state vector from the previous
time step and uses it to predict the mean of the distribution. Unlike
classical SSMs which decompose time series and estimate parameters via Kalman
filter, we learn the SSM parameters using the Transformer architecture via
attention mechanism and then combine the decomposition components (trend,
seasonality and probability terms) in a form of SSM to produce the forecasting
results.

Figure 3: Structure of SSDNet

As shown in Fig. 3, the Transformer extracts latent components $o_{t}$ from
historical time series $y_{1:T_{l}},x_{1:T_{t}}$ in Eq. (2), and the latent
components are further used to estimate the parameters of SSM. Note that the
Transformer could be replaced by any sequence model to extract $o_{t}$ and
learn the parameters of SSM. In Section VB, we evaluate the use of LSTM
instead of Transformer in SSDNet; this modified architecture is called SSDNet-
LSTM.

$o_{t}=f(y_{1:T_{l}},x_{1:T_{t}})$ (2)

Then, we employ additive time series decomposition model in the form of SSM
and generate the prediction $\hat{y_{t}}$ at step $t$ by summing up the trend
component $T_{t}$, seasonality component $S_{t}$ and probability component
$I_{t}$. We decompose the times series using SSM and generate the forecasting
result as:

$\hat{y_{t}}=z_{t}^{T}\alpha_{t}+I_{t},\quad t=1,\ldots,T_{h}$ (3)
$\alpha_{t+1}=\Gamma_{t}\alpha_{t}+c_{t}$ (4)
$I_{t}\sim\mathcal{N}(0,\sigma_{I_{t}}^{2})$ (5)

where the latent state vector $\alpha_{t}\in\Re^{s\times 1}$ contains trend
$Tr_{t}$ and seasonality $S_{t}$, $s$ is the number of seasonality and $I_{t}$
is sampled from the Gaussian distribution with zero mean and
$\sigma_{I_{t}}^{2}$ variance.

Different from the traditional SSMs, we remove the random noise term from the
state transition in Eq. (4) to prevent the state noise propagation and
estimate the probabilistic term $I_{t}$ directly as shown in Eq. (3). The SSM
part of SSDNet does not process historical series directly but instead
utilizes the latent component generated by the Transformer and encodes the
information of the time steps before the forecasting horizon in the initial
state vector $\alpha_{T_{l}}$. We also introduce the innovation term
$c_{t}\in\Re^{s\times 1}$ in the state transition Eq. (4) to allow SSDNet to
learn stochastic trends with fluctuations in time series. This is necessary as
the simple summation of linear trend and seasonality components cannot model
complex series effectively and would introduce large residuals (see Fig. 2).
While the traditional SSMs need to use Kalman filter to update the posterior
and derive the covariance matrix to optimize noise terms, SSDNet uses the
Transformer architecture resulting in less matrix computation and memory
usage.

The innovation term $c_{t}$ and the variance $\sigma_{I_{t}}^{2}$ are learnt
from the latent factor $o_{t}$ directly, as illustrated in Eq. (6) and (7).
The use of the Softplus function ensures that SSDNet always generates positive
variance, and the HardSigmoid function [22] is used for a speed-up:

$\begin{split}\sigma_{I_{t}}^{2}=g_{s}(o_{t})=&\text{Softplus}(\text{Linear}(o_{t}))\\\
=&\log(1+\exp(\text{Linear}(o_{t})))\end{split}$ (6)
$\begin{split}c_{t}=g_{c}(o_{t})=&\text{HardSigmoid}(\text{Linear}(o_{t}))-0.5\\\
=&\left\\{\begin{array}[]{ll}-0.5&\text{ if }\mathrm{x}\leq-3\\\ 0.5&\text{ if
}\mathrm{x}\geq+3\\\ \mathrm{\text{Linear}(o_{t})}/6&\text{ otherwise
}\end{array}\right.\end{split}$ (7)

The parameters state vector $\alpha_{t}$, state transition matrix $\Gamma_{t}$
and result transition matrix $z_{t}$ that determine the structure of SSM in
Eq. (3) and (4) are shown below:

$\alpha_{t}=\left(\begin{array}[]{c}Tr_{t}\\\ S_{1:s-1,t}\\\
\end{array}\right),{z}_{t}=\left(\begin{array}[]{l}1\\\ 1\\\
0_{s-2}\end{array}\right),$ (8)
$\displaystyle\Gamma_{t}=\left(\begin{array}[]{cccc}1&0_{s-2}^{\prime}&0\\\
0&-1_{s-2}^{\prime}&-1\\\ 0_{s-2}&I_{s-2}&0_{s-2}\end{array}\right)$

Note that $\Gamma_{t}$ and $z_{t}$ are non-trainable and fixed for all time
steps. In this work, we assume that the trend follows a random walk process,
and other processes such as moving average could be modelled by adjusting the
first row of the transition matrix $\Gamma_{t}$. The initial trend $Tr_{0}$
and seasonality $S_{-s+1},...,S_{-1},S_{0}$ values of the state vector
$\alpha_{0}$ are unknown and need to be generated based on the latent
components, they are learned from historical steps:

$\begin{split}\alpha_{0}&=g_{c}(o_{T_{l+1}})=\text{HardSigmoid}(\text{Linear}(o_{T_{l+1}}))-0.5\end{split}$
(9)

Alternatively, the forecasting result (3) could be written as the aggregation
of the predicted trend, seasonality and variance to provide probabilistic
result:

$\hat{y_{t}}\sim\mathcal{N}(Tr_{t}+S_{t},\,\sigma_{I_{t}}^{2})\,$ (10)

where we consider the real-world data follows a Gaussian distribution (other
distribution function could be used in Eq. (5) for higher flexibility), and
predictions are sampled from the distribution. The $\rho$-quantile output
could be generated via the inverse cumulative probability distribution:
$\hat{y_{t}}=F_{t}^{-1}(\rho)$.

At step $t$, Eq. (4) could be written in an additive form, where the trend and
seasonality evolve by adding the innovation term at each time step, as shown
in Eq. (11) and (12). The trend and seasonality models are based on the
classic random walk and dummy seasonal, and the change of trend and
seasonality are dominated by the innovation variable. The innovation variable
for trend can be different at every time step and allow SSDNet to model the
nonstationary trend with varying slope. For the extreme case when all
innovation terms are zero, the trend would be constant and the seasonality
over all periods will be the same.

$Tr_{t}=Tr_{t-1}+c_{t}[1]$ (11) $S_{t}=-\sum_{j=1}^{s-1}S_{t-j}+c_{t}[2]$ (12)

######  Theorem 1.

$\forall t>0$ and $t\in\mathbb{Z}$, the range of trend $Tr_{t}$ and
seasonality $S_{t}$ of SSDNet are bounded: $Tr_{t}\subset[-(t+1)\times
0.5,(t+1)\times 0.5]$, $S_{t}\subset[-(s-1+t)\times 0.5,(s-1+t)\times 0.5]$.

###### Proof.

Assume that innovation term $c_{t}$, trend and seasonality have maximum
values: 1) $c_{t}$ outputs its maximum value of 0.5 at all time steps, 2) the
initial trend has its maximum value $Tr_{0}=0.5$ and 3) the seasonality has
its maximum value of $S_{0}=-\sum_{j=1}^{s-1}-0.5=0.5\times(s-1)$.

Following Eq. (11) and (12), the trend and seasonality increase using the
innovation terms with maximum values of $c_{t}[1]=0.5$ and $c_{t}[2]=0.5$
respectively at each time step.

Considering the initial values and increments together, the maximum values of
the trend and seasonality are $(t+1)\times 0.5$ and $(s-1+t)\times 0.5$ at
time step $t$, correspondingly. Since $g_{c}(.)$ is a symmetric activation
function, its lower bound is the negative value of its upper bound, hence the
minimum values of trend and seasonality are the negative of their maximum
values. ∎

Theorem 1 implies that the upper and lower bound of the SSDNet trend and
seasonality components expand in the time domain and SSDNet can model time
series where the variation in the seasonal pattern increases over time. The
variation of the trend component is significantly smaller than that of the
seasonality and helps SSDNet to model a slow-moving trend and significantly
varying seasonality. The constraints of both components help SSDNet to learn
an inherent structure and allow $\alpha_{t}$ to generate meaningful waveform
(observable in our case studies) and enable the interpretability of
forecasting results.

In our case since the data is normalized to have zero mean and unit variance,
when visualizing trend and seasonality components, the bias (data mean) is
added to the trend and both the trend and seasonality are scaled (data
variance) back.

###  IV-B Loss Function

SSDNet needs to provide both accurate point and probabilistic forecasts.
Accordingly, we developed the loss function shown in Eq. (13), which considers
both the point and probabilistic forecasts by combining the MAE and NLL using
the regularization parameter $a$: higher $a$ allows the model to pay more
attention to the probabilistic forecast. The point forecast is $\hat{y}_{t}$
and it consists of trend and seasonality. The parameters of SSDNet are
optimized by minimizing this function.

$\begin{split}&L(\hat{y}_{T_{l}+1:T},\sigma_{I_{T_{l}+1:T}}^{2},y_{T_{l}+1:T},a)\\\
=&~{}a\times\operatorname{NLL}(\hat{y}_{T_{l}+1:T},\sigma_{I_{T_{l}+1:T}}^{2},y_{T_{l}+1:T})\\\
&+\operatorname{MAE}(\hat{y}_{T_{l}+1:T},y_{T_{l}+1:T})\\\
=&-\frac{a}{2T_{h}}\times\Big{(}T_{h}\log(2\pi)+\sum_{t=T_{l}+1}^{T}\log\left|\sigma_{I_{t}}^{2}\right|\\\
&+\sum_{t=T_{l}+1}^{T}(y_{t}-\hat{y}_{t})^{2}\sigma_{I_{t}}^{-2}\Big{)}+\frac{1}{T_{h}}\sum_{t=T_{l}+1}^{T}\left|y_{t}-\hat{y}_{t}\right|\end{split}$
(13)

##  V Experiments

We compare the performance of SSDNet with nine other models: six state-of-the-
art autoregressive deep learning models (DeepAR, DeepSSM, N-BEATS-G,
N-BEATS-I, LogSparse Transformer and Informer), a SSM (SARIMAX), an
interpretable regression model (Prophet) and a persistence model:

  * •

Persistence is a typical baseline in forecasting and considers the time series
of the previous day as the prediction for the next day. For the Exchange
dataset we use the last 20 steps of the input sequence.

  * •

SARIMAX [1] is an extension of the ARIMA and can handle seasonality with
exogenous factors.

  * •

Prophet [23] is an interpretable regression model predicting trend,
seasonality and holiday components.

  * •

DeepAR [4] is a widely used sequence-to-sequence probabilistic forecasting
model.

  * •

DeepSSM [5] fuses SSM with RNNs to incorporate structural assumptions and
learn complex patterns from the time series. It is the state-of-the-art deep
forecasting model that employs SSM.

  * •

N-BEATS [2] is based on backward and forward residual links and stacks of
fully connected layers. N-BEATS-G provides generic forecasting results, while
N-BEATS-I provides interpretable results by decomposing the time series into
trend and seasonality. We introduced covariates to N-BEATS at the input of
each block to facilitate multivariate series forecasting.

  * •

LogSparse Transformer [3] is a recently proposed variation of the Transformer
architecture for time series forecasting with convolutional attention and
sparse attention; it is denoted as ”LogTrans” in Table III.

  * •

Informer [7] is a Transformer-based forecasting model based on the ProbSparse
self-attention and self-attention distilling. We modified the Informer to
provide probabilistic results by using Eq. (6) as the output layer.

|  $\lambda$ |  $\delta$ |  $d_{hid}$ |  $n_{l}$ |  $d_{k\&v}$ |  $n_{h}$
---|---|---|---|---|---|---
Sanyo |  0.005 |  0 |  12 |  2 |  6 |  2
Hanergy |  0.005 |  0 |  16 |  3 |  6 |  3
Solar |  0.005 |  0.1 |  16 |  3 |  6 |  3
Electricity |  0.001 |  0.1 |  24 |  3 |  8 |  2
Exchange |  0.005 |  0 |  12 |  2 |  4 |  3
TABLE II: Hyperparameters for SSDNet |  Sanyo |  Hanery |  Solar |  Electricity |  Exchange
---|---|---|---|---|---
Persistence |  0.154/- |  0.242/- |  0.256/- |  0.069/- |  0.016/-
SARIMAX |  0.124/0.096 |  0.145/0.098 |  0.256/0.192 |  0.196/0.079 |  0.010/0.006
Prophet |  0.104/0.054 |  0.152/0.079 |  0.268/0.169 |  0.112/0.055 |  0.017/0.013
DeepAR |  0.070/0.031 |  0.092/0.045 |  0.222⋄/0.093⋄ |  0.075⋄/0.040⋄ |  0.014/0.009
DeepSSM |  0.042/0.023 |  0.070/0.053 |  0.223/0.181 |  0.083⋄/0.056⋄ |  0.014/0.012
LogTrans |  0.067/0.036 |  0.088/0.047 |  0.210⋄/0.082⋄ |  0.059⋄/0.034⋄ |  0.017/0.008
Informer |  0.046/0.022 |  0.084/0.046 |  0.215/0.115 |  0.068/0.033 |  0.014/0.009
N-BEATS-I |  0.091/ |  0.154/- |  0.215/- |  0.102/- |  0.014/-
N-BEATS-G |  0.077/- |  0.132/- |  0.212/- |  0.061/- |  0.018/-
SSDNet |  0.040/0.020 |  0.059/0.032 |  0.209/0.074 |  0.068/0.033 |  0.013/0.006
TABLE III: $\rho$0.5/$\rho$0.9-loss of data sets with various granularities.
$\diamond$ denotes results from [3].

All models were implemented with PyTorch 1.6 on Tesla T4 16GB GPU under Linux
environment. The deep learning models were optimized by mini-batch gradient
descent with the Adam optimizer and a maximum number of epochs 200. We used
Bayesian optimization for hyperparameter search for all deep learning models
with maximum number of 20 iterations. The hyperparameters with a minimum loss
on the validation set were selected. The models used for comparison were tuned
based on the authors’ recommendations. The probabilistic forecasting models
use the NLL loss, and the point forecasting model (N-BEATS) uses the mean
squared loss.

Following the experimental setup in [3], [16] and [17], we used the following
training, validation and test split: for Sanyo and Hanergy - the data from the
last year as test set, the second last year as validation set for early
stopping and the remaining data (5 years for Sanyo and 4 years for Hanergy) as
training set; for Solar and Electricity - the last week data as test set (from
25/08/2006 for Solar and 01/09/2014 for Electricity), the week before as
validation set and the remaining data as training set. For the Exchange
dataset, we used the 480 working days from 02/02/2015 is test set, the 480
working days before as validation set and the remaining data as training set.
For all data sets, the data preceding the validation set is split in the same
way into three subsets and the corresponding validation set is used to select
the best hyperparameters.

For the Transformer-based models, we used learnable position and ID (for
Solar, Electricity and Exchange sets) embedding. For SSDNet, the seasonality
$s$ was set to 20 for Sanyo, Hanergy and Exchanges and 24 for Solar and
Electricity, the loss function regularization parameter $a$ was fixed as 0.5,
the learning rate $\lambda$ was fixed, the dropout rate $\delta$ was chosen
from {0, 0.1, 0.2}, the hidden layer dimension size $d_{hid}$ and number of
layers $n_{l}$ were chosen from {8, 12, 16, 24, 32} and {2, 3, 4}, the query
and value’s dimension size $d_{k\&v}$ and number of heads $n_{h}$ were chosen
from {4, 6, 8, 12} and {2, 3, 4}. The selected best hyperparameters for SSDNet
are listed in Table II and used for the evaluation of the test set.

Following [4, 5, 3, 16, 17], we report the standard $\rho$0.5 and
$\rho$0.9-quantile losses. Note that $\rho$0.5 is equivalent to the Mean
Absolute Percentage Error (MAPE) [24]. Given the ground truth $y$ and
$\rho$-quantile of the predicted distribution $\hat{y}$, the $\rho$-quantile
loss is given by:

$\begin{split}\mathrm{QL}_{\rho}(y,\hat{y})&=\frac{2\times\sum_{t}P_{\rho}\left(y_{t},\hat{y}_{t}\right)}{\sum_{t}\left|y_{t}\right|},\\\
\quad P_{\rho}(y,\hat{y})&=\left\\{\begin{array}[]{ll}\rho(y-\hat{y})&\text{
if }y>\hat{y}\\\ (1-\rho)(\hat{y}-y)&\text{ otherwise
}\end{array}\right.\end{split}$ (14)

###  V-A Accuracy Analysis

The $\rho$0.5 and $\rho$0.9 losses are shown in Table III. As N-BEATS and
Persistence do not produce probabilistic forecasts, only the $\rho$0.5-loss is
reported for them. We can see that overall SSDNet is the most accurate model -
it outperforms all the other methods on all data sets except Electricity and
Exchange. For the Electricity data set, SSDNet has the best $\rho$0.9-loss and
second-best $\rho$0.5-loss among the probabilistic forecasting models and is
also the best among the interpretable models (N-BEATS-I, DeepSSM and Prophet).
For the Exchange data set, SSDNet is the second-best model after SARIMAX.

SSDNet performs better than DeepSSM on all datasets which indicates that the
proposed new architecture, loss function and the use of fixed SSM were
beneficial.

The two N-BEATS models perform well on Solar and Electricity (univariate data
sets) but much worse on Sanyo and Hanergy (multivariate, including the weather
features); a possible reason is that N-BEATS was designed for univariate
forecasting and our modification to handle multivariate datasets was not
sufficiently effective in this case.

We further investigated the performance on the Exchange dataset, for which the
deep neural network models, including the proposed SSDNet, performed worse
than SARIMAX. This finding is consistent with [25], and the reason is that
Exchange is the only data set without repetitive patterns and a seasonal
component. Fig. 5 and Fig. 6 show the predictions of SSDNet and SARIMAX for 2
test samples from the Exchange data set. We can see that SSDNet tends to
generate predictions with a smooth trend and some variations, while SARIMAX
tends to generate predictions that do not vary much with time. Hence, the
success of SARIMAX on the Exchange data set can be explained with the
characteristics of this data set - it is more stable and does not have strong
repetitive patterns and seasonality, as shown in Section IIA and Fig. 1.

Another interesting finding for Exchange is that when we simply use the last
time step value $\mathbf{Y}_{T_{l}}$ of the input series as the predictions
$\mathbf{Y}_{T_{l}+1:T_{l}+T_{h}}$, we obtain $\rho$0.5-loss of 0.010 which is
as good as SARIMAX and outperforms all other models. However, SSDNet and the
other deep learning models demonstrate significantly higher accuracy than
SARIMAX on the other four data sets which include repetitive patterns.

|  Sanyo |  Hanery |  Solar |  Electricity |  Exchange
---|---|---|---|---|---
DeepAR |  0.070/0.031 |  0.092/0.045 |  0.222⋄/0.093⋄ |  0.075⋄/0.040⋄ |  0.014/0.009
DeepSSM |  0.042/0.023 |  0.070/0.053 |  0.223/0.181 |  0.083⋄/0.056⋄ |  0.014/0.012
SSDNet |  0.040/0.020 |  0.059/0.032 |  0.209/0.074 |  0.068/0.033 |  0.013/0.006
SSDNet-LSTM |  0.040/0.020 |  0.066/0.037 |  0.205/0.075 |  0.071/0.037 |  0.012/0.006
TABLE IV: Ablation study - $\rho$0.5/$\rho$0.9-loss. $\diamond$ denotes results from [3]. |  Sanyo |  Hanery |  Solar |  Exchange
---|---|---|---|---
DeepSSM |  25163$\pm$201 |  9750$\pm$135 |  224343$\pm$692 |  71188$\pm$343
SSDNet |  6006$\pm$91 |  3081$\pm$91 |  61720$\pm$206 |  13817$\pm$62
SSDNet-LSTM |  5153$\pm$68 |  2079$\pm$29 |  56989$\pm$127 |  7890$\pm$99
TABLE V: Training time per epoch (milliseconds) - comparison of DeepSSM and SSDNets on Sanyo, Hanergy, Solar and Exchange data sets. |  Sanyo |  Hanery |  Solar |  Exchange
---|---|---|---|---
DeepSSM |  242$\pm$54 |  233$\pm$6 |  278$\pm$9 |  250$\pm$7
SSDNet |  411$\pm$177 |  438$\pm$132 |  573$\pm$215 |  389$\pm$125
SSDNet-LSTM |  254$\pm$4 |  257$\pm$8 |  299$\pm$12 |  270$\pm$5
TABLE VI: Testing time (milliseconds) - comparison of DeepSSM and SSDNets on
Sanyo, Hanergy, Solar and Exchange data sets.

###  V-B Ablation Analysis

To evaluate the effectiveness of the SSM part of SSDNet, we conducted
additional analysis - an ablation study. Table IV shows the performance of
SSDNet-LSTM, which is SSDNet with the Transformer part replaced by LSTM. Thus,
we keep the SSM part for decomposition and prediction the same but replace the
Transformer with LSTM for SSM parameters estimation.

SSDNet-LSTM is compared with DeepAR [4] and DeepSSM [5]. DeepAR uses LSTM to
generate predictions directly and DeepSSM uses LSTM and Kalman filter to
estimate the parameters of SSM, while SSDNet-LSTM can be considered as a
combination of DeepAR and our proposed SSM or as a variation of DeepSSM by
replacing the classic SSM with our proposed SSM without a Kalman filter.

Table IV shows that SSDNet-LSTM outperforms both DeepAR and DeepSSM,
indicating the effectiveness of the SSM part for decomposition and foresting.
SSDNet-LSTM is also competitive with SSDNet which shows the robustness and
flexibility of our fixed form SSM - its parameters can be estimated by
different architectures - the Transformer, LSTM and other neural network
models.

###  V-C Speed Analysis

We evaluate the training and testing time of the two SSDNet models (SSDNet and
SSDNet-LSTM) and compare them with DeepSSM which also employs SSM. All models
are trained on the same computer configuration; we report the average elapsed
time and the standard deviation of 10 runs. Table V and VI show the training
time per epoch and testing time in milliseconds for the Sanyo, Hanergy, Solar
and Exchange data sets. No speed test was conducted for the Electricity data
set as we report the accuracy results from their papers [5, 3] and the values
of the hyperparameters used are not known.

Comparing the training speed, we can see that both SSDNet and SSDNet-LSTM,
which do not use Kalman filter, are significantly faster than DeepSSM - the
average training speed per epoch of DeepSSM is 390.48$\%$ and 458.24$\%$ times
slower than SSDNet and SSDNet-LSTM respectively. In terms of testing speed,
DeepSSM is the fastest, very closely followed by SSDNet-LSTM (7.68$\%$ slower)
and then SSDNet. This is because the Kalman filter in DeerSSM is not required
for inference and the architecture of LSTM used by DeepSSM is simpler than the
Transformer used by SSDNet.

Overall, the training speed of SSDNet and SSDNet-LSTM is significantly
improved due to the changes in SSM and their testing speed is competitive with
DeepSSM.

Figure 4: Actual vs SSDNet predicted data with trend and seasonality
components and 95$\%$ confidence intervals: (a) and (b) - Sanyo; (c) and (d) -
Hanergy; (e) and (f) - Solar; (g) and (h) - Electricity data sets

Figure 5: Actual vs SSDNet predicted data with trend and seasonality
components and 95$\%$ confidence intervals for Exchange data set

Figure 6: Actual vs SARIMAX predicted data and 95$\%$ confidence intervals for
Exchange data set

Figure 7: SSDNet case study on Electricity data set (a) actual vs predicted
data; (b) attention patterns of the first head in the last layer and (c)
attention patterns of the second head in the last layer

###  V-D Interpretability Analysis

Fig. 4 presents the SSDNet output for a 7-day test set sample for the Sanyo,
Hanergy, Solar and Electricity data sets, and Fig. 5 presents the SSDNet
output for 2 test samples from the Exchange data set. We can see that the
results are interpretable: the trend is smooth and monotonically decreasing or
increasing, the seasonality shows regular and cyclic fluctuations (except for
the Exchange data set which doesn’t have seasonality).

SSDNet can model a rapidly changing trend - for example, in Fig. 4 (d), the
solar power generation drops rapidly on day 7 and SSDNet correctly outputs a
declining trend. SSDNet can not only learn smooth seasonal patterns but also
has the ability to learn seasonal patters with small random fluctuations (e.g.
at days 1 and 6 of Fig. 4 (a)) due to the introduction of innovation term
$c_{t}$. For the Exchange data set which doesn’t have clear seasonality, the
peak-to-peak magnitude of the learned SSDNet seasonal component is small.

Fig. 7 (a) presents the results for 8 consecutive days from the test set for
the Electricity data - 7-day past history and 1-day forecasting output of
SSDNet, and shows obvious daily seasonality. Fig. 7 (b) and (c) show the
attention patterns of two heads in the last multi-head attention layer of
SSDNet. The attention mappings show the importance of the previous time steps
for predicting the future time steps.

The first head tends to attend to the first and seventh day before the
predicted day, which is consistent with the weekly cycle of electricity data.
The second head attends to almost all days. A possible explanation is that the
first head mainly focuses on learning the seasonality, while the second head
mainly focuses on learning the trend.

##  VI Conclusion

We presented SSDNet, a novel deep learning approach, for probabilistic and
interpretable forecasting of time series data. SSDNet combines the advantages
of deep learning with the interpretability of SSM models. SSDNet employs the
Transformer architecture to learn the temporal patterns, extract latent
components and estimate the parameters of SSM. It then applies SSM to generate
the interpretable forecasting results with nonstationary trend and seasonality
components. SSDNet also applies attention mechanism to visualize the important
past time steps for the predicted future steps.

We evaluated the performance of SSDNet on five time series forecasting tasks.
The results show that in terms of accuracy SSDNet outperformed the state-of-
the-art deep learning models DeepSSM, DeepAR, LogSparse Transformer, Informer
and N-BEATS, and the statistical models SARIMAX and Prophet. It was also able
to provide interpretable results by showing clear trend and seasonality
components. Our ablation study showed the effectiveness of the SSM part and
that SSDNet with LSTM instead of Transformer also performed well. SSDNet was
also much faster to train than DeepSSM. Hence, the results show that SSDNet is
a promising method for time series forecasting.

## References

  * [1] J. Durbin and S. J. Koopman, _Time Series Analysis by State Space Methods_. Oxford University Press, 2001\.
  * [2] B. N. Oreshkin, D. Carpov, N. Chapados, and Y. Bengio, “N-BEATS: Neural basis expansion analysis for interpretable time series forecasting,” in _Proceedings of the International Conference on Learning Representations (ICLR)_ , 2020.
  * [3] S. Li, X. Jin, Y. Xuan, X. Zhou, W. Chen, Y.-X. Wang, and X. Yan, “Enhancing the locality and breaking the memory bottleneck of Transformer on time series forecasting,” in _Proceedings of the Conference on Neural Information Processing Systems (NeurIPS)_ , 2019.
  * [4] D. Salinas, V. Flunkert, J. Gasthaus, and T. Januschowski, “DeepAR: Probabilistic forecasting with autoregressive recurrent networks,” _International Journal of Forecasting_ , vol. 36, no. 3, pp. 1181 – 1191, 2020.
  * [5] S. S. Rangapuram, M. W. Seeger, J. Gasthaus, L. Stella, Y. Wang, and T. Januschowski, “Deep state space models for time series forecasting,” in _Proceedings of the Conference on Neural Information Processing Systems (NeurIPS)_ , 2018.
  * [6] Y. Wang, A. Smola, D. Maddix, J. Gasthaus, D. Foster, and T. Januschowski, “Deep factors for forecasting,” in _Proceedings of the International Conference on Machine Learning (ICML)_ , 2019.
  * [7] H. Zhou, S. Zhang, J. Peng, S. Zhang, J. Li, H. Xiong, and W. Zhang, “Informer: Beyond efficient transformer for long sequence time-series forecasting,” in _Proceedings of the Association for the Advancement of Artificial Intelligence (AAAI)_ , 2021.
  * [8] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Proceedings of the Conference on Neural Information Processing Systems (NeurIPS)_ , 2017.
  * [9] T. Guo, T. Lin, and N. Antulov-Fantulin, “Exploring interpretable LSTM neural networks over multi-variable data,” in _Proceedings of the International Conference on Machine Learning (ICML)_ , 2019.
  * [10] L. Li, J. Yan, X. Yang, and Y. Jin, “Learning interpretable deep state space model for probabilistic time series forecasting,” in _Proceedings of the International Joint Conference on Artificial Intelligence (IJCAI)_ , 2019.
  * [11] D1, “Sanyo dataset,” http://dkasolarcentre.com.au/source/alice-springs/dka-m4-b-phase, 2020.
  * [12] D2, “Hanergy dataset,” http://dkasolarcentre.com.au/source/alice-springs/dka-m16-b-phase, 2020.
  * [13] D3, “Solar dataset,” https://www.nrel.gov/grid/solar-power-data.html, 2014.
  * [14] D4, “Electricity dataset,” https://archive.ics.uci.edu/ml/datasets/
ElectricityLoadDiagrams20112014, 2015\.

  * [15] D5, “Exchange dataset,” https://github.com/laiguokun/multivariate-time-series-data, 2017.
  * [16] Y. Lin, I. Koprinska, and M. Rana, “SpringNet: Transformer and Spring DTW for time series forecasting,” in _Proceedings of the International Conference on Neural Information Processing (ICONIP)_ , 2020.
  * [17] Y. Lin, I. Koprinska, and M. Rana, “Temporal convolutional attention neural networks for time series forecasting,” in _Proceedings of the International Joint Conference on Neural Networks (IJCNN)_ , 2021.
  * [18] R. G. Krishnan, U. Shalit, and D. Sontag, “Deep Kalman filters,” _arXiv preprint arXiv: 1511.05121_ , 2015.
  * [19] P. Becker, H. Pandya, G. Gebhardt, C. Zhao, C. J. Taylor, and G. Neumann, “Recurrent Kalman networks: Factorized inference in high-dimensional deep feature spaces,” in _Proceedings of the Time Series Workshop at International Conference on Machine Learning (ICML)_ , 2019.
  * [20] S. Eleftheriadis, T. Nicholson, M. Deisenroth, and J. Hensman, “Identification of gaussian process state space models,” in _Proceedings of the Conference on Neural Information Processing Systems (NeurIPS)_ , 2017.
  * [21] R. J. Hyndman and G. Athanasopoulos, _Forecasting: principles and practice_ , 3rd ed. OTexts, 2021.
  * [22] M. Courbariaux, Y. Bengio, and J. David, “Binaryconnect: Training deep neural networks with binary weights during propagations,” in _Proceedings of the Conference on Neural Information Processing Systems (NeurIPS)_ , 2015.
  * [23] S. J. Taylor and B. Letham, “Forecasting at scale,” _The American Statistician_ , vol. 72, no. 1, pp. 37–45, 2018.
  * [24] H.-F. Yu, N. Rao, and I. S. Dhillon, “Temporal regularized matrix factorization for high-dimensional time series prediction,” in _Proceedings of the Conference on Neural Information Processing Systems (NeurIPS)_ , 2016.
  * [25] G. Lai, W.-C. Chang, Y. Yang, and H. Liu, “Modeling long- and short-term temporal patterns with deep neural networks,” in _Proceedings of the International ACM SIGIR Conference on Research and Development in Information Retrieval (SIGIR)_ , 2018.