%==============================================================
% ANALYZING FILE airpass.dat:
% International airline passenger monthly totals (in thousands),
% Jan.49-Dec.60.
% x=monthly totals
% t=1,2,...,144: time axis: January 1949 - December 1960
%
% Author: V. Vesely, The University of Malta (C) 2002
%         email: vvese01@um.edu.mt
%==============================================================
echo on
format compact
  air_pass=[crpath,'/brockwel.dat/airpass.dat'];
  [x,xinf]=getdata(air_pass);   % reading and visualizing data
  n=length(x)
  t=(1:n).';
% Graph inspection:
% 1) A clear presence of trend and seasonal periodic component
% 2) Amplitudes of oscillations growing with the trend
%    indicate heteroscedastic data with high probability.
% 3) No apparent outliers visible
% 4) No preprocessing needed.
pause;
%---------------------------------------------------------------
% Inspecting autocorrelations:
  acf(x,figure);
  xlabel(xinf(2,:));
% Slow decay => presence of trend
% Oscillations with period 12 => seasonal one-year periodicity
% of 12 months is to be expected.
  S=12
pause;
%---------------------------------------------------------------
% 1. Checking stability of variance
%---------------------------------------------------------------
  [lambdaE,y] = powtr(x,4,figure);
  xlabel(xinf(2,:));
  title('Parameter lambda estimate for the power transformation');
% Both figure and message confirm linear growth of std.
% Confidence interval for lambda:
  [lambdaE(1)-lambdaE(2),lambdaE(1)+lambdaE(2)]
% Comprising 0 but not 1 => logarithmic transformation and
% multiplicative model (MM) for the original data x is appropriate.
  pause;
%---------------------------------------------------------------------------
% Aditive Model (AM) may be used for transformed data xlog=y only:
  xlog = y;
  figure; plot(t,xlog);
  xlabel(xinf(2,:));
  title('Logarithmically transformed data');
pause;
%---------------------------------------------------------------------------
% 2. Identification of significant harmonic periodic components (periodogram)
%---------------------------------------------------------------------------
% 2.1. Preliminary detrending for the periodogram:
% Trend in xlog is clearly at least linear - we start with degree d>=1
  d=1
  [yE,trE,pE,r,s,bE] = polfs(xlog,d,[],[],[],[],figure);
% Clear residual trend present in the residual error component.
% Confidence interval for the leading coefficient of the polynomial:
  [bE(1,1)-bE(1,2),bE(1,1)+ bE(1,2)]
% Zero not included => parameter is significant
  tr = trE(:,1); % Saving estimated trend
pause;
%---------------------------------------------------------------------------
  d=d+1
  [yE,trE,pE,r,s,bE] = polfs(xlog,d,[],[],[],[],figure);
% Confidence interval for the leading coefficient of the polynomial:
  [bE(1,1)-bE(1,2),bE(1,1)+ bE(1,2)]
% Zero not included => parameter is significant
  tr = trE(:,1); % Saving estimated trend
pause;
%---------------------------------------------------------------------------
  d=d+1
  [yE,trE,pE,r,s,bE] = polfs(xlog,d,[],[],[],[],figure);
% Confidence interval for the leading coefficient of the polynomial:
  [bE(1,1)-bE(1,2),bE(1,1)+ bE(1,2)]
% Zero included => parameter not significant
% We can use the trend estimated in the previous step (degree d=2):
  xlog_tr = xlog-tr;  % Preliminary detrended data
  figure; plot(xlog_tr);
  title('Preliminary detrended data');
pause;
%---------------------------------------------------------------------------
% 2.2. Periodogram, Fisher's and Siegel's tests for periodicity:
% Check for length n being multiple of S:
  rn = rem(n,S)  % remainder after division of n by S
  N = n-rn       % length of the trailing part for periodogram
  [pdg,kF,kS] = tperiod(xlog_tr(rn+1:rn+N),figure); % Periodogram pdg
  title('PERIODOGRAM');
  xlabel(xinf(2,:));
  kF        % subscripts of significant harmonics using Fisher's test
  kS        % subscripts of significant harmonics using Siegel's test
% having periods:
  if length(kS)>length(kF)
     T=N./kS
  else
     T=N./kF
  end
% The first (most significant) harmonic component has seasonal period
% T(1)=12 as expected. Most of the remaining ones are harmonic to it
% (one half, one third, etc.)
% This confirms the dominant periodicity to be definitely
  S = T(1)
