1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
/*
* Copyright © 2005-2011 Ondra Kamenik
* Copyright © 2019 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 <http://www.gnu.org/licenses/>.
*/
// This is the mexFunction providing interface to
// DecisionRule<>::simulate(). It takes the following input
// parameters:
// order the order of approximation, needs order+1 derivatives
// nstat
// npred
// nboth
// nforw
// nexog
// ystart starting value (full vector of endogenous)
// shocks matrix of shocks (nexog x number of period)
// vcov covariance matrix of shocks (nexog x nexog)
// seed integer seed
// ysteady full vector of decision rule's steady
// dr structure containing matrices of derivatives (g_0, g_1,…)
// output:
// res simulated results
#include "dynmex.h"
#include "mex.h"
#include "decision_rule.hh"
#include "fs_tensor.hh"
#include "SylvException.hh"
#include <string>
extern "C" {
void
mexFunction(int nlhs, mxArray *plhs[],
int nhrs, const mxArray *prhs[])
{
if (nhrs != 12 || nlhs != 2)
DYN_MEX_FUNC_ERR_MSG_TXT("dynare_simul_ must have at exactly 12 input parameters and 2 output arguments.\n");
int order = static_cast<int>(mxGetScalar(prhs[0]));
int nstat = static_cast<int>(mxGetScalar(prhs[1]));
int npred = static_cast<int>(mxGetScalar(prhs[2]));
int nboth = static_cast<int>(mxGetScalar(prhs[3]));
int nforw = static_cast<int>(mxGetScalar(prhs[4]));
int nexog = static_cast<int>(mxGetScalar(prhs[5]));
const mxArray *const ystart = prhs[6];
const mxArray *const shocks = prhs[7];
const mxArray *const vcov = prhs[8];
int seed = static_cast<int>(mxGetScalar(prhs[9]));
const mxArray *const ysteady = prhs[10];
const mxArray *const dr = prhs[11];
const mwSize *const ystart_dim = mxGetDimensions(ystart);
const mwSize *const shocks_dim = mxGetDimensions(shocks);
const mwSize *const vcov_dim = mxGetDimensions(vcov);
const mwSize *const ysteady_dim = mxGetDimensions(ysteady);
int ny = nstat + npred + nboth + nforw;
if (ny != static_cast<int>(ystart_dim[0]))
DYN_MEX_FUNC_ERR_MSG_TXT("ystart has wrong number of rows.\n");
if (1 != ystart_dim[1])
DYN_MEX_FUNC_ERR_MSG_TXT("ystart has wrong number of cols.\n");
int nper = shocks_dim[1];
if (nexog != static_cast<int>(shocks_dim[0]))
DYN_MEX_FUNC_ERR_MSG_TXT("shocks has a wrong number of rows.\n");
if (nexog != static_cast<int>(vcov_dim[0]))
DYN_MEX_FUNC_ERR_MSG_TXT("vcov has a wrong number of rows.\n");
if (nexog != static_cast<int>(vcov_dim[1]))
DYN_MEX_FUNC_ERR_MSG_TXT("vcov has a wrong number of cols.\n");
if (ny != static_cast<int>(ysteady_dim[0]))
DYN_MEX_FUNC_ERR_MSG_TXT("ysteady has wrong number of rows.\n");
if (1 != ysteady_dim[1])
DYN_MEX_FUNC_ERR_MSG_TXT("ysteady has wrong number of cols.\n");
plhs[1] = mxCreateDoubleMatrix(ny, nper, mxREAL);
try
{
// initialize tensor library
TLStatic::init(order, npred+nboth+nexog);
// form the polynomial
UTensorPolynomial pol(ny, npred+nboth+nexog);
for (int dim = 0; dim <= order; dim++)
{
const mxArray *gk_m = mxGetField(dr, 0, ("g_" + std::to_string(dim)).c_str());
if (!gk_m)
DYN_MEX_FUNC_ERR_MSG_TXT(("Can't find field `g_" + std::to_string(dim) + "' in structured passed as last argument").c_str());
ConstTwoDMatrix gk(gk_m);
FFSTensor ft(ny, npred+nboth+nexog, dim);
if (ft.ncols() != gk.ncols())
DYN_MEX_FUNC_ERR_MSG_TXT(("Wrong number of columns for folded tensor: got " + std::to_string(gk.ncols()) + " but i want " + std::to_string(ft.ncols()) + '\n').c_str());
if (ft.nrows() != gk.nrows())
DYN_MEX_FUNC_ERR_MSG_TXT(("Wrong number of rows for folded tensor: got " + std::to_string(gk.nrows()) + " but i want " + std::to_string(ft.nrows()) + '\n').c_str());
ft.zeros();
ft.add(1.0, gk);
pol.insert(std::make_unique<UFSTensor>(ft));
}
// form the decision rule
UnfoldDecisionRule dr(pol, PartitionY(nstat, npred, nboth, nforw),
nexog, ConstVector{ysteady});
// form the shock realization
ConstTwoDMatrix shocks_mat(nexog, nper, ConstVector{shocks});
ConstTwoDMatrix vcov_mat(nexog, nexog, ConstVector{vcov});
GenShockRealization sr(vcov_mat, shocks_mat, seed);
// simulate and copy the results
TwoDMatrix res_mat{dr.simulate(DecisionRule::emethod::horner, nper,
ConstVector{ystart}, sr)};
TwoDMatrix res_tmp_mat{plhs[1]};
res_tmp_mat = const_cast<const TwoDMatrix &>(res_mat);
}
catch (const KordException &e)
{
DYN_MEX_FUNC_ERR_MSG_TXT("Caught Kord exception.");
}
catch (const TLException &e)
{
DYN_MEX_FUNC_ERR_MSG_TXT("Caught TL exception.");
}
catch (SylvException &e)
{
DYN_MEX_FUNC_ERR_MSG_TXT("Caught Sylv exception.");
}
plhs[0] = mxCreateDoubleScalar(0);
}
};
|
