Cross-Correlation with TimeSeries

Question

How to find hidden relationship between two TimeSeries data sets? I normally compute the Cross-Correlation function, but I cannot find an easy way of doing this on the TimeSeries objects in Mathematica.

Do you know an elegant way of computing Cross-Correlation between two TimeSeries objects?

Answer 1

The only elegant thing I've found that can be done for this question is the shifting of the time series with TimeSeriesShift to help get the answer.

{pepsi, mcds} = 
 TimeSeries[FinancialData[{"NYSE", #}, 
  "Price", {DateObject[{2012, 12, 31}], DateObject[{2015, 3, 31}], 
  "Week"}]] & /@ {"PEP", "MCD"}

Get to time series. TimeSeriesShift will shift the series the number of specified time units.

DateListPlot[{mcds, TimeSeriesShift[mcds, Quantity[-10, "Week"]]}, 
 PlotLegends -> {"before", "after shifting"}]

This can be used to calculate the cross-correlation.

crossCorrelation[tsX_, tsY_, shift_?QuantityQ] :=
 Module[{μX = Mean@tsX, σX = 
    StandardDeviation@tsX, μY = Mean@tsY, σY = 
    StandardDeviation@tsY, tsYShift, xValues, yValues, minLength, 
   pairs},
  tsYShift = TimeSeriesShift[tsY, shift];
  xValues = 
   Select[tsX["Path"], #[[1]] <= tsYShift["LastTime"] &][[All, 2]];
  yValues = 
   Select[tsYShift["Path"], #[[1]] >= tsX["FirstTime"] &][[All, 2]];
  minLength = Min[Length /@ {xValues, yValues}];
  pairs = 
   Transpose[{xValues[[;; minLength]], yValues[[;; minLength]]}];
  Last[FoldList[
     Function[{running, next}, 
      running + (next[[1]] - μX) (next[[2]] - μY)], 0, 
     pairs]]/(minLength σX σY)
  ]

You will notice that many function can take TimeSeries objects without having to strip off the time column (e.g. Mean). This makes things a little easier. Using crossCorrelation for the two time series the second series is shifted backwards and then the two series are truncated to hold the same number of points. I couldn't find a built-in function to truncate a time series.

crossCorrelation[pepsi, mcds, Quantity[-10, "Weeks"]]
(* -0.0492425 *)

This looks reasonable from the plot.

DateListPlot[{pepsi, TimeSeriesShift[mcds, Quantity[-10, "Week"]]}, 
 PlotLegends -> {"Pepsi", "McDonald's (-10 week shift)"}]

Hope this helps.

Answer 2

Using "Values" you can extract the numerical values of the time series and perform whatever statistical analysis you want on them.

For example, if the TimeSeriesare financial data:

data1 = FinancialData["SBUX", 
   "Close", {{2013, 1, 1}, {2013, 12, 1}, "Day"}];
data2 = FinancialData["SBUX", 
   "Close", {{2013, 5, 1}, {2013, 6, 1}, "Day"}];
ts1 = TimeSeries[data1[[All, 2]], {data1[[1, 1]], Automatic, "Day"}];
ts2 = TimeSeries[data2[[All, 2]], {data2[[1, 1]], Automatic, "Day"}];

Then you can calculate the correlation as

corr = ListCorrelate[ts2["Values"], ts1["Values"]]

and ListLinePlot[corr] gives

If the times of the TimeSeries are not equally spaced, you can use Interpolation (on ts["Times"]) to upsample the data at equally-spaced coordinates and preform the analysis on that.

Answer 3

Based on Edmund's answer I have written a function that will produce cross-correlation sweep. Excellent for sweeps, but not optimised for performance.

TSCorrelation[ts1_, ts2_] := Module[
  {aligned, r, z, seZ, ZtoR, up, down, error},
  aligned = Select[ 
     GatherBy[  Normal@ts1~ Join ~ Normal@ts2, First],
     Length[#] == 2 &] [[All, {1, 2}, 2]];
  (*chaeck if enough data is available*) 
  If[Length@aligned < 3, Return[
    <| "correlation" -> 0,  "error" -> 1.0 , 
     "error bounds" -> {1, 0}, "size" ->  Length@aligned |> ]];
  r = N@Correlation[   aligned[[All, 1]] , aligned[[All, 2]]  ];
  (* error estimate *)
  z = 0.5 Log[ (1 + r)/(1 - r)];
  seZ = 1/Sqrt[Length@aligned - 3];
  ZtoR[zz_] := (Exp[2 zz] - 1)/(Exp[2 zz] + 1);
  {up, down} = {ZtoR[z + seZ], ZtoR[z - seZ]};
  error = (up - down )/2;
  (* out *)
  <| "correlation" -> r,  "error" -> error , 
     "error bounds" -> {up, down}, "size" ->  Length@aligned |>
]

and

TSCrossCorrelation[ts1_, ts2_, unit_: "Day"] := Module[
  {minShift, maxShift, ccRaw},
  minShift = -QuantityMagnitude@
     DateDifference[ts1["FirstDate"], ts2["LastDate"], unit];
  maxShift = 
   QuantityMagnitude@
    DateDifference[ts2["FirstDate"], ts1["LastDate"], unit];
  ccRaw = Quiet@Table[
     Append[ TSCorrelation[TimeSeriesShift[ ts1, {shift, unit}], ts2],
       "shift" -> shift ],
     {shift, Range @@ Sort@{minShift, maxShift}}
     ];
  (*filter out high errors *)
  Select[ ccRaw, #["error"] <= 0.3 &]
]

Usage example:

data1 = FinancialData["SBUX", 
   "Close", {{2013, 1, 1}, {2013, 12, 1}, "Day"}];
data2 = FinancialData["SBUX", 
   "Close", {{2013, 5, 1}, {2013, 6, 1}, "Day"}];
ts1 = TimeSeries[data1[[All, 2]], {data1[[1, 1]], Automatic, "Day"}];
ts2 = TimeSeries[data2[[All, 2]], {data2[[1, 1]], Automatic, "Day"}];
cc = Dataset@TSCrossCorrelation[ ts1, ts2];
ListPlot @ cc[All, {"shift", "correlation"}] 

Answer 4

Thanks to rasher's answer to Unragging a matrix I have written:

PreFill := 
 Module[{t = ConstantArray["x", Length@Dates], 
    p = #[[All, 1]] /. 
      Dispatch@Thread@Rule[Dates, Range@Length@Dates]}, t[[p]] = #; 
   t] &


Filler[dt_] :=
 Module[{f, pr},
  f[x_, "x"] := x;
  f[_, x_] := x;
  pr = FoldList[f, FoldList[f, #][[-1 ;; 1 ;; -1]]][[-1 ;; 
       1 ;; -1]] & /@ dt;
  MapIndexed[(pr[[All, First@#2, 1]] = #1) &, Dates];
  pr]

to compute the correlation between timeseries of different lengths. For example, if you want to find the correlation between the stock prices of an us company and a german company you always have the problem that some dates are missing because of different bank holidays. In this case "fill" the missing value(s) with the preceding or following value.

Example:

pr = Map[FinancialData[#, "Price", {2012}] &, {"F:DAI", "F:VOW"}];

Dates = MapFirst@SelectFirst[pr, Length@# == Max@MapLength@pr &];

Correlation[First@#, Last@#] &[Map[Last, Filler@Map[PreFill, pr], {2}]]

0.92702