% For prediction purposes it is good to cancel harmonics with long
% periods which might distort the predictions
% (sometimes some experimenting with this may help).
% T(3)=[] % might be used in this case, but we shall not do that.
pause;
%---------------------------------------------------------------------------
% 2. Decomposition using all available methods
%---------------------------------------------------------------------------
% NOTE: In MATLAB with simultaneous plots the automatic color order is:
% yellow-magenta-cyan-red-green-blue-white-black
  close all;
% a) Small trend method:
%    -------------------
  MM=1; [sz1,tr1] = szsmtr(x,S,MM);          % MM with original data
  AM=0; [szlog1,trlog1] = szsmtr(xlog,S,AM); % AM with transformed data
% Simultaneous plot of both estimates:
% - periodic component:
  figure; plot(t,[sz1,exp(szlog1)]); % szlog1 of AM is to be transformed back
  title('Seasonal component using the small trend method in MM and AM');
  xlabel(xinf(2,:));
pause;
% - trend
  figure; plot(t,[tr1,exp(trlog1)]);  % trlog1 of AM is to be transformed back
  title('Trend component using the small trend method in MM and AM');
  xlabel(xinf(2,:));
pause;
%---------------------------------------------------------------------------
% b) Moving average method:
%    ----------------------
  MM=1; [sz2,tr2] = szma(x,S,MM);          % MM with original data
  AM=0; [szlog2,trlog2] = szma(xlog,S,AM); % AM with transformed data
% Simultaneous plot of both estimates:
% - periodic component:
  figure; plot(t,[sz2,exp(szlog2)]); % szlog2 of AM is to be transformed back
  title('Seasonal component using the MA method in MM and AM');
  xlabel(xinf(2,:));
pause;
% - trend
  figure; plot(t,[tr2,exp(trlog2)]);  % trlog2 of AM is to be transformed back
  title('Trend component using the MA method in MM and AM');
  xlabel(xinf(2,:));
pause;
%---------------------------------------------------------------------------
% c) Regression model I. ('polfs')
%    (hidden periods by truncated Fourier series expansion)
%    ------------------------------------------------------
% First we definitely confirm polynomial degree 2
% found with the preliminary detrending:
  close all;
  d=2;
  [yE3,trE3,pE3,r3,s,bE] = polfs(xlog,d,T,t,[],[],figure);
  Dlog3  = yE3;
  szlog3 = pE3(:,1);
  trlog3 = trE3(:,1);
  [bE(1,1)-bE(1,2),bE(1,1)+ bE(1,2)]
% Zero not included => leading parameter is significant
% in spite of being very small.
% To be sure, we can check degree 3:
pause;
%---------------------------------------------------------------------------
  [yE3,trE3,pE3,r,s,bE] = polfs(xlog,d+1,T,t,[],[],figure);
  [bE(1,1)-bE(1,2),bE(1,1)+ bE(1,2)]
% Zero included => parameter is not significant - definitely d=2
pause;
%---------------------------------------------------------------------------
% d) Regression model II. ('polper')
%    -------------------------------
  [yE4,trE4,pE4,r4] = polper(xlog,d,S,t,[],[],figure);
  szlog4 = pE4(:,1);
  trlog4 = trE4(:,1);
%---------------------------------------------------------------------------
pause;
% Simultaneous plot of estimates c) and d) after back-transformation
% - periodic components:
  figure; plot(t,exp([szlog3,szlog4]));
  title('Seasonal components using regression models I. and II.');
  xlabel(xinf(2,:));
% matching quite well.
pause;
% - trend
  figure; plot(t,exp([trlog3,trlog4]));
  title('Trend components using regression models I. and II.');
  xlabel(xinf(2,:));
% matching nearly completely.
pause;
%---------------------------------------------------------------------------
% Now simultaneous plot of all AM estimates a) to d)
% after back-transformation
% - periodic components:
  figure; plot(t,exp([szlog1,szlog2,szlog3,szlog4]));
  title('Seasonal components by all methods in AM');
  xlabel(xinf(2,:));
% cyan = 3-rd model 'polfs' is partially outlying
pause;
% - trend
  figure; plot(t,exp([trlog1,trlog2,trlog3,trlog4]));
  title('Trend components by all methods in AM');
  xlabel(xinf(2,:));
% yellow (1) and magenta (2) are not so precise in trend.
% That is why we shall exclude models 1 and 2 from further processing.
% Decision is to be done regarding 'polfs' and 'polper'.
% Let us check the residual errors for randomness:
pause;
%---------------------------------------------------------------------------
  P3=d+1+2*length(T)  % Number of parameters in 'polfs'
  [xro,xroE,mi,miE,xacf,stat] = acf(r3,figure,[],[],P3);
