% Marcet's parameterized expectations method as applied to
% ``Introduction of Derivative Securities: A General Equilibrium
% Approach'', Stock Economy: 
%
% Agent i:     max   sum delta^t c[i,t]^(1-gamma(i))/(1-gamma(i))    
%           {c[i,t]}  t 
%
%       s.t. c[i,t] + p[t]a[i,t] = a[i,t-1](p[t] + y[t])
%            y[t] = rho0 + rho1 y[t-1] + rho2 y[t-2] + eps[t]
%            a[i,t] >= ks
%            a[i,0] = 1
%
% Equilibrium conditions:
%            
%            c[1,t] + c[2,t] = 2 y[t]
%            a[1,t] + a[2,t] = 2 
%

% Ellen R. McGrattan, 11-30-90

% Reference: Ketterer and Marcet (1989).

delta=.89;
gam1=.7;
gam2= 3;
rho0=.1;
rho1=1.3;
rho2=-.4;

sige=sqrt(.0016);
ks=.935;

lambda=.5;

maxn=500;
tol=1e-5;
T=5000;

rand('normal');
rand('seed',0);
epsilon=rand(T+2,1)*sige;
y(1,1)=rho0/(1-rho1-rho2);
y(2,1)=rho0/(1-rho1-rho2);
for i=3:T+2; 
  y(i,1)=rho0+rho1*y(i-1)+rho2*y(i-2)+epsilon(i); 
end;
X=[ones(T+1,1),log(y(2:T+2)),log(y(1:T+1)),ones(T+1,1)];
[row,col]=size(X);

c=y;

beta=zeros(3,2);
beta(1,[1,2])=log(1.0/(1-delta))*ones(1,2);
tem= [..
    5.1539    5.1590
    0.0877    0.0857
   -0.0372   -0.0360
   -2.9808   -2.9791
  -45.4674  -45.4958
   -0.0527   -0.1087
    0.0282    0.0472
   47.7865   47.7553];

beta(:,1)= [ ..
    2.1885   
    0.0591  
   -0.0253];
beta(:,2)= [ ..
    2.2628  
    0.0508  
   -0.0210  ];
beta=[beta;0 0];
beta=[tem(1:4,1),tem(5:8,1)];


a(1,1)=1;
p(1,1)=delta*y(1)/(1-delta);
global zxv

for n=1:maxn;
  
  short=0;
  long=0;
  for i=2:T+1;   
    c0=c(i-1);
    diff=1;
    tem=X(i,:)*beta(:,1)-X(i,:)*beta(:,2);
    while diff>tol; 
      c1=c0+(gam2*log(2*y(i)-c0)-gam1*log(c0)-tem)/(gam1/c0+gam2/(2*y(i)-c0));
      diff=abs(c0-c1)/(1+abs(c0));
      c0=c1;
    end;
    c(i,1)=c1;
    p(i,1)=delta*exp(X(i,:)*beta(:,2))*(2*y(i)-c(i,1))^gam2; 
    a(i,1)=(a(i-1,1)*(p(i)+y(i))-c(i))/p(i);
    if (a(i)-a(i-1))<1e-10; a(i)=a(i-1); end;
    if a(i,1)<ks;
      short=short+1;
      a(i,1)=ks;
      diff=1;
      while diff>tol; 
        c1=c0-( (2*y(i)-c0)^(-gam2)*(a(i-1)*y(i)-c0)-delta* ..
           exp(X(i,:)*beta(:,2))*(ks-a(i-1)))/(gam2*(2*y(i)-c0)^(-gam2-1)* ..
           (a(i-1)*y(i)-c0)-(2*y(i)-c0)^(-gam2));
        diff=abs(c0-c1)/(1+abs(c0));
        c0=c1;
      end;
      c(i,1)=c1;
      p(i,1)=delta*exp(X(i,:)*beta(:,2))*(2*y(i)-c(i,1))^gam2; 
    end;
    if a(i,1)>2-ks;
      long=long+1;
      a(i,1)=2-ks;
      diff=1;
      while diff>tol; 
        c1=c0-(c0^(-gam1)*(a(i-1)*y(i)-c0)-delta*exp(X(i,:)*beta(:,1))* ..
           (2-ks-a(i-1)))/(-gam1*c0^(-gam1-1)*(a(i-1)*y(i)-c0)-c0^(-gam1));
        diff=abs(c0-c1)/(1+abs(c0));
        c0=c1;
      end;
      c(i,1)=c1;
      p(i,1)=delta*exp(X(i,:)*beta(:,1))*c(i,1)^gam1; 
    end;
    X(i,4)=a(i);
  end;
  sprintf('Number of times short= %g, long= %g in iteration %g\n', ..
           short,long,n)
  plot([c,y]),pause(2)
  plot(a)
  zxv(:,1)=c(3:T).^(-gam1).*(p(3:T)+y(3:T));
  zxv(:,2)=(2*y(3:T)-c(3:T)).^(-gam2).*(p(3:T)+y(3:T));
  zxv(:,3:5)=X(1:T-2,2:4);
  [betan,f,g]=uncmin(beta(:),'marcetf2');
  betan=[betan(1:4),betan(5:8)];
  if norm(betan-beta)<1e-5; break; end;
  [betan(:),beta(:)]
  beta=(1-lambda)*beta+lambda*betan;
end;
 


