function [s]=hpfilter2(y,lambda)
% This is an improved Hodrick-Prescott filter algorithm, designed to take
% advantage of sparsity. This allows for larger series to be detrended
% than can be by Ivailo Izvorski's hpfilter.m
%
% Inputs:
% y = series to be detrended
% lambda = smoothing parameter; typically 1600 for quarterly data
%
% Output:
% s = the smoothed series
%
% The original algorithm was designed by Aubhik Khan and Julia Thomas in
% FORTRAN, and is accessible at juliathomas.net
% This file was modified by Lucas M. Engelhardt, and is accessible at
% lucasmengelhardt.com
% Date last modified: 10/1/2009
% To my knowledge, this program works. However, I do not guarantee its
% accuracy. For a sufficiently short series, the results from this
% function should be compared with the results from Izvorski's hpfilter. If they
% agree, then odds are good that the answer is correct, as this algorithm
% is significantly different from that one.
if size(y,1) 2
% This allows us to solve for z's elements recursively.
z = zeros(size(y,1),1);
z(1) = y(1)/L(1,1);
z(2) = (y(2) - L(2,1)*z(1))/L(2,2);
for row = 3:N
z(row) = (y(row)-L(row,row-2)*z(row-2)-L(row,row-1)*z(row-1))/L(row,row);
end
% Using similar logic, we can now solve for the trended series. We know:
% z = U * s
% Since U's row i only has elements in columns i, i + 1, and i + 2, we can
% say:
% z(N) = U(N,N)*s(N)
% z(N-1) = U(N-1,N)*s(N) + U(N-1,N-1)*s(N-1)
% z(i) = U(i, i + 2) * s(i + 2) + U(i, i + 1) * s(i + 1) + U(i, i)*s(i)
% So, we can solve for s recursively, starting from s(N)
clear y
s=zeros(size(z,1),1);
s(N) = z(N)/U(N,N);
s(N-1) = (z(N-1) - U(N-1,N)*s(N))/U(N-1,N-1);
for row = (N-2):-1:1
s(row) = (z(row)-U(row,row+2)*s(row+2)-U(row,row+1)*s(row+1))/U(row,row);
end