Time Series fit in Mathematica

Question

Newbie in Time Series question: I need to do a very simple (I think) regression of a time serie of the type:

$Y_{t}$ = $\sum_{j=1}^{p}b_{j}Y_{t-j}$ + $\sum_{j=1}^{q}c_{j}X_{t-j}$

and I can't figure it out how do it in Mathematica. Unfortunately Mathematica documentation if not clear about this, and so I a bit lost. Does anyone know a good tutorial for Time Series in Mathematica? And how should I start?

EDIT: I also need to calculate a Granger Causality. Is there any implemented method?

EDIT: Here is some data:

Y = {{0, 2.05648}, {1, 2.05648}, {2, 2.05648}, {3, 2.05648}, {4, 2.05648}, {5, 2.05648}, {6, 2.05648}, {7, 2.05648}, {8, 2.05648}, {9, 2.05649}, {10, 2.0565}, {11, 2.05653}, {12, 2.0566}, {13, 2.05681}, {14, 2.05738}, {15, 2.05891}, {16, 2.06304}, {17, 2.07403}, {18, 2.10222}, {19, 2.16834}, {20, 2.29539}, {21, 2.46685}, {22, 2.61407}, {23, 2.69885}, {24, 2.73691}, {25, 2.75207}, {26, 2.75782}, {27, 2.75996}, {28, 2.76075}, {29, 2.76104}, {30, 2.76115}, {31, 2.76119}, {32, 2.7612}}

and

X = {{0, 18}, {1, 19}, {2, 20}, {3, 30}, {4, 19}, {5, 19}, {6, 19}, {7, 17}, {8, 17}, {9, 18}, {10, 20}, {11, 16}, {12, 16}, {13, 16}, {14, 17}, {15, 35}, {16, 30}, {17, 59}, {18, 43}, {19, 39}, {20, 39}, {21, 30}, {22, 27}, {23, 20}, {24, 23}, {25, 33}, {26, 36}, {27, 44}, {28, 61}, {29, 78}, {30, 96}, {31, 100}}

Answer 1

Here is an implementation of the regression of the type indicated in the equation.

y = {{0, 2.05648}, {1, 2.05648}, {2, 2.05648}, {3, 2.05648}, {4, 
   2.05648}, {5, 2.05648}, {6, 2.05648}, {7, 2.05648}, {8, 
   2.05648}, {9, 2.05649}, {10, 2.0565}, {11, 2.05653}, {12, 
   2.0566}, {13, 2.05681}, {14, 2.05738}, {15, 2.05891}, {16, 
   2.06304}, {17, 2.07403}, {18, 2.10222}, {19, 2.16834}, {20, 
   2.29539}, {21, 2.46685}, {22, 2.61407}, {23, 2.69885}, {24, 
   2.73691}, {25, 2.75207}, {26, 2.75782}, {27, 2.75996}, {28, 
   2.76075}, {29, 2.76104}, {30, 2.76115}, {31, 2.76119}, {32, 
   2.7612}};
x = {{0, 18}, {1, 19}, {2, 20}, {3, 30}, {4, 19}, {5, 19}, {6, 
   19}, {7, 17}, {8, 17}, {9, 18}, {10, 20}, {11, 16}, {12, 16}, {13, 
   16}, {14, 17}, {15, 35}, {16, 30}, {17, 59}, {18, 43}, {19, 
   39}, {20, 39}, {21, 30}, {22, 27}, {23, 20}, {24, 23}, {25, 
   33}, {26, 36}, {27, 44}, {28, 61}, {29, 78}, {30, 96}, {31, 100}};
xValues = x[[All, 2]];
yValues = y[[All, 2]];

The following function series[x,y,p,q] will allow me to obtain the results of the regression. That is the error between the prediction and the values of the series + the parameters of the regression.

