k_order_perturbation.cc 12.6 KB
Newer Older
george's avatar
george committed
1
/*
2
 * Copyright (C) 2008-2012 Dynare Team
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 *
 * 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 <http://www.gnu.org/licenses/>.
 */
george's avatar
george committed
19

20
21
22
23
24
25
26
27
/*
  Defines the entry point for the k-order perturbation application DLL.

  Inputs:
  1) dr
  2) M_
  3) options

28
29
30
31
  Outputs:
  - if order == 1: only g_1
  - if order == 2: g_0, g_1, g_2
  - if order == 3: g_0, g_1, g_2, g_3
32
*/
george's avatar
george committed
33

34
#include "dynamic_m.hh"
35
#include "dynamic_dll.hh"
george's avatar
george committed
36

sebastien's avatar
sebastien committed
37
38
#include <cmath>
#include <cstring>
george's avatar
george committed
39
#include <cctype>
40
#include <cassert>
george's avatar
george committed
41

42
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)  // exclude mexFunction for other applications
sebastien's avatar
sebastien committed
43

44
# include "dynmex.h"
45

sebastien's avatar
sebastien committed
46
//////////////////////////////////////////////////////
47
// Convert MATLAB Dynare endo and exo names array to a vector<string> array of string pointers
sebastien's avatar
sebastien committed
48
49
50
// Poblem is that Matlab mx function returns a long string concatenated by columns rather than rows
// hence a rather low level approach is needed
///////////////////////////////////////////////////////
51
52
void
DynareMxArrayToString(const mxArray *mxFldp, const int len, const int width, vector<string> &out)
sebastien's avatar
sebastien committed
53
54
{
  char *cNamesCharStr = mxArrayToString(mxFldp);
sebastien's avatar
sebastien committed
55

56
57
58
59
60
61
62
  out.resize(len);

  for (int i = 0; i < width; i++)
    for (int j = 0; j < len; j++)
      // Allow alphanumeric and underscores "_" only:
      if (isalnum(cNamesCharStr[j+i*len]) || (cNamesCharStr[j+i*len] == '_'))
        out[j] += cNamesCharStr[j+i*len];
sebastien's avatar
sebastien committed
63
64
}

65
66
67
68
69
70
71
72
73
74
void copy_derivatives(mxArray *destin,const Symmetry &sym,const FGSContainer *derivs,const std::string &fieldname)
{
  const TwoDMatrix* x = derivs->get(sym);
  int n = x->numRows();
  int m = x->numCols();
  mxArray *tmp = mxCreateDoubleMatrix(n, m, mxREAL);
  memcpy(mxGetPr(tmp),x->getData().base(),n*m*sizeof(double));
  mxSetField(destin,0,fieldname.c_str(),tmp);
}

