function [yE,trE,pE,r,s,bE] = polper(x,d,dp,t,tt,alpha,fig)
%         ESTIMATING POLYNOMIAL TREND tr(t) AND PERIODIC COMPONENT p(t)
% using linear regression for the data model
% x(t) = y(t) + e with random disturbances e: N(0,sigma^2*I)
% where y(t)=tr(t)+p(t), tr(t)=polynomial, p(t)=[s_1,...,s_dp,s_1,...]
% -----------------------------------------------------------------------------
% CALLING:
% [yE,trE,pE,r,s,bE] = polper(x,d,dp,t,tt,alpha,fig)
% INPUT:
% x  ... is an n x 1 vector of data observed at (time) t
% d  ... degree of the polynomial tr(t), d>=0 (integer)
% dp ... is scalar denoting the number of points in one basic period of p(t)
% t  ... is an n x 1 vector of uniform time mesh points for observations x
%        (optional: default=[]=uniform mesh (1:n)')
% tt ... N x 1 vector of time mesh for the trend tr(tt)
%        (optional: default=[]=t)
% alpha ... is a percentual confidence risk (optional: default=[]=5%)
% fig   ... figure handle for the regression plot (optional: default=[]=no plot)
%
% OUTPUT:
% yE ... =[y,Es,Ep] is an n x 3 matrix where
%        y are smoothed values y(t)=tr(t)+p(t) on mesh t
%        [y-Es,y+Es] are confidence intervals of (100-alpha) %
%                    for smoothed values y(t)
%        [y-Ep,y+Ep] are confidence intervals of (100-alpha) %
%                    for predicted values y(t)+e
% trE... =[tr,Es,Ep] is an N x 3 matrix describing trend component tr(tt)
%         analogously but on mesh tt
% pE ... =[p, Es,Ep] is an n x 3 matrix describing periodic component
%         analogously on mesh t
% r  ...  n x 1 vector of residulas x-y where y are evaluated at tt=t,
%         r is an estimate of the random error e.
% s  ...  the estimate of the standard deviation sigma (scalar)
% bE ... =[b,E] is an (d+1) x 2 matrix, where
%         b are coefficient estimates of the regression polynomial
%         tr(t)=b(1)*(t-ts)^d+...+b(d)*(t-ts)+b(d+1) where ts=(t(1)+t(n))/2
%         [b-E,b+E] are their confidence intervals of (100-alpha) %

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

error(nargchk(3,7,nargin));
if nargin<4, t = []; end
if nargin<5, tt= []; end
if nargin<6, alpha=[]; end
if nargin<7, fig  =[]; end

p = d + 1;
n = length(x);
if isempty(t), t=(1:n).'; else t=reshape(t,length(t),1); end
if n~=length(t), error('x and t must be of equal size'); end
if isempty(tt), tt=t; else tt = reshape(tt,length(tt),1); end
N = length(tt);
ts=(t(1)+t(n))./2;

D      = zeros(n,p+dp-1);
D(:,p) = ones(n,1);
C      = ones(N,p+dp-1)./dp;
C(:,p) = ones(N,1);
for k=d:-1:1
    D(:,k) = (t -ts).*D(:,k+1);
    C(:,k) = (tt-ts).*C(:,k+1);
end
m = ceil(n./dp);   % number of periods
D1= kron(ones(m,1),[zeros(1,dp-1);eye(dp-1)]);
D(:,p+1:p+dp-1) = D1(1:n,:);
C = [D;C;[zeros(dp,p),[zeros(1,dp-1);eye(dp-1)]-ones(dp,dp-1)./dp]];
[bE,rE,sE,PbE,CbE] = linregr(x,D,C,[],alpha);
   yE  = [CbE(1:n,:),PbE(1:n,2)];
   trE = [CbE(n+1:n+N,:),PbE(n+1:n+N,2)];
   pE  = [CbE(n+N+1:n+N+dp,:),PbE(n+N+1:n+N+dp,2)];
   pE  = kron(ones(m,1),pE);
   pE  = pE(1:n,:);
   r   = rE(:,1);
   s   = sE(1);
   bE  = bE(1:p,:);

if ~isempty(fig)    % plot regression curves and confidence bands
   y = yE(:,1);
   plot(t,x,'co',t,y,'w',...
     t,y-yE(:,2),'r-',t,y+yE(:,2),'r-',t,y-yE(:,3),'g-',t,y+yE(:,3),'g-');
   title('Smoothed curve');
   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,x,'co',tt,trE(:,1),'w',...
        tt,trE(:,1)-trE(:,2),'r-',tt,trE(:,1)+trE(:,2),'r-');   % ,...
%        tt,trE(:,1)-trE(:,3),'g-', tt,trE(:,1)+trE(:,3),'g-');
   title('Trend component');
   pause(2); figure; pos(1:2)=[1.4,0.6].*pos0(1:2); set(gcf,'Position',pos);
   plot(t,pE(:,1),'w',...
        t,pE(:,1)-pE(:,2),'r-',t,pE(:,1)+pE(:,2),'r-');         % ,...
%        t,pE(:,1)-pE(:,3),'g-',t, pE(:,1)+pE(:,3),'g-');
   title('Periodic component');
   pause(2); figure; pos(1:2)=[1.6,0.4].*pos0(1:2); set(gcf,'Position',pos);
   plot(t,rE(:,1),'y');
   title('Error component');
end
