%
% Single Time Series Regression for Fitting
%  Noise Models
%
%       [ a, c ] = fitnoise( y, yt, wt )
%
% where a = output parameter vector
%       c = noise parameter vector
%       y = time series
%       yt = lagged output terms (see note)
%       wt = lagged white noise terms (see note)
%
% note: "yt" and "wt" are row vectors of length
%       equal to the maximum lag of the output
%       and white noise terms respectively and which
%       contain ones and zeros indicating which
%       terms are to be fit
%       
% Bob Newell, Feb 1996
%
%-------------------------------------------------
%
function [ a, c ] = fitnoise( y, yt, wt )
%
%----------------------------------------------------------
disp( 'TWO-STEP REGRESSION WITH COLOURED NOISE' );
% size of time series
n = length( y );
%
% high-order regression ignoring noise
%
disp( 'Step One:' );
% number of terms
m = 20;
% regression matrices
yy = y(m+1:n+1-m);
xx = [];
for i = 1:m,
  xx = [ xx y(m+1-i:n+1-m-i) ];
end;
% regression
var1 = inv( xx' * xx );
sv1 = svd( xx' * xx );
p1 = var1 * xx' * yy;
resid1 = yy - xx * p1;
sos1 = resid1' * resid1;
vare1 = sos1 / ( n - 2*m );
disp( [ 'Sum of squares = ' num2str(sos1) ] );
disp( [ 'Estimated error variance = ' num2str(vare1) ] );
%--------------------------------------------------------
%
% low order regression including residuals
%
disp( 'Step Two:' );
% maximum lag
ny = length( yt );
nw = length( wt );
nl = max( ny, nw );
% regression matrices
yy = y(m+nl+1:n-m+1);
xx1 = [];
xx2 = [];
for i=1:nl,
  if i <= ny,
    if yt(i) > 0,
      xx1 = [ xx1 y(m+nl+1-i:n-m-i+1) ];
	end;
  end;
  if i <= nw,
    if wt(i) > 0,
	  xx2 = [ xx2 resid1(nl+1-i:n-2*m+1-i) ];
	end;
  end;
end;
xx = [ xx1 xx2 ];
% regression
var2 = inv( xx' * xx );
sv2 = svd( xx' * xx );
p2 = var2 * xx' * yy;
resid2 = yy - xx * p2;
sos2 = resid2' * resid2;
vare2 = sos2 / ( length(resid2) - length(p2) );
disp( 'Parameters: ' );
j = 1;
a = [];
for i = 1:ny,
  if yt(i) > 0,
    disp( [ '   a' num2str(i) ' = ' num2str(p2(j)) ] );
	a = [ a p2(j) ];
	j = j + 1;
  else
    a = [ a 0 ];
  end;
end;
c = [];
for i = 1:nw,
  if wt(i) > 0,
    disp( [ '   c' num2str(i) ' = ' num2str(p2(j)) ] );
	c = [ c p2(j) ];
	j = j + 1;
  else
    c = [ c 0 ];
  end;
end;
disp( [ 'Sum of squares = ' num2str(sos2) ] );
disp( [ 'Estimated error variance = ' num2str(vare2) ] );
disp( 'SV test:' );
disp( 100*sv2/sum(sv2) );
disp( 'Covariance Matrix =' );
d = diag( var2 );
disp( var2 ./ sqrt( d * d' ) );
%
return;
%---------------------------------------------------------
% thats all
