%%% Demo of the trendrcma program for artificial time series generated by
%%% the algorithm presented in Sections 2.2 and 2.3 of the book
%%% [ATE] C. Vamos and M. Craciun, Automatic Trend Estimation, Springer
%%% 2012.

clear all; close all;
N=1000; % the number of the time series values
P=10; % the number of the sinusoid semi-periods forming the monotonic parts of the trend
DNmin=20; %the minimum number of values in a sinusoid semi-period
          % DNmin>10 and P*DNmin<N/2
fi=0.2; % the serial correlation parameter of the AR(1) noise (0<=fi<1)
r=1; % the ratio between the amplitudes of the trend variations an the noise fluctuations

if P*DNmin>N/2
    error('The minimum number of values in a sinusoid semi-period DNmin is too small')
end

% generation of the artificial non-monotonic trend
amp=rand(1,P);
dmin=(P-1)*DNmin/(N-1-DNmin);
Np=zeros(1,P+1); % the boundaries of the sinusoid semi-periods given by Eq. (2.1) in [ATE]
f=zeros(1,N);
flagro=0;  njf=[];
while flagro==0
    d=dmin+(1-dmin)*rand(1,P);
    for p=1:P
        Np(p+1)=sum(d(1:p)); end;    
    Np=Np/Np(P+1);
    Np=1+round(Np.*(N-1));
    nv=Np(2:P+1)-Np(1:P);
    flagro=1;
    if length(find(nv<DNmin))>0 
        flagro=0; end
end
cp=0.5; Jf=1; njf(1)=1;
for p=1:P
    xant=f(Np(p));
    f(Np(p):Np(p+1))=amp(p)*sin(pi*(cp+2*cp/nv(p)*[0:Np(p+1)-Np(p)]));
    f(Np(p):Np(p+1))=f(Np(p):Np(p+1))-f(Np(p))+xant;
    if rand>0.5
        cp=-cp; % the coefficient c_p in Eq. (2.2) randomly takes
                % the values +1 and -1
        Jf=Jf+1;
        njf(Jf)=Np(p+1); end
    
end
if rand>0.5 f=-f; end
ampt=abs(max(f)-min(f));
f=f-mean(f);
%generation of an AR(1) noise with unit standard deviation (Sect. 1.3 in [ATE])
xnoise=zeros(1,N);
sigZ=sqrt((1-fi*fi));
Z=sigZ*randn(1,N);
xnoise(1)=Z(1)/sigZ;
for n=2:N, xnoise(n)=fi*xnoise(n-1)+Z(n); end
xnoise=xnoise-mean(xnoise);
ampnoise=abs(max(xnoise)-min(xnoise));
% the time series with the AR(1) noise superposed over a monotonic trend
x=f+xnoise/(ampnoise*r/ampt);

%load('serie_graf')
[trend,pest,trendacd] =trendrcma(x,1);
figure;plot(x,'.g'); hold; plot(f,'b'); plot(trend,'r')
xlabel('n'); ylabel('x_n')
legend('x_n','f_n','f_n^{est}')
title(['\deltan=',num2str(pest.dn),'  K=',num2str(pest.K),...
    '  i=',num2str(pest.nrep),'  \phi^{est}=',num2str(pest.fi,'%3.2f'),...
    '  r^{est}=',num2str(pest.r,'%3.2f')])

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% The RCMA trend can be improved if the RCMA algorithm is applied with 
% %% the reduction of dn (the minimum distance between two succesive local
% %% extrema of the RCMA trend) and the nonstationarity of the estimated
% %% residuals is reduced.
% 
% warning off all
% rez=x-trend;
% dni=floor(pest.dn/2);
% [trend2,pest2]=trendrcma(x,1,dni); 
% rez2=x-trend2;
% 
% Nbuc=floor(N/10);
% if Nbuc<100 Nbuc=100;end
% nrbuc=floor(N/Nbuc);
% for in=1:nrbuc
%     [h1v(in),pValue1(in),stat1v(in),cValue1v(in)] = kpsstest(rez((in-1)*Nbuc+1:in*Nbuc),'trend',false);
%     [h2v(in),pValue2(in),stat2v(in),cValue2v(in)] = kpsstest(rez2((in-1)*Nbuc+1:in*Nbuc),'trend',false);
% end
% %% The test of Kwiatkowski, Phillips, Schmidt and Shin (KPSS) assesses the
% %% null hypothesis that a univariate time series is trend stationary
% %% against the alternative that it is a nonstationary unit-root process.
% %% If the number of intervals on which the residuals are assessed to be
% %% nonstationary decrease and the estimated correlation parameter fi is
% %% smaller than 0.7, we consider that the RCMA trend with reduced dn is a
% %% better approximation than the first one.
% if sum(h2v)<sum(h1v)-1 && (pest.fi<0.7 || (pest.fi>0.7 && pest.r>1.5))
%     trendfinal=trend2;
%     pestf=pest2;
% % else
%     trendfinal=trend;
%     pestf=pest;
% % end;
% figure;plot(x,'.g'); hold; plot(f,'b','LineWidth',2); plot(trendfinal,'r','LineWidth',2)
% xlabel('n'); ylabel('x_n')
% legend('x_n','f_n','f_n^{est}')
% title(['dn=',num2str(pestf.dn),'  K=',num2str(pestf.K),'
% nrep=',num2str(pestf.nrep),'  fi=',num2str(pestf.fi,'%3.2f'),'  r=',num2str(pestf.r,'%3.2f')])