%ARFIT	Demonstration of modules contained in the ARfit package.
echo on
clc
%  ARfit is a collection of M-files for the modeling of multivariate 
%  time series with autoregressive (AR) models. ARfit contains modules 
%  for simulation, order selection, estimation, diagnostic checking, 
%  and spectral decomposition of autoregressive models.    

%  This demo shows their use on a bivariate AR(2) process v(k,:),
%
%       v(k,:)' = w' + A1*v(k-1,:)' + A2*v(k-2,:)' + eta(k,:)',
%
%  where the row vectors eta(k,:) are independent and identically
%  distributed Gaussian noise vectors with zero mean and covariance 
%  matrix C. The kth row v(k,:) of the 2-column  matrix v corresponds
%  to an observation of the process at instant k. The intercept
%  vector w is included to allow for a nonzero mean of the AR(p) 
%  process.

pause   % Strike any key to continue.
clc
%  Firstly, let us simulate observations from a bivariate AR(2) process, 
%  choosing the parameters 

w = [ 0.25 ;  0.1 ];

% for the intercept vector,

A1 = [ 0.4   1.2  
       0.3   0.7 ];

%  and

A2 = [ 0.35 -0.3
      -0.4  -0.5 ];

%  for the AR coefficient matrices, and 

C = [ 1.00  0.50
      0.50  1.50 ];

%  for the noise covariance matrix. The two 2x2 matrices A1 and A2 are 
%  put together in one 2x4 coefficient matrix:

A = [ A1 A2 ];

%  We use the module ARSIM to simulate 200 observations from this AR 
%  process: 

v = arsim(w,A,C,200);

pause   % Strike any key to continue.
clc
%  Suppose now that we hadn't any information about how the data set v 
%  is generated, but we want to try fitting an AR model to the time series. 
%  In this case, one first has to decide on the order p of the model. 
%  Assuming that the correct order lies somewhere between 1 and 5, we use 
%  AR to determine the optimum model order using a modified Schwarz 
%  criterion (MSC):

pmin = 1;
pmax = 5;
[west, Aest, Cest, SBC, MSC, FPE, th] = ar(v, pmin, pmax);

%  The output arguments of AR --- west, Aest, and Cest --- are the estimates 
%  for the intercept vector, for the coefficient matrix, and for the noise
%  covariance matrix, respectively. (The matrix th is needed later in
%  the computation of confidence intervals.) These parameters are estimated 
%  for a model of order popt that is chosen as the minimizer of MSC. The 
%  selected order popt in our example is:

m    = 2;                 % dimension of the state space
popt = size(Aest, 2) / m

pause   % Strike any key to continue.
clc
%  Besides the parameter estimates for the selected model, AR has
%  also returned approximations to Schwarz' Bayesian Criterion SBC, to
%  the modified Schwarz criterion MSC, and to Akaike's final prediction 
%  error estimate FPE, each for models of order pmin:pmax. In this demo, 
%  the model order was chosen as the minimizer of MSC. Here are the MSCs 
%  for the fitted models of order 1,...,5: 

disp(MSC)
 
%  To see if using the other criteria would have resulted in the same 
%  model order, take a look at the SBC:

disp(SBC)
%  and at Akaike's FPE:

disp(FPE)
%  Employing one of these other criteria, the optimum model order would 
%  have been chosen as the respective minimizer, too.

pause   % Strike any key to continue.
clc
%  Next one has to check if the model is adequate to represent the
%  time series v. A necessary condition for model adequacy is the
%  uncorrelatedness of its residuals. Invoking
 
[siglev,res] = arres(west,Aest,v);

%  we obtain the time series of residuals res as well as the
%  significance level siglev with which a modified Li-McLeod
%  portmanteau test rejects the null hypothesis of the residuals'
%  uncorrelatedness. A model passes this test if, say, siglev > 0.05. 
%  In our example, the significance level of the modified
%  Li-McLeod portmanteau statistic is  

disp(siglev);
%  (If this is persistently smaller than about 0.05, even over several
%  runs of this demo, then there most likely is a problem with the
%  random number generator that ARSIM used in the simulation of v.) 

pause   % Strike any key to continue.
clc
if exist('xcorr')    % If the Signal Processing Toolbox is installed, plot
                     % autocorrelation function of residuals ...

