% December 13,2002
%
% This MatLab code computes the learning dynamics for 
% "Are stationary hyperinflation paths learnable?"
% by K.Adam, G.Evans, and S. Honkapohja
%
% The notation used in this file correspond to that used in the Mathematica notebook (ARMA-SSE-calc3.nb)


randn('state',sum(100*clock)); % This initializes the random number generator

% The model parameters:

lambda=0;   		% the share of agents with the lagged information, change this to affect stability
a=2;					% the coefficient in front of E_t[x(t+1)]
b=-0.5;	   		% the coeffocoent in front of E_(t-1)[x(t)]
c=0;					% the constant


% Defining the variables:

sim_length=1000;				% length of the simulation

phi=randn(6,sim_length);   % phi=(alpha,beta,xi,gamma,theta,delta) is the coefficient 
									% estimates of agents PLM

phi(1,1)=-c/a+randn/4;		% initial values are the AR(1)-REE values	of phi(:,1)plus some noise 
phi(2,1)=((1-b)/a)+randn/4;
phi(3,1)=1/b+randn/4;
phi(4,1)=randn/4;
phi(5,1)=+randn/4;
phi(6,1)=+randn/4; 					

phi(:,2)=phi(:,1);			% the simulation starts with period three (due to necessary 
   								% lags) so initial parameter values are constant in the first two periods
									
z=randn(6,sim_length);		% z=(1,x(t),epsilon(t),epsilon(t-1),eta(t),eta(t-1)) is the 
                           % vector of regression variables
z(1,:)=1;						% initial values are normal random variables
                           % sunspots eta and structural shocks eta are random normal       
for i=2:sim_length
   z(4,i)=z(3,i-1);			% this is necessary for consistency of the eps(t) and eps(t-1)    
   z(6,i)=z(5,i-1);			% and eta(t) and eta(t-1)
end;


R=zeros(6,6);					% initializes the current estimate of the variance covariance 
for i=1:6 						% matrix to an identity matrix
   R(i,i)=1;
end;



% Start the simulation

r=b/a;							% Just a definition

for t=3:sim_length
   
   
   gain=1/(100+t);
   
   %Step1: calculate x(t)
   
   alpha=phi(1,t-1);			% These  are just definitions to make the code easier to read
   beta=phi(2,t-1);
   xi=phi(3,t-1);
   gamma=phi(4,t-1);
   theta=phi(5,t-1);
   delta=phi(6,t-1);
   
   A11=(c/a+(1+r)*(lambda*(1+beta)*alpha+(1-lambda)*alpha))/(1/a-(1-lambda)*beta);
   
   B_factor=1/(1/a-(1-lambda)*beta);
   
   B11=lambda*(beta^2)+r*(1-lambda)*beta;
   B12=r*lambda*(beta^2);
   B13=lambda*beta*gamma+r*(lambda*(beta*xi+gamma)+(1-lambda)*gamma);
   B14=r*lambda*beta*gamma;
   B15=lambda*beta*delta+r*(lambda*(beta*theta+delta)+(1-lambda)*delta);
   B16=r*lambda*beta*delta;
   
   C_factor=B_factor;
   
   C11=lambda*(beta*xi+gamma)+(1-lambda)*gamma+1/a;
   C12=lambda*(beta*theta+delta)+(1-lambda)*delta;
   
   % This is the implied ALM: 
   
   z(2,t)=A11+B_factor*(B11*z(2,t-1)+B12*z(2,t-2)+B13*z(3,t-1)+B14*z(3,t-2)+B15*z(5,t-1)+B16*z(5,t-2))+C_factor*(C11*z(3,t)+C12*z(5,t));
   
   %Step2: calculate new R and then the new phi(t); 
   
   new_R=R+gain*(z(:,t-1)*z(:,t-1)'-R);
   R=new_R;
   
   phi(:,t)=phi(:,t-1)+gain*inv(R)*(z(2,t-1)-phi(:,t-1)'*[1;z(2,t-2);z(3:6,t-1)])*[1;z(2,t-2);z(3:6,t-1)];
   
   
end;


%Step3: Plot the results

figure(1)
plot(phi(3:6,:)');
title('PLM-estimates of \xi, \gamma, \theta and \delta ')
xlabel('time')

figure(2)
plot(phi(1:2,:)');
title('PLM-estimates of \alpha and \beta ')
xlabel('time')

figure(3)
plot(z(2,:)');
title('Path of x(t)')
xlabel('time')

sprintf('The PLM estimates at the end of the simulation:')
sprintf('x(t)=alpha+beta*x(t-1)+xi*eps(t)+gamma*eps(t-1)+theta*eta(t)+delta*eta(t-1)')
sprintf('alpha  : %0.10g     REE-value: %0.10g',phi(1,sim_length),-c/a)
sprintf('beta   : %0.10g     REE-value: %0.10g',phi(2,sim_length),(1-b)/a)
sprintf('xi     : %0.10g     REE-value: undetermined \n',phi(3,sim_length))
sprintf('gamma  : %0.10g     REE-value: %0.10g (as implied by the xi-PLM value) \n',phi(4,sim_length),(b*phi(3,sim_length)-1)/a)
sprintf('theta  : %0.10g     REE-value: undetermined \',phi(5,sim_length))
sprintf('delta  : %0.10g     REE-value: %0.10g (as implied by the theta-PLM value) \n',phi(6,sim_length),(b*phi(5,sim_length))/a)