{'fraction_non_alphanumeric': 0.03456206058216703, 'fraction_numerical': 0.005453709179315329, 'mean_word_length': 4.5466472303206995, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modelling stochastic time delay for regression analysis

Systems with stochastic time delay between the input and output present a number of unique challenges. Time domain noise leads to irregular alignments, obfuscates relationships and attenuates inferred coefficients. To handle these challenges, we introduce a maximum likelihood regression model that regards stochastic time delay as an"error"in the time domain. For a certain subset of problems, by modelling both prediction and time errors it is possible to outperform traditional models. Through a simulated experiment of a univariate problem, we demonstrate results that significantly improve upon Ordinary Least Squares (OLS) regression.


Introduction
Consider the typical univariate regression problem where the intention is to find the relationship between x and y. When extended to time series, a number of time specific complexities arise. For the specific case where the input x affects the output y, with a random time delay in between, the estimated coefficients (or weights, predictions etc) are significantly attenuated [9]. This attenuation occurs both in traditional statistical models and machine learning models, and for certain problems can be a serious limitation. The following document proposes techniques for handling this class of time series regression.
As an example, let us take the management of blood sugar level in diabetes patients. In people with diabetes, the pancreas can't effectively regulate blood sugar levels. Therefore, these levels must be controlled by insulin injections and a special diet. The challenge for many people is that the relationship between the input (insulin) and output (blood sugar) is extremely complex [13]. The effect of the insulin injection might be observed after 15, 20 or t minutes depending on a number of factors, many of which are unknown. Due to the stochastic nature of the time delays, the actual effect can't be easily determined. It is difficult to differentiate the effect of the insulin injection from other factors and accurately determine how much one should take.
Inference for this type of problem is especially challenging. Typical regression models require a fixed alignment between cause and effect. Using standard methods, we'd need to assume that the effect occurs after some fixed time t which can be inferred from the data. However, if there is any uncertainty in the parameter t (t changes or is noisy) the resulting estimates will be significantly attenuated.
Consider the simple example in Figure 1 where the effect of the input is 1. The observed input is given by the red line, the blue line is when the effect actually occurs. The first effect happens one time point after the input. The second effect happens at the same time as the input. A fixed time delay isn't valid in this case because the time shifts differ.
If one were to model the effect using a fixed time delay and OLS, the estimate would only be half the true value because only one of the outputs is aligned. Obviously, this isn't ideal, we want the parameter estimates to be as close to the real values as possible, regardless of any noise in the lag structures. One can mitigate the problem via time aggregations, however in complicated cases with multiple factors, this is just not feasible.
Therefore, we propose a regression model which can handle stochastic time delay structures. We treat the stochastic time delay components as an 'error' in the time axis. Next, we find the maximum likelihood estimate, for a given set of parameters, considering the time error (t-axis) and the regression error (y-axis) simultaneously.
An open source Python implementation is available on GitHub.

Related Work
There is extensive existing literature on time series problems with time delay dependencies. In the statistical disciplines there a number of treatments which focus on fixed time delay dependencies. For example, Granger Causality [6] is used to determine whether one time series is useful for forecasting another, across some fixed delay. Distributed Lag [1] and Dynamic Regression models (e.g [8]) are able to handle linear and non-linear cause and effect relations that occur across multiple time lags. Nonetheless, these models assume a fixed time lag dependency.
In the machine learning literature, there are a number of models which can handle complex dependencies across time. Sequence models such as the RNN, LSTM and GRU [3] are able to generate predictions which incorporate time delays between input and the output variables. More recently, attention based models such as the Transformer have become the state-of-the-art for a variety of tasks, including time series forecasting [14]. Yet, all of these models have a tacit assumption of a fixed time lag dependency. There is also Dynamic Time Warping (DTW) [7] which tries to find the optimal alignment between time series sequences by minimizing the distance between the respective inputs. DTW is able to model varying time lag dependencies but is not purposed for regression, instead being mainly used for pattern matching and sequence alignment. The field of system identification also has considerable literature on dynamic non-linear systems. For example, the NARX [3] model can be used to identify non-linear systems with fixed time delays between input and output. At the same time, there are also treatments on systems with disturbances in the input [11], known as EIV systems. To our knowledge, these works do not deal with uncertainty in the time domain.
Finally, recent papers such as Dynamic Time Lag Regression [4] and Variable-Lag Granger Causality [2] deal with non-stationary time lag dependencies. In other words, the lag structure is assumed to evolve over time.
In this paper, we deal with the specific case of stochastic time delays i.e. time delays which vary randomly. There appears to be little existing research on this topic.

Methodology
We consider the problem as analogous to the typical error-in-variables (EIV) regression. Ordinary regression analyses (and machine learning models) define the loss function with respect to errors in the y axis only. For EIV, errors are considered in both the y-axis and the x-axis [5]. This is useful when there are measurement errors in the independent variable e.g. because the physical measurements have some degree of random error. Similarly, for this problem we assume that we have errors in the y-axis and the t-axis. That is, there are random prediction errors and random errors in the time domain.
3.1. Model Specification. Let x(t) be a discrete time series. We want to determine the functional relationship f : R −→ R to some other (target) variable y on observed data {y i } n i=1 , that is y = f (x(t)). The series x(t) must be stationary with a known value for the point at which f (x(t)) = 0. That is, we know at which point the input series has no effect on the output (when it's off). The size of the support of f (x(t)) should be small relative to number of points in the domain (we we will explore this condition further in 'Limitations'). We utilise the terminology impulse for an individual element of the support of f . Firstly, let us take the input series x(t) and decompose it into its constituent non-zero impulse components. So, if x(t) is a vector given by Given that each non-zero impulse (row) of X(t) is affected by a random time delay (denoted τ ), we then model: where X(t + T ) is the matrix of time shifted impulses, X(t) is the observed input, T is the set of single impulse time delays (τ 's), 1 is the vector (1, 1, · · · , 1) ∈ R 1×k , where k is the number of impulses and ε is a noise term.
The model definition represents the application of time shifts to the matrix rows and a subsequent reduction by summation over the columns.
We also assume that T is not a constant, but rather a random draw from some discrete distribution (e.g. discrete gaussian, poisson etc). For the discrete gaussian kernel T ∼ T (µ, σ) or for the poisson distribution T ∼ P ois(λ). Similarly one can model the errors as ε ∼ N (µ, σ ε ).

3.2.
Inference. In order to find the best estimate of f , we would like to find the function f which maximises the joint log-likelihood of the time-domain shifts and the prediction residuals. Specifically, we maximise: where θ f and θ T represent the parameters of the model f and time shifts T respectively. In other words, the L 2 term represents likelihood of time shifts and the L 1 term represents the likelihood of the prediction residuals. We maximise these terms simultaneously. For simplicity, we assume that the time shift distribution and error distributions are independent.
We will refer to this algorithm as Time Varying Stochastic (TVS) Regression.

Algorithm
Before optimisation, the values of each individual time shift τ are not known. In addition, the prediction errors ε can only be calculated if each value of τ is available (because for each time shift there is a different prediction and hence prediction error). Therefore, our algorithm finds the optimum set of time shifts in an inner optimisation (L 2 ), while iteratively searching for the optimum parameters of f in an outer optimisation loop (L 1 ). Firstly, we define a parametric form f (x; θ f ) and some initial parameters to be estimated for the function f . First, let us initialise some starting values for each of these parameters. Now, we want to find the best possible time shift τ for each input impulse in X(t). It stands to reason that the best possible time shift would be one that is not too distant from the observed impulse (i.e. has a high likelihood given some distribution) and also produces the best possible prediction. From there, the likelihood estimate is derived from the time shift and the prediction error. Finally, we iterate over a number of values of τ optimising until we maximise the likelihood for the specific impulse.
Remark 4.1. However, we must also consider that the impulses in X(t) are not independent from each other. After shifting, it's possible that two or more effects occur simultaneously. This is particularly problematic if there are multiple impulses within a short period of time or impulses have a distributed effect over multiple time points. As an example consider the series To accurately calculate the likelihood, we must optimise the time shifts simultaneously. Therefore, we treat the problem of finding the best set of time shifts as a discrete optimisation problem. For the optimisation step we utilise two assumptions.
(1) Firstly, smaller shifts are more likely than larger shifts (proportionate to the dispersion of the T distribution). The algorithm should explore the space of smaller shifts more often than larger shifts. (2) Secondly, impulses close to each other are more likely to be dependent than impulses further away. Accordingly, the optimisation procedure is: (1) Initialise the set of parameters θ f and θ T for the function f .  For the outer parameter optimisation (Step (3)), typical methods such as gradient descent, genetic algorithms and simulated annealing can be used. In our implementation, we have used the differential evolution algorithm [10] (scipiy.optimize [12]) because of its ability to handle noisy objective functions [10]. In order to improve convergence, we also standardize all input variables to the range (0,1).

Remark 4.3.
We also note that the accuracy of the final parameter estimate is relative to the ratio of the y-axis error and the effect size f (x(t)). As the ratio of noise σ ε to effect f (x(t)) increases, the time shift distribution T shrinks. In the event that the mean of θ T become 0, the model becomes a 'fixed lag' model, and the parameter estimate β tends to the standard linear regression coefficient. Therefore, the method provides no guarantee on recovering the exact time shifts, only that the coefficient estimates are equal to or better than their OLS counterparts.

Limitations
5.1. Scaling. As the length of x(t) increases, more impulses are introduced and the size of the decomposed matrix X(t) also increases. At some point, handling X(t) becomes impractical. To account for this, we assume that distant impulses do not affect each other. Concretely, if the i th and j th impulse (row) in X(t) are far away then their effect vectors are orthogonal. Hence, we can decompose the matrix X(t) into orthogonal blocks on which the inner optimisation can be run in parallel.

5.2.
Constraints. Finally, a note on problem constraints. The likelihood estimates are dependent on the inner optimisation procedure. There is no guarantee that the global maximum will be found, particularly for sequences with a high density of impulses. When the density (in time) of non-zero impulses in x(t) is high, the possible solution space of T grows exponentially. In such cases, the estimated likelihood is likely to be close to, but not exactly the same as the real value. Therefore, TVS Regression is most appropriate for sparse time series inputs.

Experiment
The following section demonstrates a simulated univariate example of TVS Regression. Figure 2 shows the simulated time series. The full code can be found via TVS Regression on GitHub. The input signal has 20 non-zero impulses which have been drawn from the standard normal distribution. The system is 'off' when x(t) = 0. In other words, when x(t) = 0, f (x(t)) = 0. The true shift distribution T is given by T ∼ P oisson(λ = 2). There is one value of τ for each of the 20 impulses. The green line in Figure 2 represents the shifted series, corresponding to the time at which the effect occurs. The blue line is the output y which includes a small amount of gaussian noise and an intercept. The values for the parameters were arbitrarily selected. The output y is defined by the following equation: y = 2(1 T X(t + T )) + 6.5 A histogram of the actual distribution of T is shown in Figure 3.
After fitting the model, we obtain the following results. Figure 4 shows the model fit (denoted TVS) and a comparison with standard OLS. Figures 5 and Figure 6 show the error distribution of the TVS fit and the convergence of the TVS model respectively.   Table 1, the estimated values for β and σ ε are significantly improved by taking into account the stochastic time delay noise. The estimated β using OLS regression is 0.5, compared to the true value 2. Even  a small amount of noise in the time axis causes attenuation limiting the usefulness of OLS for these problems. In comparison, the TVS Regression algorithm estimates β to be 2.09, much closer to the true value. While the example is based on simulated data, we believe that the experiment demonstrates clear potential for improvement in the modelling of real world systems with stochastic time delay noise.

Conclusion
We have proposed a form of regression analysis suited to the modelling of stochastic time delay problems. In addition, we have shown the feasibility of the approach and its performance on simulated data. Our approach allows for consistently improved estimation and prediction when the input is affected by noise in the time domain. To our knowledge, the method is novel for this class of problem.