dynare_simul.cc 5.08 KB
 Sébastien Villemot committed Jan 08, 2019 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 ``````// 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.hh" #include "fs_tensor.hh" #include "SylvException.hh" 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); } };``````