function [sz,tr]=szsmtr(x,d,model,fig)
%      ESTIMATING SEASONAL COMPONENT BY THE SMALL TREND METHOD
%========================================================================
% function [sz,tr]=szsmtr(x,d,model,fig)
% INPUT:
% x ... data vector n x 1 or 1 x n to be processed
% d ... period of the seasonal component (scalar)
% model ... ==0 ... additive model (default)
%           ~=0 ... multiplicative model
% fig   ... figure handle for the plot of results (optional: default=[]=no plot)
% OUTPUT:
% sz ... seasonal component of the same length as x (column vector)
% tr ... column vector of the same length as x:

% (C) 1997, V. Vesely, Masaryk University of Brno (Czech Republic)
%

error(nargchk(2,4,nargin));
if nargin<3, model=[]; end
if nargin<4, fig  =[]; end
if isempty(model), model=0; end
if min(size(x))~=1, error('data x is not a vector'); end

mult = model~=0;
if mult & any(x<=0), error('Data must be positive in multiplicative model'); end

n = length(x); m=ceil(n/d);
if n<2*d, error('x is shorter than two periods'); end

r=m*d-n;                       % length of overhead in the 1st period
x = reshape(x,n,1);
X = [zeros(r,1);x];
X = reshape(X,d,m);            % reshape periods column-by-column
tr = mean(X).';                % computing mean in each period
tr(1) = tr(1)*d/(d-r);         % correction due to overhead zeros
X = X.';                       % transpose reshaping of periods to row-by-row
corr = (1:d)<=r;               % indicator of overhead elements

% average centered periods to estimate seasonal component of one period:
if mult
   sz = sum(X./(tr*ones(1,d)))./(m-corr);
else
   X(1,1:r) = tr(1)*ones(1,r);
   sz = sum(X-tr*ones(1,d))./(m-corr);
end

% Make trend piecewise constant:
tr = kron(tr,ones(d,1));
tr = tr(r+1:n+r);

% Repeat periods of seasonal component:
sz = kron(ones(1,m),sz).';
sz = sz(r+1:n+r);
if ~isempty(fig)    % plot results
   t = 1:n;
   plot(t,x,'co',t,tr,'y');
   title('Trend component');
   pos0 = get(fig,'Position'); pos=pos0;
   pause(2); figure; pos(1:2)=[1.2,0.8].*pos0(1:2); set(gcf,'Position',pos);
   plot(t,sz,'y');,...
   title('Seasonal component');
   pause(2); figure; pos(1:2)=[1.4,0.6].*pos0(1:2); set(gcf,'Position',pos);
   if mult
      plot(t,x./(tr.*sz),'y');
   else
      plot(t,x-tr-sz,'y');
   end
   title('Error component');
end
