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

% create f(u_t| )
fu.class = 'enorm<ldmat>';
fu.rv    = u;
fu.mu    = 0;
fu.R     = 1e-1;

% create DS
DS.class = 'MpdfDS';
DS.mpdf.class  = 'mprod';
DS.mpdf.mpdfs  = {fy, epdf2mpdf(fu)};
DS.init_rv = RVtimes([y,y,y], [-1,-2,-3]);
DS.init_values = [0.1, 0.2, 0.3];

experiment.ndat=100;

%M=estimator(DS,{},experiment);
M=simulator(DS,experiment);

%%% store results
Data=[M.y'; M.u'];
drv = RVjoin([y,u]);
true_theta=[fy.A fy.const]; 
true_R=[fy.R]; 
save mpdfds_results Data drv y u true_theta true_R