function [y,E,r,S,b] = exp2filt(x,S,ACS,t,tt,alpha,fig)
%                        EXPONENTIAL SMOOTHING OF ORDER 2
% for the data model x(t) = y(t) + e with random disturbances e: N(0,sigma^2*I)
% -----------------------------------------------------------------------------
% CALLING: [y,E,r,S,b] = exp2filt(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 15 matrix recording row-by-row the history of state vectors
%        from all preceding calls of exp2filt:
%        S(m,:)=[a(0),S(0),S2(0),...,dt] ...
%               ... current initial state vector at the time t(0) where
%        a(0) ... initial estimate of smoothing parameter a: 0.7<=a<1
%        S(0) ... estimate of statistics S  at time t(0) [Cipra, eq. (7.29)]
%        S2(0)... estimate of statistics S2 at time t(0) [Cipra, eq. (7.30)]
%        dt   ... spacing of mesh points t
%        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 15 matrix with the final state vector at time t(n)
%        [a(n),S(n),S2(n),...,dt] appended to S as last row.
%        NOTE: S=[] on input => S is 2 x 15 on output (m=1) containing the
%              estimated initial state (1-st row) and the final state (2-nd row).
% b  ... =[b1,b2,...,bn] is an 2 x n matrix where
%         b# are coefficient estimates of #-th regression line
%         p#(t)=b#(1)*(t-t(#))+b#(2) estimated up to the time t(#)
%
% REFERENCE: 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

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. 68-69)
   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
   b  = polyfit((1:n0).',x(1:n0),1);  % Preliminary estimate of regr. line
   b0 = b(2); b1 = b(1);
   y0 = (b0 + b1)*ones(size(a0));
                              % One-step-prediction of x(k) [Cipra, eq. (7.34)]
   S0   = b0-b1*a0./(1-a0);     % and of S(0)
   S2_0 = b0-b1*2*a0./(1-a0);   % and of S2(0) [Cipra, eq. (7.36)]
   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;
        S0   = (1-a0)*x(k) + a0.*S0;    % Cipra, eq. (7.29)
        S2_0 = (1-a0).*S0  + a0.*S2_0;  % Cipra, eq. (7.31)
        b0   = 2*S0 - S2_0;
        b1   = (1-a0).*(S0-S2_0)./a0;   % Cipra, eq. (7.33)
        y0   = b0 + b1;      % One-step-prediction of x(k) [Cipra, eq. (7.34)]
   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 3 may be more adequate.');
   end
   a0 = a0(J(1));                 % Optimal value of initial smoothing par. a
   b  = polyfit((1:n).',x,1);     % Initial estimate of regr. line
   S0   = b(2)-b(1)*a0./(1-a0);   % and of S(0)
   S2_0 = b(2)-b(1)*2*a0./(1-a0); % and of S2(0) [Cipra, eq. (7.36)]
   S = [a0,S0,S2_0,S0,S2_0,S0,S2_0,0,0,1,zeros(1,3),0,(t(n)-t(1))/(n-1)];
                                  % Set initial a,S0,S2_0,delta,k_1,k0,Dk,Ik,dt
   clear y0 a0 SSE J
elseif any(size(S,2)~=15),
       error('The state matrix S must have 15 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];
S02 = reshape(Sc(2:7),2,3);
S0  = S02(1,:); S2 = S02(2,:); delta = Sc(8);
k_1 = Sc(9); k0 = Sc(10); Dk = Sc(11:13); Ik = Sc(14); dt = Sc(15);

b0 = 2*S0 - S2;
b1 = (1-a).*(S0-S2)./a;    % Cipra, eq. (7.33)
y0 = b0+b1;                % Initial estimates y(0) for all a(:)

r0 = zeros(n,1); r1 = zeros(n,3); b = zeros(2,n);
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=indices of optimal a(:)
     Ik = Ik + r1(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 3 ',...
                 'may be more adequate.']);
        end
        S0 = S0(J(1))*ones(1,3);
        S2 = S2(J(1))*ones(1,3);
        y0 = ( b0(J(1)) + b1(J(1)) )*ones(1,3);
        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
     S0 = (1-a)*x(k) + a.*S0;  % Cipra, eq. (7.29)
     S2 = (1-a).*S0  + a.*S2;  % Cipra, eq. (7.31)
     b0 = 2*S0 - S2;
     b1 = (1-a).*(S0-S2)./a;   % Cipra, eq. (7.33)
     b(:,k) = [b1(1);b0(1)];
     y0 = b0 + b1;             % One-step prediction [Cipra, eq. (7.34)]
     yk = b0(1);               % y(k) based on x(t), t<=k [Cipra, eq. (7.35)]
     r0(k) = x(k) - yk;        % Residual
     delta = ((kk-1)*delta + abs(r1(k,1)))/kk;  % Recursion for delta: Cipra,pp.64
     a1 = a(1);
     p = ((1-a1)/(1+a1)^3)*[2*(1-a1)^2,2*(1-a1)*(1+3*a1),1+4*a1+5*a1^2];
     p(3) = p(3) + 1;  % See Cipra, eq.(7.50)
     if t_eq_tt        % Smoothed values on mesh t==tt
        y(k) = yk;
        d0   = 1.25*sqrt(polyval(p,0)/polyval(p,1));  % Cipra, eq.(7.50)
        E(k) = d0*u*delta;      %  Confidence interval [Cipra, eq. (7.19)]
     else            % Piecewise linear smoothed values on general mesh tt
        I = (tt>t0) & (tt<=(k<n)*t(k)+(k==n)*tt(N));
        tau  = (tt(I)-t(k))./dt;
        y(I) = b0(1)+b1(1)*tau;
        dtau = 1.25*sqrt(polyval(p,tau)/polyval(p,1));  % Cipra, eq.(7.50)
        E(I) = dtau*u*delta;    %  Confidence interval [Cipra, eq. (7.19)]
        t0 = t(k);
     end
end
if k0>Sc(10), disp('Smoothing parameter has been adjusted.'); end
S02 = [S0;S2];
S = [S;[a(1),S02(:).',delta,k_1+n,k0,Dk,Ik,dt]];
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
