DYNAMIC TIME LAG REGRESSION: PREDICTING WHAT & WHEN
Year: 2020

Abstract:
This paper tackles a new regression problem, called Dynamic Time-Lag Regression (DTLR), where a cause signal drives an effect signal with an unknown time delay. The motivating application, pertaining to space weather modelling, aims to predict the near-Earth solar wind speed based on estimates of the Sun’s coronal magnetic field. DTLR differs from mainstream regression and from sequence-to-sequence learning in two respects: firstly, no ground truth (e.g., pairs of associated subsequences) is available; secondly, the cause signal contains much information irrelevant to the effect signal (the solar magnetic field governs the solar wind propagation in the heliosphere, of which the Earth’s magnetosphere is but a minuscule region). A Bayesian approach is presented to tackle the specifics of the DTLR problem, with theoretical justifications based on linear stability analysis. A proof of concept on synthetic problems is presented. Finally, the empirical results on the solar wind modelling task improve on the state of the art in solar wind forecasting.

heading: 1 INTRODUCTION
text: A significant body of work in machine learning concerns the modeling of spatio-temporal phenomena (Shi and Yeung, 2018; Rangapuram et al., 2018), including the causal analysis of time series Peters et al. (2017), with applications ranging from markets (Pennacchioli et al., 2014) to bioinformatics (Brouard et al., 2016) to climate (Nooteboom et al., 2018).
This paper focuses on the problem of modeling the temporal dependency between two spatio-temporal phenomena, where the latter one is caused by the former one (Granger, 1969; Runge, 2018) with a non-stationary time delay.
The motivating application domain is that of space weather. The sun, a perennial source of charged energetic particles, is at the origin of geomagnetic phenomena within the sun-earth system. Specifically, the sun ejects charged particles into the surrounding space in all directions and some of these particle clouds, a.k.a. solar wind, reach the Earth’s vicinity. High speed solar wind is a major threat for the modern world, causing severe damages to e.g., satellites, telecommunication infrastructures, under sea pipelines, among others.1
A key prediction task thus is to forecast the speed of the solar wind in the vicinity of the Earth (Munteanu et al., 2013; Haaland et al., 2010; Reiss et al., 2019), sufficiently early to emit an alarm and be able to prevent the damage to the best possible extent. Formally the goal is to model the dependency between heliospheric observations (available at light speed), referred to as cause series, and the solar wind speed series recorded at the Lagrangian point L1 (a point on the Sun-Earth line 1.5 million kilometers away from the Earth), referred to as effect series. The key difficulty is that the
1The adverse impact of space weather is estimated to cost 200 to 400 million USD per year, but can sporadically lead to much larger losses.
time lag between an input and its effect, the solar wind recorded at L1, varies from circa 2 to 5 days depending on, among many factors, the initial direction of emitted particles and their energy. Would the lag be constant, the solar wind prediction problem would boil down to a mainstream regression problem. The challenge here is to predict, from the solar image x(t) at time t the value y(t+ ∆t) of the solar wind speed reaching the earth at time t+ ∆t where both the value y(t+ ∆t) and the time lag ∆t depend on x(t).
Related work. Indeed, the modeling of dependencies among time series has been intensively tackled (see e.g., Zhou and Sornette (2006); Runge (2018)). When considering varying time lag, many approaches rely on dynamic time warping (DTW) (Sakoe and Chiba, 1978). For instance, DTW is used in Gaskell et al. (2015), taking a Bayesian approach to achieve the temporal alignment of both series under some restricting assumptions (considering slowly varying time lags and linear relationships between the cause and effect time series). More generally, the use of DTW in time series analysis relies on simplifying assumptions on the cause and effect series (same dimensionality and structure) and builds upon available cost matrices for the temporal alignment. Also related is sequence-to-sequence learning (Sutskever et al., 2014), primarily aimed to machine translation. While Seq2Seq modelling relaxes some former assumptions (such as the fixed or comparable sizes of the source and target series), it still relies on the known segmentation of the source series into disjoint units (the sentences), each one being mapped into a large fixed-size vector using an LSTM; and this vector is exploited by another LSTM to extract the output sequence. Attention-based mechanisms Graves (2013); Bahdanau et al. (2015) alleviate the need to encode the full source sentence into a fixed-size vector, by learning the alignment and allowing the model to search for the parts of the source sentence relevant to predict a target part. More advanced attention mechanisms (Kim et al., 2017; Vaswani et al., 2017) refine the way the source information is leveraged to produce a target part. But to our best knowledge, the end-to-end learning of the sequence-to-sequence modelling relies on the segmentation of the source and target series, and the definition of associated pairs of segments (e.g. the sentences).
Our claim is that the regression problem of predicting both what the effect is and when the effect is observed, called Dynamic TimeLag Regression (DTLR), constitutes a new ML problem: With respect to the modeling of dependencies among time series, it involves stochastic dependencies of arbitrary complexity; the relationship between the cause and the effect series can be non-linear (the what model). Furthermore, the time lag phenomenon (the when model) can be non smooth (as opposed to e.g. Zhou and Sornette (2006)). With respect to sequence-to-sequence translation, a main difference is that the end-to-end training of the model cannot rely on pairs of associated units (the sentences), adversely affecting the alignment learning. Lastly, and most importantly, in the considered DTLR problem, even if the cause series has high information content, only a small portion of it is relevant to the prediction of the effect series. On one hand, the cause series might be high dimensional (images) whereas the effect series is scalar; on the other hand, the cause series governs the solar wind speed in the whole heliosphere and not just in near-Earth space. In addition to avoiding typically one or two orders of magnitude expansion of an already large input signal dimension, inserting the time-lag inference explicitly in the model can also potentially improve its interpretability.
Organization of the paper. The Bayesian approach proposed to tackle the specifics of the DTLR regression problem is described in section 2; the associated learning equations are discussed, followed by a stability analysis and a proof of consistency (section 3). The algorithm is detailed in section 4. The experimental setting used to validate the approach is presented in section 5; the empirical validation on toy problems and on the real-world problem are discussed in section 6

