function arma21pred(a1,a2,c,lambda,N)
% function arma21pred(a1,a2,c,lambda,N)
%
% Predizione a un passo e a 2 passi di un processo ARMA(2,1)
% y(t) = a1 y(t-1) + a2 y(t-2) + w(t) + c w(t-1) , w(·)~WN(0,lambda^2)
% Confronto con il predittore banale: y^(t) = a1 y(t-1) 
% Funzioni di autocorrelazione degli errori di predizione

clc
disp('Caratteristiche processo')
if a1>0, sign_a1='-'; else sign_a1='+';end;
if a2>0, sign_a2='-'; else sign_a2='+';end;
if c<0, sign_c='-'; else sign_c='+';end;
disp(['A(z) = 1 ',sign_a1,' ',num2str(abs(a1)),'z^(-1) ',sign_a2,' ',num2str(abs(a2)),'z^(-2)'])
disp(['C(z) = 1 ',sign_c,' ',num2str(abs(c)),'z^(-1)'])
disp(' ')

if a2~=0
   poli = roots([1 -a1 -a2]);
else
   if a1~=0
      poli = roots([1 -a1]);
   else
      poli = [];
   end
end;

zero = -c;

if ~isempty(poli)
   if any(poli==zero)
      disp('A(z) e C(z) non sono coprimi: il processo corrisponde ad un AR(1)')
      if length(poli)==2
         if poli(1)==zero
            a1 = poli(2);
         else
            a1 = poli(1);
         end;
      else % 1 polo
         a1 = 0;
      end;
      a2 = 0;
      c = 0;
      if a1>0, sign_a1='-'; else sign_a1='+';end;
      disp(['A(z) = 1 ',sign_a1,' ',num2str(abs(a1)),'z^(-1)'])
      disp('C(z) = 1 ')
      poli = a1;
      disp(' ')
   end;
end;

if ~isempty(poli) & any(abs(poli)>=1)
   disp('Processo non stazionario: (A(z) instabile)')
   return
else
   disp('Processo stazionario: (A(z) stabile)')
end;
disp(' ')

if abs(c)<1
   disp('Rappresentazione canonica')
   disp(' ')
else
   disp('Rappresentazione non canonica')
   disp(' ')
   if (abs(c)==1)
      disp('Non si può applicare la predizione alla Kolmogorov-Wiener:')
      disp('C''è uno zero con modulo 1!')
      return
   else
      disp('Fattorizzazione spettrale:')
      lambda = c*lambda;
      c = 1/c;
      if a1>0, sign_a1='-'; else sign_a1='+';end;
      if a2>0, sign_a2='-'; else sign_a2='+';end;
      if c<0, sign_c='-'; else sign_c='+';end;
      disp(['A(z) = 1 ',sign_a1,' ',num2str(abs(a1)),'z^(-1) ',sign_a2,' ',num2str(abs(a2)),'z^(-2)'])
      disp(['C(z) = 1 ',sign_c,' ',num2str(abs(c)),'z^(-1)'])
      disp(' ')
   end;
end

zero = -c;
if ~isempty(poli)
   if any(poli==zero)
      disp('A(z) e C(z) non sono coprimi')
      if length(poli)==2
         if poli(1)==zero
            a1 = poli(2);
         else
            a1 = poli(1);
         end;
      else % 1 polo
         a1 = 0;
      end;
      a2 = 0;
      c = 0;
      if a1>0, sign_a1='-'; else sign_a1='+';end;
      disp(['A(z) = 1 ',sign_a1,' ',num2str(abs(a1)),'z^(-1)'])
      disp('C(z) = 1 ')
      poli = a1;
      disp(' ')
   end;
end;


eta = randn(N,1)*lambda;
v = zeros(N,1);
v(1) = eta(1);
v(2) = a1*v(1)+eta(2)+c*eta(1);
for t=3:N
   v(t) = a1*v(t-1)+a2*v(t-2)+eta(t)+c*eta(t-1);
end;
disp(['Varianza di v(t)  = ',num2str(cov(v))])

%Calcolo predittore ad 1 passo

v_pred1 = zeros(N,1);
v_pred1(1:2) = v(1:2);
for t=3:N
   v_pred1(t) = -c*v_pred1(t-1)+(c+a1)*v(t-1)+a2*v(t-2);
end;

e_pred1 = v-v_pred1;
disp(['Varianza di e1(t) = ',num2str(cov(e_pred1)),' (predittore ad 1 passo)'])


%Predittore a 2 passi

v_pred2 = zeros(N,1);
v_pred2(1:3) = v(1:3);
for t=4:N
   v_pred2(t) = -c*v_pred2(t-1)+(a2+a1*(c+a1))*v(t-2)+(a2*(c+a1))*v(t-3);
end;

e_pred2 = v-v_pred2;
disp(['Varianza di e2(t) = ',num2str(cov(e_pred2)),' (predittore a 2 passi)'])

%Predittore "stupido"
v_predb = [v(1);v(1:N-1)];
e_predb = v-v_predb;
disp(['Varianza di eb(t) = ',num2str(cov(e_predb)),' (predittore stupido)'])

figure(1)
plot(1:N,v,'r',1:N,v_pred1,'b:',1:N,v_pred2,'g:',1:N,v_predb,'m:')
title('Processo (rosso), predittore ad 1 passo (blu), predittore a 2 passi (verde), predittore stupido (magenta)')

figure(2)
iae_pred1(1)=abs(e_pred1(1));
iae_pred2(1)=abs(e_pred2(1));
iae_predb(1)=abs(e_predb(1));
for i=2:N,
   iae_pred1(i)=iae_pred1(i-1)+abs(e_pred1(i));
   iae_pred2(i)=iae_pred2(i-1)+abs(e_pred2(i));
   iae_predb(i)=iae_predb(i-1)+abs(e_predb(i));
end
%plot(1:N,abs(e_pred1),'b',1:N,abs(e_pred2),'g',1:N,abs(e_predb),'y')
plot(1:N,iae_pred1,'b',1:N,iae_pred2,'g',1:N,iae_predb,'m')
title('Int. del modulo dell''err. di pred.: pred. ad 1 passo (blu), pred. a 2 passi (verde), pred. stupido (magenta)')

figure(3)
[IR1,R1,CL1]=cra([e_pred1,e_pred1]);
[IR2,R2,CL2]=cra([e_pred2,e_pred2]);
[IRb,Rb,CLb]=cra([e_predb,e_predb]);
plot(R1(:,1),R1(:,4),'b')
hold on,plot(R2(:,1),R2(:,4),'g')
hold on,plot(Rb(:,1),Rb(:,4),'m')
hold on,plot([-20 20], [CL1 CL1],'k:')
hold on,plot([-20 20], -[CL1 CL1],'k:')
grid
hold off
title('Autocorr. dell''errore di pred.: pred. ad 1 passo (blu), pred. a 2 passi (verde), pred. stupido (magenta)')

