function [a_st,b_st]=idarx11(a,b,lambda_e,lambda_u,N,Nr)
% function [a_st,b_st]=idarx11(a,b,lambda_e,lambda_u,N,Nr)
% Proprietà asintotiche del metodo dei minimi quadrati
% Si generano Nr realizzazioni di N dati del medesimo modello ARX(1,1):
% y(t) = a*y(t-1) + b*u(t-1) + e(t), e(·)~WN(0,lambda_e^2), u(·)~WN(0,lambda_u^2)
% a_st = stima di a; b_st = stima di b
% Per ogni realizzazione si opera la stima a minimi quadrati dei parametri.
% Si confronta la distribuzione campionaria delle stime con quella prevista
% dalla teoria (distribuzione gaussiana)
%
% 1) Provare ad aumentare N (la gaussiana si dovrebbe restringere)
% 2) Modificare il file in modo che la distribuzione del disturbo non sia gaussiana


u = randn(N,1)*lambda_u;
TH = zeros(Nr,2);
for k=1:Nr
   e = randn(N,1)*lambda_e;
   y = zeros(N,1);
   y(1) = e(1);
   for t=2:N
      y(t) = a*y(t-1) + b*u(t-1) + e(t);
   end;
   TH2 = arx([y u],[1 1 1]);
   TH(k,1) = -TH2(3,1);
   TH(k,2) = TH2(3,2);
end;
PHI = [y(2:N) u(2:N)];
Mean_TH = mean(TH);
a_st = Mean_TH(1);
b_st = Mean_TH(2);
Var_TH = (lambda_e^2)*inv(PHI'*PHI);

for npar=1:2
   x=[-3*sqrt(Var_TH(npar,npar)):.001:3*sqrt(Var_TH(npar,npar))];
   figure(npar)
   plot(x+Mean_TH(npar),(1/sqrt(Var_TH(npar,npar)*2*pi))*exp(-(x.^2)/(2*Var_TH(npar,npar))),'r')
   hold on
   if npar==1
      title(['Stima di a su ' int2str(Nr) ' realizzazioni di ' num2str(N) ' dati']);
   else
      title(['Stima di b su ' int2str(Nr) ' realizzazioni di ' num2str(N) ' dati']);
   end;
   [n,h]=hist(TH(:,npar),50);
   sup = (h(2)-h(1))*Nr;
   plot(h,n/sup)
   hold off
end;
