Using the Endogenous Grid Method to solve the Euler Equation.

endogenous_grid_method.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Euler Equation Methods
% Renjie Ge

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; clc; close all; 

% Initialization

A=18.5;  alpha=0.3; beta=.9; J=100;
x1=0; xn=20; N=101; X=linspace(x1,xn,N);

F=@(x) A*x.^alpha;
Fprime=@(x) A*alpha*x.^(alpha-1);

% the euler equation iteration

%--- Method 1. endogenously backward solve for xt(xt-1)

% y=zeros(N,1)';
% uprime=@(x) 1/x;
% tic
% for j=1:J
%     for i=1:N
%         Vprime=(beta*Fprime(X(i)))/(F(X(i))-y(i));      % V'
%         NEWX(i)=((1/A)*(X(i)+1/Vprime))^(1/alpha);      % Solve non-linear equation
%     end
%     y=interp1(NEWX,X,X);                                % policy function iteration
% %     plot(NEWX,X,'o');
% end

%--- Method 2. follow the suggestion in the homework

j=0;
crit=1;
toler=1e-3;
y=zeros(N,1)';
uprime=@(x) 1./x;
tic;
while crit>toler
    Vprime=uprime(F(X)-y).*Fprime(X);
    pp=interp1(X,Vprime,'linear','pp');
    for i=1:N
        newy(i)=fsolve(@(x) uprime(F(X(i))-x)-beta*ppval(pp,x),8);  % ppval(pp,x) is the function V'(x)
    end
    crit=max(abs(newy-y));
    y=newy;
    j=j+1;
    if j>50, break, end
end


% value function

u=@(x) log(x);
Vy=zeros(N,1)';             % initial guess of Vy
c=F(X)-y;
c(find(c==0))=1e-100;      % deal with -inf of log 

for j=1:J
    Vx=u(c)+beta*Vy;
    Vy=interp1(X,Vx,y);         % update Vy given new Vx
end
time=toc

% analytical solution

E=(1/(1-beta))*(log(A*(1-alpha*beta)) + alpha*beta/(1-alpha*beta)*log(A*alpha*beta));
F=alpha/(1-alpha*beta);

X(1)=1e-100;                        % handle the -inf when taking log
vstar=E+F*log(X);                   % value functionin continuous time
ystar=alpha*beta*A*X.^alpha;        % policy function continuous time

figure
subplot(1,2,1),fplot(@(x)[E+F*log(x)],[0,20],'--');
hold on;
plot(X,Vx,'rx');
axis([0 20 31 34]);
xlabel('x');
ylabel('v(x)');
title('Continuous and Discrete Value Functions');

subplot(1,2,2),fplot(@(x)[alpha*beta*A*x^alpha],[0,20],'--');
hold on;
plot(X,y,'rx');
axis([0 20 0 12]);
xlabel('x');
ylabel('y');
title('Continuous and Discrete Policy Functions');

%----------------------- problem b compare homework 2 and 3----------------%

distance=1;
va=zeros(N,1);

% utility matrix
for i=1:N                           % Control: the uneaten part of the cake(y)
    for j=1:N                       % State: the cake left(x)
        c=A*X(j)^alpha-X(i);            
        if c>0
            util(i,j)=log(c);                % utility matrix, column vectors are different y
        else
            util(i,j)=-1000;                 % Punishment
        end
    end
end

% the value function iteration

v=zeros(N,1);
d=1;
toler=1e-6;

while d>toler
     for j=1:N                                              % for each column j choose i to maximize bellman equation
         V=@(x) -(interp1(X',util(:,j)+beta*v,x));          % generate a interpolated function using function handle
         [policy(j),tv(j)]=fminbnd(V,0,20);                 % maximize Bellman equation (choose y) at given x
     end
     d=max(abs(-tv'-v));
     v=-tv';
end

Diffv2=v'-vstar;
Diffv3=Vx-vstar;

Diffy2=policy-ystar;
Diffy3=y-ystar;

figure
subplot(1,2,1), plot(X,Diffv2,'b-o',X,Diffv3,'r-o');axis([0 20 -.0012 .0002]);
xlabel('x');
ylabel('v(x)-v*(x)');
title('Approximation Error of Value Functions');
legend('VF', 'EE');

subplot(1,2,2), plot(X,Diffy2,'b-o',X,Diffy3,'r-o');axis([0 20 -.2 .2]);
xlabel('x');
ylabel('y(x)-y*(x)');
title('Approximation Error of Policy Functions');
legend('VF', 'EE');