# Comparing multiple time series in R

If I have many time series that I'd like to compare to see if there are relationships between the variables, (I have several dependent variables and many more independent variables) how might I go about doing this (I'm working in R, just fyi)? I haven't really found too many examples seeking to compare and explore the relationships between many variables. In particular, I'd like to see if the variation in my independent time series drives variation in the dependent time series, and I'm so new to statistics (and R) that I'm really not sure how to approach this problem. Here is some sample data (I know there are 2 missing values, and I may choose to treat Y2 as an independent variable since I think Y2 and Y1 may be correlated. The column labeled "covariate" is years since last flood, because I think this might be important too):

year Y1 Y2  X1  X2  X3      X4  X5  covariate
1    40  92 0   0   20.6    91  503 3
2    54  65 0   0   21.7    33  175 4
3    59  75 1   1   22.2    34  94  5
4    68  53 8   9   22.2    24  86  6
5           5   20  20.6    5   185 7
6    76  65 8   13  22.2    32  119 8
7    76  55 16  18  23.3    0   153 9
8    82  58 18  2   24.4    19  0   1
9    60  57 28  24  23.33   0   223 2
10   58  46 18  3   22.78   0   184 3
11   49  48 2   1   23.33   0   110 4
12   28  76 0   3   22.78   0   213 5
13   56  61 0   1   22.78   12  123 6
14   105 53 56  24  23.33   0   122 7
15   99  43 28  13  24.44   0   154 8
16   119 47 46  35  23.33   0   182 9


Any guidance would be much appreciated as I'm not really sure how to proceed here.

here goes ... These are the steps that are required to form your analysis. Simply reproduce them in R or whatever tools you have available. The reason the following exercise is daunting is because the statistical problem you are asking is "daunting" and one needs to "up-armor" their solutions skills/procedures.

1. pre-whiten each of your 7 candidate regressors in order to identify an initial Tranfer Function model of the form

  MODEL COMPONENT       LAG    COEFF     STANDARD      P       T


# (BOP) ERROR VALUE VALUE

1CONSTANT -258. 330. 0.4568 -0.78

INPUT SERIES X1 Y2

2Omega (input) -Factor # 1 0 0.635 0.583 0.3079 1.09

INPUT SERIES X2 X1

3Omega (input) -Factor # 2 1 0.544 0.563 0.3623 0.97

INPUT SERIES X3 X2

4Omega (input) -Factor # 3 0 1.60 0.485 0.0110 3.29

INPUT SERIES X4 X3

5Omega (input) -Factor # 4 0 11.7 14.8 0.4493 0.80

INPUT SERIES X5 X4

6Omega (input) -Factor # 5 0 0.491 0.699 0.5020 0.70

INPUT SERIES X6 X5

7Omega (input) -Factor # 6 0 -0.151 0.164 0.3852 -0.92

INPUT SERIES X7 X6

8Omega (input) -Factor # 7 0 2.40 1.57 0.1639 1.53

Then do Interventon Detection to extract any anomalies in the data finding three

         :  NEWLY IDENTIFIED VARIABLE   X8   I~P00010       10    PULSE
:  NEWLY IDENTIFIED VARIABLE   X9   I~P00014       14    PULSE
:  NEWLY IDENTIFIED VARIABLE   X10  I~P00002        2    PULSE


This leads to an augmented model :

      MODEL COMPONENT       LAG    COEFF     STANDARD      P       T


# (BOP) ERROR VALUE VALUE

1CONSTANT                         -33.9       191.      0.8662    -0.18


INPUT SERIES X1 Y2

2Omega (input) -Factor #  1    0   1.08      0.324      0.0207     3.34


INPUT SERIES X2 X1

3Omega (input) -Factor #  2    1  0.867      0.325      0.0446     2.67


INPUT SERIES X3 X2

4Omega (input) -Factor #  3    0   2.06      0.248      0.0004     8.29


INPUT SERIES X4 X3

