1. Introduction

Multivariate time series (MTS) data arises from different fields of life including physical sciences, medicine, finance, engineering, business, etc. Modeling multivariate stationary time series data effectively is im- portant for many decision-making activities. MTS analysis refers to the modeling of one output series in terms of n input series by explaining the interactions and co-movements among the group of input series. The reliability of time series model will be increased when the output series is modeled in terms of more than one input series which seriously causes the variations in the output variable. In this paper, we consider individual time series as components of a vector time series and analyze them jointly.

The study of MTS analysis got a new momentum in 1950’s. Whittle ^[24] derived the least squares estimation equations for a non-deterministic stationary multiple process. In 1957, Quenouille ^[17] summa- rized the work of multivariate time series analysis. Akaike ^[1], Anderson ^[2], Box G.E.P and Jenkins. G.M ^[3], Brillinger. D.R ^[4], Lutkepohl ^[12], Liu. L.M and Hanssens. D.M. ^[13], Parzen. E ^[15], Richard. H. Jones ^[18], Reinsel ^[19], Tiao G.C and Box G E P ^[21], Tiao, G.C and R.S.Tsay ^[22] are some of them who made contributions to the field of MTS analysis. Priestly. M.B ^[16], Carter. C.G ^[6], Daubechies.I ^[8], Farge.M ^[9], Koopmann. L.H ^[11] are few people who studied wavelet and coherence in detail.

Time series is a single realization of the stochastic process. Multivariate time series is a series of observations where and made sequentially through time. represents observations at time t of the component time series. That is if are component time series then denote the multivariate time series vector at time t. The concept of rational approximation of spectral density function ^[10] for modeling univariate and bivariate stationary time series is developed in ^[10]. In this paper, this method is extended for the analysis of n components stationary time series. That is, a robust method is developed for estimating the model parameters of multiple input-single output transfer function noise model by applying rational approximation of spectral density function ^[10] and wavelet coherence ^[20].

This paper is organized as follows: Section 1 is the introduction. Section 2 deals with the theoretical pre-requisites needed to model MTS. The new method developed for the identification of multiple input- single output transfer function noise model is presented in Section 3, which is the main contribution of this paper. Section 4 explains the algorithm for the developed method and also shows the flow chart. Section 5 gives application of the developed method to the simulated data and also to the real world multivariate stationary time series data. Conclusion and scope for future work is presented in Section 6.

2. Theoretical Prerequisites

The rational approximation of spectral density function of a stationary time series is developed by Jessy ^[10]. Under this method, the rational form of the spectral density function of a stationary time series which is represented by an ARMA(p,q) model is obtained by applying the concept of Pade's approximation theory ^{[5, 14]}. The spectral density function , is defined as the Fourier transform of the auto-correlation function and is given by

(1)

where Since is an even function, for and So equation (1) becomes,

(2)

The rational form of the spectral density function of an ARMA(p,q) model, is

(3)

Applying Pade's approximation theory ^{[5, 14]} and rational approximation of spectral density function ^[10], obtain the estimated unique ARMA (p,q) model,

(4)

The analysis of a bivariate stationary time series is explained in ^[10]. In the bivariate case, there will be two-time series, one is the input series and other as the corresponding output series, and method is as follows. Given the stationary time series and and the transfer function noise model is

(5)

where and are polynomials in B, the back shift operator (B such that .). Then the spectral density function of (5) can be written as

(6)

where and are the spectral density functions of the univariate time series and respectively. From equation (6) we get,

(7)

The relation in equation (7) is used to estimate the transfer function A(B) explained in ^[10] and is given by,

(8)

Then the estimated transfer function noise model is

(9)

Obtain the observed noise series,

(10)

Now estimate the ARMA model for using the R-spec procedure ^[10], taking as the observed series. Let the ARMA model be Then the transfer function noise model for bivariate time series estimated is given by

(11)

3. Method

