% This code estimates Fully Modified HP filter based Trend and Cycle for a time
% sereis. 
% Reference: 
% Hanif, M.N., J. Iqbal and M.A. Choudhary (2017),
% "Fully Modified HP Filter," SBP Working paper No. 88. 
% 
load x; % actual data series (after seasonal adjusment, in case of high freqnecy data series)
z=x;
z=log(z);
siz=size(z);
t=siz(1,1);
% following is a vector for different lambda %
maxlam=10000;
lam=(1:1:maxlam)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);
aa=k'*k;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% for optimal lambda%%%%%%%%%%%%%%%%%%%
sss=zeros(maxlam,1);
for m=1:maxlam   
bb=lam(m,1)*aa;
cc=id+bb;
ddhp=cc^-1;
% tt is trended series %
tto=ddhp*z;
ss=zeros(t,1);
for i=1:t;    
ss(i,1)=(1/t)*((z(i,1)-tto(i,1))^2)*(1+(2*t)/lam(m,1));
end
sss(m,1)=sum(ss);
end
%%%%%% for end point solution%%%%%%%
% Minimum value of our statistics for MHP%
for mm=1:maxlam;
    if(sss(mm,1)==min(sss))
      pt=mm;
    end
end
lamo=lam(pt,1);% modified lambda
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%Loss function linear weight%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
table2=(0:10:2000)';
siz=size(table2);
ns=siz(1,1);
st=round(min(22, ((t/2)-5)));
result=zeros(t,ns);
for pt=1:ns
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);
aa=k'*k;
table3=zeros(t,1);
table3(st+1:t-st)=lamo;
for it=1:st
    table3(it,1)=table2(pt,1)*(st-it+1)+lamo;
end
%for it=1:st-2
   % table3(it,1)=table(pt,1)*(st-it-1)+1600;
%end
for itt=1:st
    table3(t-st+itt,1)=table2(pt,1)*(itt)+lamo;
end
bb=zeros(t,t);
for jd=1:t;
       bb(jd,:)=table3(jd,1)*aa(jd,:);
end
cc=id+bb;
dd=cc^-1;
max1=(-1*(t-1):1:(t-1));
max2=(0:1:t-1);
max3=max2';
k2=1;
thita=(0:0.1:3.1)';
function2=zeros(t,32);
for k2=1:32
gainfunc=zeros(t,1);
for i=1:t
    coss=zeros(1,t);
    for j=1:t
    coss(1,j)=dd(i,j)*cos(thita(k2,1)*(-1*max3(i,:)+max2(:,j)));
    end
    sinn=zeros(1,t);
    for j=1:t
    sinn(1,j)=dd(i,j)*sin(thita(k2,1)*(-1*max3(i,:)+max2(:,j)));
    end
    cosfunc=(sum(coss))^2;
    sinfunc=(sum(sinn))^2;
    gainfunc(i,1)=sqrt(cosfunc+sinfunc);
end
function1=zeros(t,1);
for j=1:t
    function1(j,1)=0.1*(gainfunc(j,1)-gainfunc(round(t/2),1))^2;
end
function2(:,k2)=function1;
end
function4=zeros(t,1);
for k4=1:t
    function4(k4,1)=sum(function2(k4,:));
end
result(:,pt)=function4;
end
finalresult1=sum(result);
for ss=1:ns;
    if(finalresult1(1,ss)==min(finalresult1))
      ps=ss;
    end
end
minn1=min(finalresult1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%For fully modified lambda%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);
aa=k'*k;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% for optimal lambda%%%%%%%%%%%%%%%%%%%
sss=zeros(maxlam,1);
for m=1:maxlam   
table3=zeros(t,1);
table3(st+1:t-st)=lam(m,1);
for it=1:st
    table3(it,1)=table2(ps,1)*(st-it+1)+lam(m,1);
end
for itt=1:st
    table3(t-st+itt,1)=table2(ps,1)*(itt)+lam(m,1);
end
bb=zeros(t,t);
for jd=1:t;
       bb(jd,:)=table3(jd,1)*aa(jd,:);
end
cc=id+bb;
dd=cc^-1;
% tt is trended series %
tto2=dd*z;
ss=zeros(t,1);
for i=1:t;    
ss(i,1)=(1/t)*((z(i,1)-tto2(i,1))^2)*(1+(2*t)/lam(m,1));
end
sss(m,1)=sum(ss);
end
%%%%%% for end point solution%%%%%%%
% Minimum value of our statistics for MHP%
for mm=1:maxlam;
    if(sss(mm,1)==min(sss))
      pt3=mm;
    end
