While analysing time series data, conventional linear regression methods leave us the choice with one of two assumptions: either the dependent and independent variables are synchronous or fixed time lags are assumed. However the time lag might be varying over time, or noisy. The time lag becomes an additional source of error beside the usual (measurement) noise. Ignoring the variable time lag amounts to a specification error. We use optimal transport cost as a way to account for the time lag uncertainty and make ordinary least squares (OLS) robust to this type of error or noise. Using the chain rule, this enhancement to the conventional OLS regression model of noise easily generalises to multivariate regression. Introduction We investigate the simple regression for non-synchronous time series. This type of data is common in practice in a variety of contexts. We propose a method that is relevant when all the usual assumptions for OLS are fulfilled but there is still an incorrect model specification in the sense that the time lag variation is significant in comparison to the variation of the dependent variable. This applies to situations where the dependent and independent variables are connected with a transmission link that is not instantaneous and might even exhibit shocks in the transmission delay. Current methods for this situation include OLS with multiple time lags, but with time lag averaging to zero, the best regressor is one that assumes no time lag at all. Hence we compare the proposed method with simple OLS regression with no time lag assumption as a baseline. OT distance, and optimal transportation plan are robustly calculated even for a very large number of support points [9]. While fast solvers for the optimal transport (OT) cost are available, the corresponding gradient is still challenging to compute. One convenient approach is to use some flavour of automatic differentiation [7, 8]. In practice this approach proves numerically unstable as the number of support points in the distribution increases beyond a few dozens. As an alternative, we use the transportation plan to derive an OT gradient and the corresponding chain rule. Starting with simple regression as an example and then generalising to functions only required to be continuously differentiable in all their parameters. First, we review OLS regression, Optimal Transport and the proposed Least-Square Optimal-Transport (LS-OT) regression in Section 1. Section 2 presents the experimentation methodology with synthetic data and we compare OLS and LS-OT regression results in Section 3. We conclude with a discussion of the results, limitations and further research. Section 1 1.1 OLS regression We first recall OLS basic assumptions[2| given observed data point (Y,X): - The dependent variable is modelled as linear function of the independent variables: y = slope * X - The observations (Xi,Yi) are independent and identically distributed (iid) - The non-dependent variables are linearly independent - residuals eta = Y-y are normal, uncorrelated with bounded variance The OLS [2] method aims at inferring the dependent variable y t from the contemporary independent variable xt and time lags thereof x t-1 , x t-2 , etc. The linear dependency factor slope, which is identified through an optimization process that minimises OLS cost defined as follows:𝑂𝐿𝑆𝑐𝑜𝑠𝑡= Σ𝑡(𝑌𝑡−𝑠𝑙𝑜𝑝𝑒 × 𝑋𝑡)² Regression fit can be evaluated using the coefficient of determination r-squared which is related to the ratio of unexplained over total variance:𝑅2 = 1 − Σ𝑡(𝑌𝑡−𝑦𝑡)²/ Σ𝑡(𝑌𝑡−𝐸(𝑌𝑡))² R2 = 1 in case of a perfect fit and R2 = 0 when none of the variance is explained by the model. 1.2 Optimal Transport Let’s consider the Optimal Transport cost [1]: for which the optimal plan for which T * (M)=V, M𝑂𝑇𝑐𝑜𝑠𝑡=Σ 𝐶𝑜𝑠𝑡(𝑋,𝑇(𝑋)) M and V are measures defined on the support set of points {X}. T * ( M ) is the optimal transportation plan that moves the mass of M into V at a minimal cost given the Cost matrix C. In comparison, the OLS uses the Euclidean distance as cost and sets T to be a linear function parameterized with Slope: T(X) = Slope * X. We observe that this is identical to the Wasserstein distance when the transport cost is defined by the Euclidean distance except that no error is allowed in the time axis which can be seen as an infinite cost. Section 2 2.1 Least-Squares Optimal-Transport regression Now in the OT distance minimization framework, we have the flexibility to allow a time lag and handle data that is not perfectly synchronous. We state the regression problem as an optimization problem combining the OLS and OT costs, with a trade-off parameter lambda:𝑎𝑟𝑔𝑚𝑖𝑛(𝑠𝑙𝑜𝑝𝑒) : 𝑂𝐿𝑆𝑐𝑜𝑠𝑡(𝑠𝑙𝑜𝑝𝑒) + 𝑙𝑎𝑚𝑏𝑑𝑎 × 𝑂𝑇𝑐𝑜𝑠𝑡(𝑠𝑙𝑜𝑝𝑒) Where the OT cost calculation is parameterized with an isotropic cost matrix. As will be shown in Section 1.4, the OT cost matrix can be parameterized to include the OLS cost and the trade-off between OLS and OT costs. This simplifies the formulation to:𝑎𝑟𝑔𝑚𝑖𝑛(𝑠𝑙𝑜𝑝𝑒) : 𝑐𝑠𝑡 × 𝑂𝑇𝑐𝑜𝑠𝑡(𝑠𝑙𝑜𝑝𝑒) 2.1 Gradient decomposition and chain rule We start with an univariate regression. Given an input set of points (X,Y) and initial estimated parameters for the slope and intercept also represented by a set of points (x,y). We want to find the regression parameters that minimises the OT distance between the supports (X,Y) and (x,y). We obtain the transportation plan from the backAndForth solver [source back and forth]. It is composed of 2 dimensions ( .𝑒₁,𝑒₂) We use a chain rule to back-propagate the OT distance gradient through both dimensions. Then we apply the chain rule for each dimension separately . The direction is the familiar way of understanding the gradient of a 1 argument𝑒₂ value function, that is along the y-axis. Considering a transportation plan in the𝑒₂ direction (ie using a non-isotropic transport cost where displacements along are𝑒₁ prohibitively elevated) amounts to a conventional OLS regression. Then, the direction is all about displacements along the x-axis. This requires a𝑒₁ slightly different handling as will be shown starting with the following example. Figure 1: OT plan from the estimated solution to the data to fit 2.1.1 Linear regression with parameters slope and intercept We start with 2 sets of points: the starting support points defined by the pairs (x,y) generated as per the parameters The destination support points defined by the pairs (X,Y)𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡, 𝑏𝑒𝑠𝑡. which is the given data and for which we search the true regression parameters .𝑠𝑙𝑜𝑝𝑒, 𝑏 The true relationship is defined as:𝑌 = 𝑠𝑙𝑜𝑝𝑒𝑡𝑟𝑢𝑒 × 𝑋 + 𝑏𝑡𝑟𝑢𝑒 Current estimates of the parameters give us𝑦 = 𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 × 𝑥 + 𝑏𝑒𝑠𝑡 The transportation plan is defined in 2 dimensions, along the x-axis ( ) and along the y-axis𝑒1 ( ). When looking at displacements along the x-axis, we index the y-axis according to .𝑒2𝑌𝑗 Thus defining the x-axis coordinate as the first member of the pair as per pair’s(.,𝑦𝑗)𝑒𝑠𝑡 relationship definition. And similarly when looking at displacements along the y-axis. From here on we drop the subscript for simplicity._𝑡𝑟𝑢𝑒 Positive displacements along y-axis Negative displacements along y-axis Figure 2: OT transportation plan positive and negative displacements 2.1.1.1 Parameter inference with regression Visual inspection of the displacements suggests that their amplitude can be regressed to find the gradient to the slope and to the intercept. We observe that for the displacements along y-axis ( direction ) , the terms Δ slope and Δ b can be𝑒2 obtained by linear regression by solving:Δ𝑦 = Δ𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 × 𝑥 + Δ𝑏𝑒𝑠𝑡 Where at each is given by the transportation plan.Δ𝑦𝑗𝑥𝑗 Along the x-axis ( direction ) we have:𝑒1Δ𝑥 = Δ𝑖𝑛𝑣𝑒𝑟𝑠𝑒𝑆𝑙𝑜𝑝𝑒 × 𝑦 + Δ𝑖𝑛𝑣𝑒𝑟𝑠𝑒𝐵 Here, and are given by the transportation plan and the previous estimate or Y,Δ𝑥𝑦 respectively . Again, using a conventional regression we obtain the parametersΔ𝑖𝑛𝑣𝑒𝑟𝑠𝑒𝑆𝑙𝑜𝑝𝑒, Δ𝑖𝑛𝑣𝑒𝑟𝑠𝑒𝐵 Now, using the inverse slope (obtained at previous step) , the newly1/𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 calculated optima slope is then:𝑖𝑛𝑣𝑒𝑟𝑠𝑒𝑆𝑙𝑜𝑝𝑒𝑛𝑒𝑤 𝑒𝑠𝑡 =(1/𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡+Δ𝑖𝑛𝑣𝑒𝑟𝑠𝑒𝑆𝑙𝑜𝑝𝑒) 𝑠𝑙𝑜𝑝𝑒𝑛𝑒𝑤 𝑒𝑠𝑡 =1/𝑖𝑛𝑣𝑒𝑟𝑠𝑒𝑆𝑙𝑜𝑝𝑒𝑛𝑒𝑤 𝑒𝑠𝑡 Δ𝑠𝑙𝑜𝑝𝑒 = 𝑠𝑙𝑜𝑝𝑒𝑛𝑒𝑤 𝑒𝑠𝑡−𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 = 1/(1/𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡+Δ𝑖𝑛𝑣𝑒𝑟𝑠𝑒𝑆𝑙𝑜𝑝𝑒) − 𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 And the newly calculated optimal b is thenΔ𝑏= −(1/𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡+Δ𝑖𝑛𝑣𝑒𝑟𝑠𝑒𝑆𝑙𝑜𝑝𝑒) *Δ𝑖𝑛𝑣𝑒𝑟𝑠𝑒𝐵 2.1.1.2 Parameter inference details Slope parameter Given the transportation plan in the direction is ∆Y, from y to Y𝑒₂𝑌−𝑦 = 𝑠𝑙𝑜𝑝𝑒 × 𝑋 + 𝑏 −(𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 × 𝑥 + 𝑏𝑒𝑠𝑡) 𝑌−𝑦 = (𝑠𝑙𝑜𝑝𝑒 −𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡) × 𝑋 + (𝑏 −𝑏𝑒𝑠𝑡) If we ignore the b parameter this amounts to that we estimate aΔ𝑠𝑙𝑜𝑝𝑒 = Δ𝑌 / 𝑋𝑎𝑟𝑔𝑚𝑖𝑛(∆𝑠𝑙𝑜𝑝𝑒) : 𝜮(Δ𝑠𝑙𝑜𝑝𝑒 × 𝑋− ∆𝑌)² Given the transportation plan in the direction , is ∆X, from x to X. We now have𝑒₁𝑌 = 𝑠𝑙𝑜𝑝𝑒 × 𝑋 + 𝑏 =𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 × 𝑥 + 𝑏𝑒𝑠𝑡 With the direction : y is unchanged. This part of the transportation plan describes𝑒₁ only changes along the x-axis.𝑠𝑙𝑜𝑝𝑒 =(𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 × 𝑥 + 𝑏𝑒𝑠𝑡−𝑏)/𝑋 𝑠𝑙𝑜𝑝𝑒 − 𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 =(𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 × 𝑥 + 𝑏𝑒𝑠𝑡−𝑏)/𝑋 −𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 ∆𝑠𝑙𝑜𝑝𝑒=(𝑦−𝑏)/𝑋 −𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 Which we estimate as:𝑎𝑟𝑔𝑚𝑖𝑛(∆𝑠𝑙𝑜𝑝𝑒) : 𝜮(Δ𝑠𝑙𝑜𝑝𝑒 − (𝑦−𝑏)/𝑋+𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡)² Intercept parameter For simplicity, we assume that the slope is constant or frozen For the direction :𝑒₂𝑌−𝑦 = 𝑠𝑙𝑜𝑝𝑒 × 𝑋 + 𝑏 −(𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 × 𝑥 +𝑏𝑒𝑠𝑡) =∆𝑏 Ignoring the slope parameter , we need to solve:𝑎𝑟𝑔𝑚𝑖𝑛(∆𝑏) : 𝜮(Δ𝑏 − (𝑌−𝑦))² For the direction :𝑒₁ Given the x-axis displacements and no change on the y-axisΔ𝑋 = 𝑋−𝑥𝑌 = 𝑠𝑙𝑜𝑝𝑒 × 𝑋 + 𝑏=𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 × 𝑥 +𝑏𝑒𝑠𝑡 = 𝑦 𝑏−𝑏𝑒𝑠𝑡=𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 × 𝑥−𝑠𝑙𝑜𝑝𝑒 × 𝑋 Here we use our current estimate of slope and the previous estimate𝑏−𝑏𝑒𝑠𝑡=𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 × 𝑥−𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 𝑛𝑒𝑤 × 𝑋 Δ𝑏=𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 × 𝑥−𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 𝑛𝑒𝑤 × 𝑋 Which is estimated as𝑎𝑟𝑔𝑚𝑖𝑛(∆𝑏) : 𝜮(∆𝑏−(𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 × 𝑥𝑗−𝑠𝑙𝑜𝑝𝑒𝑒𝑠𝑡 𝑛𝑒𝑤 × 𝑋𝑗))² Figure 3: left: initial y. middle: Y. right: solution for y overlaid on Y. 2.1.2 Continuously derivable function with one argument We define a 1 time continuously derivable function parameterized by θ and accepting one argument x: and its derivative wrt θ𝑦 = 𝑓(θ, 𝑥) δ𝑦/δθ = δ𝑓(θ, 𝑥)/δθ We start with input observations (X,Y) which we assume are governed by the relationship with some unknown ϑ𝑌 = 𝑓(ϑ, 𝑋) Given the transportation plan in the direction is𝑒₂∆𝑌 = 𝑌 − 𝑦Δ𝑌 = δ𝑓(θ, 𝑥)/δθ × Δθ Hence we need to solve:𝑎𝑟𝑔𝑚𝑖𝑛(∆θ) : 𝜮(Δ𝑌 − δ𝑓(θ, 𝑥)/δθ × Δθ)² Given the transportation plan in the direction , is ∆X, from x to X𝑒₁𝑌 = 𝑓(ϑ, 𝑋)=𝑓(θ, 𝑥)= 𝑦 Using the Taylor rule for𝑓(ϑ, 𝑋)=𝑓(θ, 𝑋) + δ𝑓(θ, 𝑋)/δθ × Δθ We rewrite𝑌=𝑓(θ, 𝑋) + δ𝑓(θ, 𝑋)/δθ × Δθ=𝑓(θ, 𝑥) Here we need to solve: , where𝑎𝑟𝑔𝑚𝑖𝑛(∆θ) : 𝜮(𝑓(θ, 𝑋) + δ𝑓(θ, 𝑥)/δθ × Δθ−𝑦)²𝑓(θ, 𝑥) = 𝑦 1.4.2.1 Multiple argument functions The generalisation of this chain rule for OT distance transport plan back-propagation to continuously derivable functions with several arguments (or in multiple dimensions) is straightforward if all argument variables share the same innovation in time delays. Turning to a function of multiple arguments with non-synchronous time delays: where each argument has its own time delay innovation𝑦 = 𝑓(θ, 𝑥1, 𝑥2, ..., 𝑥𝑁)𝑥𝑗 variable . For this case we use the Sinkhorn barycenter approach [4] where weℇ𝑗 calculate the gradient in turn for each of the input argument and apply the resulting gradient step until we converge to the optimal parameter , as close as possible to ϑ. Section 3 3.1 Synthetic data generator In the case of OLS the input data is assumed to be synchronous. The time lag is assumed to be fixed. We now consider a data generation process with explicit random time lag: y t = beta * X(t+eta t ) + eta Where X t is a vector: [x t , 1], beta is a constant vector: [slope, intercept] Eta is normal with variance s_e. The time lag noise eta t is normal with variance s_e t Equipped with this data generator, we proceed to estimate beta with different methods and compare results. 3.2 Implementation details For the first experiment, use python jax [3] automatic differentiation software for the OT cost gradients. The chain rule allows to backpropagate the gradient to regression model parameters, then gradient descent converges to the desired minima. We compensate for the bias in the OT distance calculation as in [4]. We run 100 simulations and report estimated slope parameters using the OLS and the Transport regression. (The complete jupyter notebook for data generation and regression experiments presented here is available at: [ https://colab.research.google.com/drive/1MH_vGrtf-cuPKtxBQbokLsDRRrXr5ero ]) In a second experiment, we use the back and forth OT solver together with the OT gradient chain rule presented in section 1.4 [6]. This solver calculates sparse plans and scales to thousands of points in the source and destination support point sets. We run 50 simulations and report the slope parameters using OLS and this second method for the gradient. (The complete jupyter notebook for data generation and regression experiments presented here is available at: [ https://colab.research.google.com/drive/1GhiKSw8g7rqstvEfQ0eOjIMMQG4sI xjM ]) In both experiments we allow only positive time lags (i.e. positive time lag noise eta t ). Section 4 4.1 Results As shown in Figure 1, automatic differentiation method, brings significant improvement on OLS alone for the regression task with a time varying lag (where the lag standard deviation is similar to several regular time steps between samples). We observe that the improvement for most samples except in the few cases where the slope estimate overshoots the true value. In most of the 100 experiments, both estimators undershoot the true value but LS-OT is robustly more accurate. Figure 4: ls shows OLS estimated Beta, LS-OT fit using automatic differentiation, and the true slope (generator parameter) The results of the second method are shown in Figure 2, using OT transportation plan based gradient. Here, the second method brings significant improvement on OLS alone for the regression task with a time varying lag (where the lag standard deviation is similar to one regular time step between samples). Figure 5: OLS estimated Beta, LS-OT fit using OT transportation plan and chain rule as gradient, and the true slope (data generator parameter) We observe that the LS-OT method is a favourable option when non-synchronous data is at hand and the preferred model when expert knowledge on the time lag variance is available. 4.2 Discussion We presented the LS-OT solver for regressions where the dependent variable gets information from the independent variable with a variable delay, modelized by a gaussian variable. This is relevant in many practical settings where multivariate analysis includes several time series of connected but disjoint phenomena in biology, mechanics, finance and economics, etc. We presented a solver that scales to thousands of points in the source and destination supports, while the automatic gradient based method exhibits numerical stability issues already with 100 points in the support. This method can easily be extended to partial least square, autoregressive models that exhibit variable time lags and vector variants thereof. Finally, we observe that in addition to informing regression parameters, the transportation plan gives cues about timing relationships among the variables, which is an avenue for further research. 4.3 Further research The automatic differentiation method is severely limited by the OT accuracy and numerical stability. It can be applied for very short times eries with 10-20 samples because of the limited accuracy of the gradient descent OT solver. The second method with OT transportation plan based gradient is much more scalable with thousands of samples and robust numerical stability and accuracy as shown in the second experiment. It is a practical implementation suitable for real life regression problems. It combines OLS parameter inference and dynamic time warping [10]. A graphical observation of the transport plan with synthetic data reveals that the dependent variable is always delayed (on the x-axis) with our synthetic data. Statistical mean and standard deviation of the components of the optimal transportation plan gives an estimate of the delay generator parameters. This indicates that the LS-OT method gives cues about time based precedence or causality among the variables. The ability of LS-OT to account for timing suggests novel analysis of lead-lag relationships. We leave these aspects for future research. In addition, we may apply the singular value decomposition (SVD) on the transportation plan to get the main (largest) eigenvalues and corresponding eigenvectors as a rough characterisation of the time lag. A single mode (of fast decaying eigenvalues) SVD would suggest that it is better to revert to OLS with the appropriate fixed time lag while multiple significant eigenvalues means that the time lag is arguably not fixed and diverse.