Comparing perturbation to VFI

Dynare help (legacy posts)
1

Hi,

I’m trying to compare the solutions to a simple neoclassical growth model using perturbation (1st and 2nd order) and value function iteration.
The model is specified as follows:
http://imageshack.us/a/img132/5375/sb5w.jpg
along with the first order/equilibrium conditions:
http://imageshack.us/a/img194/6687/yp8q.jpg

The dynare code for the perturbation method is as follows.

%
%%%PREAMBLE%%%
%
var y invest k c z;                          %Hier alle Variablen eintragen, die im Modell auftauchen: Output, Investment, Capital, Labour, InterestRate, WageRate%   
predetermined_variables k;                      %Hier nochmal alle vorherbestimmten Variablen bestimmen!!% 
varexo eps;                                        %Alle exogenen Variablen. Hier deterministischer Schock mit Epsilon=0, deshalb ist der Schock exogen. Für einen stochastischen Schock müsste man hier varexo epsilon eintragen%
parameters cbeta calpha crho ceta ctau;               %Alle Parameter. Das c vorher lieber noch dazufügen da beta selbst eine Funktion in Matlab ist%

cbeta = .99;                                     %Hier werden die Parameter kalibriert%
calpha = .33;
crho = 0.9;
ceta = 0.04;
ctau = 20;
%
%%%MODELBLOCK%%%
%
model;                                           %Hier werden die Gleichungen des Modells geschrieben%
c^(-ctau) = cbeta * c(+1)^(-ctau) * calpha * exp(z(+1)) * k(+1)^(calpha-1);    %Euler Gleichung. Das (+1) steht dafür, dass der Index der Variable t+1 ist%
k(+1) = exp(z) * k^calpha - c;
z = crho*z(-1) + eps;
y = exp(z) * k^calpha;
invest = y - c;
end;
%
%%%INITIAL VALUE BLOCK%%%
%
initval;                                        %Hier Werte nahe des Steady states angeben, Dynare findet dann den richtigen Steady state%
k = 0.2;
c = 0.77;
z = 0;
end;
steady;                                         %steady(solve_algo = 0) hier kann man Dynare sagen, welcher Algorithmus verwendet werden soll um den S.S. zu finden. Werte gehen von 0-3%
%
%%%SHOCKS%%%
%
shocks;                                         %Wir haben einen neuen Parameter csigma als Standardabweichung des Technologieschocks definiert.%
var eps = ceta^2;                               %VARIANZ hier eintragen%                                           
end;

stoch_simul(periods = 2000,order=1, irf=40);


state_range=0.08:0.01:0.25;
state_name='k';
plot_var_name='c';
plot_policy_fun(state_name,state_range,plot_var_name);

Now, the code for the value function iteration was written by Fabrice Collard (download: fabcol.free.fr/index.php?page=notes, accessed Oct 21st '13):

sigma = 20;
delta = 1;
beta = 0.95;
alpha = 0.30;
rho = 0.9;
se = 0.04;

nbk = 400;
nba = 2;
crit = 1;
epsi = 1e-6;
%
%Discretization of the shock
%
p = (1+rho)/2;
PI = [p 1-p; 1-p p];
se = 0.04;
ab = 0;
a1 = exp(-se*se/(1-rho*rho));
a2 = exp(se*se/(1-rho*rho));
A = [a1 a2];
%
%Discretization of the state space;
%
kmin = 0.08;
kmax = 0.25;
k = linspace(kmin,kmax,nbk)';
c = zeros(nbk,nba);
util = zeros(nbk,nba);
v = zeros(nbk,nba);
Tv = zeros(nbk,nba);
iter=1;
while crit>epsi;
    for i=1:nbk
        for j=1:nba
            c = A(j)*k(i).^alpha+(1-delta)*k(i)-k;
            neg = find(c<=0);
            c(neg) = NaN;
            util(:,j) = (c.^(1-sigma))/(1-sigma);
            util(neg,j) = -1e12;
        end
        [Tv(i,:),dr(i,:)] = max(util+beta*(v*PI));
    end;
    crit = max(max(abs(Tv-v)));
    v = Tv;
   iter=iter+1;
end
kp = k(dr);
for j=1:nba
    c(:,j) = A(j)*k.^alpha+(1-delta)*k-kp(:,j);
end

The notation differs a little bit from the dynare code. “se” denotes the standard error (in dynare called “ceta”) and “sigma” is the relative risk aversion parameter (in dynare called “ctau”).
Once I compare the figures of the policy functions for consumption I get really strange results:
http://imageshack.us/a/img853/9680/kz9k.jpg
In the value function iteration the AR(1) shock was discretized into two states (high/low). In my opinion the policy function curves obtained from VFI are way too high and should actually lie very close to the 2nd order perturbation solution.

I’d very much appreciate your comments on this.
Best regards,
nbt