mvar.m

时间序列分析的matlab程序

	[pc,R0] = burgv(reshape(Y',[M,1,N]),Pmax);
        [rcf,rcb,Pf,Pb] = pc2rcv(pc,R0);

        %%%%% Convert reflection coefficients RC to autoregressive parameters
        ARF = zeros(M,M*Pmax); 
        ARB = zeros(M,M*Pmax); 
        for K = 1:Pmax,
                ARF(:,K*M+(1-M:0)) = -rcf(:,:,K+1);	
                ARB(:,K*M+(1-M:0)) = -rcb(:,:,K+1);	
                for L = 1:K-1,
                        tmp                    = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
                        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));
                        ARF(:,L*M+(1-M:0))     = tmp;
                end;
        end;        
        RCF = -reshape(rcf(:,:,2:end),[M,M*Pmax]);
	PE  = reshape(Pf,[M,M*(Pmax+1)]);


        
%%%%% EXPERIMENTAL VERSIONS %%%%%

elseif Mode==83,  %%%%% de Waele 2003
	[rcf,rcb,pc,R0] = burgv_as1(reshape(Y',[M,1,N]),Pmax);
        %[ARF,ARB,Pf,Pb,RCF,RCB] = pc2parv(pc,R0);
 	RCF = -reshape(rcf(:,:,2:end),[M,M*Pmax]);
%	PE = reshape(Pf,[M,M*(Pmax+1)]);
        ARF = zeros(M,M*Pmax); 
        ARB = zeros(M,M*Pmax); 
        for K = 1:Pmax,
                ARF(:,K*M+(1-M:0)) = -rcf(:,:,K+1);	
                ARB(:,K*M+(1-M:0)) = -rcb(:,:,K+1);	
                for L = 1:K-1,
                        tmp                    = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
                        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));
                        ARF(:,L*M+(1-M:0))     = tmp;
                end;
        end;        
 

elseif Mode==84,  %%%%% de Waele 2003
	[rcf,rcb,pc,R0] = burgv_as2(reshape(Y',[M,1,N]),Pmax);
        %[ARF,ARB,Pf,Pb,RCF,RCB] = pc2parv(pc,R0);
 	RCF = -reshape(rcf(:,:,2:end),[M,M*Pmax]);
%	PE = reshape(Pf,[M,M*(Pmax+1)]);
        ARF = zeros(M,M*Pmax); 
        ARB = zeros(M,M*Pmax); 
        for K = 1:Pmax,
                ARF(:,K*M+(1-M:0)) = -rcf(:,:,K+1);	
                ARB(:,K*M+(1-M:0)) = -rcb(:,:,K+1);	
                for L = 1:K-1,
                        tmp                    = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
                        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));
                        ARF(:,L*M+(1-M:0))     = tmp;
                end;
        end;        


elseif Mode==85,  %%%%% de Waele 2003
	[rcf,rcb,pc,R0] = burgv_as1(reshape(Y',[M,1,N]),Pmax);
        %[ARF,ARB,Pf,Pb,RCF,RCB] = pc2parv(pc,R0);
 	RCF = -reshape(rcf(:,:,2:end),[M,M*Pmax]);
%	PE = reshape(Pf,[M,M*(Pmax+1)]);
        ARF = zeros(M,M*Pmax); 
        ARB = zeros(M,M*Pmax); 
        for K = 1:Pmax,
                ARF(:,K*M+(1-M:0)) = -rcf(:,:,K+1)';	
                ARB(:,K*M+(1-M:0)) = -rcb(:,:,K+1)';	
                for L = 1:K-1,
                        tmp                    = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
                        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));
                        ARF(:,L*M+(1-M:0))     = tmp;
                end;
        end;        
 

elseif Mode==86,  %%%%% de Waele 2003
	[rcf,rcb,pc,R0] = burgv_as2(reshape(Y',[M,1,N]),Pmax);
        %[ARF,ARB,Pf,Pb,RCF,RCB] = pc2parv(pc,R0);
 	RCF = -reshape(rcf(:,:,2:end),[M,M*Pmax]);
%	PE = reshape(Pf,[M,M*(Pmax+1)]);
        ARF = zeros(M,M*Pmax); 
        ARB = zeros(M,M*Pmax); 
        for K = 1:Pmax,
                ARF(:,K*M+(1-M:0)) = -rcf(:,:,K+1)';	
                ARB(:,K*M+(1-M:0)) = -rcb(:,:,K+1)';	
                for L = 1:K-1,
                        tmp                    = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
                        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));
                        ARF(:,L*M+(1-M:0))     = tmp;
                end;
        end;        

