% Listing 14.9 of the book for Example 14.7
% "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 0 0.3 1 0.5 0.2],[1 0.4],1,1,'Noisevariance',1);

% Generate input
uk = idinput(510,'prbs',[0 0.4],[-1 1]); N = length(uk);
% Simulate the process
xk = sim(proc_mod,uk);
proc_mod.Noisevariance = var(xk)/10;
yk = sim(proc_mod,uk,'Noise');
% Plot the input-output data
zk = iddata(yk,uk,1); plot(zk)
zkd = detrend(zk,0);

%% Obtain weighted LS estimate

% Estimate 4 coefficients of a FIR model using ARX routine (LS)
mod_fir = arx(zkd,[0 4 2]); 
% Compute errors and ACF
err_mod = pe(zkd,mod_fir);
acf_err = acf(err_mod.y,20);

% Construct output variables and regressors
Y = zkd.y(6:end); Neff = length(Y); ukd = zkd.u;
Phimat = [];
for k1 = 2:5,
    Phimat  = [Phimat ukd(6-k1:end-k1)];
end
% Scale the data using auto-covariances (not just ACFs)
Sigma_v = eye(Neff) + acf_err(2)*diag(ones(Neff-1,1),-1) + acf_err(2)*diag(ones(Neff-1,1),1);
Sigma_v = var(err_mod.y)*Sigma_v;

% Construct W and scale data using Cholesky factor
W = inv(Sigma_v); Cw = chol(W);
Ys = Cw*Y; Phimat_S = Cw*Phimat;
% WLS estimates of IR models
irest_wls = pinv(Phimat_S)*Ys;

%% Variance of estimates
var_irwls = inv(Phimat_S'*Phimat_S);
err_irwls = sqrt(diag(var_irwls));