% Good Day%
% This MATLAB code is for estiomation of a) Modified HP filter lambda,% 
% b) trend, and c) cyclical component of a given series (y) %
% For details see: M. Ali Choudhary, M. Nadim Hanif and Javed Iqbal (2014)%
% On smoothing macroeconomic time series using the modified HP filter, %
% Applied Economics, 46:19, 2205-2214 %
% This code is insensitive to data frequency, we can apply same code on
% annual, quterly and monthly data%
load x
% where x is original data series (monthly, quaterlly or annual, single series)%
siz=size(x);
t=siz(1,1);
maxlam=10000;
% where maxlam is maximum number of lambda to try for best fit.%
% Choice of maxlam is optional. This 10000 is based upon a working on
% variuos macroeconomic series of 93 countries %
% Following is the process to estimate Modifeid HP filter by using eq.6%
% on the paper% 
lam=(1:1:maxlam)';
% Estimate k for A =K*K' matrix, next line to eq.3 on our paper%
k=zeros(t-2,t);
for i=1:t-2;
   for j=1:t;
    if i==j 
        k(i,j)=1;
    end
    if j==i+2
        k(i,j)=1;
    end
    if j==i+1
        k(i,j)=-2;
    end
    end
end
id=eye(t,t);
a=k'*k;
sss=zeros(maxlam,1);
for m=1:maxlam   
b=lam(m,1)*a;
cc=id+b;
d=cc^-1;
% gt (eq.3, paper) is trended component %
gt=d*x;
ss=zeros(t,1);
for i=1:t; 
    % Following is eq.6 on Choudhary et.al(2014)%
ss(i,1)=((x(i,1)-gt(i,1))^2)*(1+(2*t)/lam(m,1));
end
sss(m,1)=sum(ss);
end
for mm=1:maxlam;
    if(sss(mm,1)==min(sss))
      pt=mm;
    end
end
lambda=lam(pt,1);
% Where Lambda is the smoothing parameter by modified HP filter approach%
% For trend and cyclical component from  modfied HP filter  we use lambda estimated above.%
bmhp=lambda*a;
ccmhp=id+bmhp;
dmhp=ccmhp^-1;
% Now, tcmhp and ccmhp (below) are final trend and cyclical component of %
% the series respectivly by Modified HP filter approach%
tcmhp=dmhp*x;
ccmhp=x-tcmhp;