pause;
% Once more: increasing K to P+2 to have positive number of degrees
% of freedom = 2:
  [xro,xroE,mi,miE,xacf,stat3] = acf(r3,figure,[],P3+2,P3);
  stat3
% stat3(1)>stat3(2) => WN rejected
% Except one correlation at lag 1, which contradicts white noise,
% autocorrelations seem to come from stationary process.
% We try 'trand' as a secondary criterion only
% ('trand' is not designed for residual data)
  rnd3=trand(r3)
% None of the 5 tests does reject randomness.
pause;
% The same analysis for 'polper':
  P4=d+S      % Number of parameters in 'polper'
  [xro,xroE,mi,miE,xacf,stat4] = acf(r4,figure,[],P4+2,P4);
  stat4
% Clearly not stationary.
  rnd4=trand(r4)
% 2 test reject randomness
% FINAL CONCLUSION: Only model 3 ('polfs') accepted.
% Although we could finish the analysis now, we shall try to apply
% filtration techniques for the trend.
pause;
%---------------------------------------------------------------------------
% 3. Trend estimation by moving average, alias linear filtration
%---------------------------------------------------------------------------
  close all
  P=P3;
  xlog_sz = xlog - szlog3; % Removing seasonal component
% We shall start with Spencer's weights:
  spencer=[74,67,46,21,3,-5,-6,-3];
  spencer=[fliplr(spencer(2:8)),spencer]/320
  m=length(spencer)
%
% Matrix of weights for stepwise polynomial regression of degree 2
% and length m=15:
  FIR = polyfilt(d,15);
  FIR(8,:)=spencer;        % central weights replaced by 'spencer'
% Computing moving average trend estimate using central weights by Spencer
% and boundary weights from stepwise polynomial regression:
  trsp=convfilt(FIR,xlog_sz);
  rsp =xlog_sz-trsp;          % residual noise
  [xro,xroE,mi,miE,xacf,statsp] = acf(rsp,figure,[],P+2,P);
  statsp,stat3 % statsp(1)<stat3(1)  % true => closer to not rejecting WN
  rnd3,rndsp=trand(rsp,0)
% The result is good but not better than one that from the regression:
% neither statsp(1) nor 'trand' test outcome is better.
pause;
%---------------------------------------------------------------------------
% As another attempt, let us try stepwise polynomial regression
% of degree 2 and lengths 20>=m>=5:
  d=2;
  for m=20:-1:5
      close all
      [trlog,trlogE,r]=polyfilt(d,m,xlog_sz,t);
      plot(t,[trlog3,trlog])
      [xro,xroE,mi,miE,xacf,stat] = acf(r,figure,[],P+2,P);
      m
      stat(1),stat3(1)  % stat(1)<stat(3) => closer to not rejecting WN
      rnd3,rnd=trand(r,0)
      pause;
  end
% No iteration yields better result compared with the previous ones.
pause;
%---------------------------------------------------------------------------
% We shall still keep the original result from model 3 ('polfs').
% Let us inspect residual component for the presence of cyclic
% periodicities. Period is not known - let us experiment with
% truncating:
  close all
  for k=1:n/4
      Nk=n-(k-1);
      [pdgc,kFc,kSc] = tperiod(r3(k:n)); % Periodogram pdgc
      kFc,kSc
      if length(kSc)>length(kFc)
         Tc=Nk./kSc
      else
         Tc=Nk./kFc
      end
      k
      if ~isempty(Tc)
         tperiod(r3(k:n),figure(k)); % Plot periodogram giving some Tc
      end
%      pause;
  end
pause;
%---------------------------------------------------------------------------
% Some periodicities were detected at attempt
  close all
  k=1
  [pdgc,kFc,kSc] = tperiod(r3(k:n),figure); % Periodogram pdgc
  kFc,kSc
  N1=n-(k-1);
  if length(kSc)>length(kFc)
     Tc=N1./kSc
  else
     Tc=N1./kFc
  end
% we shall therefore update the original model:
  [yE3c,trE3c,pE3c,rc,sc,bEc] = polfs(xlog,d,[T,Tc],t,[],[],figure);
  [pdgc,kFc,kSc] = tperiod(rc(rn+1:rn+N),figure); % Periodogram pdg
% No residual periodicities are present any more - all usable information
% seems to be included in the model now.
  P=d+1+2*(length(T)+ length(Tc))  % Number of parameters in 'polfs'
  [xro,xroE,mi,miE,xacf,stat] = acf(rc,figure,[],P+2,P);
  stat,stat3
  rnd=trand(rc)
