% ch8msvs.m
% Learning an msv ar(1) solution
% Case of stable ar(1) solution
% Appears in book as Figure 8.5
clear all
randn('seed',13);
alpha=2;
beta0=1.5;
beta1=-1.5;
delt=0.1;
rt=(1+4*delt*beta0*(1-beta0)^(-2))^0.5;
bre1=(beta0-1)*(1-rt)/(2*beta0);
bre2=(beta0-1)*(1+rt)/(2*beta0);
are1=alpha/(1-beta0-beta1*(1+bre1));
are2=alpha/(1-beta0-beta1*(1+bre2));
'RE values for first ar(1) solution'
are1
bre1
'Check stability conditions. Need both values <0'
beta0+beta1-1+beta1*bre1
beta0-1+2*beta1*bre1
'RE values for second ar(1) solution'
are2
bre2
'Check stability conditions. Need both values <0'
beta0+beta1-1+beta1*bre2
beta0-1+2*beta1*bre2
tmin=2;
tmax=600;
are=are2+.005
bre=bre2-.01
sig=0.02;
phi=zeros(2,tmax);
y=zeros(tmax,1);
phi0=[are,bre]';
ymean=are/(1-bre);
ylag=ymean+.001;
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)) ...
	+(delt+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-1:tmax)) 
title('Figure 8.5')
ylabel('a(t)')
subplot(3,1,2)
plot(b(tmin-1:tmax)) 
ylabel('b(t)')
subplot(3,1,3)
plot(y(tmin-1:tmax)) 
xlabel('time')
ylabel('y(t)')