function [y,E,r,S] = exp1filt(x,S,ACS,t,tt,alpha,fig)
%              FILTERING USING EXPONENTIAL SMOOTHING OF ORDER 1
% for the data model x(t) = y(t) + e with random disturbances e: N(0,sigma^2*I)
% -----------------------------------------------------------------------------
% CALLING: [y,E,r,S] = exp1filt(x,S,ACS,t,tt,alpha,fig)
% INPUT:
% x  ... is an n x 1 vector of data observed at (time) t
% S  ... is an m x 11 matrix recording row-by-row the history of state vectors
%        from all preceding calls of exp1filt:
%        S(m,:)=[a(0),y(0),...] where
%        a(0) ... initial estimate of smoothing parameter a: 0.7<=a<1
%        y(0) ... estimate of smoothed value at time t(0)
%        optional: default=[]=automatic initial estimate starting new history.
%                  Supply all available data for x in this case.
% ACS... level of Automatic Control of Smoothing parameter a, 0<=ACS<=Inf
%        (optional: default=[]=5 ... medium level)
%        increasing values suppress adjusting of a:
%        ACS=0    => maximal level of automatic control
%        ACS=Inf  => no control (a remains fixed)
% 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 smoothed values y
%        (optional: default=[]=t)
% alpha ... is a percentual confidence risk (optional: default=[]=5%)
% fig   ... figure handle for the regression plot (optional: default=[]=no plot)
%
% OUTPUT:
% y ...  is an N x 1 vector of smoothed (filtered) data
% E ...  is an N x 1 vector where
%        [y-E,y+E] are confidence intervals of (100-alpha) %
%                  for predicted values y(tt)+e
% r  ... n x 2 vector [r0,r1] of residulas
%        r0(t)=x(t)-y(t)     ... estimate error at time t
%        r1(t)=x(t)-y_t(t-1) ... one-step prediction error at time t
%        ( y_t(t-1) = y(t-1) = estimate of x(t) based on values x(k), k<=t-1 )
% S  ... is an (m+1) x 11 matrix with the final state vector at time t(n)
%        [a(n),y(n),...] appended to S as last row.
%        NOTE: S=[] on input => S is 2 x 11 on output (m=1) containing the
%              estimated initial state (1-st row) and the final state (2-nd row).
%
% REFERENCES: T. Cipra: Analyza casovych rad s aplikacemi v ekonomii,
%                       SNTL, Praha 1986 (in Czech)

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

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

if isempty(ACS), ACS = 5; end
if isempty(alpha), alpha = 5; end
alpha = 0.01*alpha;
if alpha==0.05
   u = 1.96;
else
   u = norminv(1-alpha/2);
end
Ec = u*1.25;

n = length(x);
x = reshape(x,n,1);

if isempty(t), t = (1:n).';
else, t = reshape(t,length(t),1);
end
[t,I] = sort(t);
x = x(I);
clear I;

if isempty(tt), tt=t;
else, tt = reshape(tt,length(tt),1); end
N  = length(tt);
tt = sort(tt);

t_eq_tt = (n==N);
if t_eq_tt, t_eq_tt = all(t==tt); end

if n~=length(t)
   error('The lengths of observed data ''x'' and time vector ''t'' differ ');
end

if isempty(S)            % COMPUTING THE INITIAL ESTIMATE (Cipra, pp. 61-62)
   if n<3, error('Not enough data in x to find the initial estimate of S'); end
   n0 = min(6,fix(n/2));
   a0 = 0.7:0.02:0.98;  % see Cipra: Tab. 7.3