% Residual noise still looks like stationary close to WN
% (outlying correlation lower a little bit compared with r3)
pause;
%---------------------------------------------------------------------------
% 4. CONCLUSION - FINAL PRESENTATION OF RESULTS
%---------------------------------------------------------------------------
% We definitely accept the model 3 ('polfs') updated by cyclic component
% Let us compute its final evaluation with prediction one period ahead.
% Use zooming from the menu of figure windows to inspect details.
  close all
  tt = 1:(n+S); % prolonged mesh for prediction
  [yE3c,trE3c,pE3c,rc,sc,bEc] = polfs(xlog,d,[T,Tc],t,tt,[],figure);
%
  Dlog   = yE3c(:,1);   D   = exp(Dlog);
  DlogE  = yE3c(:,2:3); DE  = exp(DlogE);
  szlog  = pE3c(:,1);   sz  = exp(szlog);
  szlogE = pE3c(:,2);   szE = exp(szlogE);
  trlog  = trE3c(:,1);  tr  = exp(trlog);
  trlogE = trE3c(:,2);  trE = exp(trlogE);
  rlog   = rc;          r   = exp(rlog);
%
pause;
%---------------------------------------------------------------------------
% Now we plot (possibly back-transformed) results along with
% point confidence regions and original data.
% -------------------------------------------
% Plotting deterministic component D_t (mean estimate) with two
% confidence regions: yellow for D_t itself,
% red for estimated values of x, i.e. of D_t+<estimated noise>
  figure;
% Case AM:
%  plot(t,x,'co',tt,D,'w-',...
%       tt,[D-DE(:,1),D+DE(:,1)],'r--',...
%       tt,[D-DE(:,2),D+DE(:,2)],'g:');
% Case MM:
  plot(t,x,'co',tt,D,'w-',...
       tt,[D./DE(:,1),D.*DE(:,1)],'r--',...
       tt,[D./DE(:,2),D.*DE(:,2)],'g:');
%
  title('Data x (o), D_t (-)');
  xlabel(xinf(2,:));
pause;
%---------------------------------------------------------------------------
% Plotting periodic component P_t with confidence region:
  figure;
% Case AM:
%  plot(tt,sz,'w-',tt,[sz-szE,sz+szE],'r--');
% Case MM:
  plot(tt,sz,'w-',tt,[sz./szE,sz.*szE],'r--');
%
  title('Data x (o), Sz_t (-)');
  xlabel(xinf(2,:));
pause;
%---------------------------------------------------------------------------
% Plotting trend component Tr_t with confidence region:
  figure
% Case AM:
%  plot(t,x,'co',tt,tr,'w-',...
%       tt,[tr-trE,tr+trE],'r--');
% Case MM:
  plot(t,x,'co',tt,tr,'w-',...
       tt,[tr.*trE,tr./trE],'r--');
%
  title('Data x (o), Tr_t (-)');
  xlabel(xinf(2,:));
pause;
%---------------------------------------------------------------------------
% Plotting residual noise component E_t:
  figure
  plot(t,r,'y-');
  title('Residual component E_t in MM');
  xlabel(xinf(2,:));
pause;
%---------------------------------------------------------------------------
% Analysis of residual noise:
  [xro,xroE,mi,miE,xacf,stat] = acf(r,figure,[],P+2,P);
% Looks like stationary noise, which is indicated by:
% 1. Fast decay of autocorrelations to zero - the trend component
%    seems to be completely described by the model.
%    Moreover: the leading correlations are close to the zero 95%
%    confidence region which signals the residual noise to be close to
%    a white noise.
% 2. Confidence region is not too wide, making the previous conclusion
%    quite reliable.
% 3. No apparent periodic behaviour of correlations - the periodic
%    component seems to be completely described by the model.
  stat
% White noise rejected by Portmanteau statistics. Unfortunately,
% because of the correction of K to P+2 (see the trailing 'acf' parameters)
% the obtained result is not much reliable.
  trand(r)
% random behaviour (nearly white noise) confirmed by all tests
% after the back transformation.
pause;
%---------------------------------------------------------------------------
% OVERALL EVALUATION: The quality of the model is very high es evidenced by:
% 1) narrow confidence bounds not opening much within the predicted
%    region => this confirms very low uncertainty of the estimates and
%    good model stability for prediction.
% 2) Residual noise stationary and close to white noise confirming that all
%    significant information is carried by the model.

%================================== END ==================================
echo off
