%%%
%%% Exercise Lecture 3.
%%%

clear all, close all
N = 1000; % Sample number
Ts = 1; % Sampling time
T = [1:N]';

randn('state',1);
e = randn(N,1); % Equation error
sigma_e = 0.5;
e = sigma_e*(e - mean(e))/std(e);

pulse_signal = 0;
square_wave = 0;

if(pulse_signal),
tp = round(N/3);
u = zeros(N,1); u(tp) = 1;
elseif(square_wave),
period = 100;
u = square(2*pi*T/period);
else,
u = idinput(N,'prbs',[0 1],[0 1]);
end;

a = -0.8;
b = 0.5;
c = 0.4;
th_star = [a b];
d = length(th_star);

armax_model = 1;
arx_model = 0;

if(armax_model),
Gz = tf([0 b],[1 a],Ts); % Deterministic ARMAX transfer function
Hz = tf([1 c],[1 a],Ts); % Stochastic ARMAX transfer function
elseif(arx_model),
Gz = tf([0 b],[1 a],Ts); % Deterministic ARX transfer function
Hz = tf([0 1],[1 a],Ts); % Stochastic ARX transfer function
end;

y = lsim([Gz Hz],[u e]); % ARX or ARMAX model simulation

P0 = 100*eye(d);
th0 = zeros(1,d);
lam = 1; % 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);

figure, subplot(211), plot(T,thm(:,1),'-',T,th_star(1)*ones(size(T)),'--')
xlabel('Samples'), ylabel('a'),
title('True (--) and tracked RLS (-) a parameters')
subplot(212), plot(T,thm(:,2),'-',T,th_star(2)*ones(size(T)),'--')
xlabel('Samples'), ylabel('b')
title('True (--) and tracked RLS (-) b parameters')

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

figure, subplot(211)
plot(T,thm(:,1),'-',T,th_star(1)*ones(size(T)),'--',T,A(2)*ones(size(T)),':',...
T,(A(2)+dA(2))*ones(size(T)),'-.',T,(A(2)-dA(2))*ones(size(T)),'-.')
xlabel('Samples'), ylabel('a')
title('True (--), tracked RLS (-) and estimated (:) a parameters')
subplot(212)
plot(T,thm(:,2),'-',T,th_star(2)*ones(size(T)),'--',T,B(2)*ones(size(T)),':',...
T,(B(2)+dB(2))*ones(size(T)),'-.',T,(B(2)-dB(2))*ones(size(T)),'-.')
xlabel('Samples'), ylabel('b')
title('True (--), tracked RLS (-) and estimated (:) b parameters')

return