%
% Time Series Seasonality Detection
%
%  - Box and Jenkins Lagged Autocorrelations
%
%        [ h, len, area ] = season( y, dt, flag )
%
% where  h    = plot handle
%        len  = vector of peak lag times
%        area = vector of peak areas
%        y    = data time series
%        dt   = time interval of data
%        flag = technique to use
%
% Bob Newell, February 1996
%
%------------------------------------------------------
%
function [ h, len, area ] = season( y, dt, flag )
%
n = length( y );
len = [];
area = [];
% number of lags
m = max( fix( n / 2 ), ( n - 100 ) );
%
if flag == 1 | flag == 3,
%------------------------------------------------------
% Box and Jenkins Lagged Autocorrelations
%
% time series properties
  ybar = mean( y );
  vary = cov( y );
% calculate for each lag
  rho = zeros( m, 1 );
  k = 0;
  area1 = 0;
  area2 = 0;
  for i = 1:m,
% calculate lagged autocorrelation (zero negatives)
    num = 0;
    for j = 1:n-i,
      num = num + ( y(j) - ybar ) * ( y(j+i) - ybar );
    end;
    rho(i) = max( num / vary, 0 );
% test print statement
%	disp( [ i*dt rho(i) ] );
% detect and calculate peak times
    if rho(i) == 0,
	  if area2 > 0,
		if k > 0,
		  len = [ len ( area1 / area2 ) ];
		  area = [ area area2 ];
	    end;
		k = k + 1;
	  end;
	  area1 = 0;
	  area2 = 0;
	else,
	  area1 = area1 + i * dt * rho(i);
	  area2 = area2 + rho(i);
	end;
  end;
% lets have a look
  h = figure;
  plot( (1:m)*dt, rho )
  title( 'CORRELOGRAM' )
  ylabel( 'Autocorrelation' )
  xlabel( 'Time Lag' )
end;
if flag == 2 | flag == 3,
%-----------------------------------------------------
% 
%
end;
%-----------------------------------------------------
% the end
