%VARMA.M
%   VARMA computes the innovations representation of a recursive
%   linear model.
%   The model occurs in the state space form
%             x(t+1) = ao*x(t) + c*w(t+1)
%               y(t) = sy*x(t) + e(t+1)
%   where e(t) is a measurement error process
%              e(t+1) = D*e(t) + ee(t+1)
%   and where ee(t+1) is a vector white noise with covariance matrix
%   R.  We assume that ee(t+1) and w(t+1) are orthogonal at all
%   leads and lags.  This specification permits the outputs of
%   our model to be measured with error.
%   The matrices ao,c,sy,D,and R must be in memory.
%   The matrices D and R must each be m x m, where [m,n]=size(sy).
%   The program uses the Kalman filter to form the innovations representation
%              xh(t+1) = ao*xh(t) + k1*u(t)
%              yy(t+1) = GG*xh(t) + u(t)
%   where GG=[sy*ao-D*sy], yy(t)=y(t+1)-D*y(t), and u(t) is the innovation
%   u(t)=y(t+1) - E[y(t+1)| y(t),y(t-1),...].  The program creates k1, GG,
%   and the covariance matrix of u, which it stores in V1.
%   The program uses evardec.m to compute a decomposition of variance for the
%   innovations system (which is equivalent with the vector autoregression).
%   The decomposition of variance is stored in tt.
   Q=c*c';
   G=sy;
   [m,n]=size(sy);
   Vm=eye(m);
   Vmm=zeros(m);
   GG=G*ao-D*G;
   Q1=Q;
   W1=Q*G';
   R1=R+G*Q*G';
   [k1,s1]=kfilter(ao,GG,Q1,R1,W1);
   V1=R1+GG*s1*GG';
   [aa,bb,cc,dd]=series(ao,k1,GG,Vm,D,Vm,Vm,Vmm);
%  [t1,t2]=evardec(aa,bb,cc,-20,-1,V1)
%disp('Variance attributable to the ith innovation')
%   for i=1:m
%   mm=m+1;
%   tt=t2(:,(i-1)*mm+1:i*mm)
%   end
