Investigate why mod-file using k_order_pert crashes Matlab
The following mod-file crashes Matlab on different Windows and Dynare versions (http://www.dynare.org/phpBB3/viewtopic.php?f=1&t=6306):
%------------------------------------------------------------------------------------------%
% VARIABLES ENDO %
%------------------------------------------------------------------------------------------%
var y c i L k khat u betateta q teta w Rk Qtilde Omega pi pireset X1 X2 mc v r rr Spr A ;
%------------------------------------------------------------------------------------------%
% VARIABLES EXO %
%------------------------------------------------------------------------------------------%
varexo
eteta % disaster probability shock
er % monetary policy shock
;
%------------------------------------------------------------------------------------------%
% PARAMETERS %
%------------------------------------------------------------------------------------------%
parameters muz, beta0, delta0, zeta, upsilon, sigxi, muxi, gamma, psi, psitilde, chi, alpha, varpi, tau,
rhor, rhoteta, phipi, phiy, sigr, sigteta, sigz, eta, yss, rss, iss, kss, piss, tetass, uss, betatetass, bk;
muz = 0.005;
delta0 = 0.02;
zeta = 0.8;
upsilon = 6;
sigxi = 0.092;
muxi = -0.007;
gamma = 3.8;
alpha = 0.33;
tau = 1.7;
rhor =0.85;
rhoteta =0.85;
phipi = 1.5;
phiy = 0.5;
sigr = 0.01;
sigteta = 0.01;
sigz = 0;
psitilde = 0.5;
bk = 0.43;
piss = 1.005;
tetass = 0.0072;
uss = 1;
qss = 1;
varpi = 2.33;
psi=1-((1-psitilde)/(1+varpi));
chi = 1- (1-gamma)/(1-psi);
beta0 = 0.999;
betatetass = beta0 * ((1 - tetass + tetass * exp((1-gamma)*log(1-bk)))^(1/(1-chi)));
%------------------------------------------------------------------------------------------%
% STEADY STATE PARAMETERS %
%------------------------------------------------------------------------------------------%
Qtildess=betatetass*exp((1-psi)*muz);
rss=exp(muz)*(piss)/Qtildess;
eta=(exp(muz)/(Qtildess*(1-tetass*bk))-1)*1/delta0 +1;
rkss=delta0*eta;
piresetss=(((piss)^(1-upsilon)-zeta)/(1-zeta))^(1/(1-upsilon));
Ass=exp(-muz)*(1 - tetass + tetass/(1-bk));
mcss=((upsilon-1)/upsilon)*(1/(piss))*((1-(zeta*Ass*Qtildess*(piss)^upsilon))/(1-(zeta*Ass*Qtildess*(piss)^(upsilon-1))))*(piresetss);
klss=(mcss*alpha/rkss)^(1/(1-alpha));
wss=mcss*(1-alpha)*(klss)^alpha;
omegass=((1-zeta)*(piresetss)^(-upsilon)*(piss)^upsilon)/(1-zeta*(piss)^upsilon);
clss=(1/omegass)*(klss)^alpha - klss*(exp(muz)/(1-tetass*bk)-1+delta0);
lss=1/(1+(varpi/wss)*clss);
kss=klss*lss;
khatss=kss;
css=clss*lss;
ylss=clss+klss*(exp(muz)/(1-tetass*bk)-1+delta0);
yss=ylss*lss;
iss=yss-css;
X1ss=yss*mcss/(1-zeta*Ass*Qtildess*(piss)^upsilon);
X2ss=yss/(1-zeta*Ass*Qtildess*(piss)^(upsilon-1));
vss=(css*(1-lss)^varpi)^(1-psi) /(1-betatetass*exp(muz*(1-psi)));
rrss = rss/piss;
sprss = rkss/rrss;
%------------------------------------------------------------------------------------------%
% MODEL %
%------------------------------------------------------------------------------------------%
model;
// exogenous variables
log(teta)=(1-rhoteta)*log(tetass)+ rhoteta*log(teta(-1))+ sigteta*eteta;
// Households
(1-L)/c = varpi / w;
pi(+1)/r*exp(muz+sigz) = Qtilde;
Qtilde = (c(+1)/c)^(-psi) * ((1-L(+1))/(1-L))^(varpi*(1-psi)) * betateta * exp((1-psi)*muz)* (exp((1-gamma)*sigz)*v(+1)^(-chi))/(exp((1-gamma)*sigz)*v(+1)^(1-chi))^(-chi/(1-chi));
(1-teta*bk)* (Qtilde*Rk(+1)*(u(+1) + 1/(delta0*eta*u(+1)^(eta-1)) * (1-delta0 * u(+1)^eta +tau*i(+1)/k * (i(+1)/k - iss/kss) - tau/2* (i(+1)/k - iss/kss)^2))) = Rk /(delta0*eta*u^(eta-1)) * exp(muz+sigz);
delta0*eta*u^(eta-1)/Rk = 1-tau*(i/k(-1) - iss/kss);
v = (c * (1-L)^varpi)^(1-psi) + betateta * exp((1-psi)*muz) * (exp((1-psi)*(1-chi)*sigz) * v(+1)^(1-chi))^(1/(1-chi));
k = (1-teta*bk)*(((1-delta0*u^eta) * k(-1) + (i/k(-1) - tau/2*((i/k(-1) - iss/kss)^2))*k(-1)))/(exp(muz+sigz));
q = 1/(1-tau*(i/k(-1) - iss/kss));
betateta = beta0 * ((1 - teta + teta * exp((1-gamma)*log(1-bk)))^(1/(1-chi)));
// Firms
y = khat^alpha * L^(1-alpha) / Omega;
w = mc * (1 - alpha) * (khat/L)^alpha;
Rk = mc * alpha * (khat/L)^(alpha-1);
khat = u*k(-1);
// price setting
A = exp(-muz) * (1 - teta + teta/(1-bk)) * exp(-sigz);
pireset = pi*upsilon/(upsilon-1)*X1/X2;
X1 = y*mc + zeta*A*Qtilde*X1(+1)*pi(+1)^upsilon;
X2 = y + zeta*A*Qtilde*X2(+1)*pi(+1)^(upsilon-1);
Omega = (1-zeta)*pireset^(-upsilon)*pi^upsilon + zeta*pi^upsilon*Omega(-1);
pi^(1-upsilon) = (1-zeta)*pireset^(1-upsilon) + zeta;
//Spreads
rr = r/pi;
Spr = Rk/rr;
// markets
y = c + i;
// Monetary policy
r= rhor*r(-1) + (1-rhor)*( phipi*(pi-piss) + phiy*(y-yss) + rss ) + sigr*er;
end;
%------------------------------------------------------------------------------------------%
% STEADY STATE INITIAL VALUES %
%------------------------------------------------------------------------------------------%
initval;
y = yss;
c = css;
i = iss;
A = Ass;
L = lss;
u = uss;
k = kss;
khat = khatss;
Omega = omegass;
pi = piss;
pireset = piresetss;
X1 = X1ss;
X2 = X2ss;
mc = mcss;
w = wss;
Rk = rkss;
teta = tetass;
r = rss;
q = qss;
Qtilde = Qtildess;
betateta = betatetass;
v = vss;
rr = rrss;
Spr = sprss;
end;
shocks;
var eteta = sigteta^2;
var er = sigr^2;
end;
resid(1);
check;
steady;
stoch_simul(order=3,k_order_solver, IRF=20) ;