george's avatar
george committed
75
extern "C" {
76
77
78
79

  void
  mexFunction(int nlhs, mxArray *plhs[],
              int nrhs, const mxArray *prhs[])
george's avatar
george committed
80
  {
81
82
    if (nrhs < 3 || nlhs < 2)
      DYN_MEX_FUNC_ERR_MSG_TXT("Must have at least 3 input parameters and takes at least 2 output parameters.");
sebastien's avatar
sebastien committed
83

84
    const mxArray *dr = prhs[0];
85
86
    const mxArray *M_ = prhs[1];
    const mxArray *options_ = prhs[2];
87
    int use_dll = (int) mxGetScalar(mxGetField(options_, 0, "use_dll"));
88
89
90

    mxArray *mFname = mxGetField(M_, 0, "fname");
    if (!mxIsChar(mFname))
91
92
      DYN_MEX_FUNC_ERR_MSG_TXT("Input must be of type char.");

sebastien's avatar
sebastien committed
93
    string fName = mxArrayToString(mFname);
94

george's avatar
george committed
95
    int kOrder;
96
    mxArray *mxFldp = mxGetField(options_, 0, "order");
george's avatar
george committed
97
    if (mxIsNumeric(mxFldp))
98
      kOrder = (int) mxGetScalar(mxFldp);
george's avatar
george committed
99
100
    else
      kOrder = 1;
sebastien's avatar
sebastien committed
101

102
103
104
105
    //if (kOrder == 1 && nlhs != 2)
    //  DYN_MEX_FUNC_ERR_MSG_TXT("k_order_perturbation at order 1 requires exactly 2 arguments in output");
    //else if (kOrder > 1 && nlhs != kOrder+2)
    //  DYN_MEX_FUNC_ERR_MSG_TXT("k_order_perturbation at order > 1 requires exactly order+2 arguments in output");
sebastien's avatar
sebastien committed
106

george's avatar
george committed
107
    double qz_criterium = 1+1e-6;
108
    mxFldp = mxGetField(options_, 0, "qz_criterium");
george's avatar
george committed
109
    if (mxIsNumeric(mxFldp))
110
111
      qz_criterium = (double) mxGetScalar(mxFldp);

112
    mxFldp = mxGetField(M_, 0, "params");
113
    double *dparams = mxGetPr(mxFldp);
114
    int npar = (int) mxGetM(mxFldp);
115
    Vector modParams(dparams, npar);
116
117
    if (!modParams.isFinite())
      DYN_MEX_FUNC_ERR_MSG_TXT("The parameters vector contains NaN or Inf");
118

119
    mxFldp = mxGetField(M_, 0, "Sigma_e");
120
    dparams = mxGetPr(mxFldp);
121
    npar = (int) mxGetN(mxFldp);
122
    TwoDMatrix vCov(npar, npar, dparams);
123
124
    if (!vCov.isFinite())
      DYN_MEX_FUNC_ERR_MSG_TXT("The covariance matrix of shocks contains NaN or Inf");
125

126
    mxFldp = mxGetField(dr, 0, "ys");  // and not in order of dr.order_var
127
    dparams = mxGetPr(mxFldp);
128
    const int nSteady = (int) mxGetM(mxFldp);
129
    Vector ySteady(dparams, nSteady);
130
131
    if (!ySteady.isFinite())
      DYN_MEX_FUNC_ERR_MSG_TXT("The steady state vector contains NaN or Inf");
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151

    mxFldp = mxGetField(dr, 0, "nstatic");
    const int nStat = (int) mxGetScalar(mxFldp);
    mxFldp = mxGetField(dr, 0, "npred");
    int nPred = (int) mxGetScalar(mxFldp);
    mxFldp = mxGetField(dr, 0, "nspred");
    const int nsPred = (int) mxGetScalar(mxFldp);
    mxFldp = mxGetField(dr, 0, "nboth");
    const int nBoth = (int) mxGetScalar(mxFldp);
    mxFldp = mxGetField(dr, 0, "nfwrd");
    const int nForw = (int) mxGetScalar(mxFldp);
    mxFldp = mxGetField(dr, 0, "nsfwrd");
    const int nsForw = (int) mxGetScalar(mxFldp);

    mxFldp = mxGetField(M_, 0, "exo_nbr");
    const int nExog = (int) mxGetScalar(mxFldp);
    mxFldp = mxGetField(M_, 0, "endo_nbr");
    const int nEndo = (int) mxGetScalar(mxFldp);
    mxFldp = mxGetField(M_, 0, "param_nbr");
    const int nPar = (int) mxGetScalar(mxFldp);
152

george's avatar
george committed
153
    nPred -= nBoth; // correct nPred for nBoth.
154

155
    mxFldp = mxGetField(dr, 0, "order_var");
156
    dparams = mxGetPr(mxFldp);
157
    npar = (int) mxGetM(mxFldp);
158
    if (npar != nEndo)
159
160
      DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect number of input var_order vars.");

161
    vector<int> var_order_vp(nEndo);
162
    for (int v = 0; v < nEndo; v++)
163
      var_order_vp[v] = (int) (*(dparams++));
164

george's avatar
george committed
165
    // the lag, current and lead blocks of the jacobian respectively
166
    mxFldp = mxGetField(M_, 0, "lead_lag_incidence");
167
    dparams = mxGetPr(mxFldp);
168
169
    npar = (int) mxGetN(mxFldp);
    int nrows = (int) mxGetM(mxFldp);
170

171
    TwoDMatrix llincidence(nrows, npar, dparams);
172
    if (npar != nEndo)
173
174
175
176
177
      {
        ostringstream strstrm;
        strstrm << "dynare:k_order_perturbation " << "Incorrect length of lead lag incidences: ncol=" << npar << " != nEndo=" << nEndo;
        DYN_MEX_FUNC_ERR_MSG_TXT(strstrm.str().c_str());
      }
sebastien's avatar
sebastien committed
178
    //get NNZH =NNZD(2) = the total number of non-zero Hessian elements
179
    mxFldp = mxGetField(M_, 0, "NNZDerivatives");
180
    dparams = mxGetPr(mxFldp);
181
    Vector NNZD(dparams, (int) mxGetM(mxFldp));
182
    if (NNZD[kOrder-1] == -1)
183
184
      DYN_MEX_FUNC_ERR_MSG_TXT("The derivatives were not computed for the required order. Make sure that you used the right order option inside the stoch_simul command");

george's avatar
george committed
185
    const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
186
187
188
189

    mxFldp = mxGetField(M_, 0, "var_order_endo_names");
    const int nendo = (int) mxGetM(mxFldp);
    const int widthEndo = (int) mxGetN(mxFldp);
190
191
    vector<string> endoNames;
    DynareMxArrayToString(mxFldp, nendo, widthEndo, endoNames);
192

193
    mxFldp = mxGetField(M_, 0, "exo_names");
194
195
    const int nexo = (int) mxGetM(mxFldp);
    const int widthExog = (int) mxGetN(mxFldp);
196
197
    vector<string> exoNames;
    DynareMxArrayToString(mxFldp, nexo, widthExog, exoNames);
198

199
    if ((nEndo != nendo) || (nExog != nexo))
200
      DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect number of input parameters.");
201

202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
    TwoDMatrix *g1m=NULL;
    TwoDMatrix *g2m=NULL;
    TwoDMatrix *g3m=NULL;
    // derivatives passed as arguments */
    if (nrhs > 3) 
      {
	const mxArray *g1 = prhs[3];
	int m = (int) mxGetM(g1);
	int n = (int) mxGetN(g1);
	g1m = new TwoDMatrix(m, n, mxGetPr(g1));
	if (nrhs > 4)
	  {
	    const mxArray *g2 = prhs[4];
	    int m = (int) mxGetM(g2);
	    int n = (int) mxGetN(g2);
	    g2m = new TwoDMatrix(m, n, mxGetPr(g2));
	    if (nrhs > 5)
	      {
		const mxArray *g3 = prhs[5];
		int m = (int) mxGetM(g3);
		int n = (int) mxGetN(g3);
		g3m = new TwoDMatrix(m, n, mxGetPr(g3));
	      }
	  }
      }
	    
228
229
230
231
232
233
234
235
236
237
238
239

    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 = 2; //params.num_threads;

    try
      {
        // make journal name and journal
        std::string jName(fName); //params.basename);
        jName += ".jnl";
        Journal journal(jName.c_str());
240
241
242

        DynamicModelAC *dynamicModelFile;
        if (use_dll == 1)
243
          dynamicModelFile = new DynamicModelDLL(fName);
244
245
        else
          dynamicModelFile = new DynamicModelMFile(fName);
246
247
248
249
250

        // intiate tensor library
        tls.init(kOrder, nStat+2*nPred+3*nBoth+2*nForw+nExog);

        // make KordpDynare object
251
        KordpDynare dynare(endoNames, nEndo, exoNames, nExog, nPar,
252
                           ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,
253
                           jcols, NNZD, nSteps, kOrder, journal, dynamicModelFile,
254
255
                           sstol, var_order_vp, llincidence, qz_criterium,
			   g1m, g2m, g3m);
256
257

        // construct main K-order approximation class
sebastien's avatar
sebastien committed
258

259
260
261
262
        Approximation app(dynare, journal,  nSteps, false, qz_criterium);
        // run stochastic steady
        app.walkStochSteady();

263
264
	app.get_rule_ders()->print();

265
266
        /* Write derivative outputs into memory map */
        map<string, ConstTwoDMatrix> mm;
sebastien's avatar
sebastien committed
267
        app.getFoldDecisionRule().writeMMap(mm, string());
268
269

        // get latest ysteady
270
        ySteady = dynare.getSteady();
271

272
        if (kOrder == 1)
273
274
          {
            /* Set the output pointer to the output matrix ysteady. */
sebastien's avatar
sebastien committed
275
276
            map<string, ConstTwoDMatrix>::const_iterator cit = mm.begin();
            ++cit;
277
            plhs[1] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL);
278
279
280
281

            // Copy Dynare++ matrix into MATLAB matrix
            const ConstVector &vec = (*cit).second.getData();
            assert(vec.skip() == 1);
282
            memcpy(mxGetPr(plhs[1]), vec.base(), vec.length() * sizeof(double));
283
          }
284
        if (kOrder >= 2)
285
          {
286
            int ii = 1;
287
288
289
            for (map<string, ConstTwoDMatrix>::const_iterator cit = mm.begin();
                 ((cit != mm.end()) && (ii < nlhs)); ++cit)
              {
290
291
292
293
294
295
296
297
298
		plhs[ii] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL);
		
		// Copy Dynare++ matrix into MATLAB matrix
		const ConstVector &vec = (*cit).second.getData();
		assert(vec.skip() == 1);
		memcpy(mxGetPr(plhs[ii]), vec.base(), vec.length() * sizeof(double));
		
		++ii;
		
299
              }
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
	    if (kOrder == 3 && nlhs > 4)
	      {
		const FGSContainer *derivs = app.get_rule_ders();
		const std::string fieldnames[] = {"gy", "gu", "gyy", "gyu", "guu", "gss", 
				      "gyyy", "gyyu", "gyuu", "guuu", "gyss", "guss"};
		// creates the char** expected by mxCreateStructMatrix()
		const char* c_fieldnames[12];
		for (int i=0; i < 12;++i)
		  c_fieldnames[i] = fieldnames[i].c_str();
		plhs[ii] = mxCreateStructMatrix(1,1,12,c_fieldnames);
		copy_derivatives(plhs[ii],Symmetry(1,0,0,0),derivs,"gy");
		copy_derivatives(plhs[ii],Symmetry(0,1,0,0),derivs,"gu");
		copy_derivatives(plhs[ii],Symmetry(2,0,0,0),derivs,"gyy");
		copy_derivatives(plhs[ii],Symmetry(0,2,0,0),derivs,"guu");
		copy_derivatives(plhs[ii],Symmetry(1,1,0,0),derivs,"gyu");
		copy_derivatives(plhs[ii],Symmetry(0,0,0,2),derivs,"gss");
		copy_derivatives(plhs[ii],Symmetry(3,0,0,0),derivs,"gyyy");
		copy_derivatives(plhs[ii],Symmetry(0,3,0,0),derivs,"guuu");
		copy_derivatives(plhs[ii],Symmetry(2,1,0,0),derivs,"gyyu");
		copy_derivatives(plhs[ii],Symmetry(1,2,0,0),derivs,"gyuu");
		copy_derivatives(plhs[ii],Symmetry(1,0,0,2),derivs,"gyss");
		copy_derivatives(plhs[ii],Symmetry(0,1,0,2),derivs,"guss");
	      }
	  }
george's avatar
george committed
324
      }
