%%%
%%% BJ Model Generation
%%%

N = 1000;
var_e = 0.1;
Ts = 1;

A = 1;
B = [0 1 0.5];
C = [1 -1 0.2];
D = [1 1.5 0.7];
F = [1 -1.5 0.7];

m0 = idpoly(A,B,C,D,F,var_e,Ts);

randn('state',0);
e = iddata([],randn(N,1));
u = iddata([],idinput(N)); 

y = sim(m0,[u e]);
z = [y u];
mi = bj(z,[2 2 2 2 1],'MaxIter',0)
m = bj(z,mi)
m.EstimationInfo
m = bj(z,m); % Continue if m.es.WhyStop shows that maxiter has 
             % been reached.
compare(z,m,mi)