mis_demo.m

uses partial correlation function and filter to implement an AR model

% mis_demo test of AR, MA or ARMA data
%
% Example data will be generated with available programs
% The example process has all poles and zeros on equal radii
% New data can be analyzed by using a row vector with observations times
% and a row vector with obserbvations as input to ARMAsel_mis
%
%  Author P.M.T. Broersen, January 2008
%  References to the literature are found in the file readme.txt
%  and also the download address of the ARMAsel software
%
%  requires ARMASA toolbox (free)
%  requires armasel_rs toolbox (free)
%  requires MATLAB optim toolbox for optimization of the likelihood
%
%  additional program called in mis_demo
%     ARMAsel_mis with four internal programs  
%  and the auxiliary program  
%     keepdata.m
%  to generate data
%  by transforming a contiguous data set into a missing data problem

clear all, close all, clc, echo off

disp('Generate a missing data problem')
disp('*******************************')

%  Variation of number of observations, missing fraction, 
%  orders of the generating process, radius and maximum AR order
%  gives some experience with the missing data algorithm.
%  Mostly, the accuracy of the selected model is hardly reduced
%  if only a small fraction of the data is missing.

disp('The number of remaining observations is :')
nobs=1000  % number of remaining observations
disp('The remaining fraction of the data is :')
gamma=.5   % fraction that remains
disp('The maximum candidate AR order is :')
ARmax= 8   % highest candidate AR order

alfa=-.8;  % radius of all poles and zeros
oAR=3;     % AR order of generating ARMA process
oMA=2;     % MA order of generating ARMA process

rc=[1 alfa alfa^2 alfa^3 alfa^4 alfa^5 alfa^6 alfa^7 alfa^8 alfa^9 ...
        alfa^10 alfa^11 alfa^12 alfa^13 alfa^14 alfa^15 alfa^16 alfa^17];
a=rc2arset(rc(1:oAR+1));
rc=[1 -alfa alfa^2 -alfa^3 alfa^4 -alfa^5 alfa^6 -alfa^7 alfa^8 -alfa^9 ...
        alfa^10 -alfa^11 alfa^12 -alfa^13 alfa^14 -alfa^15 alfa^16 -alfa^17]; 
b=rc2arset(rc(1:oMA+1));
disp('The true process parameters a and b are given by:')
disp(a)
disp(b)
disp(' ')

% True PSD and autocorrelation
[psdtrue fas]=arma2psd(a,b,500);
cortrue=arma2cor(a,b,20);

N=ceil(nobs/gamma);            % number to keep nobs remaining observations
randn('state',sum(100*clock)); % always gives new and independent random input data 
x=simuarma(a,b,N);             % generates N equidistant observations
varx=std(x)^2; 
[ti xi]=keepdata(x,gamma);     % keeps fraction gamma of the data x
mui=mean(xi);
xi=xi-mui;                     % mean is subtracted

%  Plot histogram of gaps between the data 
tdif=diff(ti);
[nt,tc]=hist(tdif,001:12);
figure
bar(tc,nt,.3,'r')
title('\fontsize{12}Histogram of time intervals between remaining observations')
xlabel('\fontsize{12}\rightarrow time [sec]')

disp(' ');
disp('ARMAsel_mis computes AR, MA and ARMA candidates')
disp('    and selects the best of all candidates.')
disp('    The parameters of AR models are used to compute')
disp('    a number of MA and ARMA candidate models.') 
disp(' ')
disp('    ******************************************************************')
disp('Some error messages from the OPTIM toolbox cannot be suppressed easily')
disp('    It is not necessary to provide gradient information')
disp('    Warnings generated by fminunc can be ignored')
disp('    ******************************************************************')
disp(' ')

[asel,bsel,sellog]=ARMAsel_mis(ti,xi,ARmax);

%  asel and bsel can be treated like parameters from ARMAsel
%
%  function [asel,bsel,sellog]=ARMAsel_mis(ti,xi,Larmax)
%
%  ti     : 1*N vector of time instants where an observation has been made
%  xi     : 1*N vector of observations at ti
%  ARmax  : highest AR order for ML estimation
%         : highest MA and ARMA orders are set automatically in ARMAsel_mis  
%
%  ti are integer numbers 1,2, ...

disp(sellog)
disp(' ')

psdar=arma2psd(sellog.ar,1,500);
psdart=arma2psd(rc2arset([1 sellog.arset{ARmax}]),1,500);
psdma=arma2psd(1,sellog.ma,500);
psdarma=arma2psd(sellog.armaar,sellog.armama,500);
psdsel=arma2psd(asel,bsel,500);