elseif Mode==90;  
	% similar to MODE=0
	%% not recommended
        C(:,1:M) = C(:,1:M)/N;
        F = Y;
        B = Y;
        PEF = PE(:,1:M);	%CHANGED%
        PEB = PE(:,1:M);	%CHANGED%
        for K=1:Pmax,
                [D,n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
	        D = D/N;
                ARF(:,K*M+(1-M:0)) = D/PEB;	
                ARB(:,K*M+(1-M:0)) = D'/PEF;	
                
                tmp        = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0))';
                B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0))';
                F(K+1:N,:) = tmp;
                
                for L = 1:K-1,
                        tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
                        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));
                        ARF(:,L*M+(1-M:0))   = tmp;
                end;
                
                RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
                RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
                
                PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
                PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
                PE(:,K*M+(1:M)) = PEF;        
        end;

        
elseif Mode==91, %%%%% Levinson-Wiggens-Robinson (LWR) algorithm using biased correlation function
	% ===== In [1,2] this algorithm is denoted "Multichannel Yule-Walker" ===== %
	% similar to MODE=1
	%% not recommended
        C(:,1:M) = C(:,1:M)/N;
        PEF = PE(:,1:M);	%CHANGED%
        PEB = PE(:,1:M);	%CHANGED%
        
        for K=1:Pmax,
                [C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
                C(:,K*M+(1:M)) = C(:,K*M+(1:M))/N;

                D = C(:,K*M+(1:M));
                for L = 1:K-1,
                        D = D - ARF(:,L*M+(1-M:0))*C(:,(K-L)*M+(1:M));
                end;
                ARF(:,K*M+(1-M:0)) = D / PEB;	
                ARB(:,K*M+(1-M:0)) = D'/ PEF;	
                for L = 1:K-1,
                        tmp                    = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
                        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));
                        ARF(:,L*M+(1-M:0))     = tmp;
                end;
                
                RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
                RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
                
                PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
                PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
                PE(:,K*M+(1:M)) = PEF;        
        end;

        
elseif Mode==96, %%%%% Levinson-Wiggens-Robinson (LWR) algorithm using unbiased correlation function
	% similar to MODE=6
	%% not recommended
        C(:,1:M) = C(:,1:M)/N;
        PEF = PE(:,1:M);	%CHANGED%
        PEB = PE(:,1:M);	%CHANGED%
        
        for K = 1:Pmax,
                [C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
                C(:,K*M+(1:M)) = C(:,K*M+(1:M))./n;

                D = C(:,K*M+(1:M));
                for L = 1:K-1,
                        D = D - ARF(:,L*M+(1-M:0))*C(:,(K-L)*M+(1:M));
                end;
                ARF(:,K*M+(1-M:0)) = D / PEB;	
                ARB(:,K*M+(1-M:0)) = D'/ PEF;	
                for L = 1:K-1,
                        tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
                        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));
                        ARF(:,L*M+(1-M:0))   = tmp;
                end;
                
                RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
                RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
                
                PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
                PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
                PE(:,K*M+(1:M)) = PEF;        
        end;

        
elseif Mode<100,
       fprintf('Warning MVAR: Mode=%i not supported\n',Mode);         

       
%%%%% IMPUTATION METHODS %%%%%
else

	Mode0 = rem(Mode,100); 	
	if ((Mode0>=10) & (Mode0<20)), 
		Mode0 = 1; 
	end;

if 0, 


