%
% load time series data
%
load evaporator_data;
t = evap(:,1);
u = evap(:,2);
y = evap(:,3);
clear evap;
%
% construct regression matrices
%
nn = 250;
tt = t(7:nn);
yy = y(7:nn);
%xx = [ y(6:nn-1) u(1:nn-6) ];
%
xx = [ y(6:nn-1) y(5:nn-2) u(2:nn-5) ];
%xx = [ y(6:nn-1) y(5:nn-2) y(4:nn-3) u(3:nn-4) ];
%
% do and evaluate regression
%
phat = xx \ yy;
yyp = xx * phat;
figure;
plot( tt, yy, '+', tt, yyp, '-' );
legend( 'Actual', 'Predicted' );
xlabel( 'Time' );
ylabel( 'Output' );
%
% residuals plots
%
e = yy - yyp;
figure;
plot( tt, e, '*' );
ax = axis;
line( ax(1:2), [0 0] );
xlabel( 'Time' );
ylabel( 'Residuals' );
figure;
plot( yy, e, '*' );
ax = axis;
line( ax(1:2), [0 0] );
xlabel( 'Output' );
ylabel( 'Residuals' );
%
% explained variation
%
et = yy - mean( yy );
varexpl = 100 * ( et' * et - e' * e ) / ( et' * et )
%
% parameter confidence limits
%
[ n, m ] = size( xx );
s = inv( xx' * xx ) * ( e' * e ) / ( n - m );
l = tdistinv( n-m, 0.05 ) * sqrt( diag( s ) );
[ phat l 100*l./abs(phat) ]
%
% parameter correlation
%
r = eye( m );
for i = 1:m-1,
  for j = i+1:m,
    r(i,j) = s(i,j) / sqrt( s(i,i) * s(j,j) );
  end;
end;
r
%
% prediction precision limits
%
s = xx * inv( xx' * xx ) * xx' * ( e' * e ) / ( n - m );
l = tdistinv( n-m, 0.05 ) * sqrt( diag( s ) );
lmean = mean( l );
[ lmean 100*lmean/mean(abs(yy)) ]
%