5Omega (input) -Factor #  4    0   1.68       8.47      0.8503     0.20


INPUT SERIES X5 X4

6Omega (input) -Factor #  5    0 -0.239      0.455      0.6218    -0.53


INPUT SERIES X6 X5

7Omega (input) -Factor #  6    0 -0.363      0.102      0.0161    -3.56


INPUT SERIES X7 X6

8Omega (input) -Factor #  7    0   3.21      0.712      0.0064     4.50


INPUT SERIES X8 I~P00010 10 PULSE

9Omega (input) -Factor #  8    0   30.4       6.52      0.0055     4.66


INPUT SERIES X9 I~P00014 14 PULSE

10Omega (input) -Factor # 9 0 14.7 6.42 0.0709 2.29

INPUT SERIES X 10 I~P00002 2 PULSE

11Omega (input) -Factor # 10 0 39.6 8.22 0.0048 4.81

which is over-specified thus we must step-down and obtain

      MODEL COMPONENT       LAG    COEFF     STANDARD      P       T


# (BOP) ERROR VALUE VALUE

1CONSTANT                         -2.32       7.25      0.7584    -0.32


INPUT SERIES X1 Y2

2Omega (input) -Factor #  1    0   1.10      0.872E-01  0.0000    12.58


INPUT SERIES X2 X1

3Omega (input) -Factor #  2    1   1.04      0.103      0.0000    10.10


INPUT SERIES X3 X2

4Omega (input) -Factor #  3    0   2.04      0.199      0.0000    10.22


INPUT SERIES X4 X5

5Omega (input) -Factor #  4    0 -0.335      0.301E-01  0.0000   -11.13


INPUT SERIES X5 X6

6Omega (input) -Factor #  5    0   2.84      0.663      0.0037     4.28


INPUT SERIES X6 I~P00010 10 PULSE

7Omega (input) -Factor #  6    0   27.8       6.48      0.0036     4.29


INPUT SERIES X7 I~P00014 14 PULSE

8Omega (input) -Factor #  7    0   21.3       6.17      0.0107     3.45


INPUT SERIES X8 I~P00002 2 PULSE

9Omega (input) -Factor #  8    0   32.3       5.99      0.0010     5.39


which now will suggest additional structure in the X's via cross-correlative tests between the current model residuals and the residuals from the pre-whitened X's .... that had previously remained unidentified.

      MODEL COMPONENT       LAG    COEFF     STANDARD      P       T


# (BOP) ERROR VALUE VALUE

1CONSTANT                         -25.8       6.12      0.0083    -4.22


INPUT SERIES X1 Y2

2Omega (input) -Factor #  1    0   1.11      0.598E-01  0.0000    18.50
3                              1 -0.211      0.516E-01  0.0095    -4.09


INPUT SERIES X2 X1

4Omega (input) -Factor #  2    1   1.05      0.692E-01  0.0000    15.23


INPUT SERIES X3 X2

5Omega (input) -Factor #  3    0   2.28      0.143      0.0000    15.99


INPUT SERIES X4 X5

6Omega (input) -Factor #  4    0 -0.311      0.193E-01  0.0000   -16.17
7                              1 -0.693E-01  0.103E-01  0.0011    -6.72


INPUT SERIES X5 X6

8Omega (input) -Factor #  5    0   2.19      0.466      0.0054     4.70


INPUT SERIES X6 I~P00010 10 PULSE

9Omega (input) -Factor #  6    0   19.8       4.31      0.0059     4.59


INPUT SERIES X7 I~P00014 14 PULSE

10Omega (input) -Factor # 7 0 18.6 4.21 0.0068 4.43

culminating in the final model

MODEL STATISTICS AND EQUATION FOR THE CURRENT EQUATION (DETAILS FOLLOW).

Estimation/Diagnostic Checking for Variable Y Y1
X1 Y2
X2 X1
X3 X2
X4 X5
X5 X6
: NEWLY IDENTIFIED VARIABLE X6 I~P00010 10 PULSE
: NEWLY IDENTIFIED VARIABLE X7 I~P00014 14 PULSE