The new method of estimating multiple input-single output transfer function noise model rely on the concept of rational approximation of spectral density function ^[10] and wavelet coherence ^[20].

Coherence ^[6] is a measure of correlation in frequency. In time series analysis, coherence is a measure of correlation between two time series at each frequency, and it is denoted by In this paper, wavelets is used to estimate coherence due to its time-frequency localized property.

A wavelet is a ‘small wave’. A wavelet grows and decays essentially in a limited time period. A real-valued function defined over the real axis and satisfying two properties,

is called a wavelet. Wavelet transform Coherence ^[6] is a method for analyzing the coherence and phase lag between two-time series as a function of both time and frequency. So wavelet coherence ^[20] can be taken as a localized correlation coefficient in time-frequency space.

The multiple input-single output transfer function noise model is of the form,

(12)

where are pairwise uncorrelated. i.e., for and are ratios of polynomials in B which are called transfer functions. Then the spectral density function of ^[10] is given by

(13)

and the cross spectral densities will be given by,

(14)

which gives;

(15)

Now to obtain the coefficients of estimate the rational form of and using the R-spec (Rational form of spectral density function) procedure ^[10] and then equate the coefficients in the equation,

(16)

Now;

(17)

is the Coherence^[6] (squared coherency) between X and Y. Equations (16) and (17) gives equation (18).

(18)

I.e.,

(19)

In matrix form, equation (19) can be written as

(20)

where each for

each for

each for

By solving equation (19), we obtain the estimates of the transfer functions But for solving equation (19), we need the squared coherences which are estimated using method of wavelet anlaysis. Now for each and

(21)

Now

and

(22)

where each and Here and are called rational form spectral density function of the time series and Now using equations (21) and (22), we get the following equations;

(23)

(24)

The above equations hold for Now applying the method of iterations ^{[5, 10]} to the system of equations (23) and (24), we get all values of and Thus, we calculated the estimated transfer functions and the estimated transfer function noise model is,

(25)

Then the noise series is obtained using,

(26)

Now estimate the ARMA model representing the noise series using the R-spec procedure ^[10]. Then the noise model will be of the form,

(27)

Then using the equations (25) and (27), we obtain the multivariate time series model as

(28)

and the variance of the white noise that we obtain along with the noise model.

4. Algorithm and Flow Chart of the Developed Method

1. Choose a multivariate time series and take as inputs and as output for and

2. Then check stationarity of and . If all of them are stationary then go to step 3. Otherwise, go to step 7.

3. Obtain the initial rational forms of the spectral density function of all the series using the R-spec procedure ^{[5, 10]}.

4. Evaluate coherence using wavelets and thus estimate the coefficients of the transfer functions which give the estimated transfer function noise model represented in equation (25).

5. Find noise series and estimate the ARMA model representing the noise series using the R-spec procedure ^{[5, 10]}.

6. Obtain the multivariate stationary time series model by combining the estimated transfer function model and the noise series.

7. If and are not stationary, then find the difference and and then check and are stationary or not. If it is stationary, then go to step 3. Otherwise go to step 8.

8. If and are not stationary, then continuously differencing and till n=5 and check stationarity at each time. If and are stationary within then go to step 3. Otherwise, go to step 9.

9. If and are not stationary then for do variance stabilizing transformation for the series and again check and are stationary or not. If and are stationary then go to step 3. Otherwise, go to step 10.

10. If and are not stationary then conclude that the given series is non-linear and non-stationary and this method is not applicable.

5. Application - Numerical Results and Comparison

The developed method is verified using simulation studies and it is applied to US hog data^[17] (US hog data studied by Quenouille from 1867 to 1948).

5.1. Simulation Study Results

The method was applied to simulated multivariate time series data and the results obtained are given in Table 5.1.

Table 5.1: For verifing the method developed, a number time series were simulated as input series and output series and estimated the parameters of multiple-input single output transfer function noise model with various parameter combinations by applying the developed method. Table 5.1 shows eleven selected examples from the simulated studies. Also estimated Mean absolute percentage error (MAPE)

