function [sz,tr,msz]=szma(x,d,model,fig)
%      ESTIMATING SEASONAL COMPONENT BY THE MOVING AVERAGE METHOD
%========================================================================
% function [sz,tr,msz]=szma(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), d=2*q or d=2*q+1
% 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:
%        tr(q+1:n-q)=preliminary trend estimate obtained from x via
%                    moving average.
%        tr(1)=...=tr(q)=tr(q+1), tr(n-q)=tr(n-q+1)=...=tr(n).
% msz... additive or multiplicative trend correction to obtain
%        mean(sz)=0 for the additive model or
%        mean(sz)=1 for the multiplicative model

% (C) 1996-7, 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); q=fix(d/2); m=fix(n/d);
if n<2*d, error('x is shorter than two periods'); end
r = n-m*d;                      % length of residual period
x = reshape(x,1,n);
w = ones(1,d)./d;               % define weights for moving average
if d==2*q
   w([1,d+1])=0.5*[w(1);w(1)];
end
nw = length(w);
sz = conv(x,w);                 % compute moving average via convolution
tr = [sz(nw)*ones(1,q),sz(nw:n),sz(n)*ones(1,q)].';
if mult
   sz = [zeros(1,q),x(q+1:n-q)./sz(nw:n),zeros(1,q)];
else
   sz = [zeros(1,q),x(q+1:n-q)-sz(nw:n),zeros(1,q)];
end
sz = reshape(sz(r+1:n),d,m).'; % reshape periods row-by-row after removing
                               % the leading residual period
corr = (1:d)<=q-r | (1:d)>d-q; % indicator of moving average boundary effects
sz = sum(sz)./(m-corr);
msz = mean(sz);
if mult
   sz = sz ./ msz;
   tr = tr .* msz;
else
   sz = sz - msz;
   tr = tr + msz;
end
sz = [sz(d-r+1:d),kron(ones(1,m),sz)].';
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