heading: 2 PROBABILISTIC DYNAMICALLY DELAYED REGRESSION

heading: 2.1 POSITION OF THE PROBLEM
text: Given two time series, the cause series x(t) (x(t) ∈ X ⊂ RD) and the observed effect series y(t), the sought model consists of a mapping f(.) which maps each input pattern x(t) to an output y(φ(t)), and a mapping g(.) which determines the time delay φ(t)− t between the input and output patterns:
y(φ(t)) = f [x(t)] (1) φ(t) = t+ g[x(t)] (2)
with f : X → R, and g : X → R+,
where t ∈ R+ represents the continuous temporal domain. The input signal x(t) is possibly high dimensional and contains the hidden cause to the effect y(t) ∈ R; y(t) is assumed to be scalar in the remainder of the paper. Function g : X → R+ represents the time delay between inputs and outputs. Vectors are written using bold fonts.
As said, Eqs 1-2 define a regression problem that differs from standard regression in two ways: Firstly, the time lag g[x(t)] is non-stationary as it depends on x(t). Secondly, g[x(t)] is unknown, i.e. it is not recorded explicitly in the training data.
Assumption. For the sake of the model identifiability and computational stability, the time warping function φ(t) = t+ g[x(t)] is assumed to be sufficiently regular w.r.t. t. Formally, φ(.) is assumed to be continuous2.

heading: 2.2 PROBABILISTIC DYNAMIC TIME-LAG REGRESSION
text: For practical reasons, cause and effect series are sampled at constant rate, and thereafter noted xt and yt with t in N . Accordingly, mapping g maps xt onto a finite set T of possible time lags, with T = {∆tmin . . . ,∆tmax} an integer interval defined from domain knowledge. The unavoidable error due to the discretization of the continuous time lag is mitigated along a probabilistic model, associating to each cause x, a set of predictors ŷ(x) = {ŷi(x), i ∈ T } and a probability distribution p̂(x) on T estimating the probability of delay of the effects of x. Overall, the DTLR solution is sought as a probability distribution conditioned on cause x, mixture of Gaussians3 centered on the predictors ŷi(x), where the mixture weights are defined from p̂(x). More formally, letting yt denote the vector of random variables {yt+i, i ∈ T }:
P [ yt|xt = x ] = ∑ {τi∈{0,1},i∈T } p̂ ( τ1, . . . , τ|T ||x)N ( ŷ(x),Σ(τ) ) (3)
with Σ = Diag(σi(τ)2) the diagonal matrix of variance parameters attached to each time-lag i ∈ T . Two simplifications are made for the sake of the analysis. Firstly, the stochastic time lag is modelled as the vector (τi), i ∈ T of binary latent variables, where τi indicates whether xt drives yt+i (τi = 1) or not (τi = 0). The assumption that every cause has a single effect is modelled by imposing: 4∑
i∈T τi = 1. (4)
From constraint (4), probability distribution p̂(x) thus is sought as the vector (p̂i(x)) for i in T , summing to 1, such that p̂i(x) stands for the probability of the effect of xt = x to occur with delay i. The second simplifying assumption is that the variance σ2i (τ) of predictor ŷi does not depend on x, by setting:
σi(τ) −2 = ( 1 + ∑ j αijτj ) σ−2,
with σ2 a default variance and αij ≥ 0 a matrix of non-negative real parameters. This particular formulation supports the tractable analysis of the posterior probability of τi (in supplementary
2 For some authors (Zhou and Sornette, 2006),the monotonicity of φ(.) is additionally required and enforced using constraints: φ(t1) ≤ φ(t2), ∀t1 ≤ t2. This assumption is not retained as one may achieve a similar effect by using regularization based smoothness penalties.
3Note that pre-processing can be used if needed to map non-Gaussian data onto Gaussian data. 4Note however that the cause-effect correspondence might be many-to-one, with an effect depending on
several causes.
material). The fact that x can influence yi through predictor ŷi(x) even when τi = 0 reflects an indirect influence due to the auto-correlation of the y series. This influence comes with a higher variance, enforced by making αij a decreasing function of |i− j|. More generally, a large value of αii compared to αij for i 6= j corresponds to a small auto-correlation time of the effect series.

heading: 2.3 LEARNING CRITERION
text: The joint distribution is classically learned by maximizing the log likelihood of the data, which can here be expressed in closed form. Let us denote respectively the dataset and parameters as {(x,y)}data and θ = (ŷ, p̂, σ, α). From Eq. (3) the conditional probability qi(x,y) def = P (τi = 1|x,y) reads:
qi(x,y) = 1
Z(x,y|θ) p̂i(x) exp
( − 1
2σ2 ∑ j∈T αji ( yj − ŷj(x) )2 + 1 2 ∑ j∈T log(1 + αji) )
(5)
with normalization constant Z(x,y|θ) = ∑ i∈T p̂i(x) exp ( − 1 2σ2 ∑ j∈T αji ( yj − ŷj(x) )2 + 1 2 ∑ j∈T log(1 + αji) ) .
The log-likelihood then reads (intermediate calculations in supplementary material, appendix A): L[{(x,y)}data|θ] = −|T | log(σ)− Edata [∑ i∈T 1 2σ2 ( yi − ŷi(x) )2 − log(Z(x,y|θ))] (6) where Edata denotes averaging over the dataset. For notational simplicity, the time index t is omitted in the following and the empirical averaging on the data is noted Edata. The hyper-parameters σ and matrix α of the model are obtained by optimizing L:
σ2
1 + αij =
Edata [( yi − ŷi(x) )2 qj(x,y) ] Edata [ qj(x,y)
] , (7) In addition the optimal ŷ and p̂ reads:
ŷi(x) = Edata
[ yi ( 1 + ∑ j∈T αijqj(x,y) )∣∣∣x] Edata [ 1 + ∑ j∈T αijqj(x,y)
∣∣∣x] (8) p̂i(x) = Edata [ qi(x,y)
∣∣∣x], (9) where the above conditional empirical averaging operates as an averaging over samples close to x.
These are self-consistent equations, since qi(x,y) depends on the parameters σ2 and αij , ŷ and p̂. The proposed algorithm detailed in section 4 implements the saddle point method defined from Eqs (7,5,8,9): alternatively, hyper-parameters σ and αij are updated from Eq. (7) based on the current ŷi and p̂i; and predictors ŷi and mixture weights p̂i are updated according to Eqs (8) and (9) respectively.

