%SPECTRAL  calculates the spectral density functions Sy(wj):
%
%          Sy(wj)=G(I-A exp(-i wj))^{-1}C C'(I-A exp(i wj))^{-1}'G'+
%                  (I-D exp(-i wj))^{-1} R (I-D exp(i wj))^{-1}'
%
%          j=0,...T-1, for the system:
%
%          x(t+1) = A x(t) + C e(t+1),  E e(t) e(t)'=I
%          y(t)   = G x(t) + v(t)
%          v(t+1) = D v(t) + u(t),      E u(t) u(t)'=R     
%
%          Thus, the program requires inputs: A,C,G,D,R,T  and returns
%      outputs: SA  and S where S(k,j)=SA(j*rg+k,k) , where
%      rg is the number of series, i.e. the number of rows of G
%      SA is an (rg*(T/2 +1)) X rg matrix.  The spectral density
%      matrix at frequency 2*pi*j/T is stored in
%      SA(j*rg+1:(j+1)*rg,:) for j=0,..., T/2.
%
%          NOTE: spectral.m is a script not a function m-file.
%
disp(' -------------------------------------------------------------- ')
disp(' You are now calculating the spectral density functions Sy(wj): ')
disp('                                                                ')
disp('  Sy(wj)=G(I-A exp(-i wj))^{-1}C C''(I-A exp(i wj))^{-1}''G''+  ')
disp('           (I-D exp(-i wj))^{-1} R (I-D exp(i wj))^{-1}''       ')
disp('                                                                ')
disp(' j=0,...T-1, for the system:                                    ')
disp('                                                                ')
disp('  x(t+1) = A x(t) + C e(t+1),  E e(t) e(t)''=I                  ')
disp('  y(t)   = G x(t) + v(t)                                        ')
disp('  v(t+1) = D v(t) + u(t),      E u(t) u(t)''=R                  ')
disp('                                                                ')
disp(' Thus, the program requires inputs: A,C,G,D,R,T  and returns    ')
disp(' outputs: SA  and S where S(k,j)=SA(j*rg+k,k) , where      ')
disp(' rg is the number of series, i.e. the number of rows of G ')
disp(' SA is an (rg*(T/2 +1)) X rg matrix.  The spectral density ')
disp(' matrix at frequency 2*pi*j/T is stored in       ')
disp(' SA(j*rg+1:(j+1)*rg,:) for j=0,..., T/2.     ')
disp(' -------------------------------------------------------------- ')
disp(' ')
[rg,cg]=size(G);
I1=eye(cg);
I2=eye(rg);
i=sqrt(-1);
CC=C*C';
SA=zeros(rg*(T/2+1),rg);S=zeros(T/2+1,rg)';
for j=0:T/2;
  ex=exp(2*pi*i*j/T);
    eval(['SA(j*rg+1:(j+1)*rg,:)=',...
        '(G/(I1-A*ex))*CC*(G/(I1-A*ex))''+((I2-D*ex)\R)*', ...
        'inv(I2-D*ex)'';'])
   for k=1:rg;
    eval(['S(k,j+1)=SA(j*rg+k,k);'])
  end;
end

t=[0:T/2]';
w=2*pi.*t/T;
