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

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

%% Data generation 

% DGP
dgp_mod = idpoly(1,[0 0 2],1,1,[1 -1.48 0.54]);
uk = idinput(2046,'prbs',[0 0.4],[-1 1]); % Full-length PRBS
yktrue = sim(dgp_mod,uk); % Noise-free simulation 

% Adjust the noise variance
dgp_mod.Noisevariance = var(yktrue)/10;
% Simulate with noise
yk = sim(dgp_mod,uk,simOptions('Addnoise',true));

% Detrend and prepare the data
dataset = iddata(yk,uk,1); Ztrain = detrend(dataset,0);

%% Estimate the IR coefficients with regularization

% Turn off regressive modeling (modeling of noise)
options = impulseestOptions; options.Advanced.AROrder = 0;
ir_regulmod = impulseest(Ztrain,options);
% Without regularization
ir_lsmod = arx(Ztrain,[0 51 0]);

% Compute IR coefficients
[ir_true,kvec] = impulse(dgp_mod,50);
[ir_regulest,kvec,~,ir_sdreg] = impulse(ir_regulmod,50);
[ir_lsest,kvec,~,ir_sdls] = impulse(ir_lsmod,50);

figure;
stem(kvec,ir_regulest,'r','Markerfacecolor','red'); hold on
stem(kvec,ir_lsest,'bd');

%% Plot significance regions for regularized estimates
ir_3sdreg = 2.58*ir_sdreg;
fillxvec = [kvec(1) kvec(1) kvec(2:end)' kvec(end:-1:2)']';
fillyvec = [-ir_3sdreg(1) ir_3sdreg(1) ir_3sdreg(2:end)' -ir_3sdreg(end:-1:2)']';
fill(fillxvec,fillyvec,[0.9 0.1 0],'FaceAlpha',0.4,'Linestyle','none');

%% Plot significance regions for the LS estimates
ir_3sdls = 2.58*ir_sdls;
fillyvec2 = [-ir_3sdls(1) ir_3sdls(1) ir_3sdls(2:end)' -ir_3sdls(end:-1:2)']';
fill(fillxvec,fillyvec2,[0.5 0.58 1],'FaceAlpha',0.2,'Linestyle','none');

%% Plot the true IR
figure
stem(kvec,ir_regulest,'r','Markerfacecolor','red'); hold on
stem(kvec,ir_true,'ks');
fill(fillxvec,fillyvec,[0.9 0.1 0],'FaceAlpha',0.4,'Linestyle','none');