function [zHats,S_1,S_2,et,MSE] = double_exp_smoothing(Z,w,S_0_1,S_0_2)
% DOUBLE_EXP_SMOOTHING - Performs double exponential smoothing
% 
% Inputs: 
%   Z  : the original data sequence
%   w  : the relaxation factor (omega)
% 
% Outputs: 
%   zHats : the double exponentially smoothed output 
%   S_1   : the values of S_n^{[1]} or first smooth
%   S_2   : the values of S_n^{[2]} or second smooth 
%   et    : the one step-ahead residuals 
%  MSE    : the mean square error   
%
% Written by:
% -- 
% John L. Weatherwax                2009-04-21
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

n = length(Z);
  
if( nargin<4 )
  % 1) compute the initial values of hat_beta_1 and hat_beta_0 ... use all of the data to determine these parameters (Equation 3.44 in the book): 
  % 
  t_vect     = (1:n) - 0.5 * ( n+1 );
  beta_hat_1 = 12 * dot( t_vect, Z ) / ( n^3 - n );
  
  z_bar      = mean(Z); 
  beta_hat_0 = z_bar - beta_hat_1 * ( (n+1)/2 ); 

  % 2) convert these initial values of hat_beta_1 and hat_beta_0 ... into initial values of our smooth (Equation 3.43 in the book): 
  %
  S_0_1 = beta_hat_0 - (     w / ( 1 - w ) ) * beta_hat_1;           % these are S_0^{[1]}, and S_0^{[2]} 
  S_0_2 = beta_hat_0 - ( 2 * w / ( 1 - w ) ) * beta_hat_1; 
end

% 3) iterate these expressions forward through all data (Equation 3.38 with l=1):
% 
et    = zeros(n,1); 
S_1   = zeros(n,1); 
S_2   = zeros(n,1); 
zHats = zeros(n,1); 

z_hat       = ( 2 + ( 1 - w ) / w ) * S_0_1    - ( 1 + ( 1 - w ) / w ) * S_0_2; 
zHats(1)    = z_hat; 
et(1)       = Z(1) - z_hat; 
S_1(1)      = ( 1 - w ) * Z(1)    + w * S_0_1;        % do the first update by hand ... 
S_2(1)      = ( 1 - w ) * S_1(1)  + w * S_0_2; 

for ii=2:length(Z),                                   % do the rest ... 
  z_hat     = ( 2 + ( 1 - w ) / w ) * S_1(ii-1) - ( 1 + ( 1 - w ) / w ) * S_2(ii-1); 
  zHats(ii) = z_hat; 
  et(ii)    = Z(ii) - z_hat; 
  S_1(ii)   = ( 1 - w ) * Z(ii)   + w * S_1(ii-1);
  S_2(ii)   = ( 1 - w ) * S_1(ii) + w * S_2(ii-1);
end

MSE = mean( et.^2 );