Update test model without "both" and without unit root variables and new Matlab test emulator of FirstOrder.cpp

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2409 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
george 2009-02-09 11:19:48 +00:00
parent b0ea237ae3
commit 24969d42a7
2 changed files with 125 additions and 8 deletions

View File

@ -0,0 +1,110 @@
function [gy]=first_order(M_, dr, jacobia)
% fd = jacobia_
% reorder jacobia_
[fd]=k_reOrderedJacobia(M_, jacobia)
%ypart=dr;
ypart.ny=M_.endo_nbr;
ypart.nyss=dr.nfwrd+dr.nboth;
ypart.nys=dr.npred;
ypart.npred=dr.npred-dr.nboth;
ypart.nboth=dr.nboth;
ypart.nforw=dr.nfwrd;
ypart.nstat =dr.nstatic
nu=M_.exo_nbr
off= 1; % = 0 in C
fyplus = fd(:,off:off+ypart.nyss-1);
off= off+ypart.nyss;
fyszero= fd(:,off:off+ypart.nstat-1);
off= off+ypart.nstat;
fypzero= fd(:,off:off+ypart.npred-1);
off= off+ypart.npred;
fybzero= fd(:,off:off+ypart.nboth-1);
off= off+ypart.nboth;
fyfzero= fd(:,off:off+ypart.nforw-1);
off= off+ypart.nforw;
fymins= fd(:,off:off+ypart.nys-1);
off= off+ypart.nys;
fuzero= fd(:,off:off+nu-1);
off=off+ nu;
n= ypart.ny+ypart.nboth;
%TwoDMatrix
matD=zeros(n,n);
% matD.place(fypzero,0,0);
matD(1:n-ypart.nboth,1:ypart.npred)= fypzero;
% matD.place(fybzero,0,ypart.npred);
matD(1:n-ypart.nboth,ypart.npred+1:ypart.npred+ypart.nboth)=fybzero;
% matD.place(fyplus,0,ypart.nys()+ypart.nstat);
matD(1:n-ypart.nboth,ypart.nys+ypart.nstat+1:ypart.nys+ypart.nstat+ypart.nyss)=fyplus;
for i=1:ypart.nboth
matD(ypart.ny()+i,ypart.npred+i)= 1.0;
end
matE=[fymins, fyszero, zeros(n-ypart.nboth,ypart.nboth), fyfzero; zeros(ypart.nboth,n)];
% matE.place(fymins;
% matE.place(fyszero,0,ypart.nys());
% matE.place(fyfzero,0,ypart.nys()+ypart.nstat+ypart.nboth);
for i= 1:ypart.nboth
matE(ypart.ny()+i,ypart.nys()+ypart.nstat+i)= -1.0;
end
matE=-matE; %matE.mult(-1.0);
% vsl=zeros(n,n);
% vsr=zeros(n,n);
% lwork= 100*n+16;
% work=zeros(1,lwork);
% bwork=zeros(1,n);
%int info;
% LAPACK_dgges("N","V","S",order_eigs,&n,matE.getData().base(),&n,
% matD.getData().base(),&n,&sdim,alphar.base(),alphai.base(),
% beta.base(),vsl.getData().base(),&n,vsr.getData().base(),&n,
% work.base(),&lwork,&(bwork[0]),&info);
[matE1,matD1,vsr,sdim,dr.eigval,info] = mjdgges(matE,matD,1);
bk_cond= (sdim==ypart.nys);
% ConstGeneralMatrix z11(vsr,0,0,ypart.nys(),ypart.nys());
z11=vsr(1:ypart.nys,1:ypart.nys);
% ConstGeneralMatrix z12(vsr,0,ypart.nys(),ypart.nys(),n-ypart.nys());
z12=vsr(1:ypart.nys(),ypart.nys+1:end);%, n-ypart.nys);
% ConstGeneralMatrix z21(vsr,ypart.nys(),0,n-ypart.nys(),ypart.nys());
z21=vsr(ypart.nys+1:end,1:ypart.nys);
% ConstGeneralMatrix z22(vsr,ypart.nys(),ypart.nys(),n-ypart.nys(),n-ypart.nys());
z22=vsr(ypart.nys+1:end,ypart.nys+1:end);
% GeneralMatrix sfder(z12,"transpose");
sfder=z12';%,"transpose");
% z22.multInvLeftTrans(sfder);
sfder=z22'\sfder;
sfder=-sfder;% .mult(-1);
%s11(matE,0,0,ypart.nys(),ypart.nys());
s11=matE1(1:ypart.nys,1:ypart.nys);
% t11=(matD1,0,0,ypart.nys(),ypart.nys());
t11=matD1(1:ypart.nys,1:ypart.nys);
dumm=(s11');%,"transpose");
%z11.multInvLeftTrans(dumm);
dumm=z11'\dumm;
preder=(dumm');%,"transpose");
%t11.multInvLeft(preder);
preder=t11\preder;
%preder.multLeft(z11);
preder= z11*preder;
% gy.place(preder,ypart.nstat,0);
% gy=(zeros(ypart.nstat,size(preder,2)) ;preder);
% sder(sfder,0,0,ypart.nstat,ypart.nys());
sder=sfder(1:ypart.nstat,1:ypart.nys);
% gy.place(sder,0,0);
% gy(1:ypart.nstat, 1:ypart.nys)=sder;
% gy=[sder;preder];
% fder(sfder,ypart.nstat+ypart.nboth,0,ypart.nforw,ypart.nys());
fder=sfder(ypart.nstat+ypart.nboth+1:ypart.nstat+ypart.nboth+ypart.nforw,1:ypart.nys);
% gy.place(fder,ypart.nstat+ypart.nys(),0);
% gy(ypart.nstat+ypart.nys,1)=fder;
gy=[sder;preder;fder];

View File

@ -17,9 +17,13 @@
// Note that there is an initial minus sign missing in equation (A1), p. S63.
//
// Michel Juillard, February 2004
// Modified for testing k_order_perturbation by GP, Jan-Feb 09
options_.usePartInfo=0;
//var m P c e W R k d n l gy_obs gp_obs Y_obs P_obs y dA P2 c2;
var m P c e W R k d n l gy_obs gp_obs y dA P2 c2;
options_.use_k_order=0;
//var m m_1 P P_1 c e W R k d n l gy_obs gp_obs Y_obs P_obs y dA P2 c2;
var m m_1 P P_1 c e W R k d n l gy_obs gp_obs y dA P2 c2;
varexo e_a e_m;
parameters alp bet gam mst rho psi del;
@ -34,7 +38,7 @@ del = 0.02;
model (use_dll);
dA = exp(gam+e_a);
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
log(m) = (1-rho)*log(mst) + rho*log(m_1(-1))+e_m;
-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c2(+1)*P2(+1)*m(+1))=0;
W = l/n;
-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
@ -46,16 +50,20 @@ m-1+d = l;
e = exp(e_a);
y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
gy_obs = dA*y/y(-1);
gp_obs = (P/P(-1))*m(-1)/dA;
gp_obs = (P/P_1(-1))*m_1(-1)/dA;
//Y_obs/Y_obs(-1) = gy_obs;
//P_obs/P_obs(-1) = gp_obs;
P2 = P(+1);
c2 = c(+1);
m_1 = m;
P_1 = P;
end;
initval;
m = mst;
m_1=mst;
P = 2.25;
P_1 = 2.25;
c = 0.45;
e = 1;
W = 4;
@ -64,11 +72,9 @@ k = 6;
d = 0.85;
n = 0.19;
l = 0.86;
y = 0.6;
gy_obs = exp(gam);
gp_obs = exp(-gam);
// Y_obs = 20000;
// P_obs = 51;
y = 0.6;
dA = exp(gam);
P2=P;
c2=c;
@ -100,12 +106,13 @@ end;
//varobs P_obs Y_obs;
varobs gp_obs gy_obs;
steady(solve_algo = 2);
//observation_trends;
//P_obs (log(mst)-gam);
//Y_obs (gam);
//end;
//options_.useAIM = 1;
estimation(datafile=fsdat,nobs=192,loglinear,mh_replic=2000,
mode_compute=4,mh_nblocks=2,mh_drop=0.45,mh_jscale=0.65);