Simulated studies from No.1 to No.5 are examples of bivariate time series models. That is one input series and one output series and its transfer function noise model are simulated by applying the method above. The model is given by

(29)

Examples from No.6 to No.8 are illustrating about simulated transfer function noise models of two inputs series and one output series and the model is given by

(30)

Examples from No.9 to No.11 pointing towards the simulated studies of three input series and one output series and its transfer function noise model is given by

(31)

5.2. Analysis of US Hog Data from 1867 to 1948

The method is applied to model US hog data ^[17]. Estimated transfer function parameters are given in Table 5.2.

US hog data consists of 5 component series namely, hog number, hog price, corn price, corn supply and wage rates, each having 82 observations. Here hog price- (Figure 2), corn price- (Figure 3), corn supply- (Figure 4) and wage rates- (Figure 5) are taken as inputs(component time series) and hog number (Figure 1) as output for the multiple-input transfer function model expressed in equation (12). Here corn supply is omitted for model building as its contribution is not significant indeed (Figure 18). Then find auto-correlations and partial auto-correlations of all inputs (Figure 10, Figure 11, Figure 14, Figure 15, Figure 19, Figure 20) and output (Figure 6 & Figure 7). First order differenced series in each of these series are found to be stationary (Figure 8, Figure 12, Figure 16, Figure 22) as the differenced autocorrelation of each of the series terminates sudden. Also calculated partial auto-correlations of all differenced series (Figure 9, Figure 13, Figure 17, Figure 22). Then applied the method above and estimated the parameters of multiple-input single-output transfer function noise model. Also estimated Mean absolute percentage error (MAPE).

Table 5.2: shows the estimated parameters of multiple-input transfer function noise model by applying the developed method. Mean Absolute Percentage Error (MAPE) is also mentioned.

The model estimated for is given by

(32)

Therefore the final model estimated for the US hog data is given by,

(33)

where is given in equation (32).

6. Computational Complexity Analysis

The Computational complexity analysis of the matlab programs of simulated studies and US hog data are given below.

7. Conclusion and Future Work

A method for modeling stationary multivariate time series by estimating the parameters of multiple-input single-output transfer function noise model using rational approximation of spectral density function ^[10] and wavelet coherence ^[20] is developed in this paper. The method is verified using simulated data and also applied to US hog data ^[17]. Applying this method to more real world time series models and exploring more mathematical properties of this model is one part of future work.

8. Appendix

8.1. Matlab Code

%R−spec Procedure for the five series Us Hog data from 1867 to 1948

tic;

Y= load('Hognumber.m');

n=length(Y);

acf= autocorr(Y);

plot(acf);

stem(acf);

title('SAMPLE AUTO−CORRELATION OF HOG NUMBER DATA');

figure,acf;

pacf=parcorr(Y);

stem(pacf);

title('SAMPLE PARTIAL AUTO CORRELATION OF HOG NUMBER DATA');

figure,pacf;

Y1=diff(Y);

n1=length(Y1);

acf1=autocorr(Y1);

plot(acf1);

stem(acf1);

title('SAMPLE AUTO−CORRELATION OF DIFFERENCED HOG NUMBER DATA');

figure,acf1;

pacf1=parcorr(Y1);

stem(pacf1);

title('SAMPLE PARTIAL AUTOCORRELATION OF DIFFERENCED HOG NUMBER DATA');

figure,pacf1;

x1 = load('Hogprice.m');

n2=length(x1);

acf2=autocorr(x1);

plot(acf2);

stem(acf2);

title('SAMPLE AUTO−CORRELATION OF HOG PRICE DATA');

figure,acf2;

pacf2=parcorr(x1);

stem(pacf2);

title('SAMPLE PARTIAL AUTO CORRELATION OF HOG PRICE DATA');

figure,pacf2;

X1=diff(x1);

n3=length(X1);

acf3=autocorr(X1);