%  Using ACF, one can also plot the autocorrelation function of, say, 
%  the first component of the time series of residuals:

acf(res(:,1));

%  95% of the lag > 0 autocorrelations should lie between the
%  dashdotted confidence bounds. 

%  Since uncorrelatedness of the residuals is only a necessary
%  condition for model adequacy, in practice further diagnostic tests
%  should be performed. However, we shall go on to the estimation of
%  confidence intervals.

pause   % Strike any key to continue.
end

clc
%  Now being reasonably confident that our model indeed provides an
%  appropriate  representation of the data, we invoke ARCONF to
%  compute confidence intervals for the estimated parameters: 

[Aest,Cest,west] = arconf(Aest,Cest,west,th);

%  The real parts of west, Aest and Cest now contain the estimated AR
%  coefficients and covariances, respectively, the imaginary parts the
%  corresponding 68.27% confidence intervals. Here's the estimated
%  intercept vector with its 68.27% confidence intervals:

disp(west)
%  And, for comparison, the `true' vector as used in the simulation:

disp(w)

pause   % Strike any key to display the other parameter estimates.
clc

%  The estimated coefficient matrix: 

disp(Aest)
%  And, for comparison, the `true' AR coefficients:

disp(A)
pause   % Strike any key to display the covariance matrix estimate.
clc
%  The estimated covariance matrix including confidence intervals is:

disp(Cest)
%  and for comparison the `true' covariances as used in the
%  simulation:

disp(C)

pause   % Strike any key to continue.

echo off
%  Compute `true' eigenmodes from model parameters:
%  Eigenvectors and eigenvalues of corresponding AR(1) coefficient
%  matrix: 
[Strue,EvTrue] = eig([A; eye(2) zeros(2)]);  
EvTrue = diag(EvTrue)';   % `true' eigenvalues
Strue = adjph(Strue);     % adjust phase of eigenmodes
clc
echo on

%  Finally, ARMODE gives us the spectral decomposition of the fitted
%  model:

[S,Sconf,per,tau,exctn] = armode(Aest,Cest,th); 

%  The columns of S are the estimated eigenmodes:

disp(S)
%  with the corresponding 68.27% confidence intervals in Sconf:

disp(Sconf)

%  Compare the estimated eigenmodes above with the modes obtained
%  from the `true' model parameters:

disp(Strue(3:4,:))
%  (Note that the estimated modes can be a permutation of the `true' 
%  modes. The sign of the modes is ambiguous, too.)

pause   % Strike any key to continue.

echo off
pertrue = 2*pi./abs(angle(EvTrue)); % `true' periods
clc

echo on
%  With these eigenmodes are associated the periods:

disp(per)
%  where the second row indicates 68.27% confidence intervals. (Note that
%  for a real eigenvalue, the period becomes infinite, resulting in the
%  indefinite expression 'NaN' for the confidence interval.) Compare
%  the above with the `true' periods:

disp(pertrue)
pause   % Strike any key to get the corresponding damping time scales.

echo off
tautrue = -1./log(abs(EvTrue)); % `true' damping time scales
clc

echo on
%  To each mode also belongs a damping time scale. Here are the
%  estimates for the damping time scales:

disp(tau)
%  with confidence intervals again in the second row. 
%  For comparison, the `true' damping time scales:

disp(tautrue)
pause   % Strike any key to get the excitation corresponding to each mode.
echo off
%  Compute `true' excitation of modes from the designed parameters:
p  = 2;              % true model order

invStr = inv(Strue); % inverse of matrix with eigenvectors as columns

% covariance matrix of corresponding decoupled AR(1) system
CovDcpld = invStr*[C zeros(2,(p-1)*2); zeros((p-1)*2, p*2)]*invStr';

% diagonal of that covariance matrix
DgCovDcpld = real(diag(CovDcpld))';

% excitation 
TrueExctn = abs(EvTrue).^(2*(p-1)) .* DgCovDcpld(1:2*p)./(1-abs(EvTrue).^2);
 
% normalize excitation 
TrueExctn = TrueExctn./sum(TrueExctn);            
clc

echo on
%  ARMODE has also returned the excitations, measures of the relative
%  dynamical importance of the modes:

disp(exctn)

%  Compare this with the `true' excitation computed from the designed
%  parameters:

disp(TrueExctn)

echo off
disp('End')

