% 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= .99;
gam2= 1;
rho0=.1;
rho1=1.3;
rho2=-.4;

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

lambda=.5;

maxn=1;
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)),.0001*rand(T+1,1)];
[row,col]=size(X);
ain=4;
c=y;

beta=zeros(col,2);
beta(1,[1,2])=log(1.0/(1-delta))*ones(1,2);

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

for n=1:maxn;
  
  short=0;
  long=0;
  for i=2:T+1;   
    c0=y(i);
    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;
  end;

  if ain; X(i,ain)=a(i); end;

  sprintf('Number of times short= %g, long= %g in iteration %g\n', ..
           short,long,n)
  plot([c,y]),pause(2)
  plot([p./y(1:length(p))]),pause(2)
  plot([a])
  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).*(p(3:T)+y(3:T))-phi1+deriv1*gamma(:,1);
    y2=(2*y(3:T)-c(3:T)).^(-gam2).*(p(3:T)+y(3:T))-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;




