% ch9msv.m
% Learning an msv ar(1) solution
% Case of stable ar(1) solution
% Appears in book as Figure 9.4
clear all
randn('seed',1);
alpha=0;
beta0=-3.53968254;
beta1=6.66666667;
beta2=-3.17460318;
delt=1;
ncon=20;
ntime=100;
tmin=2;
tmax=20000;
pmax=2000;
sig=0.02;
phi=zeros(2,tmax);
y=zeros(tmax,1);
phi2=zeros(2,tmax);
y2=zeros(tmax,1);
gam=ones(tmax,1);
gam=(1/ncon)*gam;
for j=ntime:tmax
	gam(j)=1/(ncon-ntime+j);
end
ainit=0;
binit=0.7;
phi0=[ainit,binit]';
% ymean=0 since alpha=0
ymean=0;
ylag=ymean;
z=[1,ylag]';
% We set init r equal to lim E M(z)
r=eye(2);
r(2,2)=sig^2/(1-binit^2);
y(1)=ylag;
phi(:,tmin-1)=phi0;
for t=tmin:tmax
z=[1,ylag]';
y(t)=alpha+(beta0+beta2)*phi(1,t-1) ...
	+(beta1+beta2*phi(2,t-1))*phi(1,t-1)*(1+phi(2,t-1)) ...
	+(delt+beta0*phi(2,t-1)+beta1*phi(2,t-1)^2+beta2*phi(2,t-1)^3)*ylag ...
	+sig*randn(1);
rnew=r+gam(t)*(z*z'-r);
phi(:,t)=phi(:,t-1)+gam(t)*inv(rnew)*z*(y(t)-phi(:,t-1)'*z);
ylag=y(t);
r=rnew;
end

%simulation 2
randn('seed',5);
ylag=ymean;
z=[1,ylag]';
% We set init r equal to lim E M(z)
r=eye(2);
r(2,2)=sig^2/(1-binit^2);
y2(1)=ylag;
phi2(:,tmin-1)=phi0;
for t=tmin:tmax
z=[1,ylag]';
y2(t)=alpha+(beta0+beta2)*phi2(1,t-1) ...
	+(beta1+beta2*phi2(2,t-1))*phi2(1,t-1)*(1+phi2(2,t-1)) ...
	+(delt+beta0*phi2(2,t-1)+beta1*phi2(2,t-1)^2+beta2*phi2(2,t-1)^3)*ylag ...
	+sig*randn(1);
rnew=r+gam(t)*(z*z'-r);
phi2(:,t)=phi2(:,t-1)+gam(t)*inv(rnew)*z*(y2(t)-phi2(:,t-1)'*z);
ylag=y2(t);
r=rnew;
end

a=phi(1,:)';
b=phi(2,:)';
a2=phi2(1,:)';
b2=phi2(2,:)';


figure(1)
subplot(2,1,1)
ylabel('b(t)')
plot(b(tmin-1:pmax)) 
title('Figure 9.4')
ylabel('b(t)')
subplot(2,1,2)
plot(b2(tmin-1:pmax)) 
ylabel('b(t)')
xlabel('Results from two simulations')
'value of b at 100, 1000, 2000, 10000, 20000'
b(100)
b(1000)
b(2000)
b(10000)
b(20000)

'value of b2 at 100, 1000, 2000, 10000, 20000'
b2(100)
b2(1000)
b2(2000)
b2(10000)
b2(20000)