plot(acf3);

stem(acf3);

title('SAMPLE AUTO−CORRELATION OF DIFFERENCED HOG PRICE DATA');

figure,acf3;

pacf3=parcorr(X1);

stem(pacf3);

title('SAMPLE PARTIAL AUTOCORRELATION OF DIFFERENCED HOG PRICE DATA');

figure,pacf3;

x2= load('Cornprice.m');

n4=length(x2);

acf4=autocorr(x2);

plot(acf4);

stem(acf4);

title('SAMPLE AUTO−CORRELATION OF CORN PRICE DATA');

figure,acf4;

pacf4=parcorr(x2);

stem(pacf4);

title('SAMPLE PARTIAL AUTO CORRELATION OF CORN PRICE DATA');

figure,pacf4;

X2=diff(x2);

n5=length(X2);

acf5=autocorr(X2);

plot(acf5);

stem(acf5);

title('SAMPLE AUTO−CORRELATION OF DIFFERENCED CORN PRICE DATA');

figure,acf5;

pacf5=parcorr(X2);

stem(pacf5);

title('SAMPLE PARTIAL AUTOCORRELATION OF DIFFERENCED CORN PRICE DATA');

figure,pacf5;

X3 = load('Cornsupply.m');

n6=length(X3);

acf6=autocorr(X3);

plot(acf6);

stem(acf6);

title('SAMPLE AUTO−CORRELATION OF CORN SUPPLY DATA');

figure,acf6;

pacf6=parcorr(X3);

stem(pacf6);

title('SAMPLE PARTIAL AUTO CORRELATION OF CORN SUPPLY DATA');

figure,pacf5;

x4 = load('Wagerates.m');

n7=length(x4);

acf7=autocorr(x4);

plot(acf7);

stem(acf7);

title('SAMPLE AUTO−CORRELATION OF WAGE RATES DATA');

figure,acf7;

pacf7=parcorr(x4);

stem(pacf7);

title('SAMPLE PARTIAL AUTO CORRELATION OF WAGE RATES DATA');

figure,pacf7;

X4=diff(x4);

n8=length(X4);

acf8=autocorr(X4);

plot(acf8);

stem(acf8);

title('SAMPLE AUTO−CORRELATION OF DIFFERENCED WAGE RATES DATA');

figure,acf8;

pacf8=parcorr(X4);

stem(pacf8);

title('SAMPLE PARTIAL AUTOCORRELATION OF DIFFERENCED WAGE RATES DATA');

figure,pacf8;

for k=[−1:0.1:1]

la=k*pi;

a0=0.5;

a1=acf1(2);

T01=a0+a1*cos(la);

end

for k=[−1:0.1:1]

la=k*pi;

b1=−2*acf1(2)/(1+acf1(3));

a01=0.5*(1+acf1(2)*b1);

T10=a01/(1+b1*cos(la));

end

for k=[−1:0.1:1]

la=k*pi;

b11=−2*acf1(3)/(acf1(2)+acf1(4));

a02=0.5*(1+acf1(2)*b11);

a11=0.5*(2*acf1(2)+b11*(1+acf1(3)));

T11=(a02+a11*cos(la))/(1+b11*cos(la));

end

for i=1:20

for k=[−1:0.1:1]

la=k*pi;

f(i)=(0.5+sum(acf1(i+1)*cos(i*la)));

E1(i)=f(i)−T01;

E2(i)=f(i)−T10;

E3(i)=f(i)−T11;

end

end

E11=abs(E1);

Error1=max(E11);

E12=abs(E2);

Error2=max(E12);

E13=abs(E3);

Error3=max(E13);

M=min([Error1(:);

Error2(:); Error3(:)]);

for k=[−1:0.1:1]

la=k*pi; a03=0.5;

a12=acf3(2);

t01=a03+a12*cos(la);

end

for k=[−1:0.1:1]

la=k*pi;

b12=−2*acf3(2)/(1+acf3(3));

