%
% time series fit of DO-air data
%
load ex7_8.dat
nd = input( 'no. data points to use ? ' );
t = 60 * ( ex7_8(1:nd,1) - ex7_8(1,1) );
u = ex7_8(1:nd,2);
y = ex7_8(1:nd,3);
clear ex7_8;
%----------------- now general -------------------
%
% calculate zero mean data
%
u = u - mean(u);
y = y - mean(y);
ymean = mean(y)
%
% get model structure
%
n = input( 'no. of AR (y) terms ? ' );
m = input( 'no. of MA (u) terms ? ' );
d = input( 'length of delay ? ' );
%
% construct regression matrices
%
nd = length( y ) - max( [ n d+m+1 ] );
Y = y(1:nd);
X = [];
for i = 1:n,
  X = [ X y(i+1:nd+i) ];
end;
for i = 1:m,
  X = [ X u(i+d+1:nd+d+i) ];
end;
%
% calculate parameters, predictions and residuals
%
p = X \ Y;
Yp = X * p;
e = Y - Yp;
%
% plot them
%
np = input( 'no. of plot points ? ' );
plot( t(1:np), [ Y(1:np) Yp(1:np) e(1:np) u(1:np) ] );
%
set( gca, 'fontsize', 16 );
xlabel( 'Time in minutes' );
ylabel( 'Zero-mean values' );
%
stext( 2.5, 0.75, '{\times \24 \it y_P}' );
stext( 4.2, 0.75, '{\times \24 \it y}' );
stext( 4.2, -0.07, '{\times \24 \it e}' );
stext( 4.2, -0.57, '{\times \24 \it u}' );
%
printsto -f1 -deps2 ex7_8a.eps
%
% parameter confidence limits
%
s = inv( X' * X ) * ( e' * e ) / ( nd - n - m );
stt = tdistinv( nd - n - m, 0.025 );
pcl = stt * sqrt( diag( s ) );
%
% display some results
%
disp( 'parameters and confidence limits' );
disp( [ p pcl ] );
disp( 'sum-of-squares' );
disp( e' * e );
disp( 'percentage explained' );
disp( 100 * ( cov(Y) - cov(e) ) / cov(Y) ); 
%
frspec( u, mean(diff(t)), 1, 1 );
%
set( gca, 'fontsize', 16 );
xlabel( 'Frequency in radians/minute' );
ylabel( 'Relative power' );
%
print -f2 -deps2 ex7_8b.eps
%
