MWFE
Computes least‑squares estimates of the multichannel Wiener filter coefficients for two mutually stationary multichannel time series.
Required Arguments
CXX — Array of size NCHX by NCHX by MLFIL containing the autocovariances of the input time series X. (Input)
CZX — Array of size NCHZ by NCHX by MLFIL containing the cross‑covariances between the desired output time series Z and the input time series X. (Input)
EPS — Lower bound for the normalized mean square error. (Input)
TRACE — Trace of the autocovariance matrix of the desired output time series Z at time lag zero. (Input)
LFIL — Length of the Wiener filter. (Output)
FIL — Array of size NCHZ by NCHX by MLFIL containing the multichannel Wiener filter coefficients. (Output)
ENMS — Vector of length MLFIL containing the normalized mean square error corresponding to each filter length. (Output)
Optional Arguments
NCHX — Number of input channels. (Input)
NCHX must be greater than or equal to one.
Default: NCHX = size (CXX,1).
MLFIL — Maximum length of the Wiener filter. (Input)
MLFIL must be greater than or equal to one.
Default: MLFIL = size (CXX,3).
LDCXX — Leading dimension of CXX exactly as specified in the dimension statement of the calling program. (Input)
LDCXX must be greater than or equal to NCHX.
Default: LDCXX = size (CXX,1).
MDCXX — Middle dimension of CXX exactly as specified in the dimension statement of the calling program. (Input)
MDCXX must be greater than or equal to NCHX.
Default: MDCXX = size (CXX,2).
NCHZ — Number of channels in desired output time series. (Input)
NCHZ must be greater than or equal to one.
Default: NCHZ = size (CZX,1).
LDCZX — Leading dimension of CZX exactly as specified in the dimension statement of the calling program. (Input)
LDCZX must be greater than or equal to NCHZ.
Default: LDCZX = size (CZX,1).
MDCZX — Middle dimension of CZX exactly as specified in the dimension statement of the calling program. (Input)
MDCZX must be greater than or equal to NCHX.
Default: MDCZX = size (CZX,2).
LDFIL — Leading dimension of FIL exactly as specified in the dimension statement of the calling program. (Input)
LDFIL must be greater than or equal to NCHZ.
Default: LDFIL = size (FIL,1).
MDFIL — Middle dimension of FIL exactly as specified in the dimension statement of the calling program. (Input)
MDFIL must be greater than or equal to NCHX.
Default: MDFIL = size (FIL,2).
FORTRAN 90 Interface
Generic: CALL MWFE (CXX, CZX, EPS, TRACE, LFIL, FIL, ENMS [, …])
Specific: The specific interface names are S_MWFE and D_MWFE.
FORTRAN 77 Interface
Single: CALL MWFE (NCHX, MLFIL, CXX, LDCXX, MDCXX, NCHZ, CZX, LDCZX, MDCZX, EPS, TRACE, LFIL, FIL, LDFIL, MDFIL ENMS)
Double: The double precision name is DMWFE.
Description
Routine MFFE computes least‑squares estimates of the multichannel Wiener filter coefficients for two mutually stationary multichannel time series.
Define the multichannel time series X by
X = (X1, X2, …, Xp)
where
Xj = (X1j, X2j, …, Xnj)T j = 1, 2, …, p
with p = NCHX. Similarly, define the multichannel time series Z by
Z = (Z1, Z2, …, Zq)
where
Zj = (Z1j, Z2j, …, Zmj)T j = 1, 2, …, q
with q = NCHZ. The columns of X and Z correspond to individual channels of multichannel time series and may be examined from a univariate perspective. The rows of X and Z correspond to observations of p‑variate and q‑variate time series and may be examined from a multivariate perspective. Note that an alternative characterization of a multivariate time series X considers the columns of X to be observations of a p‑variate time series with the rows of X containing univariate time series. For example, see Priestley (1981, page 692) and Fuller (1976, page 14).
Let
be the row vector containing the means of the channels of X. In particular,
where for j = 1, 2, …, p
Let
be similarly defined. In what follows, assume the channels of both X and Z have been centered about their respective means
Suppose the desired output is the multichannel time series Z defined by the model
where
and ɸk is the array of dimension p × q containing the Wiener filter coefficients
for k = 1, …, K. The array ɸk is the (k + 1)-st level of the 3‑dimensional array FIL.
The filter coefficients are computed by solving a set of normal equations. The algorithm utilizes the block Toeplitz (or Töplitz) matrix structure of these equations and is given by Robinson (1967, pages 238–246). In particular, the required input consists of the multichannel autocovariance matrices ΣX, ΣZ, and the multichannel cross‑covariance matrix ΣZX. The routine MCCF may be used to estimate these covariance matrices.
Note that successively longer filters are estimated until either the normalized mean square error is less than EPS or the filter length K = LFIL equals MLFIL. The normalized mean square error is defined by
where tr ΣZ(0) = TRACE is the trace of the multichannel autocorrelation coefficient of the desired output at lag zero. The values of Qk for the successive filters of length k = 1, 2, …, K are contained in ENMS.
Comments
1. Workspace may be explicitly provided, if desired, by use of M2FE/DM2FE. The reference is:
CALL M2FE (NCHX, MLFIL, CXX, LDCXX, MDCXX, NCHZ, CZX, LDCZX, MDCZX, EPS, TRACE, LFIL, FIL, LDFIL, MDFIL, ENMS, IWK, WK)
The additional arguments are as follows:
IWK — Work vector of length NCHX.
WK — Work vector of length NCHX * NCHX * (2 * MLFIL + 12) + NCHZ.
2. The length of the filter is determined by the arguments EPS and MLFIL. Iteration to a longer filter stops when either the normalized mean square error ENMS is less than EPS, or the filter reaches the maximum allowable length, MLFIL.
3. The routine MCCF may be used to obtain the input arguments CXX, CZX, and TRACE. For TRACE, routine MCCF may be used to obtain the autocovariances of the desired output series Z. In particular, TRACE = ZVAR(1) + … + ZVAR(NCHZ).
4. For a given lag k, the multichannel cross‑covariance coefficient between Z and X is defined as the array of size NCHZ by NCHX whose elements are the single‑channel crosscovariance coefficients CZX(i, j, k).
Example
Consider the Wölfer Sunspot Data (Box and Jenkins 1976, page 530) along with data on northern light activity and earthquake activity (Robinson 1967, page 204) to be a three‑channel time series. Routine MWFE applied to these data determines the following Wiener filter:
USE GDATA_INT
USE SCOPY_INT
USE MCCF_INT
USE MWFE_INT
USE UMACH_INT
IMPLICIT NONE
INTEGER LDCXX, LDCZX, LDFIL, LDX, LDZ, MAXLAG, MDCXX, &
MDCZX, MDFIL, MLFIL, NCHANX, NCHANZ, NOBSX, NOBSZ
REAL SSUM
PARAMETER (MLFIL=3, NCHANX=3, NCHANZ=3, NOBSX=99, NOBSZ=99, &
LDCXX=NCHANX, LDCZX=NCHANZ, LDFIL=NCHANX, LDX=NOBSX, &
LDZ=NOBSZ, MAXLAG=MLFIL-1, MDCXX=NCHANX, MDCZX=NCHANX, &
MDFIL=NCHANZ)
!
INTEGER I, J, K, LFIL, NCOL, NOUT, NROW
REAL CVXX(LDCXX,MDCXX,-MAXLAG:MAXLAG), CVXX1(3,3,3), &
CVZX(LDCZX,MDCZX,-MAXLAG:MAXLAG), CVZX1(3,3,3), &
CXX(LDCXX,MDCXX,-MAXLAG:MAXLAG), &
CZX(LDCZX,MDCZX,-MAXLAG:MAXLAG), ENMS(MLFIL), EPS, &
FIL(LDFIL,MDFIL,MLFIL), R(0:2), RDATA(100,4), &
TRACE, X(LDX,NCHANX), XMEAN(NCHANX), XVAR(NCHANX), &
YMEAN, YVAR, Z(LDZ,NCHANZ), ZMEAN(NCHANZ), &
ZVAR(NCHANZ)
!
EQUIVALENCE (CVXX(1,1,0), CVXX1(1,1,1)), (CVZX(1,1,0), CVZX1(1,1, &
1))
! Wolfer sunspot numbers
! Northern lights activity
! Earthquake activity
CALL GDATA (8, RDATA, NROW, NCOL)
!
CALL SCOPY (NOBSX, RDATA(1:,2), 1, X(1:,1), 1)
CALL SCOPY (NOBSX, RDATA(1:,3), 1, X(1:,2), 1)
CALL SCOPY (NOBSX, RDATA(1:,4), 1, X(1:,3), 1)
!
CALL SCOPY (NOBSZ, RDATA(2:,2), 1, Z(1:,1), 1)
CALL SCOPY (NOBSZ, RDATA(2:,3), 1, Z(1:,2), 1)
CALL SCOPY (NOBSZ, RDATA(2:,4), 1, Z(1:,3), 1)
! Compute multichannel ACF of Z
CALL MCCF (Z, Z, MAXLAG, CXX, XVAR=XVAR, CCV=CVXX)
! Compute TRACE
TRACE = SSUM(NCHANZ,XVAR,1)
! Compute multichannel ACF of X
CALL MCCF (X, X, MAXLAG, CXX, CCV=CVXX)
! Compute multichannel CCF of Z and X
CALL MCCF (Z, X, MAXLAG, CZX, CCV=CVZX)
! Bound normalized MSE to be positive
EPS = 0.0
! Reverse the LAG direction and scale
! to agree with Robinson (1967)
R(0) = 99.D0
R(1) = 98.D0
R(2) = 97.D0
TRACE = TRACE*R(0)
DO 10 K=0, MAXLAG
DO 10 J=1, NCHANX
DO 10 I=1, NCHANX
CVXX(I,J,K) = CVXX(I,J,-K)*R(K)
CVZX(I,J,K) = CVZX(I,J,-K)*R(K)
10 CONTINUE
! Compute multichannel Wiener filter
CALL MWFE (CVXX1, CVZX1, EPS, TRACE, LFIL, FIL, ENMS)
! Print results
CALL UMACH (2, NOUT)
WRITE (NOUT,99994) LFIL
99994 FORMAT (1X, 'Number of filter coefficients, LFIL = ', I3)
DO 30 K=1, LFIL
WRITE (NOUT,99995) K
99995 FORMAT (//, 1X, 'Wiener filter coefficient of index K = ', I3)
DO 20 I=1, NCHANX
WRITE (NOUT,99996) (FIL(I,J,K),J=1,NCHANZ)
99996 FORMAT (1X, 3F12.4)
20 CONTINUE
30 CONTINUE
WRITE (NOUT,99997)
99997 FORMAT (//, 1X, 'Normalized mean square error')
WRITE (NOUT,99998)
99998 FORMAT (1X, ' K ENMS(K)')
DO 40 K=1, LFIL
WRITE (NOUT,99999) K, ENMS(K)
99999 FORMAT (1X, I2, 5X, F12.4)
40 CONTINUE
!
END
Output
Number of filter coefficients, LFIL = 3
Wiener filter coefficient of index K = 1
1.3834 0.0348 0.0158
0.0599 0.8266 0.0629
-0.1710 -0.0332 -0.1205
Wiener filter coefficient of index K = 2
-0.7719 -0.0183 -0.0318
-0.0040 -0.2328 0.0484
-0.2170 0.1912 -0.0667
Wiener filter coefficient of index K = 3
0.0516 0.0563 -0.0138
-0.0568 0.1084 -0.1731
0.0007 0.2177 -0.0152
Normalized mean square error
K ENMS(K)
1 0.6042
2 0.5389
3 0.5174