heading: 3 THEORETICAL ANALYSIS
text: The proposed DTLR approach is shown to be consistent and analyzed in the simple case where α is a diagonal matrix (αij = αδij).

heading: 3.1 LOSS FUNCTION AND RELATED OPTIMAL PREDICTOR
text: Let us assume that the hyper-parameters of the model have been identified together with predictors ŷi(x) and weights p̂i(x). These are leveraged to achieve the prediction of the effect series. For any given input x, the sought eventual predictor is expressed as (ŷ(x), Î(x)) where Î(x) is the predicted time lag and ŷ(x) the predicted value. The associated L2 loss is:
L2(ŷ, Î) = Edata [( yÎ(x) − ŷ(x) )2] . (10)
Then it comes:
Proposition 3.1. With same notations as in Eq. (3), with αij = αδij , α > 0, the optimal composite predictor (y?, I?) is given by
y?(x) = ŷI?(x)(x) with I?(x) = arg max i
( p̂i(x) ) ,
Proof. In supplementary material, Appendix C.

heading: 3.2 LINEAR STABILITY ANALYSIS
text: The saddle point (Eqs 7, 5, 8, 9) admits among others a degenerate solution, corresponding to p̂i(x) = 1/|T |, αij = 0 for all pairs (i, j), with σ2 = σ20 . Informally the model converges toward this degenerate trivial solution when there is not enough information to build specialized predictors ŷi.
Let us denote ∆y2i (x) = ( yi − ŷi(x) )2 the square error made by predictor ŷi for x, and
σ20 = 1
|T | Edata (∑ i∈T ∆y2i (x) )
the average of MSE over the set of the predictors ŷi, i ∈ T . Let us investigate the conditions under which the degenerate solution may appear, by computing the Hessian of the log-likelihood and its eigenvalues. Under the simplifying assumption
αij = αδij ,
the model involves 2|T | functional parameters ŷ and p̂ and two hyper-parameters α and r = σ2/σ20 . After the computation of the Hessian (in supplementary material, Appendix B) the system involves three key statistical quantities, two global ones:
C1[q] = 1
σ20 Edata (∑ i∈T qi(x,y)∆y2i (x) ) , (11)
C2[q] = 1
σ40 Edata [∑ i∈T qi(x,y) ( ∆y2i (x)− ∑ j∈T qj(x,y)∆y2j (x) )2] , (12)
and a local |T |-vector of components
ui[x,q] = 1
σ20 Edata
[ qi(x,y) ( ∆y2i (x)− ∑ j∈T qj(x,y)∆y2j (x) )∣∣∣x].
Up to a constant, C1 represents the covariance between the latent variables {τi} and the normalized predictor errors. C1 smaller than one indicates a positive correlation between the latent variables and small errors; the smaller the better. For the degenerate solution, i.e. q = q0 uniform, C1[q0] = 1 and C2[q0] represents the default variability among the prediction errors. ui[x,q] informally measures the quality of predictor ŷi relatively to the other ones at x. More precisely, a negative value of ui[x,q] indicates that ŷi is doing better than average in the neighborhood of x.
At a saddle point the parameters are given by: σ2
σ20 = |T | − C1[q] |T | − 1
and α = |T | |T | − 1 1− C1[q] C1[q] .
The predictors ŷ are decoupled from the rest whenever they are centered, which we assume. So the analysis can focus on the other parameters.
If p̂ is fixed a saddle point is stable iff
C2[q] < 2C 2 1 [q] +O ( 1 |T | ) .
In particular, the degenerate solution is unstable if C2[q0] > 2 ( 1− 1 |T | ) .
Note that for ∆yi(x) iid centered with variance σ20 and relative kurtosis κ (conditionally to x) one has C2 = (2 + κ)(1 − 1/|T |). Therefore, whenever ∆y2i (x) fluctuates and the relative kurtosis is non-negative, the degenerate solution is unstable and will thus be avoided.
If p̂ is allowed to evolve (after Eq. (9)) the degenerate trivial solution becomes unstable as soon as C2[q0] is non-zero, due to the fact that the gradient points in the opposite direction to u(x) (with dp̂(x) ∝ −C2[q0]u(x)), thus rewarding the predictors with lowest errors by increasing their weights.
The system is then driven toward other solutions, among which the localized solutions of the form:
p̂i(x) = δi,I(x),
with an input dependent index I(x) ∈ T . As shown (in supplementary material, Appendix C) the maximum likelihood localized solution also minimizes the loss function (Eq. 10). The stability of such localized solutions and the existence of other (non-localized) solutions is left for further work.

heading: 4 THE DTLR ALGORITHM
text: The DTLR algorithm learns both regression models ŷ(x) and p̂(x) from series xt and yt, using alternate optimization of the model parameters and the model hyper-parameters α and σ2, after Eqs (7,5,8,9). The model search space is that of neural nets, parameterized by their weight vector θ. The inner optimization loop updates θ using mini-batch based stochastic gradient descent. At the end of each epoch, after all minibatches have been considered, the outer optimization loop computes hyper-parameters α and σ2 on the whole data.
Initialization of α and σ it←− 0 ; while it < max do
while epoch do θ ←− Optimize(L(θ, α, σ2)) ; end σ2 ←− σ20 |T |−C1[q] |T |−1 ; α←− |T ||T |−1 1−C1[q] C1[q]
; end Result: Model parameters θ = {ŷ, p̂}, hyper-parameters α, σ2
Algorithm 1: DTLR algorithm
The algorithm code is available in supplementary material and will be made public after the reviewing period. The initialization of hyper-parameters α and σ is settled using preliminary experiments (same setting for all considered problems: α ∼ U(0.75, 2); σ2 ∼ U(10−5, 5)). The neural architecture implements predictors ŷ(x) and weights p̂(x) on the top of a same feature extractor from input x. In the experiments, the architecture of the feature extractor is a 2-hidden layer fully connected network. On the top of the feature extractor are the single layer ŷ and p̂ models, each with |T | output neurons, with |T | the size of the chosen domain for the time lag.

