% script to do method # 4
% de-trend by seasonal differencing then ARMA or SARMA modelling
%
% you  must put the name of your datafile (for example suppose it is mydata.dat) 
% in the statement dataname='mydata.dat'
% data must be already present (loaded) into matlab
LOADDATA=1;
if LOADDATA

load .\stat185data\airpass.dat -ascii     
 end

ACFALPHA=.01;   % Pr[type 1 error] for threshold line to test acf for N(0,1/sqrt(n)) 
PACFALPHA=.01;  % Pr[type 1 error] for threshold line to test pacf for N(0,1/sqrt(n)) 
PERALPHA=.01;   % Pr[type 1 error] for threshold line to test periodogram for
                % periodic components
pvec=[0,1,2];  % orders of AR
%pvec=[1 2];
qvec=[0,1,2];  % orders of MA
%qvec=[0 1];
Pvec=[0 1];       % orders of seasonal AR
Qvec=[0 1];       % orders of seasonal MA

SAVEFIGS=0;
RESIDPLOTS=1; % > 0 to plot residual acfs

%dataname='data2(:,1)'; % assumed loaded
%descriptor='US Production'; % strings to be used for plotting
%xvalues='quarters';
%dindex=[1:120];	% select which of the input data to use

%dataname='brillx';     %%%%%%%%%%%%% dataname statement %%%%%%%%%%%%%%%%%%%
%descriptor='WOBBLE';
%xvalues='time';
%dindex=[1:200];
%dindex=[1:length(brillx)];

dataname='airpass';             % assumed loaded
descriptor='Air Passengers';    % strings to be used for plotting
xvalues='months';
T=12;                           % period for seasonal part
dindex=[1:144];	                % select which of the input data to use

eval(['data=',dataname,';']);   % puts the input into data
[nr,nc]=size(data);
if nc > nr
    error('data must be a column vector');
end

data=log10(data);               % log10 if desired
data=data(dindex);
nd=length(data);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     method 4 seasonal differencing 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
descriptor=[descriptor,' M4 '];
NUMDIFFS=1; % number of times to differences (like no. of derivatives to take)
TDIFF=12;    % e.g. x(t)-x(t-TDIFF) is the first difference of span TDIFF
detr1=sdiff(data,TDIFF,NUMDIFFS); % use seasonal difference with TDIFF=1
nd=length(detr1);
%%%%%%%%%%%%   plot the de-trended data
figure;
plot(detr1);
title([descriptor,' de-trended by ',int2str(T),'th difference of order ',int2str(NUMDIFFS)]);

%%%%%%%%%%%%%%%%%%%%% remove the mean; 

detr1=detr1-mean(detr1);

%%%%%%%%%%%%   compute and plot periodogram

estspec=newpgram(detr1,nd,10,floor(nd/2),8,PERALPHA,.001,1,1,descriptor);

%%%%%%%%%%%%   compute and plot periodogram of squares

newpgram(detr1.^2,nd,10,floor(nd/2),8,PERALPHA,.001,1,1,['squares of ',descriptor,' after sdiff']);
%%%%%%%%%%%%% now the acf and pacf
figure;
acfpacf(detr1,floor(nd/2),floor(nd/4),1,ACFALPHA,PACFALPHA,[descriptor,' after sdiff ']);

%%%%%%%%%%%%%%%%% now subtract the periodic mean 

% [xn,pmean,pstd]=permean1(x,per,alpha,miss,dmeanflg,normflg,plflg,datastr)
% these parameters return series with periodic mean out but variance not normalized
% NOTE !!! the following calls not needed with sdtrend because periodic parts are pretty well removed
% [detr1dm,detr1pmean,detr1pstd]=permean1(detr1,T,.01,NaN,1,0,1,descriptor);

%%%%%%%%%%%%   plot the de-trended data
% skip because permean1 plots the de-meaned data
%figure;
%plot(detr1dm);
%title([descriptor,' OLS de-trended by polynomial of order ',int2str(POLYORDER),' after permean1']);

