% TP6ex2 cas stable d'un processus autorégressif
% PROCESSUS AUTOREGRESSIF GAUSSIEN 

clear
clf
N=input('Nombre d''itérations : ');
NR=input('Nombre de réalisations : '); 
t=input('Paramètre du processus autorégréssif : '); % inférieur en norme à 1 (strictement)
v=input ('Variance associée : ');
if v~=abs(v)
   disp('Une variance est toujours positive !!');
   break
end

% Initialisation des variables
Z=[];ZZ=[];

% Debut du programme
for k=1:NR
   X=randn;
   SN=X*X;TN=0;EST=[TN];
   VN=X*X;SIG=[VN];
   for i=1:N
      % Calcul de la sortie à l'instant i
      eps=sqrt(v)*randn;Y=X;
      X=t*Y+eps;
      % Estimation des moindres carrés
      TN=TN+inv(SN)*Y*(X-TN*Y);
      VN=(i/(i+1))*VN+(1/(i+1))*(X-TN*Y)^2;
      % Mise à jour
      SN=SN+X*X;
      % Stockage
      EST=[EST TN];
      SIG=[SIG VN];
   end
   Z=[Z sqrt(N)*(TN-t)/sqrt(1-t^2)];
   ZZ=[ZZ sqrt(N)*(VN-v)/sqrt(3*v^2)];
end

% Graphiques
p=[1:N];
figure(1);
subplot(1,2,1); % 1ère ligne, 2ème colonne, 1er graphique
plot(p,EST(p),p,t*ones(N,1))
title('Estimation du paramètre')
xlabel('Loi des grands nombres')
subplot(1,2,2)
[EZ,CZ]=histo(Z,sqrt(N),0,1);
hold on
plot(CZ,dnorm(CZ),'r-')
title('Estimation du paramètre')
xlabel('Théorème limite centrale')
hold off

figure(2);
subplot(1,2,1); % 1ère ligne, 2ème colonne, 1er graphique
plot(p,SIG(p),p,v*ones(N,1))
title('Estimation de la variance')
xlabel('Loi des grands nombres')
subplot(1,2,2)
[EZZ,CZZ]=histo(ZZ,sqrt(N),0,1);
hold on
plot(CZZ,dnorm(CZZ),'r-')
title('Estimation de la variance')
xlabel('Théorème limite centrale')

hold off

% ex: nb it      1000  1000
%     nb réal      50    50
%     param       0.2   0.2
%     var           1   0.5