% script to do method # 1
% de-trend by removing seasonal polynomials and 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];       % 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);
nd=length(data);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     method 3 remove seasonal polys by OLS 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
descriptor=[descriptor,' M3 '];
POLYORDER=1;
detr1=sdtrend(data,T,POLYORDER,1,descriptor); % note de-trended in 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 sdtrend']);


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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%  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

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,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

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('M3_%0.2d',j);
        print('-deps',['.\',basename,'\',name]); % put in subdir 
    end
end
%%%%%%%%%%%%%%%%% end for now