end
lamfmhp1=lam(pt3,1);% fully modified lambda 
table3=zeros(t,1);
table3(st+1:t-st)=lamfmhp1;
for it=1:st
    table3(it,1)=table2(ps,1)*(st-it+1)+lamfmhp1;
end
for itt=1:st
    table3(t-st+itt,1)=table2(ps,1)*(itt)+lamfmhp1;
end
bb=zeros(t,t);
for jd=1:t;
       bb(jd,:)=table3(jd,1)*aa(jd,:);
end
cc2=id+bb;
dd=cc2^-1;
tt1=dd*z;% Trend from FMHP
cc1=z-tt1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%non linear part%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%Loss function linear weight%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
table2=(0:5:2000)';
siz=size(table2);
ns=siz(1,1);
st=round(min(22, ((t/2)-5)));
result=zeros(t,ns);
for pt=1:ns
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);
aa=k'*k;
table3=zeros(t,1);
table3(st+1:t-st)=lamo;
for it=1:st
    table3(it,1)=table2(pt,1)*(st-it+1)*(st-it+1)+lamo;
end
%for it=1:st-2
   % table3(it,1)=table(pt,1)*(st-it-1)+1600;
%end
for itt=1:st
    table3(t-st+itt,1)=table2(pt,1)*(itt)*(itt)+lamo;
end
bb=zeros(t,t);
for jd=1:t;
       bb(jd,:)=table3(jd,1)*aa(jd,:);
end
cc=id+bb;
dd=cc^-1;
max1=(-1*(t-1):1:(t-1));
max2=(0:1:t-1);
max3=max2';
k2=1;
thita=(0:0.1:3.1)';
function2=zeros(t,32);
for k2=1:32
gainfunc=zeros(t,1);
for i=1:t
    coss=zeros(1,t);
    for j=1:t
    coss(1,j)=dd(i,j)*cos(thita(k2,1)*(-1*max3(i,:)+max2(:,j)));
    end
    sinn=zeros(1,t);
    for j=1:t
    sinn(1,j)=dd(i,j)*sin(thita(k2,1)*(-1*max3(i,:)+max2(:,j)));
    end
    cosfunc=(sum(coss))^2;
    sinfunc=(sum(sinn))^2;
    gainfunc(i,1)=sqrt(cosfunc+sinfunc);
end
function1=zeros(t,1);
for j=1:t
    function1(j,1)=0.1*(gainfunc(j,1)-gainfunc(round(t/2),1))^2;
end
function2(:,k2)=function1;
end
function4=zeros(t,1);
for k4=1:t
    function4(k4,1)=sum(function2(k4,:));
end
result(:,pt)=function4;
end
finalresult=sum(result);
for ss=1:ns;
    if(finalresult(1,ss)==min(finalresult))
      ps2=ss;
    end
end
minn2=min(finalresult);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%For fully modified lambda%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);
aa=k'*k;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% for optimal lambda%%%%%%%%%%%%%%%%%%%
sss=zeros(maxlam,1);
for m=1:maxlam   
table3=zeros(t,1);
table3(st+1:t-st)=lam(m,1);
for it=1:st
    table3(it,1)=table2(ps2,1)*(st-it+1)*(st-it+1)+lam(m,1);
end
for itt=1:st
    table3(t-st+itt,1)=table2(ps2,1)*(itt)*(itt)+lam(m,1);
end
bb=zeros(t,t);
for jd=1:t;
       bb(jd,:)=table3(jd,1)*aa(jd,:);
end
cc=id+bb;
dd=cc^-1;
% tt is trended series %
tto2=dd*z;
ss=zeros(t,1);
for i=1:t;    
ss(i,1)=(1/t)*((z(i,1)-tto2(i,1))^2)*(1+(2*t)/lam(m,1));
end
sss(m,1)=sum(ss);
end
%%%%%% for end point solution%%%%%%%
% Minimum value of our statistics for MHP%
for mm=1:maxlam;
    if(sss(mm,1)==min(sss))
      pt3=mm;
    end
end
lamfmhp2=lam(pt3,1);% fully modified lambda 
table3=zeros(t,1);
table3(st+1:t-st)=lamfmhp2;
for it=1:st
    table3(it,1)=table2(ps2,1)*(st-it+1)*(st-it+1)+lamfmhp2;
end
for itt=1:st
    table3(t-st+itt,1)=table2(ps2,1)*(itt)*(itt)+lamfmhp2;
end
bb=zeros(t,t);
for jd=1:t;
       bb(jd,:)=table3(jd,1)*aa(jd,:);
end
cc2=id+bb;
dd=cc2^-1;
tt2=dd*z;% Trend from FMHP
cc2=z-tt2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if minn1<minn2
    ttt=tt1;
else
    ttt=tt2;
end
finaltrend=exp(ttt);
finalcycle=x-finaltrend;