%
% Input-Output Time Series Regression for Fitting
%  System and Noise Models
%
%       p = fitmodel( y, u, yt, ut, wt )
%
% where p = parameter vector
%       y = output time series
%       u = input time series
%       yt = lagged output terms (see note)
%       ut = lagged input terms (see note)
%       wt = lagged white noise terms (see note)
%
% note: "yt", "ut" and "wt" are row vectors of length
%       equal to the maximum lag of the output, input
%       and white noise terms respectively and which
%       contain ones and zeros indicating which
%       terms are to be fit
%       
% Bob Newell, Feb 1996
%
%-------------------------------------------------
%
function p2 = fitnoise( y, u, yt, ut, wt )
%
%----------------------------------------------------------
disp( 'TWO-STEP REGRESSION WITH COLOURED NOISE' );
% size of time series
n = min( length( y ), length( u ) );
%
% high-order regression ignoring noise
%
disp( 'Step One:' );
% number of terms
m = 20;
% regression matrices
yy = y(m+1:n+1-m);
xx1 = [];
xx2 = [];
for i = 1:m,
  xx1 = [ xx1 y(m+1-i:n+1-m-i) ];
  xx2 = [ xx2 u(m+1-i:n+1-m-i) ];
end;
xx = [ xx1 xx2 ];
% 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 );
nu = length( ut );
nw = length( wt );
nl = max( [ ny nu nw ] );
% regression matrices
yy = y(m+nl+1:n-m+1);
xx1 = [];
xx2 = [];
xx3 = [];
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 <= nu,
    if ut(i) > 0,
      xx2 = [ xx2 u(m+nl+1-i:n-m-i+1) ];
	end;
  end;
  if i <= nw,
    if wt(i) > 0,
	  xx3 = [ xx3 resid1(nl+1-i:n-2*m+1-i) ];
	end;
  end;
end;
xx = [ xx1 xx2 xx3 ];
% 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;
for i = 1:ny,
  if yt(i) > 0,
    disp( [ '   a' num2str(i) ' = ' num2str(p2(j)) ] );
	j = j + 1;
  end;
end;
for i = 1:nu,
  if ut(i) > 0,
    disp( [ '   b' num2str(i) ' = ' num2str(p2(j)) ] );
	j = j + 1;
  end;
end;
for i = 1:nw,
  if wt(i) > 0,
    disp( [ '   c' num2str(i) ' = ' num2str(p2(j)) ] );
	j = j + 1;
  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
