173 lines
5.3 KiB
Modula-2
173 lines
5.3 KiB
Modula-2
/*
|
||
Check the policy functions obtained by perturbation at a high approximation
|
||
order, using the Burnside (1998, JEDC) model (for which the analytical form of
|
||
the policy function is known).
|
||
|
||
As shown by Burnside, the policy function for yₜ is:
|
||
|
||
yₜ = βⁱ exp[aᵢ+bᵢ(xₜ−xₛₛ)]
|
||
|
||
|
||
where:
|
||
θ² ⎛ 2ρ 1−ρ²ⁱ⎞
|
||
— aᵢ = iθxₛₛ + σ² ─────── ⎢i − ────(1−ρⁱ) + ρ² ─────⎥
|
||
2(1−ρ)² ⎝ 1−ρ 1−ρ² ⎠
|
||
|
||
θρ
|
||
— bᵢ = ───(1−ρⁱ)
|
||
1−ρ
|
||
|
||
— xₛₛ is the steady state of x
|
||
— σ is the standard deviation of e.
|
||
|
||
With some algebra, it can be shown that the derivative of yₜ at the deterministic
|
||
steady state is equal to:
|
||
|
||
∂ᵐ⁺ⁿ⁺²ᵖ yₜ ∞ (2p)!
|
||
──────────────── = ∑ βⁱ bᵢᵐ⁺ⁿ ρᵐ ───── cᵢᵖ exp(iθxₛₛ)
|
||
∂ᵐxₜ₋₁ ∂ⁿeₜ ∂²ᵖs ⁱ⁼¹ p!
|
||
|
||
where:
|
||
— s is the stochastic scale factor
|
||
|
||
θ² ⎛ 2ρ 1−ρ²ⁱ⎞
|
||
— cᵢ = ─────── ⎢i − ────(1−ρⁱ) + ρ² ─────⎥
|
||
2(1−ρ)² ⎝ 1−ρ 1−ρ² ⎠
|
||
|
||
Note that derivatives with respect to an odd order for s (i.e. ∂²ᵖ⁺¹s) are always
|
||
equal to zero.
|
||
|
||
The policy function as returned in the oo_.dr.g_* matrices has the following properties:
|
||
— its elements are pre-multiplied by the Taylor coefficients;
|
||
— derivatives w.r.t. the stochastic scale factor have already been summed up;
|
||
— symmetric elements are folded (and they are not pre-multiplied by the number of repetitions).
|
||
|
||
As a consequence, the element gₘₙ corresponding to the m-th derivative w.r.t.
|
||
to xₜ₋₁ and the n-th derivative w.r.t. to eₜ is given by:
|
||
|
||
1 ∞ cᵢᵖ
|
||
gₘₙ = ────── ∑ ∑ βⁱ bᵢᵐ⁺ⁿ ρᵐ ──── exp(iθxₛₛ)
|
||
(m+n)! 0≤2p≤k-m-n ⁱ⁼¹ p!
|
||
|
||
where k is the order of approximation.
|
||
|
||
*/
|
||
|
||
@#define k =7
|
||
|
||
var y x;
|
||
varexo e;
|
||
|
||
parameters beta theta rho xbar;
|
||
xbar = 0.0179;
|
||
rho = -0.139;
|
||
theta = -1.5;
|
||
theta = -10;
|
||
beta = 0.95;
|
||
|
||
model;
|
||
y = beta*exp(theta*x(+1))*(1+y(+1));
|
||
x = (1-rho)*xbar + rho*x(-1)+e;
|
||
end;
|
||
|
||
shocks;
|
||
var e; stderr 0.0348;
|
||
end;
|
||
|
||
initval;
|
||
x = xbar;
|
||
y = beta*exp(theta*xbar)/(1-beta*exp(theta*xbar));
|
||
end;
|
||
|
||
steady;
|
||
|
||
stoch_simul(order=@{k},k_order_solver,irf=0,drop=0,periods=100);
|
||
|
||
% Verify that the policy function coefficients are correct
|
||
sigma2=M_.Sigma_e;
|
||
i = 1:800;
|
||
c = theta^2*sigma2/(2*(1-rho)^2)*(i-2*rho*(1-rho.^i)/(1-rho)+rho^2*(1-rho.^(2*i))/(1-rho^2));
|
||
b = theta*rho*(1-rho.^i)/(1-rho);
|
||
|
||
for ord = 0:@{k}
|
||
g = oo_.dr.(['g_' num2str(ord)])(2,:); % Retrieve computed policy function for variable y
|
||
|
||
for m = 0:ord % m is the derivation order with respect to x(-1)
|
||
v = 0;
|
||
for p = 0:floor((@{k}-ord)/2) % 2p is the derivation order with respect to s
|
||
if ord+2*p > 0 % Skip the deterministic steady state constant
|
||
v = v + sum(beta.^i.*exp(theta*xbar*i).*b.^ord.*rho^m.*c.^p)/factorial(ord)/factorial(p);
|
||
end
|
||
end
|
||
if abs(v-g(ord+1-m)) > 1e-14
|
||
error(['Error in matrix oo_.dr.g_' num2str(ord)])
|
||
end
|
||
end
|
||
end
|
||
|
||
% Verify that the simulated time series is correct
|
||
xss = oo_.steady_state(2);
|
||
xlag = xss;
|
||
|
||
for T = 1:size(oo_.endo_simul,2)
|
||
e_ = oo_.exo_simul(T);
|
||
y_ = oo_.steady_state(1);
|
||
for ord = 0:@{k}
|
||
g = oo_.dr.(['g_' num2str(ord)])(2,:); % Retrieve computed policy function for variable y
|
||
for m = 0:ord
|
||
y_ = y_ + g(ord+1-m)*(xlag-xss)^m*e_^(ord-m)*nchoosek(ord,m);
|
||
end
|
||
end
|
||
if abs(y_-oo_.endo_simul(1,T)) > 1e-14
|
||
error(['Error in dynare_simul_ DLL'])
|
||
end
|
||
xlag = oo_.endo_simul(2,T);
|
||
end
|
||
|
||
% Verify that the simulated time series is correct with the Fortran routine k_order_simul
|
||
|
||
order = options_.order;
|
||
nstat = M_.nstatic;
|
||
npred = M_.npred;
|
||
nboth = M_.nboth;
|
||
nfwrd = M_.nfwrd;
|
||
nexog = M_.exo_nbr;
|
||
ystart = oo_.dr.ys(oo_.dr.order_var,1);
|
||
ex_ = oo_.exo_simul';
|
||
ysteady = oo_.dr.ys(oo_.dr.order_var);
|
||
dr = oo_.dr;
|
||
|
||
vcov = M_.Sigma_e;
|
||
seed = options_.DynareRandomStreams;
|
||
|
||
%do 1 call to get k_order_simul into memory (dynare_simul_ was already called within stoch_simul)
|
||
fortran_endo_simul = k_order_simul(order, nstat, npred, nboth, nfwrd, nexog, ystart, ex_, ysteady, dr);
|
||
|
||
tElapsed=NaN(2,1000);
|
||
for ii=1:1000
|
||
|
||
tStart1 = tic; fortran_endo_simul = k_order_simul(order, nstat, npred, nboth, nfwrd, nexog, ystart, ex_, ysteady, dr); tElapsed(1,ii) = toc(tStart1);
|
||
tStart2 = tic; dynare_endo_simul = dynare_simul_(order, nstat, npred, nboth, nfwrd, nexog, ystart,ex_,vcov,seed, ysteady, dr); tElapsed(2,ii) = toc(tStart2);
|
||
end
|
||
|
||
|
||
if max(max(abs(oo_.endo_simul-fortran_endo_simul(oo_.dr.order_var,:)))) > 1e-10
|
||
error('Error in k_order_simul: inaccurate simulation');
|
||
end;
|
||
|
||
if max(max(abs(dynare_endo_simul-fortran_endo_simul))) > 1e-10
|
||
error('Error in k_order_simul: inaccurate simulation');
|
||
end;
|
||
|
||
tElapsed1=mean(tElapsed(1,:))
|
||
tElapsed2=mean(tElapsed(2,:))
|
||
if tElapsed1<tElapsed2
|
||
skipline()
|
||
dprintf('k_order_simul is %5.2f times faster than dynare_simul_', tElapsed2/tElapsed1)
|
||
skipline()
|
||
else
|
||
skipline()
|
||
dprintf('k_order_simul is %5.2f times slower than dynare_simul_', tElapsed1/tElapsed2)
|
||
skipline()
|
||
end
|