% ch8_1.m
% Learning a non-msv ar(1) solution
% Note fairly robust convergence to initial perturbations but slow convergence.
% Figure 8.1 in book
clear all
randn('seed',13);
alpha=2;
beta0=1.5;
beta1=-1.5;
tmin=2;
tmax=1001;
sig=0.02;
are=-alpha/beta1;
bre=(1-beta0)/beta1;
phi=zeros(2,tmax);
y=zeros(tmax,1);
phi0=[are,bre]'+[-0.01,0.005]';
ymean=alpha/(1-beta0-beta1);
ylag=ymean-.01;
z=[1,ylag]';
% We set init r equal to lim E M(z)
r=eye(2);
r(1,2)=ymean;
r(2,1)=ymean;
r(2,2)=ymean^2+sig^2/(1-bre^2);
y(1)=ylag;
phi(:,tmin-1)=phi0;
for t=tmin:tmax
z=[1,ylag]';
y(t)=alpha+beta0*phi(1,t-1)+beta1*phi(1,t-1)*(1+phi(2,t-1)) ...
	+(beta0*phi(2,t-1)+beta1*phi(2,t-1)^2)*ylag+sig*randn(1);
rnew=r+t^(-1)*(z*z'-r);
phi(:,t)=phi(:,t-1)+t^(-1)*inv(rnew)*z*(y(t)-phi(:,t-1)'*z);
ylag=y(t);
r=rnew;
end
a=phi(1,:)';
b=phi(2,:)';
a(tmin:tmax);
b(tmin:tmax);
y(tmin:tmax);
figure(1)
subplot(3,1,1)
plot(a(tmin:tmax)) 
title('Figure 8.1')
ylabel('a(t)')
subplot(3,1,2)
plot(b(tmin:tmax)) 
ylabel('b(t)')
subplot(3,1,3)
plot(y(tmin:tmax)) 
xlabel('time')
ylabel('y(t)')