a1 = -1.61;
a2 = 0.71;
a = [1 a1 a2];
lam = 1;
N = 100;
Ts = 1;
order = 2;
ar2 = idpoly(a,[],[],[],[],lam,Ts)


%% Calculate covariance matrix  using getcov in toolbox
% Average a lot of simulations

figure(1)
Nr = round(10000/N);
P = zeros(order,order);
for iter = 1:Nr
  y = sim(ar2, iddata([],zeros(N,0),Ts), simOptions('AddNoise',true));
  ar2est = ar(y,2);
  plot(ar2est.A(2),ar2est.A(3),'kx')
  hold on
  P = P + getcov(ar2est);
end
P = P/Nr
present(ar2est)

%%
% Calculate Ry(tau) by Yule-Walker equations
% Use it to calculate covariance matrix for parameter estimates

Ry = [1 a1  a2; a1 1+a2 0 ; a2 a1 1]\[lam;0;0];
Ry0 = Ry(1);
Ry1 = Ry(2);
Rbar = [Ry0 Ry1; Ry1 Ry0];
P2 = lam/N*inv(Rbar)

a1std = sqrt(P2(1,1))
a2std = sqrt(P2(2,2))
a1a2corr = P2(1,2)/(a1std*a2std)

% You should get P = P2 approximately
% Note that a1 and a2 estimates are strongly (negatively) correlated here