clear all
% name random variables
y = RV({'y'},1);
u1 = RV({'u1'},1);
u2 = RV({'u2'},1);

% create f(y_t| y_{t-3}, u_{t-1})
fy.class = 'mlnorm<ldmat>';
fy.rv    = y;
fy.rvc   = RVtimes([y,u1,u2], [-3, 0, 0]);
fy.A     = [0.5, -0.9, 0.9];
fy.const = 0;
fy.R     = 1e-2;

fu.class = 'enorm<ldmat>';
fu.rv    = RVjoin([u1,u2]);
fu.mu    = [0,0];
fu.R     = 0.1*eye(2);

f.class = 'mprod';
f.pdfs  = {fy,fu};

DS.class = 'PdfDS';
DS.pdf = f;

% create ARX estimator
A1.class = 'ARX';
A1.rv = y;
A1.rgr = RVtimes([y,u1],[-3,0]) ; % correct structure is {y,y}
A1.log_level ='logbounds,logevidence';
A1.frg = 0.95;
A1.constant = 0;

A2=A1;
A2.rgr = RVtimes([y,u2],[-3,0]) ; % correct structure is {y,y}

M=estimator(DS,{A1,A2},struct('ndat',100));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot results
ndat = size(M.DS_dt_u1,1);

true_theta1 = fy.A([1,2]);
true_theta2 = fy.A([1,3]);
true_R = fy.R;

figure(1);
subplot(2,1,1);
plot(M.DS_dt_y);
title('Output');

subplot(2,1,2);
plot(M.DS_dt_u1);
hold on
plot(M.DS_dt_u2);
title('Input');


figure(2)
hold off;
subplot(2,2,1);
hold off
plotestimates(true_theta1, ...
    M.Est0_apost_mean_theta, ...
    M.Est0_apost_lbound_theta, ...
    M.Est0_apost_ubound_theta);
set(gca,'YLim',[-1.5,1]);

subplot(2,2,2);
hold off
plotestimates(true_R, ...
    M.Est0_apost_mean_r, ...
    M.Est0_apost_lbound_r, ...
    M.Est0_apost_ubound_r);

title('Variance parameters r')

subplot(2,2,3);
hold off
plotestimates(true_theta2, ...
    M.Est1_apost_mean_theta, ...
    M.Est1_apost_lbound_theta, ...
    M.Est1_apost_ubound_theta);
set(gca,'YLim',[-1.5,1]);

subplot(2,2,4);
hold off
plotestimates(true_R, ...
    M.Est1_apost_mean_r, ...
    M.Est1_apost_lbound_r, ...
    M.Est1_apost_ubound_r);

title('Variance parameters r')