echo on
clc
%	In this example we compare several different methods of 
%	identification.

%	We start by forming some simulated data.  Here are the
%	numerator and denominator coefficients of the transfer
%	function model of a system:

B = [0 1 0.5];
A = [1 -1.5 0.7];

%	These can be put into the special "theta" structure used by
%	the toolbox:

th0 = poly2th(A,B,[1 -1 0.2]);

%	The third argument, [1 -1 0.2], gives a characterization
%	of the disturbances that act on the system.

pause	% Press any key to continue.
clc

%	Now we generate an input signal U, a disturbance signal E,
%	and simulate the response of the model to these inputs:

u = sign(randn(350,1));
e = randn(350,1);
y = idsim([u e],th0);
z = [y u]; % A matrix with the output and inputs as column vectors.

%	Plot first 100 values of the input U and the output Y.

pause, idplot(z,1:100), pause  % Press any key for plot.
clc
%	Now that we have corrupted, simulated data, we can estimate models
%	and make comparisons.  Let's start with spectral analysis.  We
%	invoke SPA which returns a spectral analysis estimate of
%	the transfer function GS and the noise spectrum NSS:

[GS,NSS] = spa(z);

pause, bodeplot(GS), pause   % Press any key for plot. 
clc

%	Now use the Least Squares method to find an ARX-model with
%	two poles, one zero, and a single delay on the input:

a2 = arx(z,[2 2 1]);

pause, present(a2), pause   % Press any key to present results.
clc

%	Now use the instrumental variable method (IV4) to find a model
%	with two poles, one zero, and a single delay on the input:

i2 = iv4(z,[2 2 1]);


pause, present(i2), pause     % Press any key to present results.
clc
%	Calculate transfer functions for the two models obtained by
%	ARX and IV4.  Show Bode plots comparing the Spectral Analysis
%	estimate with the ARX and IV4 estimates:

[Ga2,NSa2] = th2ff(a2);		% From ARX
Gi2 = th2ff(i2);			% From IV4

pause, bodeplot([GS Ga2 Gi2]), pause  % Press any key for plots.

clc
%	Calculate and plot residuals for the model obtained by IV4:

e = resid(z,i2); pause
plot(e), title('Residuals, IV4 method'), pause
clc

%	Now compute second order ARMAX model:

am2 = armax(z,[2 2 2 1]);

pause	% Press any key to continue.
clc
%	Now compute second order BOX-JENKINS model:

bj2 = bj(z,[2 2 2 2 1]);

%	Compute transfer functions and noise spectra for these models:

[Gam2,NSam2] = th2ff(am2);
[Gbj2,NSbj2] = th2ff(bj2);

%	Compare these transfer functions with those obtained by the
%	IV4 method and the spectral analysis method:

pause	% Press any key for plots.

bodeplot([Gam2 Gbj2 Gi2 GS]), pause
bodeplot([NSam2 NSbj2 NSS]), pause, clc

%	Plot residuals for the ARMAX and BOX-JENKINS models

e1 = resid(z,am2); pause
e2 = resid(z,bj2); pause
plot([e1 e2]), title('ARMAX and BOX-JENKINS residuals'), pause
clc
%	Simulate the ARMAX and BOX-JENKINS models with the real input

yham2 = idsim(u,am2);
yhbj2 = idsim(u,bj2);

pause	% Press any key for plot.

plot([y yham2 yhbj2]), ...
 title('Real output and simulated model outputs '), pause
clc
%	Now we'll compute the poles and zeros of the ARMAX and
%	BOX-JENKINS models

zpam2 = th2zp(am2);
zpbj2 = th2zp(bj2);

	% Press any key for plot.
pause	% Press keys to advance poles and zeros on plot, from the ARMAX
	% model to the Box-Jenkins model.

zpplot([zpam2 zpbj2]), pause, clc

% Clearly, AM2 and BJ2 both give good residuals and simulated outputs.
% They are also in good mutual agreement. Let's check which has the
% smaller AIC:  

[am2(2,1) bj2(2,1)]
pause	% Press any key to continue.
clc
%	Since we generated the data,  we enjoy the luxury of comparing
%	the model to the true system:

[G0,NS0] = th2ff(th0);
zp0 = th2zp(th0);

pause	% Press any key for plots of comparisons.

bodeplot([G0 Gam2]), pause
bodeplot([NS0 NSam2]), pause
zpplot([zp0 zpam2]), pause
echo off
set(gcf,'NextPlot','replace');

