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
|
// Copyright (C) 2005-2011, Ondra Kamenik
// 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
// ... order+1 matrices of derivatives
// output:
// res simulated results
#include "dynmex.h"
#include "mex.h"
#include "decision_rule.h"
#include "fs_tensor.h"
#include "SylvException.h"
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 least 12 input parameters and exactly 2 output arguments.\n");
int order = (int)mxGetScalar(prhs[0]);
if (nhrs != 12 + order)
DYN_MEX_FUNC_ERR_MSG_TXT("dynare_simul_ must have exactly 11+order input parameters.\n");
int nstat = (int)mxGetScalar(prhs[1]);
int npred = (int)mxGetScalar(prhs[2]);
int nboth = (int)mxGetScalar(prhs[3]);
int nforw = (int)mxGetScalar(prhs[4]);
int nexog = (int)mxGetScalar(prhs[5]);
const mxArray* const ystart = prhs[6];
const mxArray* const shocks = prhs[7];
const mxArray* const vcov = prhs[8];
int seed = (int)mxGetScalar(prhs[9]);
const mxArray* const ysteady = prhs[10];
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 != (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 != (int) shocks_dim[0])
DYN_MEX_FUNC_ERR_MSG_TXT("shocks has a wrong number of rows.\n");
if (nexog != (int) vcov_dim[0])
DYN_MEX_FUNC_ERR_MSG_TXT("vcov has a wrong number of rows.\n");
if (nexog != (int) vcov_dim[1])
DYN_MEX_FUNC_ERR_MSG_TXT("vcov has a wrong number of cols.\n");
if (ny != (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");
mxArray* res = mxCreateDoubleMatrix(ny, nper, mxREAL);
try {
// initialize tensor library
tls.init(order, npred+nboth+nexog);
// form the polynomial
UTensorPolynomial pol(ny, npred+nboth+nexog);
for (int dim = 0; dim <= order; dim++) {
const mxArray* gk = prhs[11+dim];
const mwSize* const gk_dim = mxGetDimensions(gk);
FFSTensor ft(ny, npred+nboth+nexog, dim);
if (ft.ncols() != (int) gk_dim[1]) {
char buf[1000];
sprintf(buf, "Wrong number of columns for folded tensor: got %d but I want %d\n",
(int) gk_dim[1], ft.ncols());
DYN_MEX_FUNC_ERR_MSG_TXT(buf);
}
if (ft.nrows() != (int) gk_dim[0]) {
char buf[1000];
sprintf(buf, "Wrong number of rows for folded tensor: got %d but I want %d\n",
(int) gk_dim[0], ft.nrows());
DYN_MEX_FUNC_ERR_MSG_TXT(buf);
}
ft.zeros();
ConstTwoDMatrix gk_mat(ft.nrows(), ft.ncols(), mxGetPr(gk));
ft.add(1.0, gk_mat);
UFSTensor* ut = new UFSTensor(ft);
pol.insert(ut);
}
// form the decision rule
UnfoldDecisionRule
dr(pol, PartitionY(nstat, npred, nboth, nforw),
nexog, ConstVector(mxGetPr(ysteady), ny));
// form the shock realization
TwoDMatrix shocks_mat(nexog, nper, (const double*)mxGetPr(shocks));
TwoDMatrix vcov_mat(nexog, nexog, (const double*)mxGetPr(vcov));
GenShockRealization sr(vcov_mat, shocks_mat, seed);
// simulate and copy the results
Vector ystart_vec((const double*)mxGetPr(ystart), ny);
TwoDMatrix* res_mat =
dr.simulate(DecisionRule::horner, nper,
ystart_vec, sr);
TwoDMatrix res_tmp_mat(ny, nper, mxGetPr(res));
res_tmp_mat = (const TwoDMatrix&)(*res_mat);
delete res_mat;
plhs[1] = res;
} catch (const KordException& e) {
DYN_MEX_FUNC_ERR_MSG_TXT("Caugth Kord exception.");
} catch (const TLException& e) {
DYN_MEX_FUNC_ERR_MSG_TXT("Caugth TL exception.");
} catch (SylvException& e) {
DYN_MEX_FUNC_ERR_MSG_TXT("Caught Sylv exception.");
}
plhs[0] = mxCreateDoubleScalar(0);
}
};
|
