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:
along with the first order/equilibrium conditions:
The dynare code for the perturbation method is as follows.
- Code: Select all
%
%%%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: http://fabcol.free.fr/index.php?page=notes, accessed Oct 21st '13):
- Code: Select all
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:
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