Number of Residuals (R) =n 15
Number of Degrees of Freedom =n-m 7
Residual Mean =Sum R / n -0.204655E-04
Sum of Squares =Sum R**2 171.679
Variance =SOS/(n) 10.7299
Standard Deviation RMSE =SQRT(Adj Var) 3.50183
Standard Error of the Mean =Standard Dev/ (n-m) 0.935902
Mean / its Standard Error =Mean/SEM -0.218671E-04
Mean Absolute Deviation =Sum(ABS(R))/n 2.50518
AIC Value ( Uses var ) =nln +2m 37.5956
SBC Value ( Uses var ) =nln +m*lnn 38.3036
BIC Value ( Uses var ) =see Wei p153 64.1437
R Square = 0.986111
Durbin-Watson Statistic =[-A(T-1)]*2/A*2 2.85557

 D-W STATISTIC IS INCONCLUSIVE.


THE DURBIN-WATSON STATISTIC IS VALID ONLY FOR MODELS THAT HAVE A WHITE NOISE ERROR TERM AND NO LAGS OF THE Y SERIES. OTHERWISE IT IS INVALID. IN THIS CASE THE TEST IS VALID.

AUTOMATICALLY REVISING MODEL

      MODEL COMPONENT       LAG    COEFF     STANDARD      P       T


# (BOP) ERROR VALUE VALUE

1CONSTANT                         -25.8       6.12      0.0029    -4.22


INPUT SERIES X1 Y2

2Omega (input) -Factor #  1    0   1.11      0.598E-01  0.0000    18.50
3                              1 -0.211      0.516E-01  0.0035    -4.09


INPUT SERIES X2 X1

4Omega (input) -Factor #  2    1   1.05      0.692E-01  0.0000    15.23


INPUT SERIES X3 X2

5Omega (input) -Factor #  3    0   2.28      0.143      0.0000    15.99


INPUT SERIES X4 X5

6Omega (input) -Factor #  4    0 -0.311      0.193E-01  0.0000   -16.17
7                              1 -0.693E-01  0.103E-01  0.0001    -6.72


INPUT SERIES X5 X6

8Omega (input) -Factor #  5    0   2.19      0.466      0.0016     4.70


INPUT SERIES X6 I~P00010 10 PULSE

9Omega (input) -Factor #  6    0   19.8       4.31      0.0018     4.59


INPUT SERIES X7 I~P00014 14 PULSE

10Omega (input) -Factor # 7 0 18.6 4.21 0.0022 4.43

Y(T) = -25.800
+[X1(T)][(+ 1.1058+ 0.211B** 1)] +[X2(T)][(+ 1.0543B** 1)] +[X3(T)][(+ 2.2848)] +[X4(T)][(- 0.311+ 0.0693B** 1)] +[X5(T)][(+ 2.1900)] +[X6(T)][(+ 19.8039)] +[X7(T)][(+ 18.6231)]

• Wow, thank you for your answer! I'm a real novice at statistics, and I do not entirely understand what you've done here. I wish I knew how to reproduce this in R so that I could examine the steps in the process, and look up the various functions being used. – Jota May 8 '11 at 21:18
• I removed your comment replying to Frank because (1) it is counterproductive, from the point of view of this site, to split discussions off to some other location and (2) it's a bad idea to post your email addresses in such a public, searchable spot. – whuber May 8 '11 at 22:38
• @whuber ok I agree . I will edit my post and eliminate those two things. – IrishStat May 8 '11 at 22:42
• Thanks so much for taking the time to help me out and to further my understanding! It was quite productive for me! – Jota May 9 '11 at 1:05
• -1 for: implying that this is the only solution ("required"), code input/output in some unspecified language, terrible formatting, poor explanation of what your model actually means. – naught101 Apr 1 '12 at 12:25