heading: 5 EXPERIMENTAL SETTING
text: The goal of experimental validation is threefold. A first issue regards the accuracy of DTLR, measured from the mean absolute error (MAE), root mean square error (RMSE) and Pearson correlation of the learned DTLR model (y?(xt), I?(x)). DTLR is compared to the natural baseline defined as the regressor with constant time lag, ŷ∆̄(xt), with ∆̄ being the average of all possible time lags in T . The predictions of DTLR and the baseline are compared with the ground truth value of the effect series. However the predicted time lag I?(xt) can only be assessed if the ground truth time-lag relationship is known. In order to do so, three synthetic problems of increasing difficulty are defined below.
Secondly, the stability and non-degeneracy of the learned model are assessed from the statistical quantities σ0 and C1 (section 3.2), compared to the degenerate solution p̂i(x) = 1/|T |. For C1 < 1, the model accurately specializes the found predictors p̂i.
Lastly, and most importantly, DTLR is assessed on the solar wind prediction problem, and compared to the best state of the art in space weather.
Synthetic Problems. Four synthetic problems of increasing difficulty are generated using Stochastic Langevin Dynamics. In all problems, the cause signal xt ∈ R10 and the effect signal yt are generated as follows (with η = 0.02, s2 = 0.7):
xt+1 = (1− η)xt +N (0, s2) (13) vt = k||xt||2 + c (14)
yt+g(xt) = f(vt), (15)
with time-lag mapping g(xt) ranges in a time interval with width 20 (except for problem I where |T | = 15). The complexity of the synthetic problems is governed by the amplitude and time-lag functions f and g (more in appendix, Table 2):
Problem f(vt) g(xt) Other I vt 5 k=10,c=0 II vt 100/vm k=1,c=10 III √ v2t+2ad ( √ v2m+2ad−v)/a k=5,a=5,d=1000,c=100 IV vt g(xt)=exp(vt)/(1+exp(vt/20)) k=10,c=40
Solar Wind Speed Prediction. The challenge of predicting solar wind speed from heliospheric data is due to the non-stationary propagation time of the solar plasma through the interplanetary medium. For the sake of a fair comparison with the best state of the art Reiss et al. (2019), the same experimental setting is used. The cause series xt includes the solar magnetic (flux tube expansion, FTE: the rate at which the magnetic flux tubes expand between the photosphere and a reference height in the corona at 2.5 solar radii) and the coronal magnetic field strength estimates produced by the Current Sheet Source Surface model (Zhao and Hoeksema, 1995; Poduval and Zhao, 2014; Poduval, 2016), exploiting the hourly magnetogram data recorded by the Global Oscillation Network Group from 2008 to 2016. The effect series, the hourly solar wind data is available from the OMNI data base from the Space Physics Data Facility 5. After domain knowledge, the time-lag ranges from 2 to 5 days, segmented in six-hour segments (thus |T | = 12). For the i-th segment, the "ground truth" solar wind yi is set to its median value over the 6 hours.
DTLR is validated using a nine fold cross-validation (Table 3 in appendix), where each fold is a continuous period corresponding to a solar rotation.6

heading: 6 EMPIRICAL VALIDATION
text: Table 1 summarizes the DTLR performance on the synthetic and solar wind problems (detailed results are provided in the appendix).

heading: 6.1 SYNTHETIC PROBLEMS
text: On the easy Problem I, DTLR predicts the correct time lag for 97.93% of the samples. The higher value of σ0 in problems I and II compared to the other problems is explained from the higher variance in the effect series y(t).
5https://omniweb.gsfc.nasa.gov 6 The Sun completes a rotation (or Carrington rotation) in approximately 27 days.
On Problem II, DTLR accurately learns the inverse relationship between xt, g(xt) and yt on average. The time lag is overestimated in the regions with low time lag (with high velocity), which is blamed on the low sample density in this region, due to the data generation process. Interestingly, Problems III and IV are better handled by DTLR, despite a more complex dynamic time lag relationship. In both latter cases however, the model tends to under-estimate the time lag in the high time lag region and conversely to over-estimate it in the low time lag region.

heading: 6.2 THE SOLAR WIND PROBLEM
text: DTLR finds an operational solar wind model (Table 1), though the significantly higher difficulty of the solar wind problem is witnessed by the C1 value close to the degenerate value 1. The detailed comparison with the state of the art Reiss et al. (2019) (Fig. 1, Left) shows that DTLR improves on the current best state of the art (on all variants including ensemble approaches, and noting that median models are notoriously hard to beat). (Fig. 1, Right) shows the good correlation between the predicted solar wind7 and the measured solar wind.
Model M.A.E R.M.S.E WS 74.09 85.27 DCHB 83.83 103.43 WSA 68.54 82.62 Ensemble Median (WS) 71.52 83.36 Ensemble Median (DCHB) 78.27 100.04 Ensemble Median (WSA) 62.24 74.86 Persistence (4 days) 130.48 161.99 Persistence (27 days) 66.54 78.86 Fixed Lag Baseline 67.33 80.39 DTLR 60.19 72.64
(a) Comparative assessment on the Solar Wind problem compared to the state of the art Reiss et al. (2019, Table 1)
300
400
500
600
300 400 500 600 700
v (km/s)
v̂ (k
m /s
)
(b) Scatter Chart (9 fold CV)
Figure 1: DTLR on the solar wind problem. Left: comparative quantitative assessment w.r.t. the state of the art (Carrington rotation 2077). Right: qualitative assessment of the prediction.

