Correlating noisy time series with a phase difference (lag) I am trying to find the correlation (or any other indicator of "similarity") between a real-world time series (example: monthly sales of tractors - seasonal over the year) and some market index, like the S&P GSCI Agriculture Index (link). I am trying to figure out how to think about the problem and find ways to get real insight from both time series. Both time series would be expected to have autocorrelation, seasonal components, and systematic noise.
Would you smooth both with some moving average then check for correlation for different lags? or decompose both time series and look at the trend correlation?
What is your method of doing something like this?
Thanks!
 A: Here's some Python code I wrote once to test out the max correlation method. The idea here is to perform correlations between the two signals, and then whichever phase angle produces the maximum such correlation is the winner. Naturally, more noise would make this process more difficult: I would recommend a median filter for smoothing data before trying this, although there are loads of options.
import numpy as np


def float_equal(x, y) -> bool:
    """
    This function measures whether two floating-point numbers are
    equal, to within a certain tolerance
    :param x:
    :param y:
    :return: Whether the two numbers are equal, to within machine
                epsilon.
    """

    tol = 0.0000001
    return np.abs(x - y) < tol


# Define the array length, and the begin and end of the x location
# array.
array_length = 10000
begin_x_location = 0
end_x_location = 2 * np.pi

# x array of values from 0 to 2Pi.
signal1_x = np.arange(begin_x_location,
                      end_x_location,
                      end_x_location / array_length)

# Need a basic cosine.
signal1_y = np.cos(signal1_x)

# A shifted cosine.
signal2_y = np.cos(signal1_x - np.pi / 4)

# Add some noise to each y signal above, so that the correlations aren't
# perfect. Want the magnitude to be significantly less than the signal.
full_signal1 = signal1_y + np.random.normal(0, 0.01, array_length)
full_signal2 = signal2_y + np.random.normal(0, 0.01, array_length)

# Debugging how to do the correlations.
print(str(len(full_signal1[:array_length-1250])))
print(str(len(full_signal2[1250:])))

# Note that 1250 corresponds to pi/4, the shift. We expect this to be
# the maximum correlation.
corr = np.corrcoef(full_signal1[:array_length-1250],
                   full_signal2[1250:])

# The correlation at the known phase shift.
print(str(corr[0][1]))

print("Now we find max correlation coefficient:")

# Create a vector of correlation coefficients. Each coefficient is
# generated by shifting one array relative to the other. We need the
# two arrays to be the same length, so we simply trim off the unwanted
# numbers to get two equal-sized arrays.
corrs = np.array([np.corrcoef(full_signal1[:array_length-i],
                              full_signal2[i:])[0][1]
                  for i in range(array_length-100)])

# This is the max correlation coefficient.
print('Maximum correlation: ' + str(np.max(corrs)))

# This is where that maximum occurred. Use the float_equal function.
index_location = np.where(float_equal(corrs, np.max(corrs)))[0][0]
print('Maximum occurred at index: ' + str(index_location))

# Now convert index to x location and print out where that happened.
# Note that the answer is in radians.
x_location = (end_x_location - begin_x_location) * index_location \
             / array_length + begin_x_location
print('Maximum occurred at x value: ' + str(x_location))
```