a04=0.5*(1+acf3(2)*b12);

t10=a04/(1+b12*cos(la));

end

for k=[−1:0.1:1]

la=k*pi;

b13=−2*acf3(3)/(acf3(2)+acf3(4));

a05=0.5*(1+acf3(2)*b13);

a13=0.5*(2*acf3(2)+b13*(1+acf3(3)));

t11=(a05+a13*cos(la))/(1+b13*cos(la));

end

for i=1:20

for k=[−1:0.1:1]

la=k*pi;

f(i)=(0.5+sum(acf3(i+1)*cos(i*la)));

e1(i)=f(i)−t01;

e2(i)=f(i)−t10;

e3(i)=f(i)−t11;

end

end

e11=abs(e1);

Error11=max(e11);

e12=abs(e2);

Error21=max(e12);

e13=abs(e3);

Error31=max(e13);

M1=min([Error11(:); Error21(:); Error31(:)]);

for k=[−1:0.1:1]

la=k*pi;

a014=0.5;

a111=acf5(2);

t012=a014+a111*cos(la);

end

for k=[−1:0.1:1]

la=k*pi;

b110=−2*acf5(2)/(1+acf5(3));

a015=0.5*(1+acf5(2)*b110);

t102=a015/(1+b110*cos(la));

end

for k=[−1:0.1:1]

la=k*pi;

b111=−2*acf5(3)/(acf5(2)+acf5(4));

a016=0.5*(1+acf5(2)*b111);

a112=0.5*(2*acf5(2)+b111*(1+acf5(3)));

t112=(a016+a112*cos(la))/(1+b111*cos(la));

end

for i=1:20

for k=[−1:0.1:1]

la=k*pi;

f(i)=(0.5+sum(acf5(i+1)*cos(i*la)));

ec1(i)=f(i)−t012;

ec2(i)=f(i)−t102;

ec3(i)=f(i)−t112;

end

end

ec11=abs(ec1);

Error12=max(ec11);

ec22=abs(ec2);

Error22=max(ec22);

ec33=abs(ec3);

Error32=max(ec33);

M2=min([Error12(:); Error22(:); Error32(:)]);

for k=[−1:0.1:1]

la=k*pi;

a017=0.5;

a113=acf6(2);

t013=a017+a113*cos(la);

end

for k=[−1:0.1:1]

la=k*pi;

b112=−2*acf6(2)/(1+acf6(3));

a018=0.5*(1+acf6(2)*b112);

t103=a018/(1+b112*cos(la));

end

for k=[−1:0.1:1]

la=k*pi;

b113=−2*acf6(3)/(acf6(2)+acf6(4));

a019=0.5*(1+acf6(2)*b113);

a114=0.5*(2*acf6(2)+b113*(1+acf6(3)));

t113=(a019+a114*cos(la))/(1+b113*cos(la));

end

for i=1:20

for k=[−1:0.1:1]

la=k*pi;

f(i)=(0.5+sum(acf6(i+1)*cos(i*la)));

es1(i)=f(i)−t013;

es2(i)=f(i)−t103;

es3(i)=f(i)−t113;

end

end

es11=abs(es1);

Error13=max(es11);

es22=abs(es2);

Error23=max(es22);

es33=abs(es3);

Error33=max(es33);

M3=min([Error13(:); Error23(:); Error33(:)]);

for k=[−1:0.1:1]

la=k*pi;

a020=0.5;

a115=acf8(2);

t014=a020+a115*cos(la);

end

for k=[−1:0.1:1]

la=k*pi;

b114=−2*acf8(2)/(1+acf8(3));

a021=0.5*(1+acf8(2)*b114);

t104=a021/(1+b114*cos(la));

end

for k=[−1:0.1:1]

la=k*pi;

b115=−2*acf8(3)/(acf8(2)+acf8(4));

a022=0.5*(1+acf8(2)*b115);

a116=0.5*(2*acf8(2)+b115*(1+acf8(3)));

