% ch8sun.m
% Learning an arma(1,1) sunspot solution
% Figure 8.3
clear all
randn('seed',0);
alpha=2;
beta0=1.5;
beta1=-1.5;
N=20;
tmin=2;
tmax=10001;
pmax=601;
sig=0.02;
are=-alpha/beta1;
bre=(1-beta0)/beta1;
phi=zeros(4,tmax);
y=zeros(tmax,1);
cre=0.50;
dre=0.50;
gam=ones(tmax,1);
gam=(1/N)*gam;
for j=N:tmax
	gam(j)=1/j;
end
phi0=[are,bre,cre,dre]'+[0,0.1,0,0]';
ymean=alpha/(1-beta0-beta1);
ylag=ymean;
vlag=0;
% We set init r equal to lim E M(z)
r=eye(4);
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;
y(1)=ylag;
phi(:,tmin-1)=phi0;
for t=tmin:tmax
elag=randn(1);
z=[1,ylag,vlag,elag]';
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)*ylag ...
	+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);
ylag=y(t);
vlag=v;
r=rnew;
end
a=phi(1,:)';
b=phi(2,:)';
c=phi(3,:)';
d=phi(4,:)';
figure(1)
subplot(4,1,1)
plot(a(tmin:pmax)) 
title('Figure 8.3')
ylabel('a(t)')
subplot(4,1,2)
plot(b(tmin:pmax)) 
ylabel('b(t)')
subplot(4,1,3)
plot(c(tmin:pmax)) 
ylabel('c(t)')
subplot(4,1,4)
plot(d(tmin:pmax)) 
ylabel('d(t)')
%figure(2)
%plot(y(tmin:pmax)) 
xlabel('time')
%ylabel('y(t)')
'parameters a,b,c,d at t=100'
phi(:,100)
'parameters a,b,c,d at t=1000'
phi(:,1000)
'parameters a,b,c,d at t=10000'
phi(:,10000)