% SPECTRS.M
%  Uses the Tiao-Grupe formula to calculate the spectral density matrix
%  for an (n X n) proces y(t) that has hidden periodicity p.
%  The spectral density of y is given by
%   Sy(z) = Q(z)*H*(I - Ah*z^p)^(-1)*Ch*Ch'*(I - Ah*z^(-p))'^(-1)H'Q(z^(-1))
%  where z=exp(2*pi*j*i/T) for j=0,1, ..., T-1, and
%  where Ah, Ch, and H are from the stacked system
%
%       X(t+1) = Ah * X(t) + Ch* W(t+1)
%        Y(t)  = H * X(t)  ,
%   where X(t)' = [x(tp-p+1)', x(tp-p+2)', ..., x(tp)'],
%         Y(t)' = [y(tp-p+1)', y(tp-p+2)', ..., y(tp)']
%         W(t)' = [w(tp-p+1)', w(tp-p+2)', ..., w(tp)']
%   The program returns the spectral density matrices for frequencies
%   wj = 2*pi*j/T for j =0,1, ... ,T-1 in the matrices Sy0, Sy1, ... , SyT-1.
%   The spectral densities of the individual series are returned in
%   the matrix S.
%
%   SOLVES.M and SIMULS.M must be run first, and their outputs in memory.
%   nnc (the index of the constant) must be specified.  Edit the
%   file to specify T and the variables to be included in y(t).
%   SEE ALSO:  SPECTRAL.M, SOLVES.M, SIMULS.M, SEASLA10.M.
t0=clock
										
CC=Ch*Ch';
i=sqrt(-1);
T=2^7;                 %  edit this line to set T
										
for j=1:p
nj=int2str(j);
eval(['G',nj,'=[sc',nj,';si',nj,'(1,:)];'])  % edit this line to change observable
end
disp('delete rows and columns corresponding to constants')
AA=Ah;
										
for j=nnc:n-1:p*(n-1)
AA=cold(AA,j);
AA=rowd(AA,j);
CC=rowd(CC,j);
CC=cold(CC,j);
end
m1=max(size(AA));
I1=eye(m1);
for j=1:p
nj=int2str(j);
eval(['G',nj,'=cold(G',nj,',nnc);'])
end
										
   end
ns=2;
rg=ns;       %  ns is the number of series
ST=zeros(rg*(T/2+1),rg);
S=zeros(rg,T/2+1);
for j=0:T/2
   nj=int2str(j);
   ex=exp(2*pi*i*j/T);
   ns=2;                    %  ns is the number of series
   H=zeros(ns,(n-1)*p)';        %  edit to add more series
     nn=n-1;     %  Adjustment because constants have been removed
   rg=min(size(H));
 for h=1:p
   nh=int2str(h);
   eval(['H((h-1)*nn+1:h*nn,1:ns)=p^(-.5)*ex^(h-1)*G',nh,''';'])
  end
    eval(['ST(j*rg+1:(j+1)*rg,:)=',...
    '(H.''/(I1-AA*ex^p))*CC*(H.''/(I1-AA*ex^p))'';'])
  for k=1:ns
   S(k,j+1)=ST(j*rg+k,k);
   end

 end

et=etime(clock,t0)


%  Now plot the spectral densities
f=[0:T/2].*2*pi/T;
semilogy(f',real(S(1,1:T/2+1))),title('spectr of series 1'),pause
semilogy(f',real(S(2,1:T/2+1))),title('spectr of series 2'),pause