t114=(a022+a116*cos(la))/(1+b115*cos(la));

end

for i=1:20

for k=[−1:0.1:1]

la=k*pi;

f(i)=(0.5+sum(acf8(i+1)*cos(i*la)));

ew1(i)=f(i)−t014;

ew2(i)=f(i)−t104;

ew3(i)=f(i)−t114;

end

end

ew11=abs(ew1);

Error14=max(ew11);

ew22=abs(ew2);

Error24=max(ew22);

ew33=abs(ew3);

Error34=max(ew33);

M4=min([Error14(:);Error24(:); Error34(:)]);

A01=a0;

A11=a1;

A21=(a1*b13);

B01=a03;

B11=a12;

syms W01 W11

[W01,W11]=solve(W01ˆ2+W11ˆ2−K11*A01==0,2*W01*W11−K11*A11==0);

w01=double(W01(1));

w11=double(W11(1));

A02=a0;

A12=(a0*b111+a1);

A22=(a1*b111);

B02=a016; B12=(a112);

syms W02 W12 W22 [W02,W12,W22]=solve(W02ˆ2+W12ˆ2+W22ˆ2−K12*A02==0,2*W02*W12+ 2*W12*W22−K12*A12==0,2*W02*W22−K12*A22==0);

w02=double(W02(1));

w12=double(W12(1));

w22=double(W22(1));

A04=(a0);

A14=(a0*b114+a1);

A24=(a1*b114);

B04=a021;

syms W04 W14 W24 [W04,W14,W24]=solve(W04ˆ2+W14ˆ2+W24ˆ2−K14*A04==0,2*W04*W14+ 2*W14*W24−K14*A14==0,2*W04*W24−K14*A24==0);

w04=double(W04(1));

w14=double(W14(1));

w24=double(W24(1))

for t=10:80

for n=0:1:8

for i=0:1:6

ycap(t)=(−1)ˆn*w01*delta11ˆi*X1(t−i)+(−1)ˆn*w11*delta11ˆi*X1(t−(i+1))

+(−1)ˆn*w02*delta12ˆi*X2(t−i)+(−1)ˆn*w12*delta12ˆi*X2(t−(i+1))+

(−1)ˆn*w22*delta12ˆi*X2(t−(i+2))+(−1)ˆn*w04*delta14ˆi*X4(t−i)+

(−1)ˆn*w14*delta14ˆi*X4(t−(i+1))+(−1)ˆn*w24*delta14ˆi*X4(t−(i+2));

end

end

end

for t=10:80

noise(t)=Y1(t)−ycap(t);

end

plot(noise);

stem(noise)

title('NOISE OF MULTIVARIATE DATA');

figure,noise;

acf11=autocorr(noise);

plot(acf11);

stem(acf11);

title('SAMPLE AUTO−CORRELATION OF NOISE SERIES');

figure,acf11;

n12=length(noise);

pacf11=parcorr(noise);

plot(pacf11);

stem(pacf11);

title('SAMPLE PARTIAL AUTO−CORRELATION OF NOISE SERIES');

figure,pacf11;

fprintf('\n\n PARTICAL AUTOCORRELATION OF NOISE SERIES \n');

for t=3:80

D11(t)=(ycap(t)−Y1(t))ˆ2;

D12=sum(D11(t));

RMSE=sqrt(D12/80);

end

plot(RMSE);

stem(RMSE);

title('ROOT MEAN SQUARE ERROR');

fprintf('\n\n The value of RMSE=%f',RMSE);

for t=3:80

D(t)=abs((ycap(t)−Y1(t))/Y1(t));

D1=sum(D(t));

MAPE=(D1/n1)*100;

end

plot(MAPE);

stem(MAPE);

title('MEAN ABSOLUTE PERCENTAGE ERROR');

fprintf('\n\n The value of MAPE=%f',MAPE);

toc;

References