%   a0 = 0.7:0.001:0.999; % Admissible range for values of smoothing parameter
   y0 = mean(x(1:n0))*ones(size(a0));  % Preliminary estimate of y(0)
   SSE = zeros(size(a0));
   for  k=1:n  % Compute sum of squares of one-step prediction error for all a0
        SSE = SSE  + (x(k)-y0).^2;
        y0  = (1-a0)*x(k) + a0.*y0; % Cipra, eq. (7.8)
   end
   m = min(SSE); J = find(m==SSE); % Find indices of smoothing parameters
                                   % minimizing SSE
   if J(1)==1 | J(length(J))==n    % If minimum achieved at the boundary
      disp(['WARNING: Initial estimate of the optimal ',...
                      'smoothing constant out of range [0.7,1).']);
      disp('Exponential smoothing of order 2 may be more adequate.');
   end
   y0 = mean(x);
   S = [a0(J(1)),y0,y0,y0,0,0,1,zeros(1,3),0];
                                  % Set initial a,y(0),delta,k_1,k0,Dk,Ik
   clear y0 a0 SSE J
elseif any(size(S,2)~=11),
       error('The state matrix S must have 11 columns');
end

% Get initial values a,a-0.05,a+0.05,
% the latter two are close alternatives for optimal adjusting of a:
Sc = S(size(S,1),:);
a  = [Sc(1),Sc(1)-0.05,Sc(1)+0.05];

y0 = Sc(2:4);         % Initial estimates y(0) for all a(:)
delta = Sc(5);        % Initial delta for confidence halfwidth: E(0)=Ec*delta
k_1 = Sc(6); k0 = Sc(7); Dk = Sc(8:10); Ik = Sc(11);

r0 = zeros(n,1); r1 = zeros(n,3);
y  = zeros(N,1); E = y; t0 = -Inf;
for  k=1:n
     kk = k_1+k;
     r1(k,:) = x(k) - y0;             % One-step prediction errors for all a(:)
     Dk = ( (kk-k0)*Dk+abs(r1(k,:)) )./(kk-k0+1);
                                      % Their mean values from time k0
                                      % at which a(:) was last adjusted
                                      % used as measure of optimality of a(:)
     m = min(Dk); J = find(m==Dk);    % J=indeces of optimal a(:)
     Ik = sum(r1(k0:k,1));
     if abs(Ik/Dk(1))>ACS & J(1)>1    % Adjust a(:) to more optimal values
        a = [a(J(1)),a(J(1))-0.05,a(J(1))+0.05];
        if a(1)<0.7 | a(1)>=1
           disp(['WARNING: Adjusted optimal value ',num2str(a(1)),...
                 ' of smoothing constant out of range [0.7,1).']);
           disp(['Higher ACS or exponential smoothing of order 2 ',...
                 'may be more adequate.']);
        end
        y0= y0(J(1))*ones(1,3); % Adjusted initial estimate
        r1(k,:) = x(k) - y0;    % and the one-step prediction error
        Dk = abs(r1(k,:));
        Ik = r1(k,1);
        k0 = kk;                % Save time of adjusting
     end
     y0    = (1-a)*x(k) + a.*y0;% Recursion for y(k) based on its prediction
                                % at time k-1: Cipra, eq. (7.8)
     r0(k) = x(k) - y0(1);      % Residual
     delta = ((kk-1)*delta + abs(r1(k,1)))/kk; % Recursion for delta: Cipra,pp.64
     if t_eq_tt      % Smoothed values on mesh t==tt
        y(k) = y0(1);
        E(k) = delta;
     else            % Piecewise constant smoothed values on general mesh tt
        I = (tt>t0) & (tt<=(k<n)*t(k)+(k==n)*tt(N));
        m = sum(I);
        y(I) = y0(1)*ones(m,1);
        E(I) = delta*ones(m,1);
        t0 = t(k);
     end
end
if k0>Sc(7), disp('Smoothing parameter has been adjusted.'); end
E = Ec*E;
S = [S;[a(1),y0(1),y0(1),y0(1),delta,k_1+n,k0,Dk,Ik]];
r = [r0,r1(:,1)];

if ~isempty(fig)    % plot regression curves and confidence bands
   plot(t,x,'go',tt,y,'w',tt,[y-E,y+E],'y:');
end