heading: 7 DISCUSSION AND PERSPECTIVES
text: The contribution of the paper is twofold. A new ML setting, Dynamic Time Lag Regression has been defined, aimed at the modelling of varying time-lag dependency between time series. The introduction of this new setting is motivated by an important scientific and practical problem from the domain of space weather, an open problem for over two decades.
Secondly, a Bayesian formalization has been proposed to tackle the DTLR problem, relying on a saddle point optimization process. A closed form analysis of the training procedure stability under simplifying assumptions has been conducted, yielding a practical alternate optimization formulation, implemented in the DTLR algorithm. This algorithm has been successfully validated on synthetic and real-world problems, although some bias toward the mean has been detected in some cases.
On the methodological side, this work opens a short term perspective (handling the bias) and a longer term perspective, extending the proposed nested inference procedure and integrating the model selection step within the inference architecture. The challenge is to provide the algorithm with the means of assessing online the stability and/or the degeneracy of the learning trajectory.
Regarding the motivating solar wind prediction application, a next step consists of enriching the data sources and the description of the cause series xt, typically by directly using the solar images. Another perspective is to consider other applications of the general DTLR setting, e.g. considering fine-grained modelling of diffusion phenomena.
7The predicted values, every 6 hours, are interpolated for comparison with the hourly measured solar wind.

heading: APPENDIX A LOG LIKELIHOOD OF THE LATENT MODEL (3)
text: A.1 DIRECT COMPUTATION
Due to the single effect constraint (4) the mixture model (3) can be expressed simply as P (y|x) = (∑ i∈T p̂i(x) ∏ j∈T √ 1 + αji 2πσ2 e− 1 2σ2 (1+αji) ( yj−ŷj(x) )2)
= (∑ i∈T p̂i(x) ∏ j∈T √ 1 + αji 2πσ2 e− 1 2σ2 αji ( yj−ŷj(x) )2) exp ( − 1 2σ2 ∑ j∈T ( yj − ŷj(x) )2)
Let θ def= (ŷ, p̂, σ, α) denote the parameters of the model and consider the probability that predictor ŷi is the good one conditionally to a pair of observation (x,y):
qi(x,y) = P (τi = 1|x,y)
= 1
Z(x,y|θ) p̂i(x) exp
( − 1
2σ2 ∑ j∈T αji ( yj − ŷj(x) )2 + 1 2 ∑ j∈T log(1 + αji) )
with
Z(x,y|θ) = ∑ i∈T p̂i(x) exp ( − 1 2σ2 ∑ j∈T αji ( yj − ŷj(x) )2 + 1 2 ∑ j∈T log(1 + αji) ) .
This gives immediately L[{(x,y)}data|θ] = −|T | log(σ)− Edata [∑ i∈T 1 2σ2 ( yi − ŷi(x) )2 − log(Z(x,y|θ))] A.2 LARGE DEVIATION ARGUMENT
Even though the log likelihood can be obtained by direct summation, for sake of generality we show how this can result from a large deviation principle. Assume that the number of learning samples tends to infinity, and so that in a small volume dv = dxdy around a given joint configuration (x,y), the number of data Nx,y becomes large. Restricting the likelihood to this subset of the data yields the following:
Lx,y = Nx,y∏ m=1 ∑ {τ(m)} p̂(τ (m)|x)∏ i∈T √ 2π σi(τ (m)) exp ( −1 2 ∑ i∈T ( yi − ŷi(x) )2 σi(τ (m))2 ) .
Upon introducing the relative frequencies:
qi(x,y) = 1
Nx,y Nx,y∑ m=1 τ (m) i satisfying
∑ i∈T qi(x,y) = 1,
the sum over the τ (m)i is replaced by a sum over these new variables, with the summand obeying a large deviation principle
Lx,y ∑ q exp ( −Nx,yFx,y [ q ])
where the rate function reads Fx,y [ q ] = |T | log(σ) + ∑ i∈T [( yi − ŷi(x) )2 1 +∑j∈T αijqj 2σ2 − 1 2 qi ∑ j∈T log(1 + αji) + qi log qi p̂i ] .
Taking the saddle point for qi yield as a function of (x,y) expression (7). Inserting this into F and taking the average over the data set yields the log likelihood (5) with opposite sign:
L[{(x,y)}data|θ] = −Edata [ Fx,y [ q(x,y) ]] .
A.3 SADDLE POINT EQUATIONS
Now we turn to the self-consistent equations relating the parameters θ of the model at a saddle point of the log likelihood function. First, the optimization of the predictors ŷ yields:
∂L ∂ŷi(x) = 1 σ2 Edata
[( yi − ŷi(x) )( 1 + ∑ j∈T αijqj(x,y) )∣∣∣x].
Then the optimization of p̂ gives: ∂L
∂p̂i(x) = Edata [qi(x,y) p̂i(x) − λ(x) ∣∣∣x],
= 1
p̂i(x) Edata
[ qi(x,y) ∣∣∣x]− λ(x)
with λ(x) a Lagrange multiplier to insure that ∑ i p̂i(x) = 1 This gives
p̂i(x) = 1
λ(x) Edata
[ qi(x,y) ∣∣∣x] Hence ∑
i∈T p̂i(x) =
1
λ(x) = 1 ∀x
in order to fulfill the normalization constraint, yielding finally expression (9).
Finally the optimization of α reads:
∂L ∂αij = 1 2(1 + αij) Edata
[ qj(x,y) ] − 1
2σ2 Edata
[( yi − ŷi(x) )2 qj(x,y) ] .

heading: APPENDIX B PROOF OF PROPOSITION 3.1
text: Given I(x) a candidate index function we associate the point-like measure
pi(x) = δi,I(x).
Written in terms of p the loss function reads L2(ŷ, p) = Ex,y [∑ i∈T pi(x) ( yi − ŷ(x) )2] .
Under (3) (with αij = αδij) the loss is equal to
L2(ŷ, p) = Ex [∑ i∈T pi(x) (( ŷi(x)− ŷ(x) )2 − p̂i(x) ασ2 1 + α )] + σ2
The minimization w.r.t. ŷ yields ŷ(x) = ∑ i∈T pi(x)ŷi(x). (16)
In turn, as a function of pi the loss being a convex combination, its minimization yields
pi(x) = δi,I(x), (17)
I(x) = arg min i∈T
(( ŷi(x)− ŷ(x) )2 − p̂i(x) ασ2 1 + α ) . (18)
Combining these equations (16,17,18) we get
I(x) = arg max i∈T
( p̂i(x) ) ,
which concludes the proof.

