% script to do method # 1
% de-trend by first differencing and then remove periodic terms by permean1,
% 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


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];       % orders of seasonal MA


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);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     method 2 remove trend by first differences 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
descriptor=[descriptor,' M2 '];
NUMDIFFS=1; % number of times to differences (like no. of derivatives to take)
TDIFF=1;    % 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 difference of order ',int2str(NUMDIFFS)]);



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

%%%%%%%%%%%%% now the acf and pacf
figure;
acfpacf(detr1,floor(nd/2),floor(nd/4),1,ACFALPHA,PACFALPHA,[descriptor,' de-trended ']);

%%%%%%%%%%%%%%%%% 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
[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']);

%%%%%%%%%%%%   compute and plot periodogram of squares
figure;

%[xt,xxt]=sqplot(x,np1,np2,fftlen,plsw,dbsw,alpha,halflen,dataname)
pgram(detr1dm.^2,nd,[2:floor(nd/2)],8,PERALPHA,['squares of ',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(detr1dm,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=pgram(detr1dm,nd,[],0,.01,descriptor);
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);

Cp=fliplr(C); % convert to matlab poly form for use in polyval

Ap=fliplr(A);

thetaspec=polyval(Cp,e);
phispec=polyval(Ap,e);
resids=pe(mbest,detr1dm);
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,modspec,'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
figure;

%[xt,xxt]=sqplot(x,np1,np2,fftlen,plsw,dbsw,alpha,halflen,dataname)
pgram(resids.^2,nd,[2:floor(nd/2)],8,PERALPHA,['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

nfig=gcf;
basename=descriptor;
mkdir(basename);
for j=1:nfig
    figure(j);
    name=sprintf('M2_%0.2d',j);
    print('-deps',['.\',basename,'\',name]); % put in subdir 
end

%%%%%%%%%%%%%%%%% end for now