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.