[1]	Akaike H. Fitting Autoregressive Models for Prediction. Annals of the Institute of Statistical Math- ematics, 21 (1969): 243-247.
	In article		View Article
[2]	Anderson T W. An Introduction to Multivariate Statistical Analysis, 2nd ed. New York, 1984, John Wiley & Sons Inc.
	In article
[3]	Box G E P and Jenkins G M. Time Series Analysis, Forecasting and Control, San Francisco, Holden-Day (1970).
	In article
[4]	Brillinger D R. Time Series data analysis & Theory, Extended edition, Holden-Day, San Francisco, 1981.
	In article		PubMed
[5]	Blum E K. Numerical Analysis and Computation Theory and Practice, Menlo Park, California, Addison Wesley Publishing Company, Inc, (1972).
	In article
[6]	Carter C G. Coherence and time-delay estimation. Proc. IEEE, 75, 1987: 236-255.
	In article		View Article
[7]	Capon J. High resolution frequency-wavenumber spectrum anlysis, Proc. of the IEEE, vol. 57, pp. 1408-1418, Aug, 1969.
	In article		View Article
[8]	Daubechies I. The wavelet transforms time-frequency localization and signal analysis. IEEE Trans. Inform. Theory, Ten Lectures on Wavelets. Society for Industrial and Applied Mathematics, 357 pp, (1990).
	In article
[9]	Farge M. Wavelet Transforms and their applications to turbulence, Annu. Rev.Fluid Mech., 24, 1992: 395-457.
	In article		View Article
[10]	Jessy J C. Arma Modelling of Time Series Based On Rational Approximation of Spectral Density Functions, Department of Mathematics and Statistics,University of Cochin, 1985.
	In article		PubMed
[11]	Koopmans L H. On the coefficient of coherence for weakly stationary stochastic processes. Ann. Math. Statist. 35, 1964: 532-549.
	In article		View Article
[12]	Lutkepohl H. Introduction to Multiple Time Series. Berlin, 1991: Springer.
	In article		View Article
[13]	Liu L M and Hanssens D.M. Identification of Multiple-Input Transfer Function Models, Communi-cations in Statistics - Theory and Methods, 11(3), 1982: 297-314.
	In article		View Article
[14]	Padé H. Sur la répresentation approchée d’une fonction par des fractions rationelles, Thesis, [Ann.E´ cole Nor. (3), 9, 1892, pp. 1:93 supplement.
	In article
[15]	Parzen E. On empirical multiple time series analysis, Technical Report No.13, Stanford University, July 13, 1965.
	In article
[16]	Priestly M B. Spectral Analysis and Time Series, vol 1, Univariate Series, 653 pp., Academic, San Diego, Calif., 1981 a.
	In article
[17]	Quenouille M H. The Analysis of Multiple Time Series, London: Griffin (1957).
	In article
[18]	Richard H Jones. Prediction of Multivariate time series. Department of Statistics, Johns Hopkins University, Baltimore Md (1964).
	In article		PubMed
[19]	Reinsel G C (1997). Elements of Multivariate Time Series Analysis, 2nd Ed. New York: Springer.
	In article		View Article
[20]	Torrence C and G P Combo. A practical guide to wavelet analysis, Bull. Am. Meterol. Soc., 79, 1998: 61-7.
	In article		View Article
[21]	Tiao G.C and Box G E P. Modelling Multiple Time Series with Applications, Journal of the American Statistical Association, Vol. 76. 1981: 802-816.
	In article
[22]	Tiao G C and Tsay R S. Model Specification in Multivariate Time Series, J.R. Statist. Soc.B (1989) 51, No. 2, pp 157-213.
	In article
[23]	Wei, William W.S (1989). Time Series Analysis Univariate and Multivariate Methods, Addison- Wesley, New York.
	In article
[24]	Whittle.P. The Analysis of Multiple Stationary Time Series, Journal of the Royal Statistical Society. Series B (Methodological), vol 15, no.1, 1953: 125-139.
	In article