octave_packages/tsa-4.2.4/mvar.m

   1 function [ARF,RCF,PE,DC,varargout] = mvar(Y, Pmax, Mode);
   2 % MVAR estimates parameters of the Multi-Variate AutoRegressive model
   3 %
   4 %    Y(t) = Y(t-1) * A1 + ... + Y(t-p) * Ap + X(t);
   5 % whereas
   6 %    Y(t) is a row vecter with M elements Y(t) = y(t,1:M)
   7 %
   8 % Several estimation algorithms are implemented, all estimators
   9 % can handle data with missing values encoded as NaNs.
  10 %
  11 %       [AR,RC,PE] = mvar(Y, p);
  12 %       [AR,RC,PE] = mvar(Y, p, Mode);
  13 %
  14 % with
  15 %       AR = [A1, ..., Ap];
  16 %
  17 % INPUT:
  18 %  Y     Multivariate data series
  19 %  p     Model order
  20 %  Mode  determines estimation algorithm
  21 %
  22 % OUTPUT:
  23 %  AR    multivariate autoregressive model parameter
  24 %  RC    reflection coefficients (= -PARCOR coefficients)
  25 %  PE    remaining error variances for increasing model order
  26 %          PE(:,p*M+[1:M]) is the residual variance for model order p
  27 %
  28 % All input and output parameters are organized in columns, one column
  29 % corresponds to the parameters of one channel.
  30 %
  31 % Mode determines estimation algorithm.
  32 %  1:  Correlation Function Estimation method using biased correlation function estimation method
  33 %               also called the 'multichannel Yule-Walker' [1,2]
  34 %  6:  Correlation Function Estimation method using unbiased correlation function estimation method
  35 %
  36 %  2:  Partial Correlation Estimation: Vieira-Morf [2] using unbiased covariance estimates.
  37 %               In [1] this mode was used and (incorrectly) denominated as Nutall-Strand.
  38 %               Its the DEFAULT mode; according to [1] it provides the most accurate estimates.
  39 %  5:  Partial Correlation Estimation: Vieira-Morf [2] using biased covariance estimates.
  40 %               Yields similar results than Mode=2;
  41 %
  42 %  3:  buggy: Partial Correlation Estimation: Nutall-Strand [2] (biased correlation function)
  43 %  9:  Partial Correlation Estimation: Nutall-Strand [2] (biased correlation function)
  44 %  7:  Partial Correlation Estimation: Nutall-Strand [2] (unbiased correlation function)
  45 %  8:  Least Squares w/o nans in Y(t):Y(t-p)
  46 % 10:  ARFIT [3]
  47 % 11:  BURGV [4]
  48 % 13:  modified BURGV -
  49 % 14:  modified BURGV [4]
  50 % 15:  Least Squares based on correlation matrix
  51 % 22: Modified Partial Correlation Estimation: Vieira-Morf [2,5] using unbiased covariance estimates.
  52 % 25: Modified Partial Correlation Estimation: Vieira-Morf [2,5] using biased covariance estimates.
  53 %
  54 % 90,91,96: Experimental versions
  55 %
  56 %    Imputation methods:
  57 %  100+Mode:
  58 %  200+Mode: forward, past missing values are assumed zero
  59 %  300+Mode: backward, past missing values are assumed zero
  60 %  400+Mode: forward+backward, past missing values are assumed zero
  61 % 1200+Mode: forward, past missing values are replaced with predicted value
  62 % 1300+Mode: backward, 'past' missing values are replaced with predicted value
  63 % 1400+Mode: forward+backward, 'past' missing values are replaced with predicted value
  64 % 2200+Mode: forward, past missing values are replaced with predicted value + noise to prevent decaying
  65 % 2300+Mode: backward, past missing values are replaced with predicted value + noise to prevent decaying
  66 % 2400+Mode: forward+backward, past missing values are replaced with predicted value + noise to prevent decaying
  67 %
  68 %
  69 %
  70
  71 % REFERENCES:
  72 %  [1] A. Schloegl, Comparison of Multivariate Autoregressive Estimators.
  73 %       Signal processing, Elsevier B.V. (in press).
  74 %       available at http://dx.doi.org/doi:10.1016/j.sigpro.2005.11.007
  75 %  [2] S.L. Marple 'Digital Spectral Analysis with Applications' Prentice Hall, 1987.
  76 %  [3] Schneider and Neumaier)
  77 %  [4] Stijn de Waele, 2003
  78 %  [5]  M. Morf, et al. Recursive Multichannel Maximum Entropy Spectral Estimation,
  79 %      IEEE trans. GeoSci. Elec., 1978, Vol.GE-16, No.2, pp85-94.
  80 %
  81 %
  82 % A multivariate inverse filter can be realized with
  83 %   [AR,RC,PE] = mvar(Y,P);
  84 %   e = mvfilter([eye(size(AR,1)),-AR],eye(size(AR,1)),Y);
  85 %
  86 % see also: MVFILTER, MVFREQZ, COVM, SUMSKIPNAN, ARFIT2
  87
  88 %       $Id: mvar.m 9609 2012-02-10 10:18:00Z schloegl $
  89 %       Copyright (C) 1996-2006,2011,2012 by Alois Schloegl <alois.schloegl@ist.ac.at>
  90 %       This is part of the TSA-toolbox. See also
  91 %       http://pub.ist.ac.at/~schloegl/matlab/tsa/
  92 %       http://octave.sourceforge.net/
  93 %       http://biosig.sourceforge.net/
  94 %
  95 %    This program is free software: you can redistribute it and/or modify
  96 %    it under the terms of the GNU General Public License as published by
  97 %    the Free Software Foundation, either version 3 of the License, or
  98 %    (at your option) any later version.
  99 %
 100 %    This program is distributed in the hope that it will be useful,
 101 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
 102 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 103 %    GNU General Public License for more details.
 104 %
 105 %    You should have received a copy of the GNU General Public License
 106 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 107
 108
 109 % Inititialization
 110 [N,M] = size(Y);
 111
 112 if nargin<2,
 113         Pmax=max([N,M])-1;
 114 end;
 115
 116 if iscell(Y)
 117         Pmax = min(max(N ,M ),Pmax);
 118         C    = Y;
 119 end;
 120 if nargin<3,
 121         % according to [1], and other internal validations this is in many cases the best algorithm
 122         Mode=2;
 123 end;
 124
 125 [C(:,1:M),n] = covm(Y,'M');
 126 PE(:,1:M)  = C(:,1:M)./n;
 127 if 0,
 128 elseif Mode==0;  % this method is broken
 129         fprintf('Warning MVAR: Mode=0 is broken.\n')
 130         C(:,1:M) = C(:,1:M)/N;
 131         F = Y;
 132         B = Y;
 133         PEF = C(:,1:M);  %?% PEF = PE(:,1:M);
 134         PEB = C(:,1:M);
 135         for K=1:Pmax,
 136                 [D,n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
 137                 D = D/N;
 138                 ARF(:,K*M+(1-M:0)) = D/PEB;
 139                 ARB(:,K*M+(1-M:0)) = D'/PEF;
 140
 141                 tmp        = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0))';
 142                 B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0))';
 143                 F(K+1:N,:) = tmp;
 144
 145                 for L = 1:K-1,
 146                         tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 147                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 148                         ARF(:,L*M+(1-M:0))   = tmp;
 149                 end;
 150
 151                 RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
 152                 RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
 153
 154                 PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
 155                 PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
 156                 PE(:,K*M+(1:M)) = PEF;
 157         end;
 158
 159
 160 elseif Mode==1, %%%%% Levinson-Wiggens-Robinson (LWR) algorithm using biased correlation function
 161         % ===== In [1,2] this algorithm is denoted 'Multichannel Yule-Walker' ===== %
 162         C(:,1:M) = C(:,1:M)/N;
 163         PEF = C(:,1:M);
 164         PEB = C(:,1:M);
 165
 166         for K=1:Pmax,
 167                 [C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
 168                 C(:,K*M+(1:M)) = C(:,K*M+(1:M))/N;
 169
 170                 D = C(:,K*M+(1:M));
 171                 for L = 1:K-1,
 172                         D = D - ARF(:,L*M+(1-M:0))*C(:,(K-L)*M+(1:M));
 173                 end;
 174                 ARF(:,K*M+(1-M:0)) = D / PEB;
 175                 ARB(:,K*M+(1-M:0)) = D'/ PEF;
 176                 for L = 1:K-1,
 177                         tmp                    = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 178                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 179                         ARF(:,L*M+(1-M:0))     = tmp;
 180                 end;
 181
 182                 RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
 183                 RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
 184
 185                 PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
 186                 PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
 187                 PE(:,K*M+(1:M)) = PEF;
 188         end;
 189
 190
 191 elseif Mode==6, %%%%% Levinson-Wiggens-Robinson (LWR) algorithm using unbiased correlation function
 192         C(:,1:M) = C(:,1:M)/N;
 193         PEF = C(:,1:M);  %?% PEF = PE(:,1:M);
 194         PEB = C(:,1:M);
 195
 196         for K = 1:Pmax,
 197                 [C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
 198                 C(:,K*M+(1:M)) = C(:,K*M+(1:M))./n;
 199                 %C{K+1} = C{K+1}/N;
 200
 201                 D = C(:,K*M+(1:M));
 202                 for L = 1:K-1,
 203                         D = D - ARF(:,L*M+(1-M:0))*C(:,(K-L)*M+(1:M));
 204                 end;
 205                 ARF(:,K*M+(1-M:0)) = D / PEB;
 206                 ARB(:,K*M+(1-M:0)) = D'/ PEF;
 207                 for L = 1:K-1,
 208                         tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 209                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 210                         ARF(:,L*M+(1-M:0))   = tmp;
 211                 end;
 212
 213                 RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
 214                 RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
 215
 216                 PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
 217                 PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
 218                 PE(:,K*M+(1:M)) = PEF;
 219         end;
 220
 221
 222 elseif Mode==2 %%%%% Partial Correlation Estimation: Vieira-Morf Method [2] with unbiased covariance estimation
 223         %===== In [1] this algorithm is denoted 'Nutall-Strand with unbiased covariance' =====%
 224         %C(:,1:M) = C(:,1:M)/N;
 225         F = Y;
 226         B = Y;
 227         %PEF = C(:,1:M);
 228         %PEB = C(:,1:M);
 229         PEF = PE(:,1:M);
 230         PEB = PE(:,1:M);
 231         for K = 1:Pmax,
 232                 [D,n] = covm(F(K+1:N,:),B(1:N-K,:),'M');
 233                 D = D./n;
 234
 235                 ARF(:,K*M+(1-M:0)) = D / PEB;
 236                 ARB(:,K*M+(1-M:0)) = D'/ PEF;
 237
 238                 tmp        = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0)).';
 239                 B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0)).';
 240                 F(K+1:N,:) = tmp;
 241
 242                 for L = 1:K-1,
 243                         tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 244                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 245                         ARF(:,L*M+(1-M:0))   = tmp;
 246                 end;
 247
 248                 RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
 249                 RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
 250
 251                 [PEF,n] = covm(F(K+1:N,:),F(K+1:N,:),'M');
 252                 PEF = PEF./n;
 253
 254                 [PEB,n] = covm(B(1:N-K,:),B(1:N-K,:),'M');
 255                 PEB = PEB./n;
 256
 257                 PE(:,K*M+(1:M)) = PEF;
 258         end;
 259
 260
 261 elseif Mode==5 %%%%% Partial Correlation Estimation: Vieira-Morf Method [2] with biased covariance estimation
 262         %===== In [1] this algorithm is denoted 'Nutall-Strand with biased covariance' ===== %
 263
 264         F = Y;
 265         B = Y;
 266         PEF = C(:,1:M);
 267         PEB = C(:,1:M);
 268         for K=1:Pmax,
 269                 [D,n]  = covm(F(K+1:N,:),B(1:N-K,:),'M');
 270                 %D=D/N;
 271
 272                 ARF(:,K*M+(1-M:0)) = D / PEB;
 273                 ARB(:,K*M+(1-M:0)) = D'/ PEF;
 274
 275                 tmp        = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0)).';
 276                 B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0)).';
 277                 F(K+1:N,:) = tmp;
 278
 279                 for L = 1:K-1,
 280                         tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 281                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 282                         ARF(:,L*M+(1-M:0))   = tmp;
 283                 end;
 284
 285                 RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
 286                 RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
 287
 288                 [PEB,n] = covm(B(1:N-K,:),B(1:N-K,:),'M');
 289                 %PEB = D/N;
 290
 291                 [PEF,n] = covm(F(K+1:N,:),F(K+1:N,:),'M');
 292                 %PEF = D/N;
 293
 294                 %PE(:,K*M+(1:M)) = PEF;
 295                 PE(:,K*M+(1:M)) = PEF./n;
 296         end;
 297
 298
 299 elseif 0, Mode==3 %%%%% Partial Correlation Estimation: Nutall-Strand Method [2] with biased covariance estimation
 300         %%% OBSOLETE because its buggy, use Mode=9 instead.
 301         % C(:,1:M) = C(:,1:M)/N;
 302         F = Y;
 303         B = Y;
 304         PEF = C(:,1:M);
 305         PEB = C(:,1:M);
 306         for K=1:Pmax,
 307                 [D, n]  = covm(F(K+1:N,:),B(1:N-K,:),'M');
 308                 D = D./N;
 309
 310                 ARF(:,K*M+(1-M:0)) = 2*D / (PEB+PEF);
 311                 ARB(:,K*M+(1-M:0)) = 2*D'/ (PEF+PEB);
 312
 313                 tmp        = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0)).';
 314                 B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0)).';
 315                 F(K+1:N,:) = tmp;
 316
 317                 for L = 1:K-1,
 318                         tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 319                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 320                         ARF(:,L*M+(1-M:0))   = tmp;
 321                 end;
 322
 323                 %RCF{K} = ARF{K};
 324                 RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
 325
 326                 [PEF,n] = covm(F(K+1:N,:),F(K+1:N,:),'M');
 327                 PEF = PEF./N;
 328
 329                 [PEB,n] = covm(B(1:N-K,:),B(1:N-K,:),'M');
 330                 PEB = PEB./N;
 331
 332                 %PE(:,K*M+(1:M)) = PEF;
 333                 PE(:,K*M+(1:M)) = PEF./n;
 334         end;
 335
 336
 337 elseif Mode==7 %%%%% Partial Correlation Estimation: Nutall-Strand Method [2] with unbiased covariance estimation
 338
 339         F = Y;
 340         B = Y;
 341         PEF = PE(:,1:M);
 342         PEB = PE(:,1:M);
 343         for K = 1:Pmax,
 344                 [D]  = covm(F(K+1:N,:),B(1:N-K,:),'M');
 345                 %D = D./n;
 346
 347                 ARF(:,K*M+(1-M:0)) = 2*D / (PEB+PEF);
 348                 ARB(:,K*M+(1-M:0)) = 2*D'/ (PEF+PEB);
 349
 350                 tmp        = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0)).';
 351                 B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0)).';
 352                 F(K+1:N,:) = tmp;
 353
 354                 for L = 1:K-1,
 355                         tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 356                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 357                         ARF(:,L*M+(1-M:0))   = tmp;
 358                 end;
 359
 360                 %RCF{K} = ARF{K};
 361                 RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
 362
 363                 [PEF] = covm(F(K+1:N,:),F(K+1:N,:),'M');
 364                 %PEF = PEF./n;
 365
 366                 [PEB] = covm(B(1:N-K,:),B(1:N-K,:),'M');
 367                 %PEB = PEB./n;
 368
 369                 %PE{K+1} = PEF;
 370                 PE(:,K*M+(1:M)) = PEF;
 371         end;
 372
 373
 374 elseif any(Mode==[3,9]) %%%%% Partial Correlation Estimation: Nutall-Strand Method [2] with biased covariance estimation
 375         %% same as 3 but with fixed normalization
 376         F = Y;
 377         B = Y;
 378         PEF = C(:,1:M);
 379         PEB = C(:,1:M);
 380         for K=1:Pmax,
 381                 [D, n]  = covm(F(K+1:N,:),B(1:N-K,:),'M');
 382
 383                 ARF(:,K*M+(1-M:0)) = 2*D / (PEB+PEF);
 384                 ARB(:,K*M+(1-M:0)) = 2*D'/ (PEF+PEB);
 385
 386                 tmp        = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0)).';
 387                 B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0)).';
 388                 F(K+1:N,:) = tmp;
 389
 390                 for L = 1:K-1,
 391                         tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 392                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 393                         ARF(:,L*M+(1-M:0))   = tmp;
 394                 end;
 395
 396                 %RCF{K} = ARF{K};
 397                 RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
 398
 399                 [PEF,nf] = covm(F(K+1:N,:),F(K+1:N,:),'M');
 400                 [PEB,nb] = covm(B(1:N-K,:),B(1:N-K,:),'M');
 401
 402                 %PE(:,K*M+(1:M)) = PEF;
 403                 PE(:,K*M+(1:M)) = PEF./nf;
 404         end;
 405
 406
 407 elseif Mode==4,  %%%%% Kay, not fixed yet.
 408         fprintf('Warning MVAR: Mode=4 is broken.\n')
 409
 410         C(:,1:M) = C(:,1:M)/N;
 411         K = 1;
 412         [C(:,M+(1:M)),n] = covm(Y(2:N,:),Y(1:N-1,:));
 413         C(:,M+(1:M)) = C(:,M+(1:M))/N;  % biased estimate
 414
 415         D = C(:,M+(1:M));
 416         ARF(:,1:M) = C(:,1:M)\D;
 417         ARB(:,1:M) = C(:,1:M)\D';
 418         RCF(:,1:M) = ARF(:,1:M);
 419         RCB(:,1:M) = ARB(:,1:M);
 420         PEF = C(:,1:M)*[eye(M) - ARB(:,1:M)*ARF(:,1:M)];
 421         PEB = C(:,1:M)*[eye(M) - ARF(:,1:M)*ARB(:,1:M)];
 422
 423         for K=2:Pmax,
 424                 [C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
 425                 C(:,K*M+(1:M)) = C(:,K*M+(1:M)) / N; % biased estimate
 426
 427                 D = C(:,K*M+(1:M));
 428                 for L = 1:K-1,
 429                         D = D - C(:,(K-L)*M+(1:M))*ARF(:,L*M+(1-M:0));
 430                 end;
 431
 432                 ARF(:,K*M+(1-M:0)) = PEB \ D;
 433                 ARB(:,K*M+(1-M:0)) = PEF \ D';
 434                 for L = 1:K-1,
 435                         tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 436                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 437                         ARF(:,L*M+(1-M:0)) = tmp;
 438                 end;
 439                 RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0)) ;
 440                 RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0)) ;
 441
 442                 PEF = PEF*[eye(M) - ARB(:,K*M+(1-M:0)) *ARF(:,K*M+(1-M:0)) ];
 443                 PEB = PEB*[eye(M) - ARF(:,K*M+(1-M:0)) *ARB(:,K*M+(1-M:0)) ];
 444                 PE(:,K*M+(1:M))  = PEF;
 445         end;
 446
 447
 448 elseif Mode==8,  %%%%% Least Squares
 449         ix  = Pmax+1:N;
 450         y   = repmat(NaN,N-Pmax,M*Pmax);
 451         for k = 1:Pmax,
 452                 y(:,k*M+[1-M:0]) = Y(ix-k,:);
 453         end;
 454         ix2 = ~any(isnan([Y(ix,:),y]),2);
 455         ARF = Y(ix(ix2),:)'/y(ix2,:)';
 456         PE  = covm(Y(ix,:) - y*ARF');
 457         RCF = zeros(M,M*Pmax);
 458
 459
 460 elseif Mode==10,  %%%%% ARFIT
 461         try
 462                 RCF = [];
 463                 [w, ARF, PE] = arfit(Y, Pmax, Pmax, 'sbc', 'zero');
 464         catch
 465                 ARF = zeros(M,M*Pmax);
 466                 RCF = ARF;
 467         end;
 468
 469
 470 elseif Mode==11,  %%%%% de Waele 2003
 471         %%% OBSOLETE
 472         warning('MVAR: mode=11 is obsolete use Mode 13 or 14!');
 473         [pc,R0] = burgv(reshape(Y',[M,1,N]),Pmax);
 474         try,
 475                 [ARF,ARB,Pf,Pb,RCF,RCB] = pc2parv(pc,R0);
 476         catch
 477                 [RCF,RCB,Pf,Pb] = pc2rcv(pc,R0);
 478                 [ARF,ARB,Pf,Pb] = pc2parv(pc,R0);
 479         end;
 480         ARF = -reshape(ARF(:,:,2:end),[M,M*Pmax]);
 481         RCF = -reshape(RCF(:,:,2:end),[M,M*Pmax]);
 482         PE = reshape(Pf,[M,M*(Pmax+1)]);
 483
 484
 485 elseif Mode==12,
 486         %%% OBSOLETE
 487         warning('MVAR: mode=12 is obsolete use Mode 13 or 14!');
 488         % this is equivalent to Mode==11
 489         [pc,R0] = burgv(reshape(Y',[M,1,N]),Pmax);
 490         [rcf,rcb,Pf,Pb] = pc2rcv(pc,R0);
 491
 492         %%%%% Convert reflection coefficients RC to autoregressive parameters
 493         ARF = zeros(M,M*Pmax);
 494         ARB = zeros(M,M*Pmax);
 495         for K = 1:Pmax,
 496                 ARF(:,K*M+(1-M:0)) = -rcf(:,:,K+1);
 497                 ARB(:,K*M+(1-M:0)) = -rcb(:,:,K+1);
 498                 for L = 1:K-1,
 499                         tmp                    = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 500                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 501                         ARF(:,L*M+(1-M:0))     = tmp;
 502                 end;
 503         end;
 504         RCF = -reshape(rcf(:,:,2:end),[M,M*Pmax]);
 505         PE  = reshape(Pf,[M,M*(Pmax+1)]);
 506
 507
 508 elseif Mode==13,
 509         % this is equivalent to Mode==11 but can deal with missing values
 510         %%%%%%%%%%% [pc,R0] = burgv_nan(reshape(Y',[M,1,N]),Pmax,1);
 511         % Copyright S. de Waele, March 2003 - modified Alois Schloegl, 2009
 512         I = eye(M);
 513
 514         sz = [M,M,Pmax+1];
 515         pc= zeros(sz); pc(:,:,1) =I;
 516         K = zeros(sz); K(:,:,1)  =I;
 517         Kb= zeros(sz); Kb(:,:,1) =I;
 518         P = zeros(sz);
 519         Pb= zeros(sz);
 520
 521         %[P(:,:,1)]= covm(Y);
 522         [P(:,:,1)]= PE(:,1:M);  % normalized
 523         Pb(:,:,1)= P(:,:,1);
 524         f = Y;
 525         b = Y;
 526
 527         %the recursive algorithm
 528         for i = 1:Pmax,
 529                 v = f(2:end,:);
 530                 w = b(1:end-1,:);
 531
 532                 %% normalized, unbiased
 533                 Rvv = covm(v); %Pfhat
 534                 Rww = covm(w); %Pbhat
 535                 Rvw = covm(v,w); %Pfbhat
 536                 Rwv = covm(w,v); % = Rvw', written out for symmetry
 537                 delta = lyap(Rvv*inv(P(:,:,i)),inv(Pb(:,:,i))*Rww,-2*Rvw);
 538
 539                 TsqrtS = chol( P(:,:,i))'; %square root M defined by: M=Tsqrt(M)*Tsqrt(M)'
 540                 TsqrtSb= chol(Pb(:,:,i))';
 541                 pc(:,:,i+1) = inv(TsqrtS)*delta*inv(TsqrtSb)';
 542
 543                 %The forward and backward reflection coefficient
 544                 K(:,:,i+1) = -TsqrtS *pc(:,:,i+1) *inv(TsqrtSb);
 545                 Kb(:,:,i+1)= -TsqrtSb*pc(:,:,i+1)'*inv(TsqrtS);
 546
 547                 %filtering the reflection coefficient out:
 548                 f = (v'+ K(:,:,i+1)*w')';
 549                 b = (w'+Kb(:,:,i+1)*v')';
 550
 551                 %The new R and Rb:
 552                 %residual matrices
 553                 P(:,:,i+1)  = (I-TsqrtS *pc(:,:,i+1) *pc(:,:,i+1)'*inv(TsqrtS ))*P(:,:,i);
 554                 Pb(:,:,i+1) = (I-TsqrtSb*pc(:,:,i+1)'*pc(:,:,i+1) *inv(TsqrtSb))*Pb(:,:,i);
 555         end %for i = 1:Pmax,
 556         R0 = PE(:,1:M);
 557
 558         %% [rcf,rcb,Pf,Pb] = pc2rcv(pc,R0);
 559         rcf  = zeros(sz); rcf(:,:,1)  = I;
 560         Pf   = zeros(sz); Pf(:,:,1)   = R0;
 561         rcb  = zeros(sz); rcb(:,:,1) = I;
 562         Pb   = zeros(sz); Pb(:,:,1)  = R0;
 563
 564         for p = 1:Pmax,
 565                 TsqrtPf = chol( Pf(:,:,p))'; %square root M defined by: M=Tsqrt(M)*Tsqrt(M)'
 566                 TsqrtPb= chol(Pb(:,:,p))';
 567                 %reflection coefficients
 568                 rcf(:,:,p+1) = -TsqrtPf *pc(:,:,p+1) *inv(TsqrtPb);
 569                 rcb(:,:,p+1)= -TsqrtPb*pc(:,:,p+1)'*inv(TsqrtPf );
 570                 %residual matrices
 571                 Pf(:,:,p+1)  = (I-TsqrtPf *pc(:,:,p+1) *pc(:,:,p+1)'*inv(TsqrtPf ))*Pf(:,:,p);
 572                 Pb(:,:,p+1) = (I-TsqrtPb*pc(:,:,p+1)'*pc(:,:,p+1) *inv(TsqrtPb))*Pb(:,:,p);
 573         end %for p = 2:order,
 574         %%%%%%%%%%%%%% end %%%%%%
 575
 576
 577         %%%%% Convert reflection coefficients RC to autoregressive parameters
 578         ARF = zeros(M,M*Pmax);
 579         ARB = zeros(M,M*Pmax);
 580         for K = 1:Pmax,
 581                 ARF(:,K*M+(1-M:0)) = -rcf(:,:,K+1);
 582                 ARB(:,K*M+(1-M:0)) = -rcb(:,:,K+1);
 583                 for L = 1:K-1,
 584                         tmp                    = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 585                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 586                         ARF(:,L*M+(1-M:0))     = tmp;
 587                 end;
 588         end;
 589         RCF = -reshape(rcf(:,:,2:end),[M,M*Pmax]);
 590         PE  = reshape(Pf,[M,M*(Pmax+1)]);
 591
 592
 593 elseif Mode==14,
 594         % this is equivalent to Mode==11 but can deal with missing values
 595         %%%%%%%%%%%%%% [pc,R0] = burgv_nan(reshape(Y',[M,1,N]),Pmax,2);
 596         % Copyright S. de Waele, March 2003 - modified Alois Schloegl, 2009
 597         I = eye(M);
 598
 599         sz = [M,M,Pmax+1];
 600         pc= zeros(sz); pc(:,:,1) =I;
 601         K = zeros(sz); K(:,:,1)  =I;
 602         Kb= zeros(sz); Kb(:,:,1) =I;
 603         P = zeros(sz);
 604         Pb= zeros(sz);
 605
 606         %[P(:,:,1),nn]= covm(Y);
 607         [P(:,:,1)]= C(:,1:M);   %% no normalization
 608         Pb(:,:,1)= P(:,:,1);
 609         f = Y;
 610         b = Y;
 611
 612         %the recursive algorithm
 613         for i = 1:Pmax,
 614                 v = f(2:end,:);
 615                 w = b(1:end-1,:);
 616
 617                 %% no normalization
 618                 [Rvv,nn] = covm(v); %Pfhat
 619                 [Rww,nn] = covm(w); %Pbhat
 620                 [Rvw,nn] = covm(v,w); %Pfbhat
 621                 [Rwv,nn] = covm(w,v); % = Rvw', written out for symmetry
 622                 delta = lyap(Rvv*inv(P(:,:,i)),inv(Pb(:,:,i))*Rww,-2*Rvw);
 623
 624                 TsqrtS = chol( P(:,:,i))'; %square root M defined by: M=Tsqrt(M)*Tsqrt(M)'
 625                 TsqrtSb= chol(Pb(:,:,i))';
 626                 pc(:,:,i+1) = inv(TsqrtS)*delta*inv(TsqrtSb)';
 627
 628                 %The forward and backward reflection coefficient
 629                 K(:,:,i+1) = -TsqrtS *pc(:,:,i+1) *inv(TsqrtSb);
 630                 Kb(:,:,i+1)= -TsqrtSb*pc(:,:,i+1)'*inv(TsqrtS);
 631
 632                 %filtering the reflection coefficient out:
 633                 f = (v'+ K(:,:,i+1)*w')';
 634                 b = (w'+Kb(:,:,i+1)*v')';
 635
 636                 %The new R and Rb:
 637                 %residual matrices
 638                 P(:,:,i+1)  = (I-TsqrtS *pc(:,:,i+1) *pc(:,:,i+1)'*inv(TsqrtS ))*P(:,:,i);
 639                 Pb(:,:,i+1) = (I-TsqrtSb*pc(:,:,i+1)'*pc(:,:,i+1) *inv(TsqrtSb))*Pb(:,:,i);
 640         end %for i = 1:Pmax,
 641         R0 = PE(:,1:M);
 642
 643         %%%%%%%%%% [rcf,rcb,Pf,Pb] = pc2rcv(pc,R0);
 644         sz = [M,M,Pmax+1];
 645         rcf  = zeros(sz); rcf(:,:,1)  = I;
 646         Pf   = zeros(sz); Pf(:,:,1)   = R0;
 647         rcb  = zeros(sz); rcb(:,:,1) = I;
 648         Pb   = zeros(sz); Pb(:,:,1)  = R0;
 649         for p = 1:Pmax,
 650                 TsqrtPf = chol( Pf(:,:,p))'; %square root M defined by: M=Tsqrt(M)*Tsqrt(M)'
 651                 TsqrtPb = chol(Pb(:,:,p))';
 652                 %reflection coefficients
 653                 rcf(:,:,p+1)= -TsqrtPf *pc(:,:,p+1) *inv(TsqrtPb);
 654                 rcb(:,:,p+1)= -TsqrtPb *pc(:,:,p+1)'*inv(TsqrtPf );
 655                 %residual matrices
 656                 Pf(:,:,p+1) = (I-TsqrtPf *pc(:,:,p+1) *pc(:,:,p+1)'*inv(TsqrtPf))*Pf(:,:,p);
 657                 Pb(:,:,p+1) = (I-TsqrtPb *pc(:,:,p+1)'*pc(:,:,p+1) *inv(TsqrtPb))*Pb(:,:,p);
 658         end %for p = 2:order,
 659         %%%%%%%%%%%%%% end %%%%%%
 660
 661         %%%%% Convert reflection coefficients RC to autoregressive parameters
 662         ARF = zeros(M,M*Pmax);
 663         ARB = zeros(M,M*Pmax);
 664         for K = 1:Pmax,
 665                 ARF(:,K*M+(1-M:0)) = -rcf(:,:,K+1);
 666                 ARB(:,K*M+(1-M:0)) = -rcb(:,:,K+1);
 667                 for L = 1:K-1,
 668                         tmp                    = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 669                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 670                         ARF(:,L*M+(1-M:0))     = tmp;
 671                 end;
 672         end;
 673         RCF = -reshape(rcf(:,:,2:end),[M,M*Pmax]);
 674         PE  = reshape(Pf,[M,M*(Pmax+1)]);
 675
 676
 677 elseif Mode==15,  %%%%% Least Squares
 678         %% FIXME
 679         for K=1:Pmax,
 680                 [C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
 681 %                C(K*M+(1:M),:) = C(K*M+(1:M),:)/N;
 682         end;
 683         ARF = C(:,1:M)'/C(:,M+1:end)';
 684         PE  = covm(C(:,1:M)' - ARF * C(:,M+1:end)');
 685         RCF = zeros(M, M*Pmax);
 686
 687
 688 elseif 0,
 689         %%%%% ARFIT and handling missing values
 690         % compute QR factorization for model of order pmax
 691         [R, scale]   = arqr(Y, Pmax, 0);
 692          % select order of model
 693         popt         = Pmax; % estimated optimum order
 694         np           = M*Pmax; % number of parameter vectors of length m
 695         % decompose R for the optimal model order popt according to
 696         %
 697         %   | R11  R12 |
 698         % R=|          |
 699         %   | 0    R22 |
 700         %
 701         R11   = R(1:np, 1:np);
 702         R12   = R(1:np, Pmax+1:Pmax+M);
 703         R22   = R(M*Pmax+1:Pmax+M, Pmax+1:Pmax+M);
 704         A     = (R11\R12)';
 705         % return covariance matrix
 706         dof   = N-Pmax-M*Pmax;                % number of block degrees of freedom
 707         PE    = R22'*R22./dof;        % bias-corrected estimate of covariance matrix
 708
 709         try
 710                 RCF = [];
 711                 [w, ARF, PE] = arfit(Y, Pmax, Pmax, 'sbc', 'zero');
 712         catch
 713                 ARF = zeros(M,M*Pmax);
 714                 RCF = ARF;
 715         end;
 716
 717
 718 elseif Mode==22 %%%%% Modified Partial Correlation Estimation: Vieira-Morf [2,5] using unbiased covariance estimates.
 719         %--------Initialization----------
 720         F = Y;
 721         B = Y;
 722         [PEF, n] = covm(Y(2:N,:),'M');
 723         PEF = PEF./n;
 724         [PEB, n] = covm(Y(1:N-1,:),'M');
 725         PEB = PEB./n;
 726         %---------------------------------
 727         for K = 1:Pmax,
 728             %---Update the estimated error covariances(15.89) in [2]---
 729             [PEFhat,n] = covm(F(K+1:N,:),'M');
 730             PEFhat = PEFhat./n;
 731
 732             [PEBhat,n] = covm(B(1:N-K,:),'M');
 733             PEBhat = PEBhat./n;
 734
 735             [PEFBhat,n] = covm(F(K+1:N,:),B(1:N-K,:),'M');
 736             PEFBhat = PEFBhat./n;
 737
 738             %---Compute estimated normalized partial correlation matrix(15.88)in [2]---
 739             Rho = inv(chol(PEFhat)')*PEFBhat*inv(chol(PEBhat));
 740
 741             %---Update forward and backward reflection coefficients (15.82) and (15.83) in [2]---
 742             ARF(:,K*M+(1-M:0)) = chol(PEF)'*Rho*inv(chol(PEB)');
 743             ARB(:,K*M+(1-M:0)) = chol(PEB)'*Rho'*inv(chol(PEF)');
 744
 745             %---------------------------------
 746             tmp        = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0)).';
 747             B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0)).';
 748             F(K+1:N,:) = tmp;
 749
 750             for L = 1:K-1,
 751                     tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 752                     ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 753                     ARF(:,L*M+(1-M:0))   = tmp;
 754             end;
 755
 756             RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
 757             RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
 758
 759             %---Update forward and backward error covariances (15.75) and (15.76) in [2]---
 760             PEF = (eye(M)-ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0)))*PEF;
 761             PEB = (eye(M)-ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0)))*PEB;
 762
 763             PE(:,K*M+(1:M)) = PEF;
 764         end
 765
 766
 767 elseif Mode==25 %%%%Modified Partial Correlation Estimation: Vieira-Morf [2,5] using biased covariance estimates.
 768         %--------Initialization----------
 769         F = Y;
 770         B = Y;
 771         [PEF, n] = covm(Y(2:N,:),'M');
 772         PEF = PEF./N;
 773         [PEB, n] = covm(Y(1:N-1,:),'M');
 774         PEB = PEB./N;
 775         %---------------------------------
 776         for K = 1:Pmax,
 777             %---Update the estimated error covariances(15.89) in [2]---
 778             [PEFhat,n] = covm(F(K+1:N,:),'M');
 779             PEFhat = PEFhat./N;
 780
 781             [PEBhat,n] = covm(B(1:N-K,:),'M');
 782             PEBhat = PEBhat./N;
 783
 784             [PEFBhat,n] = covm(F(K+1:N,:),B(1:N-K,:),'M');
 785             PEFBhat = PEFBhat./N;
 786
 787             %---Compute estimated normalized partial correlation matrix(15.88)in [2]---
 788             Rho = inv(chol(PEFhat)')*PEFBhat*inv(chol(PEBhat));
 789
 790             %---Update forward and backward reflection coefficients (15.82) and (15.83) in [2]---
 791             ARF(:,K*M+(1-M:0)) = chol(PEF)'*Rho*inv(chol(PEB)');
 792             ARB(:,K*M+(1-M:0)) = chol(PEB)'*Rho'*inv(chol(PEF)');
 793
 794             %---------------------------------
 795             tmp        = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0)).';
 796             B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0)).';
 797             F(K+1:N,:) = tmp;
 798
 799             for L = 1:K-1,
 800                     tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 801                     ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 802                     ARF(:,L*M+(1-M:0))   = tmp;
 803             end;
 804
 805             RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
 806             RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
 807
 808             %---Update forward and backward error covariances (15.75) and (15.76) in [2]---
 809             PEF = (eye(M)-ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0)))*PEF;
 810             PEB = (eye(M)-ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0)))*PEB;
 811
 812             PE(:,K*M+(1:M)) = PEF;
 813         end
 814
 815
 816
 817
 818 %%%%% EXPERIMENTAL VERSIONS %%%%%
 819
 820 elseif Mode==90;
 821         % similar to MODE=0
 822         %% not recommended
 823         C(:,1:M) = C(:,1:M)/N;
 824         F = Y;
 825         B = Y;
 826         PEF = PE(:,1:M);        %CHANGED%
 827         PEB = PE(:,1:M);        %CHANGED%
 828         for K=1:Pmax,
 829                 [D,n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
 830                 D = D/N;
 831                 ARF(:,K*M+(1-M:0)) = D/PEB;
 832                 ARB(:,K*M+(1-M:0)) = D'/PEF;
 833
 834                 tmp        = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0))';
 835                 B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0))';
 836                 F(K+1:N,:) = tmp;
 837
 838                 for L = 1:K-1,
 839                         tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 840                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 841                         ARF(:,L*M+(1-M:0))   = tmp;
 842                 end;
 843
 844                 RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
 845                 RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
 846
 847                 PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
 848                 PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
 849                 PE(:,K*M+(1:M)) = PEF;
 850         end;
 851
 852
 853 elseif Mode==91, %%%%% Levinson-Wiggens-Robinson (LWR) algorithm using biased correlation function
 854         % ===== In [1,2] this algorithm is denoted 'Multichannel Yule-Walker' ===== %
 855         % similar to MODE=1
 856         %% not recommended
 857         C(:,1:M) = C(:,1:M)/N;
 858         PEF = PE(:,1:M);        %CHANGED%
 859         PEB = PE(:,1:M);        %CHANGED%
 860
 861         for K=1:Pmax,
 862                 [C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
 863                 C(:,K*M+(1:M)) = C(:,K*M+(1:M))/N;
 864
 865                 D = C(:,K*M+(1:M));
 866                 for L = 1:K-1,
 867                         D = D - ARF(:,L*M+(1-M:0))*C(:,(K-L)*M+(1:M));
 868                 end;
 869                 ARF(:,K*M+(1-M:0)) = D / PEB;
 870                 ARB(:,K*M+(1-M:0)) = D'/ PEF;
 871                 for L = 1:K-1,
 872                         tmp                    = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 873                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 874                         ARF(:,L*M+(1-M:0))     = tmp;
 875                 end;
 876
 877                 RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
 878                 RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
 879
 880                 PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
 881                 PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
 882                 PE(:,K*M+(1:M)) = PEF;
 883         end;
 884
 885
 886 elseif Mode==96, %%%%% Levinson-Wiggens-Robinson (LWR) algorithm using unbiased correlation function
 887         % similar to MODE=6
 888         %% not recommended
 889         C(:,1:M) = C(:,1:M)/N;
 890         PEF = PE(:,1:M);        %CHANGED%
 891         PEB = PE(:,1:M);        %CHANGED%
 892
 893         for K = 1:Pmax,
 894                 [C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
 895                 C(:,K*M+(1:M)) = C(:,K*M+(1:M))./n;
 896
 897                 D = C(:,K*M+(1:M));
 898                 for L = 1:K-1,
 899                         D = D - ARF(:,L*M+(1-M:0))*C(:,(K-L)*M+(1:M));
 900                 end;
 901                 ARF(:,K*M+(1-M:0)) = D / PEB;
 902                 ARB(:,K*M+(1-M:0)) = D'/ PEF;
 903                 for L = 1:K-1,
 904                         tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
 905                         ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
 906                         ARF(:,L*M+(1-M:0))   = tmp;
 907                 end;
 908
 909                 RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
 910                 RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
 911
 912                 PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
 913                 PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
 914                 PE(:,K*M+(1:M)) = PEF;
 915         end;
 916
 917 elseif Mode<100,
 918        fprintf('Warning MVAR: Mode=%i not supported\n',Mode);
 919
 920 %%%%% IMPUTATION METHODS %%%%%
 921 else
 922
 923         Mode0 = rem(Mode,100);
 924         if ((Mode0 >= 10) && (Mode0 < 20)),
 925                 Mode0 = 1;
 926         end;
 927
 928 if 0,
 929
 930
 931 elseif Mode>=2400,  % forward and backward
 932 % assuming that past missing values are already IMPUTED with the prediction value + innovation process
 933 % no decaying
 934
 935         [ARF,RCF,PE2] = mvar(Y, Pmax, Mode0);
 936         c = chol( PE2 (:, M*Pmax+(1:M)));
 937
 938         Y1 = Y;
 939         Y1(1,isnan(Y1(1,:))) = 0;
 940         z  = [];
 941         for k = 2:size(Y,1)
 942                 [z1,z] = mvfilter(ARF,eye(M),Y1(k-1,:)',z);
 943                 ix = isnan(Y1(k,:));
 944                 z1 = z1 + (randn(1,M)*c)';
 945                 Y1(k,ix) = z1(ix);
 946         end;
 947
 948         Y2 = flipud(Y);
 949         [ARB,RCF,PE] = mvar(Y2, Pmax, Mode0);
 950         Y2(1,isnan(Y2(1,:))) = 0;
 951         z  = [];
 952         for k = 2:size(Y2,1)
 953                 [z2,z] = mvfilter(ARB,eye(M),Y2(k-1,:)',z);
 954                 ix = isnan(Y(size(Y,1)-k+1,:));
 955                 z2 = z2 + (randn(1,M)*c)';
 956                 Y2(k,ix) = (z2(ix)' + Y2(k,ix))/2;
 957         end;
 958         Y2 = flipud(Y2);
 959
 960         Z = (Y2+Y1)/2;
 961         Y(isnan(Y)) = Z(isnan(Y));
 962
 963         [ARF,RCF,PE] = mvar(Y, Pmax, rem(Mode,100));
 964
 965
 966 elseif Mode>=2300,  % backward prediction
 967 % assuming that past missing values are already IMPUTED with the prediction value + innovation process
 968 % no decaying
 969
 970         Y  = flipud(Y);
 971         [ARB,RCF,PE] = mvar(Y, Pmax, Mode0);
 972         c = chol(PE(:,M*Pmax+(1:M)));
 973         Y1 = Y;
 974         Y1(1,isnan(Y1(1,:))) = 0;
 975         z  = [];
 976         for k = 2:size(Y,1)
 977                 [z1,z] = mvfilter(ARB,eye(M),Y1(k-1,:)',z);
 978                 ix = isnan(Y1(k,:));
 979                 z1 = z1 + (randn(1,M)*c)';
 980                 Y1(k,ix) = z1(ix);
 981         end;
 982         Y1 = flipud(Y1);
 983         [ARF,RCF,PE] = mvar(Y1, Pmax, rem(Mode,100));
 984
 985
 986 elseif Mode>=2200,  % forward predictions,
 987 % assuming that past missing values are already IMPUTED with the prediction value + innovation process
 988 % no decaying
 989         [ARF,RCF,PE] = mvar(Y, Pmax, Mode0);
 990         c = chol(PE(:,M*Pmax+(1:M)));
 991         Y1 = Y;
 992         Y1(1,isnan(Y1(1,:))) = 0;
 993         z  = [];
 994         for k = 2:size(Y,1)
 995                 [z1,z] = mvfilter(ARF,eye(M),Y1(k-1,:)',z);
 996                 ix = isnan(Y1(k,:));
 997                 z1 = z1 + (randn(1,M)*c)';
 998                 Y1(k,ix) = z1(ix);
 999         end;
1000         [ARF,RCF,PE] = mvar(Y1, Pmax, rem(Mode,100));
1001
1002
1003 elseif Mode>=1400,  % forward and backward
1004 %assuming that past missing values are already IMPUTED with the prediction value
1005         [ARF,RCF,PE] = mvar(Y, Pmax, Mode0);
1006         Y1 = Y;
1007         Y1(1,isnan(Y1(1,:))) = 0;
1008         z  = [];
1009         for k = 2:size(Y,1)
1010                 [z1,z] = mvfilter(ARF,eye(M),Y1(k-1,:)',z);
1011                 ix = isnan(Y1(k,:));
1012                 Y1(k,ix) = z1(ix);
1013         end;
1014
1015         Y2 = flipud(Y);
1016         [ARB,RCF,PE] = mvar(Y2, Pmax, Mode0);
1017         Y2(1,isnan(Y2(1,:))) = 0;
1018         z  = [];
1019         for k = 2:size(Y2,1)
1020                 [z2,z] = mvfilter(ARB,eye(M),Y2(k-1,:)',z);
1021                 ix = isnan(Y2(k,:));
1022                 Y2(k,ix) = z2(ix)';
1023         end;
1024         Y2 = flipud(Y2);
1025
1026         Z = (Y2+Y1)/2;
1027         Y(isnan(Y)) = Z(isnan(Y));
1028
1029         [ARF,RCF,PE] = mvar(Y, Pmax, rem(Mode,100));
1030
1031
1032 elseif Mode>=1300,  % backward prediction
1033         Y  = flipud(Y);
1034         [ARB,RCF,PE] = mvar(Y, Pmax, Mode0);
1035         Y2 = Y;
1036         Y2(1,isnan(Y2(1,:))) = 0;
1037         z  = [];
1038         for k = 2:size(Y,1)
1039                 [z2,z] = mvfilter(ARB,eye(M),Y2(k-1,:)',z);
1040                 ix = isnan(Y2(k,:));
1041                 Y2(k,ix) = z2(ix);
1042         end;
1043         Y2 = flipud(Y2);
1044         [ARF,RCF,PE] = mvar(Y2, Pmax, rem(Mode,100));
1045
1046
1047 elseif Mode>=1200,  % forward predictions,
1048 %assuming that past missing values are already IMPUTED with the prediction value
1049         [ARF,RCF,PE] = mvar(Y, Pmax, Mode0);
1050         Y1 = Y;
1051         Y1(1,isnan(Y1(1,:))) = 0;
1052         z  = [];
1053         for k = 2:size(Y,1)
1054                 [z1,z] = mvfilter(ARF,eye(M),Y1(k-1,:)',z);
1055                 ix = isnan(Y1(k,:));
1056                 Y1(k,ix) = z1(ix);
1057         end;
1058         [ARF,RCF,PE] = mvar(Y1, Pmax, rem(Mode,100));
1059
1060
1061 elseif Mode>=400,  % forward and backward
1062 % assuming that 'past' missing values are ZERO
1063         [ARF,RCF,PE] = mvar(Y, Pmax, Mode0);
1064         Y1 = Y;
1065         Y1(isnan(Y)) = 0;
1066         Z1 = mvfilter([zeros(M),ARF],eye(M),Y1')';
1067         Y1(isnan(Y)) = Z1(isnan(Y));
1068
1069         Y  = flipud(Y);
1070         [ARB,RCF,PE] = mvar(Y, Pmax, Mode0);
1071         Y2 = Y;
1072         Y2(isnan(Y)) = 0;
1073         Z2 = mvfilter([zeros(M),ARB],eye(M),Y2')';
1074         Y2(isnan(Y)) = Z2(isnan(Y));
1075         Y2 = flipud(Y2);
1076
1077         [ARF,RCF,PE] = mvar((Y1+Y2)/2, Pmax, rem(Mode,100));
1078
1079
1080 elseif Mode>=300,  % backward prediction
1081 % assuming that 'past' missing values are ZERO
1082         Y  = flipud(Y);
1083         [ARB,RCF,PE] = mvar(Y, Pmax, Mode0);
1084         Y2 = Y;
1085         Y2(isnan(Y)) = 0;
1086         Z2 = mvfilter([zeros(M),ARB],eye(M),Y2')';
1087         Y2(isnan(Y)) = Z2(isnan(Y));
1088         Y2 = flipud(Y2);
1089
1090         [ARF,RCF,PE] = mvar(Y2, Pmax, rem(Mode,100));
1091
1092
1093 elseif Mode>=200,
1094 % forward predictions, assuming that past missing values are ZERO
1095         [ARF,RCF,PE] = mvar(Y, Pmax, Mode0);
1096         Y1 = Y;
1097         Y1(isnan(Y)) = 0;
1098         Z1 = mvfilter([zeros(M),ARF],eye(M),Y1')';
1099         Y1(isnan(Y)) = Z1(isnan(Y));
1100
1101         [ARF,RCF,PE] = mvar(Y1, Pmax, rem(Mode,100));
1102
1103
1104 elseif Mode>=100,
1105         [ARF,RCF,PE] = mvar(Y, Pmax, Mode0);
1106         Z1 = mvfilter(ARF,eye(M),Y')';
1107         Z1 = [zeros(1,M); Z1(1:end-1,:)];
1108         Y(isnan(Y)) = Z1(isnan(Y));
1109         [ARF,RCF,PE] = mvar(Y, Pmax, rem(Mode,100));
1110
1111
1112 end;
1113 end;
1114
1115
1116 if any(ARF(:)==inf),
1117 % Test for matrix division bug.
1118 % This bug was observed in LNX86-ML5.3, 6.1 and 6.5, but not in SGI-ML6.5, LNX86-ML6.5, Octave 2.1.35-40; Other platforms unknown.
1119 p = 3;
1120 FLAG_MATRIX_DIVISION_ERROR = ~all(all(isnan(repmat(0,p)/repmat(0,p)))) | ~all(all(isnan(repmat(nan,p)/repmat(nan,p))));
1121
1122 if FLAG_MATRIX_DIVISION_ERROR,
1123         %fprintf(2,'%%% Warning MVAR: Bug in Matrix-Division 0/0 and NaN/NaN yields INF instead of NAN.  Workaround is applied.\n');
1124         warning('MVAR: bug in Matrix-Division 0/0 and NaN/NaN yields INF instead of NAN.  Workaround is applied.');
1125
1126         %%%%% Workaround
1127         ARF(ARF==inf)=NaN;
1128         RCF(RCF==inf)=NaN;
1129 end;
1130 end;
1131
1132 %MAR   = zeros(M,M*Pmax);
1133 DC     = zeros(M);
1134 for K  = 1:Pmax,
1135 %       VAR{K+1} = -ARF(:,K*M+(1-M:0))';
1136         DC  = DC + ARF(:,K*M+(1-M:0))'.^2; %DC meausure [3]
1137 end;
1138