% script to do method # 1
% de-trend by OLS 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
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,2];
Pvec=[0 1];       % orders of seasonal AR
Qvec=[0];       % 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 1 remove poly by OLS 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
descriptor=[descriptor,' M1 '];
POLYORDER=2;
[detr1,polystuff]=polyres(data,POLYORDER); % note de-trended in detr1

figure;
y=polyval(polystuff,dindex)'; % make it a column
plot(dindex',data(dindex),'-',dindex',y,'--r');
title(['OLS polyfit of order ',int2str(POLYORDER),' to ',descriptor]);
xlabel(xvalues);
%%%%%%%%%%%%   plot the de-trended data
figure;
plot(detr1);
title([descriptor,' OLS de-trended by polynomial of order ',int2str(POLYORDER)]);

%%%%%%%%%%%%   compute and plot periodogram
%figure;
%xt=newpgram(x,nfft,np1,np2,halflen,alpha,rejalpha,plsw,logsw,dataname)

newpgram(detr1,nd,10,floor(nd/2),8,PERALPHA,.001,1,1,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;
%%% old pgram(detr1dm,nd,[2:floor(nd/2)],8,PERALPHA,[descriptor,' after permean1']);
newpgram(detr1dm,nd,10,floor(nd/2),8,PERALPHA,.001,1,1,[descriptor,' after permean1']);
%%%%%%%%%%%%   compute and plot periodogram of squares
%figure;

% old pgram(detr1dm.^2,nd,[2:floor(nd/2)],8,PERALPHA,['squares of ',descriptor,' after permean1']);
newpgram(detr1dm.^2,nd,10,floor(nd/2),8,PERALPHA,.001,1,1,['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,estspec,'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;

%%% old pgram(resids.^2,nd,[2:floor(nd/2)],8,PERALPHA,['squares of ',descriptor,' SARMA residuals']);
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=descriptor;
    mkdir(basename);
    for j=1:nfig
        figure(j);
        name=sprintf('M1_%0.2d',j);
        print('-deps',['.\',basename,'\',name]); % put in subdir 
    end
end
%%%%%%%%%%%%%%%%% end for now