% Listing 21.1 for estimating ARX model
% "Principles of System Identification: Theory and Practice"

% Arun K. Tangirala
% December 06, 2015
% Freely distributed under no warranty

%% Data generation

% Create the plant and noise model objects 
proc_arx = idpoly([1 -0.5],[0 0 0.6 -0.2],1,1,1,'Noisevariance',0.05);
% Create input sequence 
uk = idinput(2555,'prbs',[0 0.2],[-1 1]); 
% Simulate the process 
yk = sim(proc_arx,uk,simOptions('Addnoise',true));

% Build iddata objects and remove means
z = iddata(yk,uk,1); zd = detrend(z,0); 

%% Estimating the ARX model 

% Compute IR for time-delay estimation 
mod_fir = impulseest(zd);
figure; impulseplot(mod_fir,'sd',3)
% Time-delay = 2 samples
% Estimate ARX model (assume known orders ) 
na = 1; nb = 2; nk = 2; 
mod_arx = arx(zd,[na nb nk];

%% Residual analysis

% Check the residual plot	 
figure ; resid(zd, mod_arx);
% Analyze model errors
present(mod_arx) 