%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Spectral analysis of a contrived time series
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Start by constructing a time series we'd like to evaluate
%    We are designing these so x has unit variance.

ntim = 1000;         %  The time series will have 1500 points

%  Use one of the following two routines - generate an AR1 OR an AR2
%  process.

%  AR1:
if false;
  x = zeros(ntim,1);   %  Initialize x
  x(1) = randn(1);     %  Initial condition - random

  rho1 = 0.8;           %  One lag autocorrelation
  for i = 2:ntim;
    x(i) = rho1*x(i-1)+sqrt(1-rho1^2)*randn(1);
  end
end

%  AR2:
if true;
  x = zeros(ntim,1);   %  Initialize x
  x(1) = randn(1);     %  Initial condition - random
  
  rho1 = 0.5;          %  One lag autocorrelation
  rho2 = -0.35;        %  Two lag autocorrelation

    %  Get parameters alpha1 and alpha2 via Yule-Walker equations
    fi = [rho1; rho2];
    Big_Fi = [1 rho1; rho1 1];    %  X has unit variance
  
    alph = Big_Fi \ fi;           %  Solve Yule-Walker equations
    alpha1 = alph(1); alpha2 = alph(2);
  
  %  Now, generate time series
  noiseamp = sqrt(1 - alpha1*rho1 - alpha2*rho2);

  for i = 3:ntim;
    x(i) = alpha1*x(i-1)+alpha2*x(i-2)+noiseamp*randn(1);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Let's analyze this using overlapping windows
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nfft = 128;              %  Length of window
noverlap = nfft/2;       %  Amount of overlap per segment

%window = boxcar(nfft);   %  Window we'll use ...
window = hanning(nfft);  %  Or this one
%window = hamming(nfft);  %  Or this one

Fs = 1/1;  %  Sampling rate:  observations per day (in the case of pentad
           %  data, we might use the following):
           %    Fs = 0.2 = (1 observation)/(5 days);
           %  Or, for monthly data:
           %    Fs = 12 = (12 observations)/(1 year);

DOF = 2*ntim/nfft;       %  Degrees of Freedom per spectral estimate

%  Calculate the spectrum
[p,f] = spectrum(x, nfft, noverlap, window, Fs);

%  Generate null hypotheses:
h = [0:(nfft/2)];

%  AR1:
if false;
  pAR = (1-rho1^2)./(1+rho1^2-2*rho1*cos(2*pi*h/nfft));
end

%  AR2:
if true;
  g = alpha1*(1-alpha2)*cos(2*pi*h/nfft)+alpha2*cos(4*pi*h/nfft);
  pAR = noiseamp^2./(1+alpha1^2+alpha2^2-2*g);
end

%  Normalize to make sure everything has equal variance (noting that the
%  first two spectral estimates are typically strongly affected by the
%  window).
pAR = pAR * sum(p(3:end,1))/sum(pAR(3:end));

%  Get f-value:
fsig = finv(0.95, DOF, 100000);  %  If you do not have the stats toolbox,
                                 %  then just look this up in a table.

%  Plot the spectrum in a number of ways:
figure(1); orient landscape; clf;

%  Linear
subplot(2,2,1);
  h = plot(f, p(:,1), '-k', f, pAR, '-r', f, pAR*fsig, '--r');
  set(h(1), 'linewidth', 2);  %  Make the spectrum a thicker line
  axis([0 0.5 -.1 14.1]);
  set(gca, 'XTick', [0:.1:0.5], 'YTick', 0:2:14);
  grid on;
  t1 = title('Spectrum:   Linear plot', 'fontsize', 12);
  xl = xlabel('Frequency (cycles per day)', 'fontsize', 10);
  yl = ylabel('Power', 'fontsize', 10);
  
%  Log-Linear
subplot(2,2,2);
  h = semilogy(f, p(:,1), '-k', f, pAR, '-r', f, pAR*fsig, '--r');
  set(h(1), 'linewidth', 2);  %  Make the spectrum a thicker line
  axis([0 0.5 -.1 14.1]);
  set(gca, 'XTick', [0:.1:0.5]);
  grid on;
  t1 = title('Spectrum:  Log-linear plot', 'fontsize', 12);
  xl = xlabel('Frequency (cycles per day)', 'fontsize', 10);
  yl = ylabel('log_1_0(Power)', 'fontsize', 10);
  
%  Log-log
subplot(2,2,3);
  h = loglog(f, p(:,1), '-k', f, pAR, '-r', f, pAR*fsig, '--r');
  set(h(1), 'linewidth', 2);  %  Make the spectrum a thicker line
  axis([0 0.5 -.1 14.1]);
  grid on;
  set(gca, 'XTick', 2.^[-6:-1], 'XTickLabel', 2.^[6:-1:1]);
  set(gca, 'XMinorGrid', 'off');
  t1 = title('Spectrum:  Log-log plot', 'fontsize', 12);
  xl = xlabel('Period (days per cycle)', 'fontsize', 10);
  yl = ylabel('log_1_0(Power)', 'fontsize', 10);

%  f*log(P) vs. log(f)
subplot(2,2,4);
  h = plot(log(f), f.*log(p(:,1)), '-k', ...
           log(f), f.*log(pAR'), '-r', ...
	   log(f), f.*log(pAR'*fsig), '--r');
  set(h(1), 'linewidth', 2);  %  Make the spectrum a thicker line
  axis([log(f(2)) log(0.5) -1.75 0.5]);
  grid on;
  set(gca, 'XTick', log(2.^[-6:-1]), 'XTickLabel', 2.^[6:-1:1]);
  set(gca, 'XMinorGrid', 'off');
  t1 = title('Spectrum:  f*ln(P) vs ln(f)', 'fontsize', 12);
  xl = xlabel('Period f (days per cycle)', 'fontsize', 10);
  yl = ylabel('f*ln(P)', 'fontsize', 10);