% Listing 21.2 for estimating ARMAX 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 
p_armax = idpoly([1 -0.5],[0 0 0.6 -0.2],[1 -0.3],1,1,'Noisevariance',0.1);

% Create input sequence 
uk = idinput(2555,'prbs',[0 0.2],[-1 1]); 

% Simulate the process 
yk = sim(p_armax,uk,simOptions('Addnoise',true));

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

%% Estimate ARMAX model

% Compute IR for time-delay estimation 
mod_fir = impulseest(Ztrain);
figure; impulseplot(mod_fir,'sd',3);
% Time-delay = 2 samples
% Estimate ARMAX model 
na = 1; nb = 2; nc = 1; nk = 2; 
mod_armax = armax(Ztrain,[na nb nc nk]) 

%% Residual analysis

% Check the residual plot	 
figure ; resid(Ztrain, mod_armax);

% Analyze the model
present(mod_armax) 