%%% Demo of the trendacd 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=3; % the number of the sinusoid semi-periods forming the monotonic trend
DNmin=20; %the minimum number of values in a sinusoid semi-period
          % DNmin>10 and P*DNmin<N/2
fi=0.3; % the serial correlation parameter of the AR(1) noise (0<=fi<1)
r=1; % the ratio between the amplitudes of the trend variations and 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 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;
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; % the coefficient c_p in Eq. (2.2) is constant
        % such that the artificial trend should be monotonic
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;
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);

trend=trendacd(x);
figure;plot(x,'.g'); hold; plot(f,'b'); plot(trend,'r')
xlabel('n'); ylabel('x_n')
legend('x_n','f_n','f_n^{est}')