#!/usr/bin/octave
% parametry AR(2)
fi = [0.75 -0.5];
%
G = 1./roots([-fi(2) -fi(1) 1]);
% czy to jest proces stacjonarny?
if (
  abs(G)<1 
  && fi(1)+fi(2) < 1 
  && fi(2)+fi(1) < 1
  && -1 < fi(2) 
  && fi(2) < 1
  )
  printf "proces stacjonarny\n"
end
% funkcje autokorelacji
ro(1)=1;
ro(2)=fi(1)/(1-fi(2));
for i =3:19
  ro(i)=fi(1)*ro(i-1)+fi(2)*ro(i-2);
end
% liczba próbek
t = [0:200];
% biały szum
a = randn(1,length(t));
% wartości procesu
z(1:2) = a(1:2);
for i = 3:length(t)
  z (i) = fi(1) * z(i-1) + fi(2)*z(i-2)+a(i);
end
sigma = sqrt(1/(1-dot(ro(2:3),fi(1:2))));
f0 = abs(atan2(imag(G(1)),real(G(1)))/(2*pi));
T  = 1/f0;
d = abs(G(1));
F = atan2( (1+d^2)*tan(2*pi*f0) , 1-d^2);
for i =1:length(ro)
  ro2(i) =              \
  sign(fi(1))           \
  * (d^(i-1) / sin (F)) \
  * (sin(2*pi*f0*(i-1) + F)) ;
end
if(ro-ro2>3*eps) exit end
% czas relaksacji
tau = (-log(d))^(-1)
f=0:0.01:0.5;
pzz=2./(1+fi(1)^2+fi(2)^2-2*fi(1)*(1-fi(2))* cos(2*pi*f)-2*fi(2)*cos(4*pi*f));
% szereg czasowy
subplot(3,1,1)
plot(t,z)
% autokorelacje
subplot(3,1,2)
bar(0:length(ro)-1,ro)
axis([-0.5 length(ro)], "xaxis")
% gęstość widmowa
subplot(3,1,3)
plot(f,pzz)
% pause
print -dsvg "ar.svg"