invert_blocktri.m

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

function Ainv= invert_blocktri(A, B, C)
% function Ainv= invert_blocktri(A, B, C)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Huang/McColl, J. Phys. A: Math. Gen. 30 (1997) 7919--7933

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
M= 10;
N= 10;

A= [];
B= [];
C= [];

%for n= 1:N
%   A= [A -n*eye(M)];
%   B= [B eye(M)];
%   C= [C n*eye(M)];
%%   C= [C n*zeros(M)];
%end;

for n= 1:N
   A= [A randn(M)];
   B= [B randn(M)];
   C= [C randn(M)];
end;

A(:, 1:M)= zeros(M);
C(:, (N-1)*M+1:end)= zeros(M);

X= zeros(M*N);
for i= 1:N
   for j= 1:N
      if(i-1==j)
         X((i-1)*M+1:(i)*M, (j-1)*M+1:(j)*M)= A(:, (i-1)*M+1:(i)*M);
      end;
      if(i==j)
         X((i-1)*M+1:(i)*M, (j-1)*M+1:(j)*M)= B(:, (j-1)*M+1:(j)*M);
      end;
      if(j-1==i)
         X((i-1)*M+1:(i)*M, (j-1)*M+1:(j)*M)= C(:, (j-2)*M+1:(j-1)*M);
      end;
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M= size(A, 1);
N= size(A, 2)/M;

% Compute Z
Z= [eye(M) inv(C(:, 1:M)) * B(:, 1:M)];

for i= 2:N-1
   Z= [Z inv(C(:, (i-1)*M+1:(i)*M)) * ( B(:, (i-1)*M+1:(i)*M)*Z(:, (i-1)*M+1:(i)*M) - A(:, (i-1)*M+1:(i)*M)*Z(:, (i-2)*M+1:(i-1)*M) )];
end;


% Compute Y
Y= [inv(A(:, (N-1)*M+1:end)) * B(:, (N-1)*M+1:end) eye(M)];
for j= N-1:-1:2
   Y= [inv(A(:, (j-1)*M+1:(j)*M)) * ( B(:, (j-1)*M+1:(j)*M)*Y(:, 1:M) - C(:, (j-1)*M+1:(j)*M)*Y(:, M+1:2*M) ) Y];
end;

% Main Block Diag
Ainv= zeros(M*N);
Ainv(1:M, 1:M)= inv( B(:, 1:M) - C(:, 1:M)*Y(:, 1*M+1:2*M)*inv(Y(:, 1:M)) );
for j= 2:N-1
   Ainv((j-1)*M+1:j*M, (j-1)*M+1:j*M)= inv( B(:, (j-1)*M+1:(j)*M) - A(:, (j-1)*M+1:(j)*M)*Z(:, (j-2)*M+1:(j-1)*M)*inv(Z(:, (j-1)*M+1:(j)*M)) - C(:, (j-1)*M+1:(j)*M)*Y(:, (j)*M+1:(j+1)*M)*inv(Y(:, (j-1)*M+1:(j)*M)));
end;
Ainv((N-1)*M+1:end, (N-1)*M+1:end)= inv( B(:, (N-1)*M+1:end) - A(:, (N-1)*M+1:end)*Z(:, (N-2)*M+1:(N-1)*M)*inv(Z(:, (N-1)*M+1:N*M)) );

% The Rest
for i= 1:N
   for j= 1:N
      if(j<i)
         Ainv((i-1)*M+1:(i)*M, (j-1)*M+1:(j)*M)= ...
            -Y(:, (i-1)*M+1:(i)*M)*inv(Y(:, (i-2)*M+1:(i-1)*M))*Ainv((i-2)*M+1:(i-1)*M, (j-1)*M+1:(j)*M);
      end;
   end;
end;
for i= N:-1:1
   for j= N:-1:1
      if(i<j)
         Ainv((i-1)*M+1:(i)*M, (j-1)*M+1:(j)*M)= ...
            -Z(:, (i-1)*M+1:(i)*M)*inv(Z(:, (i  )*M+1:(i+1)*M))*Ainv((i  )*M+1:(i+1)*M, (j-1)*M+1:(j)*M);
      end;
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
norm(Ainv-inv(X))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%