Remove unused and unmaintained unit test for k-order
Sébastien Villemot
parent
6f38dcd1d5
commit
911437378c
|
|
@ -1,129 +0,0 @@
|
|||
function [gy]=first_order(M_, dr, jacobia)
|
||||
% Emulation of Dynare++ c++ first_order.cpp for testing pruposes
|
||||
|
||||
% Copyright © 2009 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
% Dynare is free software: you can redistribute it and/or modify
|
||||
% it under the terms of the GNU General Public License as published by
|
||||
% the Free Software Foundation, either version 3 of the License, or
|
||||
% (at your option) any later version.
|
||||
%
|
||||
% Dynare is distributed in the hope that it will be useful,
|
||||
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
% GNU General Public License for more details.
|
||||
%
|
||||
% You should have received a copy of the GNU General Public License
|
||||
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
||||
|
||||
% 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];
|
||||
|
|
@ -1,119 +0,0 @@
|
|||
// This file replicates the estimation of the CIA model from
|
||||
// Frank Schorfheide (2000) "Loss function-based evaluation of DSGE models"
|
||||
// Journal of Applied Econometrics, 15, 645-670.
|
||||
// the data are the ones provided on Schorfheide's web site with the programs.
|
||||
// http://www.econ.upenn.edu/~schorf/programs/dsgesel.ZIP
|
||||
// You need to have fsdat.m in the same directory as this file.
|
||||
// This file replicates:
|
||||
// -the posterior mode as computed by Frank's Gauss programs
|
||||
// -the parameter mean posterior estimates reported in the paper
|
||||
// -the model probability (harmonic mean) reported in the paper
|
||||
// This file was tested with dyn_mat_test_0218.zip
|
||||
// the smooth shocks are probably stil buggy
|
||||
//
|
||||
// The equations are taken from J. Nason and T. Cogley (1994)
|
||||
// "Testing the implications of long-run neutrality for monetary business
|
||||
// cycle models" Journal of Applied Econometrics, 9, S37-S70.
|
||||
// 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;
|
||||
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;
|
||||
|
||||
alp = 0.33;
|
||||
bet = 0.99;
|
||||
gam = 0.003;
|
||||
mst = 1.011;
|
||||
rho = 0.7;
|
||||
psi = 0.787;
|
||||
del = 0.02;
|
||||
|
||||
model (use_dll);
|
||||
dA = exp(gam+e_a);
|
||||
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;
|
||||
R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
|
||||
1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
|
||||
c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
|
||||
P*c = m;
|
||||
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(-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;
|
||||
R = 1.02;
|
||||
k = 6;
|
||||
d = 0.85;
|
||||
n = 0.19;
|
||||
l = 0.86;
|
||||
y = 0.6;
|
||||
gy_obs = exp(gam);
|
||||
gp_obs = exp(-gam);
|
||||
dA = exp(gam);
|
||||
P2=P;
|
||||
c2=c;
|
||||
end;
|
||||
|
||||
shocks;
|
||||
var e_a; stderr 0.014;
|
||||
var e_m; stderr 0.005;
|
||||
end;
|
||||
|
||||
//unit_root_vars P_obs Y_obs;
|
||||
|
||||
steady(solve_algo = 2);
|
||||
|
||||
check;
|
||||
|
||||
estimated_params;
|
||||
alp, beta_pdf, 0.356, 0.02;
|
||||
bet, beta_pdf, 0.993, 0.002;
|
||||
gam, normal_pdf, 0.0085, 0.003;
|
||||
mst, normal_pdf, 1.0002, 0.007;
|
||||
rho, beta_pdf, 0.129, 0.223;
|
||||
psi, beta_pdf, 0.65, 0.05;
|
||||
del, beta_pdf, 0.01, 0.005;
|
||||
stderr e_a, inv_gamma_pdf, 0.035449, inf;
|
||||
stderr e_m, inv_gamma_pdf, 0.008862, inf;
|
||||
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);
|
||||
|
||||
|
|
@ -1,328 +0,0 @@
|
|||
/*
|
||||
* Copyright © 2008-2011 Dynare Team
|
||||
*
|
||||
* This file is part of Dynare.
|
||||
*
|
||||
* Dynare is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* Dynare is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License
|
||||
* along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
/*************************************
|
||||
* This main() is for testing k_order DLL entry point by linking to
|
||||
* the k_ord library statically and passing its hard-coded data:
|
||||
* parameters, covar, ysteady and the variable names from fs2000a.mod model
|
||||
* The main has been derived from mxFunction used for K-Order DLL
|
||||
***************************************/
|
||||
|
||||
// #include "stdafx.h"
|
||||
#include "dynamic_dll.h"
|
||||
#include "k_ord_dynare.h"
|
||||
|
||||
int
|
||||
main(int argc, char* argv[])
|
||||
{
|
||||
|
||||
double qz_criterium = 1 + 1e-6;
|
||||
const int check_flag = 0;
|
||||
const char* fName = "./fs2000k"; // mxArrayToString(mFname);
|
||||
const char* dfExt = ".mexa64"; // Dynamic file extension, e.g.".dll";
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_order_perturbation: check_flag = %d , fName = %s .\n", check_flag, fName);
|
||||
#endif
|
||||
int kOrder = 2;
|
||||
int npar = 7; //(int)mxGetM(mxFldp);
|
||||
double dparams[7] = {0.3300, 0.9900, 0.0030, 1.0110, 0.7000, 0.7870, 0.0200};
|
||||
Vector* modParams = new Vector(dparams, npar);
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_ord_perturbation: nParams=%d .\n", npar);
|
||||
for (int i = 0; i < npar; i++)
|
||||
{
|
||||
mexPrintf("k_ord_perturbation: dParams[%d]= %g.\n", i, dparams + i * (sizeof(double)));
|
||||
}
|
||||
for (int i = 0; i < npar; i++)
|
||||
{
|
||||
mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]);
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
double d2Dparams[4] = {//(double *) mxGetData(mxFldp);
|
||||
0.1960e-3, 0.0, 0.0, 0.0250e-3};
|
||||
npar = 2; //(int)mxGetN(mxFldp);
|
||||
TwoDMatrix* vCov = new TwoDMatrix(npar, npar, (d2Dparams));
|
||||
double dYSparams[16] = {
|
||||
// 27 mxGetData(mxFldp);
|
||||
// 1.0110, 2.2582, 5.8012, 0.5808,
|
||||
1.0110, 2.2582, 0.4477, 1.0000, 4.5959, 1.0212, 5.8012, 0.8494, 0.1872,
|
||||
0.8604, 1.0030, 1.0080, 0.5808, 1.0030, 2.2582, 0.4477
|
||||
//, 1.0110, 2.2582, 0.4477, 1.0000, 0.1872, 2.2582, 0.4477
|
||||
};
|
||||
const int nSteady = 16; // 27 //31;//29, 16 (int)mxGetM(mxFldp);
|
||||
Vector* ySteady = new Vector(dYSparams, nSteady);
|
||||
|
||||
double nnzd[3] = {77, 217, 0};
|
||||
const Vector* NNZD = new Vector(nnzd, 3);
|
||||
|
||||
// mxFldp = mxGetField(dr, 0,"nstatic" );
|
||||
const int nStat = 7; //(int)mxGetScalar(mxFldp);
|
||||
// mxFldp = mxGetField(dr, 0,"npred" );
|
||||
const int nPred = 2; // 6 - nBoth (int)mxGetScalar(mxFldp);
|
||||
// mxFldp = mxGetField(dr, 0,"nspred" );
|
||||
const int nsPred = 4; //(int)mxGetScalar(mxFldp);
|
||||
// mxFldp = mxGetField(dr, 0,"nboth" );
|
||||
const int nBoth = 2; // (int)mxGetScalar(mxFldp);
|
||||
// mxFldp = mxGetField(dr, 0,"nfwrd" );
|
||||
const int nForw = 5; // 3 (int)mxGetScalar(mxFldp);
|
||||
// mxFldp = mxGetField(dr, 0,"nsfwrd" );
|
||||
const int nsForw = 7; //(int)mxGetScalar(mxFldp);
|
||||
|
||||
// mxFldp = mxGetField(M_, 0,"exo_nbr" );
|
||||
const int nExog = 2; // (int)mxGetScalar(mxFldp);
|
||||
// mxFldp = mxGetField(M_, 0,"endo_nbr" );
|
||||
const int nEndo = 16; // 16(int)mxGetScalar(mxFldp);
|
||||
// mxFldp = mxGetField(M_, 0,"param_nbr" );
|
||||
const int nPar = 7; //(int)mxGetScalar(mxFldp);
|
||||
// it_ should be set to M_.maximum_lag
|
||||
// mxFldp = mxGetField(M_, 0,"maximum_lag" );
|
||||
const int nMax_lag = 1; //(int)mxGetScalar(mxFldp);
|
||||
|
||||
int var_order[] //[18]
|
||||
= {
|
||||
5, 6, 8, 10, 11, 12, 14, 7, 13,
|
||||
1, 2, 3, 4, 9, 15, 16
|
||||
// 5, 6, 8, 10, 11, 12, 16, 7, 13, 14, 15, 1, 2, 3, 4, 9, 17,
|
||||
// 18
|
||||
};
|
||||
// Vector * varOrder = new Vector(var_order, nEndo);
|
||||
vector<int>* var_order_vp = new vector<int>(nEndo); // nEndo));
|
||||
for (int v = 0; v < nEndo; v++)
|
||||
(*var_order_vp)[v] = var_order[v];
|
||||
|
||||
const double ll_incidence[] //[3][18]
|
||||
= {1, 5, 21, 2, 6, 22, 0, 7, 23, 0, 8, 24, 0, 9, 0, 0, 10, 0, 3, 11, 0, 0, 12, 0,
|
||||
0, 13, 25, 0, 14, 0, 0, 15, 0, 0, 16, 0, 4, 17, 0, 0, 18, 0, 0, 19, 26, 0, 20, 27};
|
||||
TwoDMatrix* llincidence = new TwoDMatrix(3, nEndo, ll_incidence);
|
||||
|
||||
const int jcols = nExog + nEndo + nsPred + nsForw; // Num of Jacobian columns
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_order_perturbation: jcols= %d .\n", jcols);
|
||||
#endif
|
||||
// mxFldp= mxGetField(M_, 0,"endo_names" );
|
||||
const int nendo = 16; // 16(int)mxGetM(mxFldp);
|
||||
const int widthEndo = 6; // (int)mxGetN(mxFldp);
|
||||
const char* cNamesCharStr = "mPceWRkdnlggydPc yp A22 __ oo "
|
||||
" bb ss ";
|
||||
// const char** endoNamesMX= DynareMxArrayToString( mxFldp,nendo,widthEndo);
|
||||
const char** endoNamesMX = DynareMxArrayToString(cNamesCharStr, nendo, widthEndo);
|
||||
#ifdef DEBUG
|
||||
for (int i = 0; i < nEndo; i++)
|
||||
{
|
||||
mexPrintf("k_ord_perturbation: EndoNameList[%d][0]= %s.\n", i, endoNamesMX[i]);
|
||||
}
|
||||
#endif
|
||||
// mxFldp = mxGetField(M_, 0,"exo_names" );
|
||||
const int nexo = 2; //(int)mxGetM(mxFldp);
|
||||
const int widthExog = 3; //(int)mxGetN(mxFldp);
|
||||
// const char** exoNamesMX= DynareMxArrayToString( mxFldp,nexo,widthExog);
|
||||
const char* cExoNamesCharStr = "ee__am";
|
||||
const char** exoNamesMX = DynareMxArrayToString(cExoNamesCharStr, nexo, widthExog);
|
||||
#ifdef DEBUG
|
||||
for (int i = 0; i < nexo; i++)
|
||||
{
|
||||
mexPrintf("k_ord_perturbation: ExoNameList[%d][0]= %s.\n", i, exoNamesMX[i]);
|
||||
}
|
||||
#endif
|
||||
if ((nEndo != nendo) || (nExog != nexo)) //(nPar != npar)
|
||||
{
|
||||
mexPrintf("Incorrect number of input parameters.\n");
|
||||
return;
|
||||
}
|
||||
|
||||
#ifdef DEBUG
|
||||
for (int i = 0; i < nEndo; i++)
|
||||
{
|
||||
mexPrintf("k_ord_perturbation: EndoNameList[%d]= %s.\n", i, endoNamesMX[i]);
|
||||
}
|
||||
// for (int i = 0; i < nPar; i++) {
|
||||
//, , mexPrintf("k_ord_perturbation: params_vec[%d]= %g.\n", i, params_vec[i] ); }
|
||||
for (int i = 0; i < nPar; i++)
|
||||
{
|
||||
mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]);
|
||||
}
|
||||
for (int i = 0; i < nSteady; i++)
|
||||
{
|
||||
mexPrintf("k_ord_perturbation: ysteady[%d]= %g.\n", i, (*ySteady)[i]);
|
||||
}
|
||||
|
||||
mexPrintf("k_order_perturbation: nEndo = %d , nExo = %d .\n", nEndo, nExog);
|
||||
#endif
|
||||
/* Fetch time index */
|
||||
// int it_ = (int) mxGetScalar(prhs[3]) - 1;
|
||||
|
||||
const int nSteps
|
||||
= 0; // Dynare++ solving steps, for time being default to 0 = deterministic steady state
|
||||
const double sstol = 1.e-13; // NL solver tolerance from
|
||||
|
||||
THREAD_GROUP::max_parallel_threads = 1; // 2 params.num_threads;
|
||||
|
||||
try
|
||||
{
|
||||
// make journal name and journal
|
||||
std::string jName(fName); // params.basename);
|
||||
jName += ".jnl";
|
||||
Journal journal(jName.c_str());
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_order_perturbation: Call tls init\n");
|
||||
#endif
|
||||
|
||||
tls.init(kOrder, (nStat + 2 * nPred + 3 * nBoth + 2 * nForw + nExog));
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_order_perturbation: Calling dynamicDLL constructor.\n");
|
||||
#endif
|
||||
// DynamicFn * pDynamicFn = loadModelDynamicDLL (fname);
|
||||
DynamicModelDLL dynamicDLL(fName, nEndo, jcols, nMax_lag, nExog, dfExt);
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_order_perturbation: Calling dynare constructor.\n");
|
||||
#endif
|
||||
// make KordpDynare object
|
||||
KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar, // paramNames,
|
||||
ySteady, vCov, modParams, nStat, nPred, nForw, nBoth, jcols, NNZD, nSteps,
|
||||
kOrder, journal, dynamicDLL, sstol, var_order_vp, // var_order
|
||||
llincidence, qz_criterium);
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_order_perturbation: Call Approximation constructor \n");
|
||||
#endif
|
||||
Approximation app(dynare, journal, nSteps, false, qz_criterium);
|
||||
// run stochastic steady
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_order_perturbation: Calling walkStochSteady.\n");
|
||||
#endif
|
||||
app.walkStochSteady();
|
||||
|
||||
// open mat file
|
||||
std::string matfile(fName); //(params.basename);
|
||||
matfile += ".mat";
|
||||
FILE* matfd = NULL;
|
||||
if (NULL == (matfd = fopen(matfile.c_str(), "wb")))
|
||||
{
|
||||
fprintf(stderr, "Couldn't open %s for writing.\n", matfile.c_str());
|
||||
exit(1);
|
||||
}
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_order_perturbation: Filling Mat file outputs.\n");
|
||||
#endif
|
||||
std::string ss_matrix_name(fName); // params.prefix);
|
||||
ss_matrix_name += "_steady_states";
|
||||
ConstTwoDMatrix(app.getSS()).writeMat4(matfd, ss_matrix_name.c_str());
|
||||
|
||||
// write the folded decision rule to the Mat-4 file
|
||||
app.getFoldDecisionRule().writeMat4(matfd, fName); // params.prefix);
|
||||
|
||||
fclose(matfd);
|
||||
|
||||
map<string, ConstTwoDMatrix> mm;
|
||||
app.getFoldDecisionRule().writeMMap(mm, string());
|
||||
#ifdef DEBUG
|
||||
app.getFoldDecisionRule().print();
|
||||
mexPrintf("k_order_perturbation: Map print: \n");
|
||||
for (map<string, ConstTwoDMatrix>::const_iterator cit = mm.begin(); cit != mm.end(); ++cit)
|
||||
{
|
||||
mexPrintf("k_order_perturbation: Map print: string: %s , g:\n", (*cit).first.c_str());
|
||||
(*cit).second.print();
|
||||
}
|
||||
#endif
|
||||
|
||||
// get latest ysteady
|
||||
double* dYsteady = (dynare.getSteady().base());
|
||||
ySteady = (Vector*)(&dynare.getSteady());
|
||||
}
|
||||
catch (const KordException& e)
|
||||
{
|
||||
printf("Caugth Kord exception: ");
|
||||
e.print();
|
||||
return 1; // e.code();
|
||||
}
|
||||
catch (const TLException& e)
|
||||
{
|
||||
printf("Caugth TL exception: ");
|
||||
e.print();
|
||||
return 2; // 255;
|
||||
}
|
||||
catch (SylvException& e)
|
||||
{
|
||||
printf("Caught Sylv exception: ");
|
||||
e.printMessage();
|
||||
return 3; // 255;
|
||||
}
|
||||
catch (const DynareException& e)
|
||||
{
|
||||
printf("Caught KordpDynare exception: %s\n", e.message());
|
||||
return 4; // 255;
|
||||
}
|
||||
catch (const ogu::Exception& e)
|
||||
{
|
||||
printf("Caught ogu::Exception: ");
|
||||
e.print();
|
||||
return 5; // 255;
|
||||
}
|
||||
|
||||
// bones for future developement of the Matlab output.
|
||||
const int nrhs = 5;
|
||||
const int nlhs = 2;
|
||||
|
||||
mxArray* prhs[nrhs];
|
||||
mxArray* plhs[nlhs];
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_order_perturbation: Filling MATLAB outputs.\n");
|
||||
#endif
|
||||
|
||||
double *dgy, *dgu, *ysteady;
|
||||
int nb_row_x;
|
||||
|
||||
ysteady = NULL;
|
||||
if (nlhs >= 1)
|
||||
{
|
||||
/* Set the output pointer to the output matrix ysteady. */
|
||||
plhs[0] = mxCreateDoubleMatrix(nEndo, 1, mxREAL);
|
||||
/* Create a C pointer to a copy of the output ysteady. */
|
||||
ysteady = mxGetPr(plhs[0]);
|
||||
}
|
||||
|
||||
dgy = NULL;
|
||||
if (nlhs >= 2)
|
||||
{
|
||||
/* Set the output pointer to the output matrix gy. */
|
||||
plhs[1] = mxCreateDoubleMatrix(nEndo, jcols, mxREAL);
|
||||
// plhs[1] = (double*)(gy->getData())->base();
|
||||
/* Create a C pointer to a copy of the output matrix gy. */
|
||||
dgy = mxGetPr(plhs[1]);
|
||||
}
|
||||
|
||||
dgu = NULL;
|
||||
if (nlhs >= 3)
|
||||
{
|
||||
/* Set the output pointer to the output matrix gu. */
|
||||
plhs[2] = mxCreateDoubleMatrix(nEndo, nExog, mxREAL);
|
||||
// plhs[2] = (double*)((gu->getData())->base());
|
||||
/* Create a C pointer to a copy of the output matrix gu. */
|
||||
dgu = mxGetPr(plhs[2]);
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
Loading…
Reference in New Issue