function trend = hpfit(y,lambda) ;
%
% Hodrick-Prescott filter with parameter lambda, where lambda is the
%  value of the Lagrange multiplier on the constraint.
% Literally, the routine minimizes
%  sum((y-trend).^2) + lambda*sum(trend(1:T-2) - 2*trend(2:T-1) + trend(3:T));
%  which is the definition given in Prescott, Quart. Rev. Minn. Fed., 1986.
% The time series will be taken to be in columns.  If there is more than
%  one column of data, the routine will HP-filter each column separately.
%
% credits: written by Eric T. Swanson, 1996. (with companion program hpfilter.m)
%

if (nargin<2); lambda=1600; end ; % Prescott used 1600, it's the default.
if (isempty(y)); trend = y; return; end ;

[nrows,ncols] = size(y) ;

% Use the Matlab diag routine to construct the matrix A quickly:
first = lambda*ones(nrows-2,1) ;
second = -4*lambda*ones(nrows-1,1) ;
third = (1+6*lambda)*ones(nrows,1) ;

A = diag(third) + diag(second,1) + diag(second,-1) ...
      + diag(first,2) + diag(first,-2) ;

A(1,1:3) = [1+lambda, -2*lambda, lambda] ;
A(2,1:4) = [-2*lambda, 1+5*lambda, -4*lambda, lambda] ;
A(nrows-1,nrows-3:nrows) = [lambda, -4*lambda, 1+5*lambda, -2*lambda] ; 
A(nrows,nrows-2:nrows) = [lambda, -2*lambda, 1+lambda] ;

% If y isn't missing any data, the filtering is easy:
if (all(isfinite(y))) ;
  trend = A \ y ;
else ; % otherwise, zero out the ~isfinite data and ignore them in the calcs:
  for i=1:ncols ;
    good = isfinite(y(:,i)) ;
    Aloc = A - diag(~good) ;
    trend(:,i) = Aloc \ nantozip(y(:,i)) ;
  end ;
end ;

