## Contents

clear all


## Generate data from a "true" system with parameters w_0

select some values for the parameters w_0

M = 4;
w_0 = [1.5 -0.8 2.3 0.5]';
% generate the input u(t) as an iid Gaussian sequence
Nd = 200;
Nw = 100;
u = randn(Nd,1);
e = randn(Nd,1);
y  = zeros(Nd,1);
% the output of the sistem w_0 is y(t)= u(t)^T*w_0 + e(t)
for i = (M+1):Nd
y(i) = w_0' * u(i:-1:(i-M+1)) + e(i);
end
% remove the warming up period 1:Nw
u(1:Nw) = [];
y(1:Nw) = [];
N  = Nd-Nw;


## Estimate the parameters w based on the data u and y

Select the data matrix A by using the Covariance method

Ac=[];
for j = M:-1:1
Ac = [Ac u(j+(0:(N-M)))];
end
yc = y(M:N);
% Estimate the LS parameters
hat_w = (Ac'*Ac)\(Ac'*yc);
[w_0 hat_w]

% use the compact formulations to get the LS solution with i2 = N
i1 = M; i2 = N;
Phi = zeros(M,M); psi = zeros(M,1);
for i = i1:i2
vect_u = u(i:-1:(i-M+1));
Phi = Phi + vect_u* vect_u';
psi = psi + vect_u*y(i);
end
hat_w2 =Phi\psi;
Pn = inv(Phi);
[w_0 hat_w hat_w2]

% use the compact formulations to get the LS solution with i2 = N-1
i1 = M; i2 = N-1;
Phi = zeros(M,M); psi = zeros(M,1);
for i = i1:i2
vect_u = u(i:-1:(i-M+1));
Phi = Phi + vect_u* vect_u';
psi = psi + vect_u*y(i);
end
hat_w3 =Phi\psi;
Pn1 = inv(Phi);
param_true_estN_estN1 = [w_0 hat_w2 hat_w3]

% check the recursive equation for Pn
i= N;
vect_u = u(i:-1:(i-M+1));
T = Pn1*vect_u;
Pnrec = Pn1- (Pn1*vect_u*vect_u'*Pn1')/(1+vect_u'*Pn1*vect_u);
diffP = Pnrec-Pn
% check the recursive equation for w
alpha_n = y(N)-hat_w3'*vect_u;
w_rec = hat_w3 + Pnrec*vect_u*alpha_n;
diff_wn = w_rec-hat_w2

ans =

1.5000    1.4270
-0.8000   -0.9892
2.3000    2.3528
0.5000    0.5748

ans =

1.5000    1.4270    1.4270
-0.8000   -0.9892   -0.9892
2.3000    2.3528    2.3528
0.5000    0.5748    0.5748

param_true_estN_estN1 =

1.5000    1.4270    1.4342
-0.8000   -0.9892   -0.9931
2.3000    2.3528    2.3600
0.5000    0.5748    0.5679

diffP =

1.0e-17 *

-0.3469    0.0108    0.0651    0.0054
0.0108    0.6939    0.0217    0.0867
0.0651    0.0217    0.3469    0.0054
0.0054    0.0867    0.0054    0.3469

diff_wn =

1.0e-15 *

0.4441
-0.4441
0
0.1110



## Use the RLS algorithm for recursively computing w for i = 1...N

Initialize

w = zeros(M,1);
delta = 100; lambda = 1;
P = delta*eye(M);
W = [];
for i = (M+1):N
vect_u = u(i:-1:(i-M+1));
piv = vect_u'*P;
gamma = lambda + piv*vect_u;
k = piv'/gamma;
alpha = y(i) - vect_u'*w;
w = w + k*alpha;
Pprime = k*piv;
P = (P-Pprime)/lambda;
% store the parameters
W(:,i) = w;
end
param_true_estN_estN1 = [w_0 hat_w2 hat_w3 w]
diffw = hat_w2-w

param_true_estN_estN1 =

1.5000    1.4270    1.4342    1.4137
-0.8000   -0.9892   -0.9931   -0.9856
2.3000    2.3528    2.3600    2.3184
0.5000    0.5748    0.5679    0.5499

diffw =

0.0133
-0.0037
0.0344
0.0249



## Plot the norm of w-w_0

param_err = [];
for i = 1:N
w = W(:,i);
param_err(i) = sqrt( (w-w_0)'*(w-w_0));
end
plot(param_err)
xlabel('Iteration i')
ylabel('|| w-w_0||')