Title: DYNAMIC TIME LAG REGRESSION: PREDICTING WHAT & WHEN

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.

1. INTRODUCTION:

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. 1A 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 L 1 (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 time lag between an input and its effect, the solar wind recorded at L 1 , 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 2 PROBABILISTIC DYNAMICALLY DELAYED REGRESSION

2.1. POSITION OF THE PROBLEM:

Given two time series, the cause series x(t) (x(t) ∈ X ⊂ R D ) 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 .

2.2. PROBABILISTIC DYNAMIC TIME-LAG REGRESSION:

For practical reasons, cause and effect series are sampled at constant rate, and thereafter noted x t and y t with t in N . Accordingly, mapping g maps x t onto a finite set T of possible time lags, with T = {∆t min . . . , ∆t max } 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 ŷ(x) = {ŷ 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 ŷi (x), where the mixture weights are defined from p(x). More formally, letting y t denote the vector of random variables {y t+i , i ∈ T }:

EQUATION

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 x t drives y t+i (τ i = 1) or not (τ i = 0). The assumption that every cause has a single effect is modelled by imposing:4 

EQUATION

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 pi (x) stands for the probability of the effect of x t = x to occur with delay i.

The second simplifying assumption is that the variance σ 2 i (τ ) of predictor ŷ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 material). The fact that x can influence y i through predictor ŷ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 = j corresponds to a small auto-correlation time of the effect series.

2.3. LEARNING CRITERION:

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 θ = (ŷ, p, σ, α). From Eq. (3) the conditional probability q i (x, y) def = P (τ i = 1|x, y) reads:

EQUATION

with normalization constant

Z(x, y|θ) = i∈T pi (x) exp - 1 2σ 2 j∈T α ji y j -ŷj (x) 2 + 1 2 j∈T log(1 + α ji ) .

The log-likelihood then reads (intermediate calculations in supplementary material, appendix A):

EQUATION

where E data 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 E data . The hyper-parameters σ and matrix α of the model are obtained by optimizing L:

EQUATION

In addition the optimal ŷ and p reads:

EQUATION

where the above conditional empirical averaging operates as an averaging over samples close to x.

These are self-consistent equations, since q i (x, y) depends on the parameters σ 2 and α ij , ŷ 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 ŷi and pi ; and predictors ŷi and mixture weights pi are updated according to Eqs ( 8) and ( 9) respectively.

3. THEORETICAL ANALYSIS:

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

3.1. LOSS FUNCTION AND RELATED OPTIMAL PREDICTOR:

Let us assume that the hyper-parameters of the model have been identified together with predictors ŷi (x) and weights pi (x). These are leveraged to achieve the prediction of the effect series. For any given input x, the sought eventual predictor is expressed as (ŷ(x), Î(x)) where Î(x) is the predicted time lag and ŷ(x) the predicted value. The associated L 2 loss is:

EQUATION

)

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) = ŷI (x) (x) with I (x) = arg max i pi (x) ,

Proof. In supplementary material, Appendix C.

3.2. LINEAR STABILITY ANALYSIS:

The saddle point (Eqs 7, 5, 8, 9) admits among others a degenerate solution, corresponding to pi (x) = 1/|T |, α ij = 0 for all pairs (i, j), with σ 2 = σ 2 0 . Informally the model converges toward this degenerate trivial solution when there is not enough information to build specialized predictors ŷi .

Let us denote ∆y:

2 i (x) = y i -ŷi (x)

2 the square error made by predictor ŷi for x, and

σ 2 0 = 1 |T | E data i∈T ∆y 2 i (x)

the average of MSE over the set of the predictors ŷ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 ŷ and p and two hyper-parameters α and r = σ 2 /σ 2 0 . After the computation of the Hessian (in supplementary material, Appendix B) the system involves three key statistical quantities, two global ones:

EQUATION

EQUATION

)

and a local |T |-vector of components

u i [x, q] = 1 σ 2 0 E data q i (x, y) ∆y 2 i (x) - j∈T q j (x, y)∆y 2 j (x) x .

Up to a constant, C 1 represents the covariance between the latent variables {τ i } and the normalized predictor errors. C 1 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 = q 0 uniform, C 1 [q 0 ] = 1 and C 2 [q 0 ] represents the default variability among the prediction errors. u i [x, q] informally measures the quality of predictor ŷi relatively to the other ones at x. More precisely, a negative value of u i [x, q] indicates that ŷi is doing better than average in the neighborhood of x.

At a saddle point the parameters are given by:

σ 2 σ 2 0 = |T | -C 1 [q] |T | -1 and α = |T | |T | -1 1 -C 1 [q] C 1 [q] .

The predictors ŷ 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

C 2 [q] < 2C 2 1 [q] + O 1 |T | .

In particular, the degenerate solution is unstable if

C 2 [q 0 ] > 2 1 - 1 |T | .

Note that for ∆y i (x) iid centered with variance σ 2 0 and relative kurtosis κ (conditionally to x) one has C 2 = (2 + κ)(1 -1/|T |). Therefore, whenever ∆y 2 i (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 C 2 [q 0 ] is non-zero, due to the fact that the gradient points in the opposite direction to u(x) (with dp(x) ∝ -C 2 [q 0 ]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:

pi (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.

4. THE DTLR ALGORITHM:

The DTLR algorithm learns both regression models ŷ(x) and p(x) from series x t and y t , 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 ←-σ 2 0 |T |-C1[q] |T |-1 ; α ←-|T | |T |-1 1-C1[q]

C1[q] ; end Result: Model parameters θ = {ŷ, 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 ŷ(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 ŷ and p models, each with |T | output neurons, with |T | the size of the chosen domain for the time lag.

5. EXPERIMENTAL SETTING:

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 (x t ), I (x)). DTLR is compared to the natural baseline defined as the regressor with constant time lag, ŷ ∆(x t ), 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 (x t ) 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 C 1 (section 3.2), compared to the degenerate solution pi (x) = 1/|T |. For C 1 < 1, the model accurately specializes the found predictors pi .

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 x t ∈ R 10 and the effect signal y t are generated as follows (with η = 0.02, s 2 = 0.7):

EQUATION

with time-lag mapping g(x t ) 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 √ v 2 t +2ad ( √ v 2 m +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 x t 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 Facility5 . 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 y i 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 

6. EMPIRICAL VALIDATION:

Table 1 summarizes the DTLR performance on the synthetic and solar wind problems (detailed results are provided in the appendix). 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).

On Problem II, DTLR accurately learns the inverse relationship between x t , g(x t ) and y t 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.

6.2. THE SOLAR WIND PROBLEM:

DTLR finds an operational solar wind model ( 

7. DISCUSSION AND PERSPECTIVES:

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 x t , 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.

Let θ def = (ŷ, p, σ, α) denote the parameters of the model and consider the probability that predictor ŷi is the good one conditionally to a pair of observation (x, y):

q i (x, y) = P (τ i = 1|x, y) = 1 Z(x, y|θ) pi (x) exp - 1 2σ 2 j∈T α ji y j -ŷj (x) 2 + 1 2 j∈T log(1 + α ji ) with Z(x, y|θ) = i∈T pi (x) exp - 1 2σ 2 j∈T α ji y j -ŷj (x) 2 + 1 2 j∈T log(1 + α ji ) .

This gives immediately

L[{(x, y)} data |θ] = -|T | log(σ) -E data i∈T 1 2σ 2 y i -ŷ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 N x,y becomes large. Restricting the likelihood to this subset of the data yields the following:

L x,y = Nx,y m=1 {τ (m) } p(τ (m) |x) i∈T √ 2π σ i (τ (m) ) exp - 1 2 i∈T y i -ŷi (x) 2 σ i (τ (m) ) 2 .

Upon introducing the relative frequencies:

q i (x, y) = 1 N x,y Nx,y m=1 τ (m) i satisfying i∈T q i (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 L x,y q exp -N x,y F x,y q where the rate function reads

F x,y q = |T | log(σ) + i∈T y i -ŷi (x) 2 1 + j∈T α ij q j 2σ 2 - 1 2 q i j∈T log(1 + α ji ) + q i log q i pi .

Taking the saddle point for q i 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 |θ] = -E data F x,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 ŷ yields:

∂L ∂ ŷi (x) = 1 σ 2 E data y i -ŷi (x) 1 + j∈T α ij q j (x, y) x .

Then the optimization of p gives:

∂L ∂ pi (x) = E data q i (x, y) pi (x) -λ(x) x , = 1 pi (x) E data q i (x, y) x -λ(x)

with λ(x) a Lagrange multiplier to insure that i pi (x) = 1 This gives

pi (x) = 1 λ(x) E data q i (x, y) x Hence i∈T pi (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 ) E data q j (x, y) - 1 2σ 2 E data y i -ŷi (x) 2 q j (x, y) .

APPENDIX B PROOF OF PROPOSITION 3.1

Given I(x) a candidate index function we associate the point-like measure

p i (x) = δ i,I(x) .

Written in terms of p the loss function reads

L 2 (ŷ, p) = E x,y i∈T p i (x) y i -ŷ(x) 2 . Under (3) (with α ij = αδ ij ) the loss is equal to L 2 (ŷ, p) = E x i∈T p i (x) ŷi (x) -ŷ(x) 2 -pi (x) ασ 2 1 + α + σ 2

The minimization w.r.t. ŷ yields ŷ

EQUATION

)

In turn, as a function of p i the loss being a convex combination, its minimization yields

EQUATION

EQUATION

)

Combining these equations (16,17,18) we get

I(x) = arg max i∈T pi (x) ,

which concludes the proof.

APPENDIX C STABILITY ANALYSIS:

The analysis is restricted for simplicity to the case α ij = αδ ij . The log likelihood as a function of r = σ 2 /σ 2 0 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β) + E data log(Z) -λ(x) i∈T pi (x) with Z = i pi (x) exp - β 2σ 2 0 ∆y 2 i (x) ,

and where λ(x) is a Lagrange multiplier which has been added to impose the normalization of p.

The gradient reads

∂L ∂r = 1 2r 2 |T |(1 -r) + βr 2 1 + βr , ∂L ∂β = r 2(1 + rβ) - 1 2 C 1 [q], ∂L ∂ ŷi (x) = 1 σ 2 E data y i -ŷi (x) 1 + αq i (x, y) x . ∂L ∂ pi (x) = E data q i (x, y)|x pi (x) -λ(x), with C 1 [q] = 1 σ 2 0 E data i∈T q i (x, y)∆y 2 i (x) ,

This leads to the following relation at the saddle point:

r = |T | -C 1 [q] |T | -1 , α = |T | |T | -1 1 -C 1 [q] C 1 [q] , ŷi (x) = E data y i 1 + αq i (x, y) x E data 1 + αq i (x, y) x pi (x) = E data q i (x, y)|x .

Let us now compute the Hessian. It is easy to see that the block corresponding to the predictors ŷ decouples from the rest as soon as these predictors are centered.

Denoting C 2 [q] = 1 σ 4 0 E data i∈T q i (x, y) ∆y 2 i (x) - n j=1 q j (x, y)∆y 2 j (x) 2 ,

we have

∂ 2 L ∂r 2 = 1 2r 2 -|T | + 2 |T | |T | -1 C 1 [q] -1 -β 2 C 2 1 [q] ∂ 2 L ∂r∂β = 1 2r 2 C 2 1 [q] ∂ 2 L ∂β 2 = 1 4 C 2 [q] -2C 2 1 [q] ∂ 2 L ∂ pi (x)∂ pj (x) = - E data q i (x, y)q j (x, y)|x pi (x)p j (x) ∂ 2 L ∂r∂ pi (x) = 0 ∂ 2 L ∂β∂ pi (x) = - u i [x, q] 2p i (x) , where u i [x, q] def = 1 σ 2 0 E data q i (x, y) ∆y 2 i (x) - j∈T q j (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 pi . The stability of the first one depends on the sign of C 2 [q] -2C 2 1 [q] for |T | large while the second block is always stable as being an average of the exterior product of the vector (q 1 (x, y)/p 1 (x), . . . , q |T | (x, y)/p |T | (x)) by itself. At the degenerate point α = 0, r = 1, pi = 1/|T | the Hessian simplifies as follows. Denote

dη = dre 1 + dβe 2 + dx i∈T dp i (x)e i+2 (x)

a given vector of perturbations, decomposed onto a set of unit tangent vectors, {e 1 and e 2 } being respectively associated to r and β, while e i (x) associated to pi (x) for all i ∈ T and x ∈ X . Denote

u = i∈T dxu i [x]e i (x) v(x) = i∈T e i (x) with C 2 = 1 |T |σ 4 0 E data i∈T ∆y 2 i (x) - 1 |T | j∈T ∆y 2 j (x) 2 . u i [x] = 1 σ 2 0 E data ∆y 2 i (x) -σ 2 0 |x .

With these notations the Hessian reads:

H = 1 2 -|T |e 1 e t 1 + e 1 e t 2 + e 2 e t 1 + C 2 2 -1 e 2 e t 2 -ue t 2 -e 2 u t -dxv(x)v t (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 η = η 1 e 1 + η 2 e 2 + η 3 u.

In this subspace the Hessian reads

H = 1 2       -|T | 1 0 1 C 2 2 -1 -M |T |C 2 0 -M |T |C 2 0      

, where M is the number of data points, resulting from the fact that

i∈T dxu i [x] 2 = M σ 4 0 E data i∈T ∆y 2 i (x) -σ 2 0 2 , = M C 2 ,

because E data (•|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 |(C 2 /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η ∝ Λe 1 + Λ(|T | + Λ)e 2 -M |T |(|T | + Λ)C 2 u,

which corresponds to increasing r and α while decreasing for each x the pi having the highest mean relative error

u i [x].

Concerning solutions for which pi

(x) = δ i Î(x)

is concentrated on some index Î(x), the analysis is more complex. In that case C 2 [p] = 0 and C 1 [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 E data i∈T ∆y 2 i (x) - α 2σ 2 E data ∆y 2 I(x) (x)

so we get the following optimal solution

σ 2 = 1 |T | E data i∈T ∆y 2 i (x) , 1 1 + α = E data ∆y 2 I(x) (x) σ 2 , I(x) = arg min i∈T E data ∆y 2 i (x)|x .

Figure 2 : Geometry of the CSSS Model (Zhao and Hoeksema, 1995) . The spherical surface with radius R cp is the cusp surface, representing the locus of cusp points of helmet streamers and the location of the inner corona. The source surface (radius R ss ) is free to be placed anywhere outside the cusp surface. Figure reproduced from Zhao and Hoeksema (1995) .

APPENDIX D EXPERIMENTS: ADDITIONAL DETAILS:

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; Odstrčil et al., 1996; Odstrčil and Pizzo, 1999a; b; Odstrčil, 2003; Odstrčil et al., 2004) of NOAA/SWPC. Mathematically, FTE is defined as:

EQUATION

where, B r (phot) and B r (ss) are the radial component of magnetic fields at the photosphere and the source surface, and R and R ss 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, 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.2 SOLAR WIND PREDICTION q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 300 400 500 600 0 200 400 600 t (hours) km/s q q

Actual Prediction (a) Hourly forecasts for period 2008-11-20 07:00 to 2008-12-17 14:00 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 300 400 500 600 700 0 200 400 600 t (hours) km/s q q

Actual Prediction (b) Hourly forecasts for period 2016-11-16 17:00 to 2016-12-14 01:00 

:

The adverse impact of space weather is estimated to cost

to 400 million USD per year, but can sporadically lead to much larger losses.

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.

Note that pre-processing can be used if needed to map non-Gaussian data onto Gaussian data.

Note however that the cause-effect correspondence might be many-to-one, with an effect depending on several causes.

https://omniweb.gsfc.nasa.gov

The Sun completes a rotation (or Carrington rotation) in approximately 27 days.

The predicted values, every 6 hours, are interpolated for comparison with the hourly measured solar wind.

In the solar wind problem, the training and test data sizes correspond to one cross-validation split

Figures and Tables:

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.

Figure 3: Problem II, Results

Figure 6: Solar Wind Prediction: reconstructed time series predictions

DTLR performance: accuracy (MAE and RMSE, the lower the better; Pearson, the higher the better) and stability σ 0 and C 1 (the lower the better). For each indicator, is reported the DTLR value (9-fold CV), the baseline value and the time-lag error.

Table content: <table><tr><td>Problem</td><td>M.A.E</td><td>R.M.S.E</td><td>Pearson Corr.</td><td>σ 0</td><td>C 1</td></tr><tr><td>I</td><td colspan="3">8.82 / 21.79 / 0.021 12.35 / 28.79 / 0.26 0.98 / 0.87 / -</td><td>29.8</td><td>0.14</td></tr><tr><td>II</td><td>10.15 / 27.40 / 0.4</td><td colspan="4">13.70 / 35.11 / 0.67 0.95 / 0.73 / 0.70 26.83 0.16</td></tr><tr><td>III</td><td>3.17 / 11.01 / 0.17</td><td>4.63 / 14.99 / 0.42</td><td colspan="3">0.98 / 0.79 / 0.84 11.84 0.09</td></tr><tr><td>IV</td><td>3.88 / 12.28 / 0.34</td><td>5.33 / 15.89 / 0.64</td><td>0.98 /0.79/ 0.81</td><td colspan="2">12.18 0.13</td></tr><tr><td colspan="2">Solar Wind 56.35 / 66.45 / -</td><td>74.20 / 84.53 / -</td><td>0.6 / 0.41 / -</td><td colspan="2">76.46 0.89</td></tr><tr><td colspan="2">6.1 SYNTHETIC PROBLEMS</td><td/><td/><td/><td/></tr></table>

, though the significantly higher difficulty of the solar wind problem is witnessed by the C 1 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 wind 7 and the measured solar wind.

Table content: <table><tr><td>Model</td><td colspan="2">M.A.E R.M.S.E</td><td/><td/><td/></tr><tr><td>WS</td><td>74.09</td><td>85.27</td><td/><td>600</td><td/></tr><tr><td>DCHB</td><td>83.83</td><td>103.43</td><td/><td/><td/></tr><tr><td>WSA</td><td>68.54</td><td>82.62</td><td/><td/><td/></tr><tr><td colspan="2">Ensemble Median (WS) Ensemble Median (DCHB) 78.27 71.52</td><td>83.36 100.04</td><td>v ^ (km/s)</td><td>500</td><td/></tr><tr><td>Ensemble Median (WSA)</td><td>62.24</td><td>74.86</td><td/><td/><td/></tr><tr><td>Persistence (4 days)</td><td colspan="2">130.48 161.99</td><td/><td>400</td><td/></tr><tr><td>Persistence (27 days)</td><td>66.54</td><td>78.86</td><td/><td/><td/></tr><tr><td>Fixed Lag Baseline</td><td>67.33</td><td>80.39</td><td/><td/><td/></tr><tr><td>DTLR</td><td>60.19</td><td>72.64</td><td/><td>300</td><td/></tr><tr><td colspan="2">(a) Comparative assessment on the Solar</td><td/><td/><td>300</td><td>400</td><td>500 v (km/s)</td><td>600</td><td>700</td></tr><tr><td colspan="2">Wind problem compared to the state of the</td><td/><td/><td/><td/></tr><tr><td>art Reiss et al. (2019, Table 1)</td><td/><td/><td/><td colspan="3">(b) Scatter Chart (9 fold CV)</td></tr></table>

Synthetic and Real-World Problems

Table content: <table><tr><td>Problem # train</td><td># test</td><td>d</td><td>|T |</td></tr><tr><td colspan="3">I 10, 000 2, 000 10</td><td>15</td></tr><tr><td colspan="3">II 10, 000 2, 000 10</td><td>20</td></tr><tr><td colspan="3">III 10, 000 2, 000 10</td><td>20</td></tr><tr><td colspan="3">IV 10, 000 2, 000 10</td><td>20</td></tr><tr><td colspan="4">Solar Wind 77, 367 2, 205 374 12</td></tr></table>

Cross validation splits used to evaluate DTLR on the solar wind forecasting task

Table content: <table><tr><td>Split Id Carrington Rotation Start</td><td>End</td></tr><tr><td>1</td><td/></tr></table>

Peer review (iclr2020):
This paper proposes a bayesian method, Dynamic Time Lag Regression (DTLR), to predict effect-data y from a cause x. While both, cause and effect, are time series (x = x(t) and y = y(t)), the cause is assumed to have influence on the effect via a time lag, y(x(t), t) = y(phi(t)) = y(t + g(x(t))). The novelity of the problem posed is that the time lag g(x(t)) is allowed to be non-stationary and supposed to be unknown. The authors outlie a a Bayesian approach. The likelihood of the cause-effect data is conidered used to optimize the regressors for a given, finite set of time lags and the probability of these time lags. This discussion is very brief and detailled derivations are given in the appendix. An accompanying stability analysis is performed. In the implementation a neural network-based approach is used, where in turns the parameters of the model are optimised by a typical gradient descent algorithm and afterwards the hyperparemeters are set to values determined by the parameters. The theoretical deviation is accompanied by practical expirements on synthetically generated data and (real) space weather data. The performance of the model is compared to the one of a prediction with a fixed time lag, which it in fact outperforms. On the space weather (solar wind) data the model is shown to outperform several state-of-the-art models. A codebase for an implementation done in the Scala language is given. The paper presented provides a meaningful contribution to the field of time series analysis and Bayesian methods. The novelty of the problem imposed as well as the accompanying solution is outlined very clearly. The problem is well justified and the motivation is clearly given. I recommend to accept it, even though a few details could have been described more precisely. 1. The definition of the time lag g(x(t)) could have been made more clearly, providing a few examples at this point could help to clarify these definitions here. 2. For the synthetic data generation a justification would be helpfull in order to motivate the choice of the functions f and v. 3. A chart with sample time series data (y) against the time would be helpful to better visualize the problem.