Created several other functions to facilitate the computation of graphs, etc..

    series[x_, y_, p_, q_] := 
     Module[{b = Symbol["b" <> ToString[#]] & /@ Range[p], 
       c = Symbol["c" <> ToString[#]] & /@ Range[q], l = Length@x, 
       min = Min[p, q], fc, sc}, 
      fc = Table[Sum[b[[j]] y[[t - j]], {j, 1, p}], {t, p + 1, l}]; 
      sc = Table[Sum[c[[j]] x[[t - j]], {j, 1, q}], {t, q + 1, l}]; 
      If[p > q, fc = Join[ConstantArray[0, p - q], fc], 
       sc = Join[ConstantArray[0, q - p], sc]]; 
      NMinimize[
       Total[(#[[1]] - #[[2]])^2 & /@ 
         Transpose[{Total /@ Transpose[{fc, sc}], Drop[y, min]}]], 
       Join[b, c], Reals]]
seriesValues[x_, y_, p_, q_] := 
 Module[{b = Symbol["b" <> ToString[#]] & /@ Range[p], 
   c = Symbol["c" <> ToString[#]] & /@ Range[q], l = Length@x, 
   min = Min[p, q], fc, sc}, 
  fc = Table[Sum[b[[j]] y[[t - j]], {j, 1, p}], {t, p + 1, l}]; 
  sc = Table[Sum[c[[j]] x[[t - j]], {j, 1, q}], {t, q + 1, l}]; 
  If[p > q, fc = Join[ConstantArray[0, p - q], fc], 
   sc = Join[ConstantArray[0, q - p], sc]]; 
  Last@NMinimize[
    Total[(#[[1]] - #[[2]])^2 & /@ 
      Transpose[{Total /@ Transpose[{fc, sc}], Drop[y, min]}]], 
    Join[b, c], Reals]]
seriesPrediction[x_, y_, p_, q_] := 
 Module[{b = Symbol["b" <> ToString[#]] & /@ Range[p], 
   c = Symbol["c" <> ToString[#]] & /@ Range[q], l = Length@x, 
   min = Min[p, q], fc, sc}, 
  fc = Table[Sum[b[[j]] y[[t - j]], {j, 1, p}], {t, p + 1, l}]; 
  sc = Table[Sum[c[[j]] x[[t - j]], {j, 1, q}], {t, q + 1, l}]; 
  If[p > q, fc = Join[ConstantArray[0, p - q], fc], 
   sc = Join[ConstantArray[0, q - p], sc]]; 
  Total /@ Transpose[{fc, sc}]]

Explore different alternatives

Flatten[Table[{series[xValues, Most@yValues, i, j], i, j}, {i, 1, 
    8}, {j, 1, 8}], 1] // Sort
(*{{{0.000245998, {b1 -> 3.84299, b2 -> -6.85701, b3 -> 8.38577, 
    b4 -> -8.01898, b5 -> 6.28504, b6 -> -3.97163, b7 -> 1.67097, 
    b8 -> -0.324682, c1 -> -0.000407285, c2 -> 0.000414386, 
    c3 -> 0.000725565, c4 -> -0.000600485, c5 -> -0.000876084, 
    c6 -> -0.000100851, c7 -> -0.000505129, c8 -> 0.000108438}}, 8, 
  8}, {{0.000288402, {b1 -> 3.48744, b2 -> -5.66612, b3 -> 6.38801, 
    b4 -> -5.55299, b5 -> 3.8124, b6 -> -1.97005, b7 -> 0.514932, 
    c1 -> -0.0003576, c2 -> 0.000292707, c3 -> 0.000799256, 
    c4 -> -0.000233652, c5 -> -0.000873094, c6 -> -0.00036767, 
    c7 -> -0.000579723}}, 7, 7},...8*)

best = %[[All, 2 ;;]]
g = ListPlot[
     Evaluate[{seriesPrediction[xValues, 
         Most@yValues, #[[1]], #[[2]]] /. 
        seriesValues[xValues, Most@yValues, #[[1]], #[[2]]], 
       Drop[yValues, Min[#[[1]], #[[2]]] + 1]}], 
     Joined -> {True, False}, 
     PlotLabel -> 
      StringJoin["p=", ToString[#[[1]]], ", q=", 
       ToString[#[[2]]]]] & /@ best;
Grid[Partition[g, 4]]

You can implement the Granger Causality following @belisarius lead.