% ch2.m
% Muth model
% Figure 2.1
randn('seed',9);
alpha=-0.5;
mu=5;
gamma=1;
tmax=200;
sig=0.5;
phi=zeros(2,tmax);
p=zeros(tmax,1);
phinaught=[1,2]';
w=randn(1);
z=[1,w]';
% Note here limit M(z)=eye(2)
p(1)=(mu+alpha*phinaught(1))+(gamma+alpha*phinaught(2))*w+sig*randn(1);
r=eye(2);
phi(:,1)=phinaught+inv(r)*z*(p(1)-phinaught'*z);
for t=2:tmax
w=randn(1);
z=[1,w]';
p(t)=(mu+alpha*phi(1,t-1))+(gamma+alpha*phi(2,t-1))*w+sig*randn(1);
rnew=r+t^(-1)*(z*z'-r);
phi(:,t)=phi(:,t-1)+t^(-1)*inv(rnew)*z*(p(t)-phi(:,t-1)'*z);
r=rnew;
end
a=phi(1,:)';
b=phi(2,:)';
figure(1)
subplot(3,1,1)
plot(a) 
title('Figure 2.1')
ylabel('a(t)')
subplot(3,1,2)
plot(b) 
ylabel('b(t)')
subplot(3,1,3)
plot(p) 
xlabel('time')
ylabel('p(t)')