% esmf.m

function g=esmf(y,lam)

%  program to apply the "Exponential Smoothing" filter to time series
%  using a matrix inversion technique. The data is assumed to be
%  organized into a matrix "y" with the rows of y
%  being the observations and the columns being the series.

%  The value of the smoothing parameter is the second argument.

%  ESM filtering involves the solution of the difference equation:

%  y(t) =  - lam*g(t+1) + (1+2*lam)*g(t) -lam*g(t-1) 

%  with 
%  y(1) = (1+lambda)*g(1) - lam*g(2)
%  y(T) = (1+lambda)*g(T) - lam*g(T-1)

% If there are a smaller number of observations than series, then
% there is most likely a mistake: convert the vector/matrix of inputs:

oy = size(y);
ny = max(oy);  % length of time series

if (oy(1)<oy(2))
 y=y';
end

disp('Computing Exponential Smoothing Filtered Time Series with Matrix Inversion')
disp('Growth Component is Returned as g')

%  Strategy: Structure difference equation as a matrix equation:
%               M g = y
%            and then invert M.

M  = zeros(ny,ny);

d1=ones(ny-1,1);
d1=d1*lam;

d2=ones(ny,1);
d2=d2*(1+2*lam);
d2(1)=1+lam;
d2(ny)=1+lam;

M = diag(d1,1)+diag(d2)+diag(d1,-1)
M 
pause

g = inv(M)*y;

% convert if necessary

if (oy(1)<oy(2))
 g=g';
end

