% ARdemo.m - Time series prediction demonstration
% 
% call tsgenf.m, autocorr.m
% then perform prediction using
% (a) auto-regressive method
% (b) TDNN method
%
% copyright (C) 2000 by Yu Hen Hu
% created: 4/12/2000
% modified: 10/28/2001 add time seris analysis, 
%                      AR model is updated
% modified: 10/7/2005 add titles to figure 2, no pause

clear all,clf
% generate time series data
nmax=1000; %   length of time series
nr=950;    %   length of this sequence used for training
nstep=1;   %   n step prediction
nlag=3;    %   # of past samples used for prediction 
disp('****************************')
disp(['Generate a time series of ' int2str(nmax) ' points.'])
disp(['Use the first ' int2str(nr) ' points for training, and'])
disp(['last ' int2str(nmax-nr) ' points for testing'])
disp(['Use ' int2str(nstep) '-step prediction based on past ' int2str(nlag) ' samples.'])
disp('****************************')
disp('0: AR(1) time series (default)')
disp('1: logistic time series')
disp('2: rounding time series')
tstype=input('Enter your choice: ');
if isempty(tstype), tstype=0; end

[x0,train,test]=tsgenf(nmax,nr,tstype,nstep,nlag);
% x0 is the time series
% train is nr-nlag-nstep+1 by nlag+1,  test is nmax-nr by nlag + 1
%
% test= [x(nmax-nlag-nstep+1)     x(nmax-nstep) |        x(nmax)]
%       [                     ...               |      x(nmax-1)]
%       [x(nr-nlag-nstep+2)   ... x(nr-nstep+1) |        x(nr+1)]
%
% train=[x(nr-nlag-nstep+1)   ... x(nr-nstep)   |        x(nr)  ]
%       [x(nr-nlag-nstep)     ... x(nr-nstep-1) |       x(nr-1) ]
%       [                     ...               |               ]
%       [x(1)      x(2)       ... x(nlag)       |  x(nlag+nstep)]
%

% **************************
% analyzing the time series
% **************************

xmean=mean(x0);
x=x0-xmean; % make x zero mean
rxx=autocorr(x,50); % compute the first 50 auto-correlation lags
figure(1),clf,subplot(211),plot([1:nmax],x)
switch tstype
case 0, title('AR(1) time series')
case 1, title('Logistic time series')
case 2, title('Rounding time series')
end
subplot(212),stem([0:49],rxx),title('auto-correlation lags')

% disp('Press any key to continue ...');pause

% **************************
% auto-regressive method   *
% **************************
a0=levinson(rxx(1:nlag+1),nlag+1);  % a row vector, a(1)=1
a=fliplr(-a0(2:nlag+1)); % filter coefficients a=[-a(nlag) -a(nlag-1) ... -a(1)]
% the reverse of the order is because the train file is organized 
% that each row is [x(n-nlag) x(n-nlag+1) ... x(n-1)]

dap=train(:,1:nlag)*a';     % AR model training set fitting result
tap=test(:,1:nlag)*a';      % AR model nstep prediction results
lag=test(1,1:nlag);   % last nlag of training data
for i=1:length(tap), % recursive prediction using only last nlag training data!
   tr(i)=lag*a';  % i-th prediction
   lag=[lag(2:nlag) tr(i)];
end
figure(2),subplot(211),plot([nlag+1:nr],train(:,nlag+1),'g-',...
   [nlag+1:nr],dap,'b:')
legend('original','training out'), title('training set input vs output')
figure(2),subplot(212),plot([nr+1:nmax],test(:,nlag+1),'g:d',...
   [nr+1:nmax],tap,'r-.',...
   [nr+1:nmax],tr,'b-.'), 
legend('target output','one-step prediction','recursive prediction')
title('testing set input vs output');