325
    catch (const KordException &e)
george's avatar
george committed
326
      {
327
        e.print();
328
329
330
        ostringstream strstrm;
        strstrm << "dynare:k_order_perturbation: Caught Kord exception: " << e.get_message();
        DYN_MEX_FUNC_ERR_MSG_TXT(strstrm.str().c_str());
george's avatar
george committed
331
      }
332
333
334
    catch (const TLException &e)
      {
        e.print();
335
        DYN_MEX_FUNC_ERR_MSG_TXT("dynare:k_order_perturbation: Caught TL exception");
george's avatar
george committed
336
      }
337
338
339
    catch (SylvException &e)
      {
        e.printMessage();
340
        DYN_MEX_FUNC_ERR_MSG_TXT("dynare:k_order_perturbation: Caught Sylv exception");
george's avatar
george committed
341
      }
342
343
    catch (const DynareException &e)
      {
344
345
346
        ostringstream strstrm;
        strstrm << "dynare:k_order_perturbation: Caught KordDynare exception: " << e.message();
        DYN_MEX_FUNC_ERR_MSG_TXT(strstrm.str().c_str());
george's avatar
george committed
347
      }
348
349
    catch (const ogu::Exception &e)
      {
350
351
352
        ostringstream strstrm;
        strstrm << "dynare:k_order_perturbation: Caught general exception: " << e.message();
        DYN_MEX_FUNC_ERR_MSG_TXT(strstrm.str().c_str());
sebastien's avatar
sebastien committed
353
      }
354
    plhs[0] = mxCreateDoubleScalar(0);
sebastien's avatar
sebastien committed
355
356
  } // end of mexFunction()
} // end of extern C
357
#endif // ifdef MATLAB_MEX_FILE  to exclude mexFunction for other applications