;+ ; NAME: ; C_TIMECORRELATE ; ; PURPOSE: ; This function computes the "time cross correlation" Pxy(L) or ; the "time cross covariance" between 2 arrays (this is some ; kind of c_correlate but for multidimenstionals arrays) as a ; function of the lag (L). ; ; CATEGORY: ; Statistics. ; ; CALLING SEQUENCE: ; Result = c_timecorrelate(X, Y, Lag) ; ; INPUTS: ; X: an Array which last dimension is the time dimension of ; size n, float or double. ; ; Y: an Array which last dimension is the time dimension of ; size n, float or double. ; ; LAG: A scalar or n-element vector, in the interval [-(n-2), (n-2)], ; of type integer that specifies the absolute distance(s) between ; indexed elements of X. ; ; KEYWORD PARAMETERS: ; COVARIANCE: If set to a non-zero value, the sample cross ; covariance is computed. ; ; DOUBLE: If set to a non-zero value, computations are done in ; double precision arithmetic. ; ; EXAMPLE ; Define two n-element sample populations. ; x = [3.73, 3.67, 3.77, 3.83, 4.67, 5.87, 6.70, 6.97, 6.40, 5.57] ; y = [2.31, 2.76, 3.02, 3.13, 3.72, 3.88, 3.97, 4.39, 4.34, 3.95] ; ; Compute the cross correlation of X and Y for LAG = -5, 0, 1, 5, 6, 7 ; lag = [-5, 0, 1, 5, 6, 7] ; result = c_timecorrelate(x, y, lag) ; ; The result should be: ; [-0.428246, 0.914755, 0.674547, -0.405140, -0.403100, -0.339685] ; ; PROCEDURE: ; ; ; FOR L>=0 ; ; n-L-1 ; sigma (X[k]-Xmean)(Y[k+L]-Ymean) ; k=0 ; correlation(X,Y,L)=------------------------------------------------------ ; n-1 n-1 ; sqrt( (sigma (X[k]-Xmean)^2)*(sigma (Y[k]-Ymean)^2)) ; k=0 k=0 ; ; ; ; n-L-1 ; sigma (X[k]-Xmean)(Y[k+L]-Ymean) ; k=0 ; covariance(X,Y,L)=------------------------------------------------------ ; n ; ; FOR L<0 ; ; ; n-L-1 ; sigma (X[k+L]-Xmean)(Y[k]-Ymean) ; k=0 ; correlation(X,Y,L)=------------------------------------------------------ ; n-1 n-1 ; sqrt( (sigma (X[k]-Xmean)^2)*(sigma (Y[k]-Ymean)^2)) ; k=0 k=0 ; ; ; ; n-L-1 ; sigma (X[k+L]-Xmean)(Y[k]-Ymean) ; k=0 ; covariance(X,Y,L)=------------------------------------------------------ ; n ; ; Where Xmean and Ymean are the time means of the sample populations ; x=(x[t=0],x[t=1],...,x[t=n-1]) and y=(y[t=0],y[t=1],...,y[t=n-1]), ; respectively. ; ; ; ; REFERENCE: ; INTRODUCTION TO STATISTICAL TIME SERIES ; Wayne A. Fuller ; ISBN 0-471-28715-6 ; ; MODIFICATION HISTORY: ; 01/03/2000 Sebastien Masson (smasson@lodyc.jussieu.fr) ; Based on the C_CORRELATE procedure of IDL ;- FUNCTION TimeCross_Cov, X, Y, M, nT, Double = Double, ZERO2NAN = zero2nan ; if double then one = 1.0d ELSE one = 1.0 ; Sample cross covariance function. TimeDim = size(X, /n_dimensions) Xmean = TOTAL(X, TimeDim, Double = Double) / nT Xmean = Xmean[*]#replicate(one, nT - M) Ymean = TOTAL(Y, TimeDim, Double = Double) / nT Ymean = Ymean[*]#replicate(one, nT - M) ; case TimeDim of 1:res = TOTAL((X[0:nT - M - 1L] - Xmean) * (Y[M:nT - 1L] - Ymean) $ , Double = Double) 2:res = TOTAL((X[*, 0:nT - M - 1L] - Xmean) * (Y[*, M:nT - 1L] - Ymean) $ , TimeDim, Double = Double) 3:res = TOTAL((X[*, *, 0:nT - M - 1L] - Xmean) * (Y[*, *, M:nT - 1L] - Ymean) $ , TimeDim, Double = Double) 4:res = TOTAL((X[*, *, *, 0:nT - M - 1L] - Xmean) * (Y[*, *, *, M:nT - 1L] - Ymean) $ , TimeDim, Double = Double) ENDCASE if keyword_set(zero2nan) then begin zero = where(res EQ 0) if zero[0] NE -1 then res[zero] = !values.f_nan ENDIF ; RETURN, res END FUNCTION C_Timecorrelate, X, Y, Lag, Covariance = Covariance, Double = Double ;Compute the sample cross correlation or cross covariance of ;(Xt, Xt+l) and (Yt, Yt+l) as a function of the lag (l). ON_ERROR, 2 xsize = SIZE(X) ysize = SIZE(Y) nt = float(xsize[xsize[0]]) NDim = xsize[0] if total(xsize[0:xsize[0]] NE ysize[0:ysize[0]]) NE 0 then $ MESSAGE, "X and Y arrays must have the same size and the same dimensions" ;Check length. if nt lt 2 then $ MESSAGE, "Time dimension of X and Y arrays must contain 2 or more elements." ;If the DOUBLE keyword is not set then the internal precision and ;result are identical to the type of input. if N_ELEMENTS(Double) eq 0 then $ Double = (Xsize[Xsize[0]+1] eq 5 or ysize[ysize[0]+1] eq 5) if n_elements(lag) EQ 0 then lag = 0 nLag = N_ELEMENTS(Lag) if nLag eq 1 then Lag = [Lag] ;Create a 1-element vector. case NDim of 1:if Double eq 0 then Correl = FLTARR(nLag) else Correl = DBLARR(nLag) 2:if Double eq 0 then Correl = FLTARR(Xsize[1], nLag) else Correl = DBLARR(Xsize[1], nLag) 3:if Double eq 0 then Correl = FLTARR(Xsize[1], Xsize[2], nLag) $ else Correl = DBLARR(Xsize[1], Xsize[2], nLag) 4:if Double eq 0 then Correl = FLTARR(Xsize[1], Xsize[2], Xsize[3], nLag) $ else Correl = DBLARR(Xsize[1], Xsize[2], Xsize[3], nLag) endcase if KEYWORD_SET(Covariance) eq 0 then begin ;Compute Cross Correlation. for k = 0, nLag-1 do begin if Lag[k] ge 0 then BEGIN case NDim of 1:Correl[k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / $ sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $ TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan)) 2:Correl[*, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / $ sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $ TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan)) 3:Correl[*, *, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / $ sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $ TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan)) 4:Correl[*, *, *, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / $ sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $ TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan)) endcase ENDIF else BEGIN case NDim of 1:Correl[k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / $ sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $ TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan)) 2:Correl[*, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / $ sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $ TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan)) 3:Correl[*, *, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / $ sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $ TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan)) 4:Correl[*, *, *, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / $ sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $ TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan)) endcase ENDELSE endfor endif else begin ;Compute Cross Covariance. for k = 0, nLag-1 do begin if Lag[k] ge 0 then BEGIN case NDim of 1:Correl[k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / nT 2:Correl[*, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / nT 3:Correl[*, *, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / nT 4:Correl[*, *, *, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / nT ENDCASE ENDIF else BEGIN case NDim of 1:Correl[k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / nT 2:Correl[*, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / nT 3:Correl[*, *, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / nT 4:Correl[*, *, *, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / nT ENDCASE ENDELSE endfor endelse if Double eq 0 then RETURN, FLOAT(Correl) else $ RETURN, Correl END