elseif Mode>=2400,  % forward and backward
% assuming that past missing values are already IMPUTED with the prediction value + innovation process
% no decaying 

 	[ARF,RCF,PE2] = mvar(Y, Pmax, Mode0);	
 	c = chol( PE2 (:, M*Pmax+(1:M)));

	Y1 = Y; 
	Y1(1,isnan(Y1(1,:))) = 0; 
	z  = [];
	for k = 2:size(Y,1)
		[z1,z] = mvfilter(ARF,eye(M),Y1(k-1,:)',z);
		ix = isnan(Y1(k,:)); 
		z1 = z1 + (randn(1,M)*c)'; 
		Y1(k,ix) = z1(ix); 
	end; 

	Y2 = flipud(Y); 
 	[ARB,RCF,PE] = mvar(Y2, Pmax, Mode0);	
	Y2(1,isnan(Y2(1,:))) = 0; 
	z  = [];
	for k = 2:size(Y2,1)
		[z2,z] = mvfilter(ARB,eye(M),Y2(k-1,:)',z);
		ix = isnan(Y(size(Y,1)-k+1,:)); 
		z2 = z2 + (randn(1,M)*c)'; 
		Y2(k,ix) = (z2(ix)' + Y2(k,ix))/2; 
	end; 
	Y2 = flipud(Y2); 
	
	Z = (Y2+Y1)/2;
	Y(isnan(Y)) = Z(isnan(Y));
	
 	[ARF,RCF,PE] = mvar(Y, Pmax, rem(Mode,100));	


elseif Mode>=2300,  % backward prediction
% assuming that past missing values are already IMPUTED with the prediction value + innovation process
% no decaying 

	Y  = flipud(Y);
 	[ARB,RCF,PE] = mvar(Y, Pmax, Mode0);	
 	c = chol(PE(:,M*Pmax+(1:M))); 
	Y1 = Y; 
	Y1(1,isnan(Y1(1,:))) = 0; 
	z  = [];
	for k = 2:size(Y,1)
		[z1,z] = mvfilter(ARB,eye(M),Y1(k-1,:)',z);
		ix = isnan(Y1(k,:)); 
		z1 = z1 + (randn(1,M)*c)'; 
		Y1(k,ix) = z1(ix); 
	end; 	
	Y1 = flipud(Y1);
 	[ARF,RCF,PE] = mvar(Y1, Pmax, rem(Mode,100));	


elseif Mode>=2200,  % forward predictions, 
% assuming that past missing values are already IMPUTED with the prediction value + innovation process
% no decaying 
 	[ARF,RCF,PE] = mvar(Y, Pmax, Mode0);	
 	c = chol(PE(:,M*Pmax+(1:M))); 
	Y1 = Y; 
	Y1(1,isnan(Y1(1,:))) = 0; 
	z  = [];
	for k = 2:size(Y,1)
		[z1,z] = mvfilter(ARF,eye(M),Y1(k-1,:)',z);
		ix = isnan(Y1(k,:)); 
		z1 = z1 + (randn(1,M)*c)'; 
		Y1(k,ix) = z1(ix); 
	end; 	
 	[ARF,RCF,PE] = mvar(Y1, Pmax, rem(Mode,100));	


elseif Mode>=1400,  % forward and backward
%assuming that past missing values are already IMPUTED with the prediction value
 	[ARF,RCF,PE] = mvar(Y, Pmax, Mode0);	
	Y1 = Y; 
	Y1(1,isnan(Y1(1,:))) = 0; 
	z  = [];
	for k = 2:size(Y,1)
		[z1,z] = mvfilter(ARF,eye(M),Y1(k-1,:)',z);
		ix = isnan(Y1(k,:)); 
		Y1(k,ix) = z1(ix); 
	end; 

	Y2 = flipud(Y); 
 	[ARB,RCF,PE] = mvar(Y2, Pmax, Mode0);	
	Y2(1,isnan(Y2(1,:))) = 0; 
	z  = [];
	for k = 2:size(Y2,1)
		[z2,z] = mvfilter(ARB,eye(M),Y2(k-1,:)',z);
		ix = isnan(Y2(k,:)); 
		Y2(k,ix) = z2(ix)'; 
	end; 
	Y2 = flipud(Y2); 
	
	Z = (Y2+Y1)/2;
	Y(isnan(Y)) = Z(isnan(Y));
	
 	[ARF,RCF,PE] = mvar(Y, Pmax, rem(Mode,100));	


elseif Mode>=1300,  % backward prediction
	Y  = flipud(Y);
 	[ARB,RCF,PE] = mvar(Y, Pmax, Mode0);	
	Y2 = Y; 
	Y2(1,isnan(Y2(1,:))) = 0; 
	z  = [];
	for k = 2:size(Y,1)
		[z2,z] = mvfilter(ARB,eye(M),Y2(k-1,:)',z);
		ix = isnan(Y2(k,:)); 
		Y2(k,ix) = z2(ix); 
	end; 	
	Y2 = flipud(Y2);
 	[ARF,RCF,PE] = mvar(Y2, Pmax, rem(Mode,100));	


elseif Mode>=1200,  % forward predictions, 
%assuming that past missing values are already IMPUTED with the prediction value
 	[ARF,RCF,PE] = mvar(Y, Pmax, Mode0);	
	Y1 = Y; 
	Y1(1,isnan(Y1(1,:))) = 0; 
	z  = [];
	for k = 2:size(Y,1)
		[z1,z] = mvfilter(ARF,eye(M),Y1(k-1,:)',z);
		ix = isnan(Y1(k,:)); 
		Y1(k,ix) = z1(ix); 
	end; 	
 	[ARF,RCF,PE] = mvar(Y1, Pmax, rem(Mode,100));	


elseif Mode>=400,  % forward and backward
% assuming that "past" missing values are ZERO
 	[ARF,RCF,PE] = mvar(Y, Pmax, Mode0);	
	Y1 = Y; 
	Y1(isnan(Y)) = 0; 
	Z1 = mvfilter([zeros(M),ARF],eye(M),Y1')';
	Y1(isnan(Y)) = Z1(isnan(Y));

	Y  = flipud(Y);
 	[ARB,RCF,PE] = mvar(Y, Pmax, Mode0);	
	Y2 = Y; 
	Y2(isnan(Y)) = 0; 
	Z2 = mvfilter([zeros(M),ARB],eye(M),Y2')';
	Y2(isnan(Y)) = Z2(isnan(Y));
	Y2 = flipud(Y2);

 	[ARF,RCF,PE] = mvar((Y1+Y2)/2, Pmax, rem(Mode,100));	


elseif Mode>=300,  % backward prediction
% assuming that "past" missing values are ZERO
	Y  = flipud(Y);
 	[ARB,RCF,PE] = mvar(Y, Pmax, Mode0);	
	Y2 = Y; 
	Y2(isnan(Y)) = 0; 
	Z2 = mvfilter([zeros(M),ARB],eye(M),Y2')';
	Y2(isnan(Y)) = Z2(isnan(Y));
	Y2 = flipud(Y2);

 	[ARF,RCF,PE] = mvar(Y2, Pmax, rem(Mode,100));	


elseif Mode>=200,  
% forward predictions, assuming that past missing values are ZERO
 	[ARF,RCF,PE] = mvar(Y, Pmax, Mode0);	
	Y1 = Y; 
	Y1(isnan(Y)) = 0; 
	Z1 = mvfilter([zeros(M),ARF],eye(M),Y1')';
	Y1(isnan(Y)) = Z1(isnan(Y));

 	[ARF,RCF,PE] = mvar(Y1, Pmax, rem(Mode,100));	


elseif Mode>=100,  
 	[ARF,RCF,PE] = mvar(Y, Pmax, Mode0);	
	Z1 = mvfilter(ARF,eye(M),Y')';
	Z1 = [zeros(1,M); Z1(1:end-1,:)];
	Y(isnan(Y)) = Z1(isnan(Y)); 
 	[ARF,RCF,PE] = mvar(Y, Pmax, rem(Mode,100));	


end;
end;


if any(ARF(:)==inf),
% Test for matrix division bug. 
% 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.
p = 3;
FLAG_MATRIX_DIVISION_ERROR = ~all(all(isnan(repmat(0,p)/repmat(0,p)))) | ~all(all(isnan(repmat(nan,p)/repmat(nan,p))));

if FLAG_MATRIX_DIVISION_ERROR, 
	%fprintf(2,'### Warning MVAR: Bug in Matrix-Division 0/0 and NaN/NaN yields INF instead of NAN.  Workaround is applied.\n');
	warning('MVAR: bug in Matrix-Division 0/0 and NaN/NaN yields INF instead of NAN.  Workaround is applied.');

	%%%%% Workaround 
	ARF(ARF==inf)=NaN;
	RCF(RCF==inf)=NaN;
end;
end;	

%MAR   = zeros(M,M*Pmax);
DC     = zeros(M);
for K  = 1:Pmax,
%       VAR{K+1} = -ARF(:,K*M+(1-M:0))';
        DC  = DC + ARF(:,K*M+(1-M:0))'.^2; %DC meausure [3]
end;