simple_mis_demo.m

uses partial correlation function and filter to implement an AR model

% simple_mis_demo 
%
% most simple demo of a missing data problem with only AR candidate models
%
% Author Piet M.T. Broersen, January 2008
%
% The three statements :
%           *********************************** 
%           *  asel=ARMAsel_mis(ti,xi,ARmax); *
%           *  psdsel=arma2psd(asel,1);       *
%           *  corsel=arma2cor(asel,1,20);    *
%           *********************************** 
%  compute the selected estimated spectrum and autocorrelation function.  
%

clear all, close all, clc, echo off

disp('The true AR process parameters a are given by:')
a = [1 1.2 .8]
disp('The maximum candidate AR order is:')
ARmax= 5   % highest candidate AR order

N=100;
disp('100 equidistant observations are generated')
x=simuarma(a,1,N);               % generates N equidistant observations
disp('Afterward, 50 observations are discarded at arbitrary places')
disp('The remaining observations are a missing data problem, with 50 % missing')
[ti xi]=keepdata(x,.5);          % keeps fraction 0.5 of the data x
mui=mean(xi);
xi=xi-mui;                       % mean is subtracted
plot(1:100,x,'.-',ti,xi,'r+')
title('50 remaining observations, 50 % missing')
xlabel('\rightarrow time axis')
legend('all 100','remaining 50')

disp(' ')
disp('*****')
disp('Some warning messages from the OPTIM toolbox cannot be suppressed easily')
disp('It is not necessary to provide gradient information')
disp('Warnings generated by the MATLAB OPTIM routine fminunc can be ignored')
disp('*****')
disp(' ')

asel=ARMAsel_mis(ti,xi,ARmax);

% the most simple input of ARMAsel_mis is 
%       a row vector ti with the remaining observation times
%       a row vector xi with the remainig observations
%       the highest AR order to be computed
% the most simple output is an estimated AR parameter vector

psdsel=arma2psd(asel,1);
corsel=arma2cor(asel,1,20);

disp(' ')
disp(['The program selected from AR candidate orders 0, 1, ..,',int2str(ARmax)])
disp('In this example, the true process parameters are known') 
disp('and the true PSD and autocorrelation are plotted as a reference')

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

figure
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','estimated ARMAsel-mis');
xlabel('\fontsize{12}\rightarrow  time lag')
axis tight

figure
fig=semilogy(fas,psdtrue,fas,psdsel);
title('\fontsize{12}True and estimated spectrum')
set(fig,'Linewidth',2)
legend('true','estimated ARMAsel-mis',2);
xlabel('\fontsize{12}\rightarrow  frequency')
axis tight