function [zHats,S_1,S_2,S_3,et,MSE] = triple_exp_smoothing(Z,w,S_0_1,S_0_2,S_0_3)
% TRIPLE_EXP_SMOOTHING - Performs triple exponential smoothing
% 
% Inputs: 
%   Z     : the original data sequence
%   w     : the relaxation factor (omega)
%   S_0_1 : the initial value for S_n^{[1]} or first smooth (optional)
%   S_0_2 : the initial value for S_n^{[2]} or second smooth (optional)
%   S_0_3 : the initial value for S_n^{[3]} or the third smooth (optional)
% 
% 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 
%   S_3   : the values of S_n^{[3]} or the third 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.
% 
%-----

alpha = 1 - w; 
n     = length(Z); 

if( nargin < 5 )
  % 1) compute the initial values of hat_beta_0, hat_beta_1, and hat_beta_2 ... 
  %    we will use all of the data to determine these parameters (Equation 3.44 in the book): 
  % 
  X          = [ ones(1,n); 1:n; (1:n).^2 ].';
  beta_hat   = (X' * X) \ ( X' * Z ); 
  beta_hat_0 = beta_hat(1); 
  beta_hat_1 = beta_hat(2);
  beta_hat_2 = beta_hat(3); 
  
  % 2) convert these initial values of hat_beta_0, hat_beta_1, and hat_beta_2 ... 
  %    into initial values of our smooths these are S_0^{[1]}, S_0^{[2]}, and S_0^{[3]} 
  %
  S_0_1 = beta_hat_0 - (     w / alpha ) * beta_hat_1 + (     w * ( 2 -     alpha ) / ( 2 * alpha^2 ) ) * beta_hat_2; 
  S_0_2 = beta_hat_0 - ( 2 * w / alpha ) * beta_hat_1 + ( 2 * w * ( 3 - 2 * alpha ) / ( 2 * alpha^2 ) ) * beta_hat_2; 
  S_0_3 = beta_hat_0 - ( 3 * w / alpha ) * beta_hat_1 + ( 3 * w * ( 4 - 3 * alpha ) / ( 2 * alpha^2 ) ) * beta_hat_2; 
end
  
% 3) iterate these expressions forward through all data (Equation 3.47 / 3.48 with l=1):
% 
et    = zeros(n,1); 
S_1   = zeros(n,1); 
S_2   = zeros(n,1); 
S_3   = zeros(n,1); 
zHats = zeros(n,1); 

z_hat       = ( 3 + ( alpha * ( 6 - 5 * alpha ) / ( 2 * w^2 ) ) + ( alpha^2 / (2 * w^2 ) ) ) * S_0_1  ...
            - ( 3 + ( alpha * ( 5 - 4 * alpha ) / ( 2 * w^2 ) ) + ( alpha^2 / (    w^2 ) ) ) * S_0_2  ...
            + ( 1 + ( alpha * ( 4 - 3 * alpha ) / ( 2 * w^2 ) ) + ( alpha^2 / (2 * w^2 ) ) ) * S_0_3;
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; 
S_3(1)      = ( 1 - w ) * S_2(1)  + w * S_0_3; 

for ii=2:length(Z),                                   % do the rest ... 
  z_hat     = ( 3 + ( alpha * ( 6 - 5 * alpha ) / ( 2 * w^2 ) ) + ( alpha^2 / (2 * w^2 ) ) ) * S_1(ii-1)  ...
            - ( 3 + ( alpha * ( 5 - 4 * alpha ) / ( 2 * w^2 ) ) + ( alpha^2 / (    w^2 ) ) ) * S_2(ii-1)  ...
            + ( 1 + ( alpha * ( 4 - 3 * alpha ) / ( 2 * w^2 ) ) + ( alpha^2 / (2 * w^2 ) ) ) * S_3(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);
  S_3(ii)   = ( 1 - w ) * S_2(ii) + w * S_3(ii-1);
end

MSE = mean( et.^2 );