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

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

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

% create ARX estimator
A1.class = 'ARX';
A1.yrv = y;
A1.rgr = RVtimes([y,u],[-3,-1]) ; % correct structure is {y,y}
A1.log_level ='logbounds,logevidence';

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

C2 = C1;
C2.windsurfer = 1;

N = 10;
losses_plain=zeros(1,N);
losses_wind =zeros(1,N);

for i=1:10
    M= controlloop(DS,{C1});
    losses_plain(i) = (M.DS_dt_y - C1.yreq)'*(M.DS_dt_y - C1.yreq);

    M= controlloop(DS,{C2});
    losses_wind(i) = (M.DS_dt_y - C1.yreq)'*(M.DS_dt_y - C1.yreq);
    if losses_wind(i)>100
      %  keyboard;
    end
end
[mean(losses_plain), mean(losses_wind)]
[std(losses_plain), std(losses_wind)]