%=============================================================================
% AR(2) process: 
% y(t) = a1 y(t-1) + a2 y(t-2) + e(t)
% (with conjugate poles in the transfer function)
%
% Optimal one-step-ahead predictor:
% yhat(t) = a1 y(t-1) + a2 y(t-2)
%=============================================================================

N = 100;

%-------------------------------------
% (z - p)(z - conj(p))
%  = z^2 - (2 Re p) z + |p|^2
% := z^2   - a1 z      - a2
pole_mod = 0.9;
pole_phase = 0.5;
a1 = 2*pole_mod*cos(pole_phase);
a2 = -pole_mod^2;
%-------------------------------------

% process noise
e = randn(N, 1);

% process
y = filter([1 0 0], [1 -a1 -a2], e);

% output of the best predictor
yhat = filter([a1 a2], [1 0 0], y);
% residual variance
variance1 = sum((y(2:N) - yhat(1:N-1)).^2)/(N-1);   % UNBIASED, because we KNOW that the mean is 0
disp(sprintf('Residual variance of optimal predictor: %6.3f', variance1))

% output of a slightly different predictor, for comparison
yhat2 = filter([(a1+0.2) (a2-0.2)], [1 0 0], y);
variance2 = sum((y(2:N) - yhat2(1:N-1)).^2)/(N-1);
disp(sprintf('Residual variance of slightly different predictor: %6.3f', variance2))

% output of a completely wrong predictor, for comparison
yhat3 = filter([-0.8 0.4], [1 0 0], y);
variance3 = sum((y(2:N) - yhat3(1:N-1)).^2)/(N-1);
disp(sprintf('Residual variance of completely wrong predictor: %6.3f', variance3))

%-------------------------------------
AXIS = [0 100 -10 10];
figure(1)
clf;
hold on;
axis(AXIS)
plot(1:N, y, 'k');
title('AR(2) process with resonance')
pause;
plot(2:N, yhat(1:N-1), 'b');
title(sprintf('AR(2) process; blue: OPTIMAL ESTIMATOR; residual variance: %6.3f', variance1))
pause;
plot(2:N, yhat2(1:N-1), 'r');
title(sprintf('AR(2) process; red: SUB-OPTIMAL ESTIMATOR; residual variance: %6.3f', variance2))
pause;

clf;
hold on;
axis(AXIS)
plot(1:N, y, 'k');
plot(2:N, yhat(1:N-1), 'b');
plot(2:N, yhat3(1:N-1), 'r');
title(sprintf('AR(2) process; red: WRONG ESTIMATOR; residual variance: %6.3f', variance3))
pause;
%-------------------------------------