%%%%%%%%%%%%   compute and plot periodogram
% figure;
% pgram(detr1dm,nd,[2:floor(nd/2)],8,PERALPHA,[descriptor,' after permean1']);

%%%%%%%%%%%%% now the acf and pacf with periodic mean out
% figure;
% acfpacf(detr1dm,floor(nd/2),floor(nd/4),1,ACFALPHA,PACFALPHA,[descriptor,' after permean1 ']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%  now do SARMA or ARMA %%%%%%%%%%%%%%%%%%%%%%
%
% Here use SARMA because it can be set for regular ARMA by making Pvec=[0] and Qvec=[0]
% for method 1 
d=0;
D=0;

[mbest,minaic,pbest,qbest,Pbest,Qbest]=sarmabat(detr1,pvec,qvec,Pvec,Qvec,T);

disp(['SUMMARY OF SARMA MODELS FOR DATA : ',descriptor]);
disp(sprintf('AIC min at (p,q)(P,Q)=(%d,%d)(%d,%d) aic=%g',pbest,qbest,Pbest,Qbest,minaic));
[A,B,C,D,F]=polydata(mbest);
disp(['\phi coefficients : ',num2str(A)]);
disp(['\theta coefficients : ',num2str(C)]);

% now plot periodogram and estimated spectra from model
%%%%%%%%%%%%%% compute periodogram


estspec=estspec(1:floor(nd/2));


%%%%%%%%%%%%%%  plot |W(\lambda)|^2
NLAM=length(estspec);	% no. pts to plot in  [0,2\pi]
dlam=pi/NLAM;
lams=[0:dlam:pi-dlam];
clear i;                % avoids any mix-up
e=exp(-i*lams);


thetaspec=polyval(C,e); % flipping evidently not needed
phispec=polyval(A,e);

resids=pe(mbest,detr1);
sigr=std(resids);

W=thetaspec./phispec;

%estspec=avspec(ddata',NLAM*2,NLAM,0)/(2*pi);
modspec=((sigr^2)/(2*pi))*(abs(W).*abs(W));

figure;
subplot(211);
semilogy(lams(2:end),estspec(2:end),'r',lams(2:end),modspec(2:end),'b-.');
xlabel(' 0 \leq \lambda < \pi');
v=axis;
text(v(1)+.2*(v(2)-v(1)),v(3)+.8*(v(4)-v(3)),['theta:',...
      sprintf(' %4.2f',C),sprintf('\n phi:'),sprintf(' %4.2f',A)]);
title('spectral density estimates - log scale');
legend('periodogram','estimated model');


%%%%%%%%%%%%%%  plot residuals 

subplot(212);
plot([1:nd],resids,'-');
title(['Residuals from best SARMA(',sprintf('%d,%d)(%d,%d)',pbest,qbest,Pbest,Qbest),' for data: ',descriptor]);
nres=length(resids);
msres=(resids'*resids)/nres;
legend(['ms resids :',num2str(msres)]);

%%%%%%%%%%%%   compute and plot periodogram of squares

newpgram(resids.^2,nd,10,floor(nd/2),8,PERALPHA,.001,1,1,['squares of ',descriptor,' SARMA residuals']);
%%%%%%%%%%%%%   acf and pacf of residuals

figure;

acfpacf(resids,floor(nres/2),floor(nres/4),1,ACFALPHA,PACFALPHA,[descriptor,' SARMA resids ']);


%%%%%%%%%%%%%%%%%  save figures in files
if SAVEFIGS
    nfig=gcf;
    basename=deblank(descriptor);
    mkdir(basename);
    for j=1:nfig
        figure(j);
        name=sprintf('M4_%0.2d',j);
        print('-deps',['.\',basename,'\',name]); % put in subdir 
    end
end
%%%%%%%%%%%%%%%%% end for now