heading: APPENDIX C STABILITY ANALYSIS
text: The analysis is restricted for simplicity to the case αij = αδij . The log likelihood as a function of r = σ2/σ20 and β = α/r after inserting the optimal q = q(x,y) reads in that case
L(r, β) = −|T | 2 log(r)− |T | 2r + 1 2 log(1 + rβ) + Edata
[ log(Z)− λ(x) ∑ i∈T p̂i(x) ]
with
Z = ∑ i p̂i(x) exp ( − β 2σ20 ∆y2i (x) ) ,
and where λ(x) is a Lagrange multiplier which has been added to impose the normalization of p̂. The gradient reads
∂L ∂r = 1 2r2
( |T |(1− r) + βr 2
1 + βr
) ,
∂L ∂β = r 2(1 + rβ) − 1 2 C1[q],
∂L ∂ŷi(x) = 1 σ2 Edata
[( yi − ŷi(x) )( 1 + αqi(x,y) )∣∣∣x]. ∂L
∂p̂i(x) =
Edata [ qi(x,y)|x ] p̂i(x) − λ(x),
with
C1[q] = 1
σ20 Edata (∑ i∈T qi(x,y)∆y 2 i (x) ) ,
This leads to the following relation at the saddle point:
r = |T | − C1[q] |T | − 1 ,
α = |T | |T | − 1 1− C1[q] C1[q] ,
ŷi(x) = Edata
[ yi ( 1 + αqi(x,y) )∣∣∣x] Edata [ 1 + αqi(x,y)
∣∣∣x] p̂i(x) = Edata [ qi(x,y)|x ] .
Let us now compute the Hessian. It is easy to see that the block corresponding to the predictors ŷ decouples from the rest as soon as these predictors are centered.
Denoting
C2[q] = 1
σ40 Edata [∑ i∈T qi(x,y) ( ∆y2i (x)− n∑ j=1 qj(x,y)∆y 2 j (x) )2] ,
we have
∂2L ∂r2 = 1 2r2
( −|T |+ 2 |T | |T | − 1 ( C1[q]− 1 ) − β2C21 [q] ) ∂2L ∂r∂β = 1 2r2 C21 [q]
∂2L ∂β2 = 1 4
( C2[q]− 2C21 [q] ) ∂2L
∂p̂i(x)∂p̂j(x) = −
Edata [ qi(x,y)qj(x,y)|x ] p̂i(x)p̂j(x)
∂2L ∂r∂p̂i(x) = 0
∂2L ∂β∂p̂i(x) = −ui[x,q] 2p̂i(x) ,
where ui[x,q] def = 1
σ20 Edata
[ qi(x,y) ( ∆y2i (x)− ∑ j∈T qj(x,y)∆y 2 j (x) ) |x ] .
There are two blocks in this Hessian, the one corresponding to r and β and the one corresponding to derivatives with respect to p̂i. The stability of the first one depends on the sign of C2[q]− 2C21 [q] for |T | large while the second block is always stable as being an average of the exterior product of the vector (q1(x,y)/p̂1(x), . . . , q|T |(x,y)/p̂|T |(x)) by itself. At the degenerate point α = 0, r = 1, p̂i = 1/|T | the Hessian simplifies as follows. Denote
dη = dre1 + dβe2 + ∫ dx ∑ i∈T dp̂i(x)ei+2(x)
a given vector of perturbations, decomposed onto a set of unit tangent vectors, {e1 and e2} being respectively associated to r and β, while ei(x) associated to p̂i(x) for all i ∈ T and x ∈ X . Denote
u = ∑ i∈T ∫ dxui[x]ei(x)
v(x) = ∑ i∈T ei(x)
with
C2 = 1
|T |σ40 Edata [∑ i∈T ( ∆y2i (x)− 1 |T | ∑ j∈T ∆y2j (x) )2] .
ui[x] = 1
σ20 Edata
[ ∆y2i (x)− σ20 |x ] .
With these notations the Hessian reads:
H = 1
2
( −|T |e1et1 + e1et2 + e2et1 + (C2 2 − 1 ) e2e t 2 − uet2 − e2ut − ∫ dxv(x)vt(x) ) .
In fact we are interested in the eigenvalues of H in the subspace of deformations which conserve the norm of p̂, i.e. orthogonal to v(x), thereby given by
η = η1e1 + η2e2 + η3u.
In this subspace the Hessian reads
H = 1
2  −|T | 1 0 1 C2 2 − 1 −M |T |C2
0 −M |T |C2 0  , where M is the number of data points, resulting from the fact that∑
i∈T
∫ dxui[x] 2 = M
σ40 Edata [∑ i∈T ( ∆y2i (x)− σ20 )2] ,
= MC2,
because Edata(·|x) as a function of x is actually a point-wise function on the data. If |u|2 > 0 or if |u| = 0 and 1 + |T |(C2/2 − 1) > 0 there is at least one positive eigenvalue. Let Λ be such an eigenvalue. After eliminating dr and dβ from the eigenvalue equations in dη, the deformation along this mode verifies dη ∝ Λe1 + Λ(|T |+ Λ)e2 −M |T |(|T |+ Λ)C2u, which corresponds to increasing r and α while decreasing for each x the p̂i having the highest mean relative error ui[x].
Concerning solutions for which p̂i(x) = δiÎ(x)
is concentrated on some index Î(x), the analysis is more complex. In that case C2[p] = 0 and C1[p] > 0. The (r, β) sector has 2 negative eigenvalues, while the p̂ block is (−) a covariance matrix, so it has as well negative eigenvalues. The coupling between these two blocks could however in principle generate in some cases some instabilities.
Still, the log likelihood of such solutions reads
L = −|T | 2 log(σ2) + 1 2 log(1 + α)− 1 2σ2 Edata [∑ i∈T ∆y2i (x) ] − α 2σ2 Edata [ ∆y2I(x)(x) ] so we get the following optimal solution
σ2 = 1
|T | Edata [∑ i∈T ∆y2i (x) ] ,
1
1 + α =
Edata [ ∆y2I(x)(x) ] σ2 ,
I(x) = arg min i∈T
Edata [ ∆y2i (x)|x ] .

