% 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] + pb[t]b[i,t] = b[i,t-1] + y[t]
%            y[t] = rho0 + rho1 y[t-1] + rho2 y[t-2] + eps[t]
%            b[i,t] >= kb
%            b[i,0] = 0
%
% Equilibrium conditions:
%            
%            c[1,t] + c[2,t] = 2 y[t]
%            b[1,t] + b[2,t] = 0 
%

% Ellen R. McGrattan, 11-30-90

% Reference: Ketterer and Marcet (1989).

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

sige=sqrt(.0016);
kb=-1;

lambda=.1;

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

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)),.00001*rand(T+1,1)];
[row,col]=size(X);

c=y;

beta=zeros(col,2);
beta(1,[1,2])=log(1.0/(1-delta))*ones(1,2);
tem=[ ..
  -2.9740e-04   1.7508e-01
  -1.2585e+00  -1.1612e+00
   4.0345e-01   3.6488e-01
  -3.8506e-02  -5.1163e+00
  -2.9740e-04   1.7508e-01
  -1.2585e+00  -1.1612e+00
   4.0345e-01   3.6488e-01
  -3.8506e-02  -5.1163e+00];
beta(:,1)=tem(1:4,2);
beta(:,2)=tem(5:8,2);


b(1,1)=0;
pb(1,1)=delta;

for n=1:maxn;
  
  short=0;
  long=0;
  for i=2:T+1;   
    c0=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;
    pb(i,1)=delta*exp(X(i,:)*beta(:,2))*(2*y(i)-c(i,1))^gam2; 
    b(i,1)=(b(i-1,1)+y(i)-c(i))/pb(i);
    if b(i,1)<kb;
      short=short+1;
      b(i,1)=kb;
      diff=1;
      p0=pb(i-1);
      while diff>tol; 
        p1=p0-( (b(i-1)+y(i)-p0*kb)^(-gam1)*p0-delta*exp(X(i,:)*beta(:,1)))/ ..
           ( gam1*(b(i-1)+y(i)-p0*kb)^(-gam1-1)*p0*kb+(b(i-1)+y(i)-p0*kb)^(-gam1));
        diff=abs(p0-p1)/(1+abs(p0));
        p0=p1;
      end;
p1
      pb(i,1)=p1;
      c(i,1)=b(i)+y(i)-pb(i)*kb;
    end;
    if b(i,1)>-kb;
      long=long+1;
      b(i,1)=-kb;
      diff=1;
      p0=pb(i-1);
      while diff>tol; 
        p1=p0-( (b(i-1)+y(i)+p0*kb)^(-gam1)*p0-delta*exp(X(i,:)*beta(:,1)))/ ..
           ( -gam1*(b(i-1)+y(i)+p0*kb)^(-gam1-1)*p0*kb+(b(i-1)+y(i)+p0*kb)^(-gam1));
        diff=abs(p0-p1)/(1+abs(p0));
        p0=p1;
      end;
      pb(i,1)=p1;
      c(i,1)=b(i)+y(i)-pb(i)*kb;
    end;
    X(i,4)=b(i);
    if (c(i)==NaN); return; end;
  end;
  sprintf('Number of times short= %g, long= %g in iteration %g\n', ..
           short,long,n)
  plot([c,y]),pause(2)
  plot([pb./y(1:length(pb))]),pause(2)
  plot([b])
  gamma=beta;
  diff=1;
  while diff>tol
    phi1=exp(X(1:T-2,:)*gamma(:,1));
    phi2=exp(X(1:T-2,:)*gamma(:,2));
    deriv1=[phi1*ones(1,col)].*X(1:T-2,:);
    deriv2=[phi2*ones(1,col)].*X(1:T-2,:);
    y1=c(3:T).^(-gam1)-phi1+deriv1*gamma(:,1);
    y2=(2*y(3:T)-c(3:T)).^(-gam2)-phi2+deriv2*gamma(:,2);
    gammn(:,1)=(deriv1'*deriv1)\(deriv1'*y1);
    gammn(:,2)=(deriv2'*deriv2)\(deriv2'*y2);
    [gammn,gamma]
    diff=norm(gammn-gamma);
    gamma=gammn;
  end;
  if norm(gamma-beta)<tol; break; end;
  [gamma(:),beta(:)]
  beta=(1-lambda)*beta+lambda*gamma;
end;
 