figure
fig=semilogy(fas,psdar,fas,psdart,':',fas,psdma,'--',fas,psdarma,'-.');
title('\fontsize{12}Estimated spectra for missing data problem')
set(fig,'Linewidth',2)
legend('AR(select)',['AR(',int2str(ARmax),') (max)'],'MA(select)','ARMA(select)',3);
xlabel('\fontsize{12}\rightarrow  frequency')
axis tight

corar =arma2cor(sellog.ar,1,20);
corart=arma2cor(rc2arset([1 sellog.arset{ARmax}]),1,20);
corma=arma2cor(1,sellog.ma,20);
corarma=arma2cor(sellog.armaar,sellog.armama,20);
corsel=arma2cor(asel,bsel,20);
 
figure
fig=plot(0:20,corar,0:20,corart,':',0:20,corma,'--',0:20,corarma,'-.');
title('\fontsize{12}Estimated autocorrelations for missing data problem')
set(fig,'Linewidth',2)
legend('AR(select)',['AR(',int2str(ARmax),') (max)'],'MA(select)','ARMA(select)');
xlabel('\fontsize{12}\rightarrow  time lag')
axis tight

figure
subplot(221)
fig=semilogy(fas,psdtrue,fas,psdsel,'r');
title('\fontsize{12}True and selected spectra for missing data problem')
set(fig,'Linewidth',2)
legend('true','ARMAsel-mis',3);
xlabel('\fontsize{12}\rightarrow  frequency')
axis tight
subplot(222)
fig=plot(0:20,cortrue,0:20,corsel,'r');
title('\fontsize{12}True and selected autocorrelations for missing data problem')
set(fig,'Linewidth',2)
legend('true','ARMAsel-mis');
xlabel('\fontsize{12}\rightarrow  time lag')
axis tight
subplot(223)
fig=semilogy(fas,psdar,fas,psdart,':',fas,psdma,'--',fas,psdarma,'-.');
title('\fontsize{12}Selected spectra for different model types')
set(fig,'Linewidth',2)
legend('AR(select)',['AR(',int2str(ARmax),') (max)'],'MA(select)','ARMA(select)',3);
xlabel('\fontsize{12}\rightarrow  frequency')
axis tight

subplot(224)
fig=plot(0:20,corar,0:20,corart,':',0:20,corma,'--',0:20,corarma,'-.');
title('\fontsize{12}Selected autocorrelations for different model types')
set(fig,'Linewidth',2)
legend('AR(select)',['AR(',int2str(ARmax),') (max)'],'MA(select)','ARMA(select)');
xlabel('\fontsize{12}\rightarrow  time lag')
axis tight
subplot

figure
semilogy(fas,psdtrue,'r')
hold on
for io=1:ARmax
    arpar=rc2arset([1 sellog.arset{io}]);
    psdt=arma2psd(arpar,1,500);
    semilogy(fas,psdt*2^(-io))
    if io==ARmax
        title('\fontsize{12}Shifted spectra of all estimated AR models')
        xlabel('\fontsize{12}\rightarrow frequency axis')
        ylabel('\fontsize{12}\rightarrow Shifted logarithm of spectra')
        legend('True','AR(1)','AR(k) shifted',3)
    end
    hold on
end
hold off

disp('The model order and type with the smallest GIC value will be selected.')
disp('That is the single ARMAsel_mis model from automatic selection.')
disp(' ')
disp('The GIC values of all estimated AR, MA and ARMA candidate models are:')
disp(sellog.GIC)
disp(' ')
disp('Model accuracies with the model error ME as measure')
disp(' ')
disp('The basic ME of the, nothing explaining, white noise model for the data is:')
ME0=moderr(1,1,a,b,nobs);
disp(ME0)
disp('The ME values of the AR models for the data are:')
for j=1:ARmax
    ME_AR(j)=moderr(rc2arset([1 sellog.arset{j}]),1,a,b,nobs);
end
disp(ME_AR)
disp('The ME value of the selected model from only AR candidates :')
ME_ARsel=moderr(sellog.ar,1,a,b,nobs);
disp(ME_ARsel)
disp('The ME value of the selected model from only MA candidates :')
ME_MAsel=moderr(1,sellog.ma,a,b,nobs);
disp(ME_MAsel)
disp('The ME value of the selected model from only ARMA(r,r-1) candidates :')
ME_ARMAsel=moderr(sellog.armaar,sellog.armama,a,b,nobs);
disp(ME_ARMAsel)
disp('The ME value of the model automatically selected from all previous candidates :')
MEsel=moderr(asel,bsel,a,b,nobs);
disp(MEsel)