下面是我写的一个小程序,希望有用。
function yf=fun(x,lamda)
%Hodrick-Prescott filter
%argin:x: the data need to be filtered, time dimension in the row
% lamda: define the weight for the trend. default value 1600
%argout: yf: the HP filtered data
%
%by Huabin Wu. Dec 10, 2007
[t,nvar]=size(x);
if t<nvar
x=x';
end
[t,nvar]=size(x);
if nargin<2
lamda=1600;
elseif nargin>2
error('Too many arguments!');
end
if t>4
d=[lamda -4*lamda 6*lamda+1];
d=ones(t,1)*d;
m=diag(d(:,3))+diag(d(1:t-1,2),1)+diag(d(1:t-1,2),-1)...
+diag(d(1:t-2,1),2)+diag(d(1:t-2,1),-2);
m(1,1)=1+lamda; m(1,2)=-2*lamda;
m(2,1)=-2*lamda; m(2,2)=5*lamda+1;
m(t-1,t-1)=5*lamda+1; m(t-1,t)=-2*lamda;
m(t,t-1)=-2*lamda; m(t,t)=1+lamda;
yf=inv(m)*x;
yf=x-yf;
else
error('The data is insufficient to be filtered!')
end
[此贴子已经被作者于2008-12-3 21:49:00编辑过]