heading: APPENDIX D EXPERIMENTS: ADDITIONAL DETAILS
text: Here we provide some additional details and context to the experimental validation of the DTLR methodology described in section 5. Table 2 provides some information about the datasets used in the synthetic and solar wind prediction problems8. Sections D.1.1 and D.1.2 give additional plots for evaluating the experimental results.
For the solar wind prediction task, the solar wind data was mapped into standardized Gaussian space using a quantile-quantile and inverse probit mapping. Nine fold cross-validation was performed using splits as specified in table 3. To compare the DTLR results with the state of the art solar wind forecasting, we used results from Reiss et al. (2019, Table 1). Since Reiss et al. (2019) compared the various forecasting methods on only one solar rotation (first row of table 3), comparing these results with DTLR can be considered as a preliminary examination. Nevertheless, the results presented in table 1a show encouraging signs for the competitiveness and usefulness of the DTLR method.
It is well-established that the solar wind speed observed near the Earth’s orbit (at the L1 Lagrange point about a million kilometers upstream solar wind) bears an inverse correlation with the rate (FTE) at which the magnetic flux tubes expand between the photosphere and a reference height called the source surface, typically placed at 2.5 solar radii in many coronal models (Wang and Sheeley, 1990). That is, regions of the Sun with smaller FTE give rise to fast wind and regions with larger values of FTE emanate slow solar wind. This inverse correlation forms the basis of the current solar wind prediction techniques (Poduval and Zhao, 2014; Poduval, 2016) including the WSA-ENLIL model (Wang and Sheeley, 1990; Odstrčil et al., 1996; Odstrčil and Pizzo, 1999a;b; Odstrčil, 2003; Odstrčil et al., 2004) of NOAA/SWPC. Mathematically, FTE is defined as:
FTE = ( R Rss )2 Br(phot) Br(ss) (19)
where, Br(phot) and Br(ss) are the radial component of magnetic fields at the photosphere and the source surface, and R and Rss are the respective radii. The coronal magnetic field is computed using CSSS model, the geometry of which is shown in figure 2, by extrapolating the photospheric magnetic field measured by various ground based observatories and spacecraft. For the present work,
8In the solar wind problem, the training and test data sizes correspond to one cross-validation split
we used the photospheric field measured using the Global Oscillation Network Group (GONG), a network of ground based solar observatories over different places on the globe.
D.1 SUPPLEMENTARY PLOTS
D.1.1 SYNTHETIC PROBLEMS
200
300
400
200 300 400 Actual Output
P re
di ct
ed O
ut pu
t
(a) Problem II, Goodness of fit, Output y(x)
2
4
6
8
10
2 3 4 5 6 7 Actual Time Lag
P re
di ct
ed T
im e
La g
(b) Problem II, Goodness of fit, Time lag τ(t)
0.0
0.4
0.8
1.2
0 2 4 6 8 Error: Time Lag
de ns
ity
(c) Problem II, Error of time lag prediction
2
3
4
5
6
200 300 400 Output
T im
e La
g
Type predicted actual
(d) Problem II, Output vs Time Lag Relationship
Figure 3: Problem II, Results
150
175
200
225
250
275
150 175 200 225 250 275 Actual Output
P re
di ct
ed O
ut pu
t
(a) Problem III, Goodness of fit, Output y(x)
4
5
6
7
8
9
3 4 5 6 7 Actual Time Lag
P re
di ct
ed T
im e
La g
(b) Problem III, Goodness of fit, Time lag τ(t)
50
100
150
50 100 150 Actual Output
P re
di ct
ed O
ut pu
t
(a) Problem IV, Goodness of fit, Output y(x)
2.5
5.0
7.5
2 3 4 5 6 7 Actual Time Lag
P re
di ct
ed T
im e
La g
(b) Problem IV, Goodness of fit, Time lag τ(t)
D.1.2 SOLAR WIND PREDICTION
●●●● ●●●●
●● ●●●●●●●●
●
●● ●● ●●
●●
●●●
●●●●
●●●
●● ● ● ● ● ● ●● ●●●●
● ●● ●●
●● ●● ●● ● ●●● ●
● ●
● ● ● ● ● ● ● ● ●●●
● ● ●●●● ●●●● ● ●● ● ● ●
● ● ● ● ● ● ● ● ● ● ●● ●●●●●
● ●●●●●●●●●●●●●●●●●●●●● ●●●●●● ●●●●●● ●● ●●●●●
●●●●●●●●● ●● ●
●● ●● ●●● ● ●●● ●●● ●
●●●●●● ●●●
● ● ●● ● ● ●●
● ●
● ●●● ● ● ● ● ●● ●● ●● ●●●● ●
●●●●● ● ● ● ●●● ● ●●●
●●●● ●● ● ● ● ●●●
● ● ●●●● ●●●●
● ●●● ● ●●●
●●●●●●●●●●●●● ●●●
●●●●●●● ●● ●● ●●●●
●●● ●● ● ●
●● ●● ●● ●● ●● ●
● ●
● ●
●● ● ●●● ● ●●●
●● ● ●● ● ● ● ● ●● ● ● ● ● ●●● ●
●●● ●●●●●●●●●
● ●●●●●●●●●●●● ●●●●
● ● ● ●
●● ● ●
● ● ●● ● ● ●● ●
● ● ● ● ● ● ● ● ●● ●
● ● ●● ●● ● ● ● ●● ● ●●●● ● ● ●● ●● ● ● ●● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ●●● ●●● ● ●● ● ●●●● ●● ● ● ● ●●●
● ●● ●● ● ●
●
● ● ●
●
●
●
●
●
●
●
● ●
●
● ● ● ●
●
●
● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
● ●
● ● ● ● ● ● ● ● ● ● ● ● ●● ●
● ● ● ● ● ● ● ●
●
●
● ● ● ● ● ● ● ●
● ● ● ● ● ●● ● ● ● ● ●
● ● ●
● ● ● ●● ● ● ●● ● ● ● ● ●
● ● ● ● ●●● ●● ●● ●
● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●
● ● ● ● ●● ●●● ● ●● ●● ● ● ● ●
●● ● ●●● ● ●●
● ● ●● ● ● ●●●●● ● ● ●
●
●● ● ●● ●● ● ●
●
●
●
●
●
● ● ● ● ●
●
● ● ● ●● ● ● ● ● ● ● ● ●
● ● ● ●●● ● ● ● ● ●
● ●
● ● ●
●
● ● ● ● ● ● ● ● ● ●
●
●●
● ● ● ● ● ● ●
● ● ●● ● ● ● ● ● ●
● ● ●
● ●
● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ●●●
● ● ● ●● ●●● ● ● ● ● ●
● ● ● ●● ● ● ● ●● ●
● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ●
● ●
●
● ● ● ● ● ● ●● ●
● ● ● ●
● ● ● ● ● ● ● ●
●
●
● ●
● ●● ● ● ● ●●●
●●
●
●
●
●
●
● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●
● ●
● ●
● ● ● ● ● ●● ● ●
●● ●● ● ● ● ●●● ●● ● ●
● ● ●● ● ●●● ● ●● ● ● ● ●
● ●● ● ● ● ● ● ●
●
● ● ● ● ● ● ● ●● ●● ● ●
●● ●● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ● ● ●
● ● ●● ● ● ● ● ●● ● ●
● ● ● ● ● ● ● ●
● ● ●● ● ● ● ● ● ●● ● ● ● ● ●●
●● ● ●
● ●● ● ● ●
300
400
500
600
0 200 400 600
t (hours)
km /s
● ●Actual Prediction
(a) Hourly forecasts for period 2008-11-20 07:00 to 2008-12-17 14:00
● ●● ● ● ● ●● ● ● ● ● ● ● ●●● ●●●● ●● ●●
●● ●● ●● ●● ● ●●
● ●●● ● ●●● ●● ●●● ● ● ●●
● ● ●● ● ●●● ● ● ● ●● ● ● ●● ● ● ● ●
●
●
● ●
●
● ●
● ●
● ●
● ●
●
● ●
● ●
●
● ● ●● ● ●
●●● ●● ● ● ●● ●● ● ● ● ● ●
● ● ● ● ●● ●● ●●
● ● ● ●
● ●● ● ● ●●● ●
● ● ●● ● ● ●● ● ● ●● ● ● ● ●
● ●● ● ● ●●● ● ●●●
● ● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ●● ● ● ● ● ● ● ● ●●● ● ● ● ●● ●●●
● ●● ● ● ● ●● ●●
● ● ●● ●● ●●●● ● ● ●
●
● ●
● ● ●●● ● ● ● ●
● ●
● ●
●
● ●
● ● ● ●
● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ● ●●● ● ●●
●● ● ●
● ● ●●
●● ● ●●
● ●● ● ● ● ● ●●
● ● ●
● ● ●
● ● ● ● ● ●● ●● ●● ●●
● ●● ●
● ● ● ● ●● ● ● ● ●●
●
● ●
●
●
● ●
● ●●
●● ●● ● ●●●
● ●●
● ● ● ● ● ● ● ● ● ●● ● ● ●●●
● ●
● ●● ● ● ● ● ● ● ● ● ●
● ● ●● ● ● ●●● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ●
● ●
● ● ● ● ●● ●● ●
● ● ●● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ●
● ● ● ●
● ● ● ● ● ● ● ●●
● ● ●● ● ● ● ● ●● ● ● ● ●
● ●
●● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ●
● ● ● ●● ● ●
● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ●
● ●
●
●
●
● ●
● ● ● ● ● ●
● ● ●
● ● ●
● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ●
●
● ● ● ●
● ● ● ● ● ● ● ●
● ● ● ● ● ●
● ●
● ● ● ● ● ● ● ●
● ● ● ●
● ● ● ●
● ●
●● ● ●
●
● ●
● ●
● ●
● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ●
●
● ● ● ● ● ●
● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ●
● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ●●● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ●● ● ● ●● ● ●
● ● ● ● ● ● ● ● ●
● ● ● ● ●
● ● ● ●● ● ● ●
● ● ● ●
● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ●
●
●
●
● ●
● ● ● ● ●
● ●
●
● ● ● ●
● ● ●
●
●
●
● ●
● ● ● ● ● ● ●
●
● ●
●
●
● ● ● ● ● ●
●
● ● ●
● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ●
● ● ●
● ● ● ● ● ● ● ● ● ●●
● ● ● ● ● ● ● ● ● ● ●
●
● ● ● ● ● ● ● ●● ● ● ● ●
● ● ● ●
● ● ● ● ●
● ● ● ●● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ●
●● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●
●● ● ● ● ● ● ● ● ● ●
● ●
● ● ● ● ● ● ● ● ● ●
● ● ●
● ●
300
400
500
600
700
0 200 400 600
t (hours)
km /s
● ●Actual Prediction
(b) Hourly forecasts for period 2016-11-16 17:00 to 2016-12-14 01:00
Figure 6: Solar Wind Prediction: reconstructed time series predictions

heading: APPENDIX E NEURAL NETWORK ARCHITECTURE DETAILS

Review summaries:
The paper is motivated by a specific application about the prediction of solar wind. It is a major threat to the modern world. Authors propose a new regression model called Dynamic Time-Lag Regression which is a special case for spatio-temporal data. For solving this regression problem, the authors use a standard Bayesian approach to fit the model. The analysis of saddle point in Section 3.2 is quite impressive and meaningful. Correspondingly, the practical DTLR algorithm iteratively updates the parameters and hyper-parameters. In experiments, it reveals that the proposed method improves the performance for this particular problem significantly. The application solved by authors is very important, and it is improved by the proposed method. However, the novelty of the work is not presented clearly.