param_kroots.m

用于模拟时变非平稳的ARMA过程

function [CKmn, CK0n]= param_kroots(Cmn, initval)
% function [CKmn, CK0n]= param_kroots(Cmn, initval)
%   This file is part of the TFPM toolbox v0.9 (c)
%   [e-mail] and underlies the GPL.
% 
% Computes the Kamen roots (Kamen, LINALG, vol98, pp263, 1988) of
% a monic TVAR/MA system given through Cmn. The switch initval
% determines the starting values of the recursion:
% 's': stationary assumption,
% 'c': cyclic assumption,
% or a (M-1) vector of specified values. 
% TAKE CARE: ERROR PROPAGATION!! CHANGE!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;%tfpm;
MAR  =   2;
LAR  =   0;
MMA  =   3;
LMA  =   0;
N    =   8;
re_im= 'r';
mo_no= 'n';
tfpm_file_gen;
%-------------
Cml= Bml;
Cmn= Bmn;
initval= 'c';
initval= 's';
initval= [1:33];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(nargin==1)
   initval= 'c';
end;

% Dimensions:
NM= size(Cmn);
N= NM(1);
M= NM(2)-1;

CK0n= [];
for n= 0:N-1
   ck0= Cmn(n+1, 1);
   Cmn(n+1, :)= Cmn(n+1, :)/ck0;
   CK0n= [CK0n; ck0];
end;

if(~M)%% Multiplier has no roots!
   CKmn= zeros(N, M);
else
   if(M==1)%% Order-1 system has very simple roots!
      CKmn= -Cmn(:, 2);
   else

% Set initval value (STATIONARY!!):
      if(initval=='s')
         Initval= roots(Cmn(1, :));
%         Initval= Initval(1);
         Initval= Initval(M)*ones(M-1, 1);
      else

% Set initval value (CYCLIC!!):
         if(initval=='c')
            Initval= [];
            for mm= 1:M-1
               Init= roots(Cmn(N+1-mm, :));
%               Init= Init(1);
               Init= Init(M);
               Initval= [Init; Initval];
            end;
         else
            initval= initval(:);
            Initval= initval(1:M-1);
         end;%if(initval=='c')
      end;%if(initval=='s')

      if(M>2)%% The recursion
         PM= Initval;
         Emn= [];
         for n= 0:N-1
            eMM= 0;
            for MM= M-1:-1:1
               eMM= [eMM; double(( eMM(M-MM)-Cmn(n+1, 2+MM) )/PM(n+M-MM))];
            end;
            Emn= [Emn; [1; flipud(eMM(2:end))].'];
            PM= [PM; eMM(M)-Cmn(n+1, 2)];
         end;
         PM= PM(M:end);
         CKmn= [param_kroots(Emn, initval) PM];
      else%% The final step
         CKmn= [];
         for n= 0:N-1
            p1= Cmn(n+1, 3)/Initval;
            p2= -p1-Cmn(n+1, 2);
            Initval= p2;
            CKmn= [CKmn; [p1 p2]];
         end;
      end;%if(M>2)%% The recursion
   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(99);clf
hold on
plot(CKmn+j*eps*ones(size(CKmn)), '.');
n= -50:50;
theta= 2*n/100*pi;
plot(sin(theta), cos(theta)) 
axis([-1.425 1.425 -1.1 1.1])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%