% Listing 15.2 of the book for Example 15.1
% "Principles of System Identification: Theory and Practice"

% Arun K. Tangirala
% September 04, 2015
% Freely distributed under no warranty

%% Generate data

% Process
proc_mod = idpoly([1 -0.7],[0 2],1,1,1,'Noisevariance',1);
uk = idinput(2046,'prbs',[0 0.4],[-1 1]);
% Measurement with SNR approx 10
ykstar = sim(proc_mod,uk);
proc_mod.Noisevariance = (1 - 0.49)*var(ykstar)/10;
yk = sim(proc_mod,uk,'Noise');
Z = [yk uk];

%%  ML estimation

% Assume known delay and order
th_mle = fminsearch(@(pars) neglogl_arx(Z,pars),[-0.4 1 0.4]');

%% LS estimates
Phi = [-yk(1:end-1) uk(1:end-1)]; yvec = yk(2:end);
th_ls = pinv(Phi)*yvec;