function [vdsdata,vseascomp] = gpdeseas(vdata,lag) 
% GPDESEAS, General Purpose DESEASonalization
%     For time series data, and integer lags,
%        this subtracts out "seasonal means"
%        i.e. mean of every "lag-th" observation.
%     For noninteger "lags", this uses smoothing to 
%        compute the seasonal component, from the scatterplot
%        where the x-coord. has been reduced modulo lag.
%     This is most sensible for a "detrended" series, so may
%        want to subtract a polynomial (linear?) fit, or a
%        large bandwidth smooth before applying this.
%
% Inputs:
%     vdata  -  n x 1 column vector of time series data
%     lag    -  number of time points in each cycle
%                    (e.g. 12 for monthly data, and annal effect)
%
% Outputs:
%     vdsdata    -  n x 1 column vector of deseasonalized data 
%                      (vdata with seasonal component subtracted out)
%     vseascomp  -  n x 1 column vector containing the 
%                      seasonal component  (optional)
%                      Note:  use parts 1:lag to one cycle of this

%    Copyright (c) J. S. Marron 1998



%  set initial parameters
%
n = length(vdata) ;
ncycle = floor(n / lag) ;
nleftover = n - (lag * ncycle) ;



if ncycle == 0 ;    %  then not enough data to deseasonalize at this lag,
                    %  so do trivial return
  vseascomp = zeros(n,1) ;
  vdsdata = vdata ;

else ;       %  then do real deseasonalization


  if abs(floor(lag) - lag) < 10^(-12) ;   %   if lag is an integer,
                                           %   then work with means

    if ncycle == 1 ;    %  then don't need to take first blockwise mean
      vseascomp = vdata(1:lag)' ;
    else ;     %  get first blockwise mean
      mdata = reshape(vdata(1:(lag*ncycle)),lag,ncycle)' ;
            %  note reshape fills out columns first
      vseascomp = mean(mdata)  ;
    end ;


    %  handle leftovers, if needed
    %
    if nleftover > 0 ;
      w = ncycle / (ncycle + 1) ;
            %  weight for old mean in updated mean
      vseascomp(1:nleftover) = w * vseascomp(1:nleftover) + ...
                                   (1 - w) * vdata((n-nleftover+1):n)' ;
            %  updated mean
    end ;


    %  create vector seasonal component
    %
    vseascomp = kron(ones(ncycle,1),vseascomp') ;
    if nleftover > 0 ;
      vseascomp = [vseascomp; vseascomp(1:nleftover)] ;
    end ;



  else ;    % then lag is not an integer, so use scatterplot smoothing

  
    vt = (1:n)' ;
           %  variable representing "time"
    vt = vt / lag ;
           %  map to [0,1)
    vt = vt - floor(vt) ;
           %  remainder after division
    vt = vt * lag ;
           %  map back to [0,lag)
    bevt = [(vt - lag); vt; (vt + lag)] ;
    bevdata = [vdata; vdata; vdata] ;
           %  boundary extended versions


    [vseasgrid, xgrid] = gpnpr([bevt, bevdata],0,[-lag; (2 * lag)]) ;
           %  seasonal component on boundary extended grid
           %  0 for RSW smooth


    vseascomp = interp1(xgrid,vseasgrid,vt,'spline') ;
           %  interpolate back to get seasonal component

  end ;



  vdsdata = vdata - vseascomp ;

end ;


