%%%
%%% Exercise 2 for Lecture 3.
%%%

clear all, close all

N = 5000; % Sample number
Ts = 1; % Sampling time
T = [1:N]';

randn('state',1);
e = randn(N,1); % Equation error
sigma_e = 0.1;

e = sigma_e*(e - mean(e))/std(e);
u = idinput(N,'prbs',[0 1],[0 1]); % Identification deterministic input
tc = round(N/2);

u1 = u(1:tc,1);
u2 = u(tc+1:N,1);
e1 = e(1:tc,1);
e2 = e(tc+1:N,1);

a1 = 0.8;
b1 = 0.5;
a2 = 0.2;
b2 = 0.1;

th1_star = [a1 b1];
th2_star = [a2 b2];
d = length(th1_star);

G1z = tf([0 b1],[1 a1],Ts); % Deterministic ARX transfer function
H1z = tf([0 1],[1 a1],Ts); % Stochastic ARX transfer function
G2z = tf([0 b2],[1 a2],Ts); % Deterministic ARX transfer function
H2z = tf([0 1],[1 a2],Ts); % Stochastic ARX transfer function

y1 = lsim([G1z H1z],[u1 e1]); % ARX 1st model simulation
y2 = lsim([G2z H2z],[u2 e2]); % ARX 2nd model simulation
y = [y1;y2];

P0 = 100*eye(d);
th0 = zeros(1,d);
lam = 1.5; % 0.98

thm = rarx([y u],[1 1 1],'ff',lam,th0,P0);

m_arx = arx([y u],[1 1 1]); [A,B,dA,dB] = arxdata(m_arx);
am = [a1*ones(size(u1)); a2*ones(size(u2))];
bm = [b1*ones(size(u1)); b2*ones(size(u2))];

figure, subplot(211), plot(T,thm(:,1),'-',T,am,'--')
xlabel('Samples'), ylabel('a'),
title('True (--) and tracked RLS (-) a parameter')
subplot(212), plot(T,thm(:,2),'-',T,bm,'--'), xlabel('Samples'), ylabel('b')
title('True (--) and tracked RLS (-) b parameters')

figure, subplot(211), plot(T,thm(:,1),'-',T,am,'--',T,A(2)*ones(size(T)),':')
xlabel('Samples'), ylabel('a')
title('True (--), tracked RLS (-) and LS estimated (:) a parameter')
subplot(212), plot(T,thm(:,2),'-',T,bm,'--',T,B(2)*ones(size(T)),':')
xlabel('Samples'), ylabel('b')
title('True (--), tracked RLS (-) and LS estimated (:) b parameter')

return