% ch9sunu.m
% Learning an arma(1,1) sunspot solution. 
% Adds y(-2) for strong E-stability check.
% Illustrates instability with y(-2) included
% Use as Figure 9.2 in book
% Uses fixed gain for j=1,...,N. Then decreasing gain
clear all
randn('seed',0);
alpha=2;
beta0=1.5;
beta1=-1.5;
N=20;
tmin=2;
tmax=52;
pmax=51;
sig=0.02;
are=-alpha/beta1;
bre=(1-beta0)/beta1;
phi=zeros(5,tmax);
gam=ones(tmax,1);
gam=(1/N)*gam;
for j=N:tmax
	gam(j)=1/j;
end
y=zeros(tmax,1);
cre=0.25;
dre=0.25;
fre=0;
phi0=[are,bre,cre,dre,fre]'+[0,0.1,0,0,0.05]';
ymean=alpha/(1-beta0-beta1);
y2lag=ymean;
ylag=ymean;
vlag=0;
% We set init r equal to lim E M(z)
r=eye(5);
r(1,2)=ymean;
r(2,1)=ymean;
r(2,2)=ymean^2+(dre^2+(1+cre^2)*sig^2)/(1-bre^2);
r(2,3)=sig^2;
r(3,2)=sig^2;
r(3,3)=sig^2;
r(1,5)=ymean;
r(5,1)=ymean;
r(2,5)=are*ymean+bre*r(2,2)+cre*sig^2;
r(5,2)=are*ymean+bre*r(2,2)+cre*sig^2;
r(5,5)=r(2,2);
y(1)=ylag;
phi(:,tmin-1)=phi0;
for t=tmin:tmax
elag=randn(1);
z=[1,ylag,vlag,elag,y2lag]';
v=sig*randn(1);
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+phi(5,t-1)))*ylag ...
	+(beta0+beta1*phi(2,t-1))*phi(5,t-1)*y2lag ...
	+phi(3,t-1)*(beta0+beta1*phi(2,t-1))*vlag ...
	+phi(4,t-1)*(beta0+beta1*phi(2,t-1))*elag+v;
rnew=r+gam(t)*(z*z'-r);
phi(:,t)=phi(:,t-1)+gam(t)*inv(rnew)*z*(y(t)-phi(:,t-1)'*z);
y2lag=ylag;
ylag=y(t);
vlag=v;
r=rnew;
end
a=phi(1,:)';
b=phi(2,:)';
c=phi(3,:)';
d=phi(4,:)';
f=phi(5,:);
figure(1)
subplot(5,1,1)
plot(a(tmin:pmax)) 
title('Figure 9.2')
ylabel('a(t)')
subplot(5,1,2)
plot(b(tmin:pmax)) 
ylabel('b1(t)')
subplot(5,1,3)
plot(c(tmin:pmax)) 
ylabel('c1(t)')
subplot(5,1,4)
plot(d(tmin:pmax)) 
ylabel('d1(t)')
subplot(5,1,5)
plot(f(tmin:pmax)) 
ylabel('b2(t)')
xlabel('time')