% 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;

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

% create ARX estimator
A1.class = 'ARX';
A1.rv = y;
A1.rgr = RVtimes([y,u1,u2],[-3,0,0]) ; % correct structure is {y,y}
A1.options ='logbounds,logll';

C1.class = 'LQG_ARX';
C1.ARX = A1;
C1.Qu = 1*eye(2);
C1.Qy = eye(1);
C1.yreq = 1;
C1.horizon = 100;

M= controlloop(DS,{C1});


%%%%%%%%%%%%%%%%%%%%%
%  PLOT

figure;
subplot(2,1,1);
plot(M.DS_y);
subplot(2,1,2);
hold off 
plot(M.DS_u1);
hold on
plot(M.DS_u2);