k_ord_dynare.cc 12.2 KB
Newer Older
george's avatar
george committed
1
/*
2
 * Copyright (C) 2008-2017 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

// GP, based on work by O.Kamenik

#include <vector>
23
#include "first_order.hh"
24
#include "dynamic_abstract_class.hh"
george's avatar
george committed
25

26
#include <cmath>
sebastien's avatar
sebastien committed
27
#include <sstream>
george's avatar
george committed
28

29 30 31
#include <iostream>
#include <fstream>

george's avatar
george committed
32 33 34 35
/**************************************************************************************/
/*       Dynare DynamicModel class                                                                 */
/**************************************************************************************/

36 37
KordpDynare::KordpDynare(const std::vector<std::string> &endo, int num_endo,
                         const std::vector<std::string> &exo, int nexog, int npar,
38 39
                         Vector &ysteady, TwoDMatrix &vcov, Vector &inParams, int nstat,
                         int npred, int nforw, int nboth, const int jcols, const Vector &nnzd,
40
                         const int nsteps, int norder,
41 42
                         Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg, double sstol,
                         const std::vector<int> &var_order, const TwoDMatrix &llincidence, double criterium) noexcept(false) :
43 44 45 46 47 48
  nStat{nstat}, nBoth{nboth}, nPred{npred}, nForw{nforw}, nExog{nexog}, nPar{npar},
  nYs{npred + nboth}, nYss{nboth + nforw}, nY{num_endo}, nJcols{jcols}, NNZD{nnzd}, nSteps{nsteps},
  nOrder{norder}, journal{jr}, ySteady{ysteady}, params{inParams}, vCov{vcov},
  md{1}, dnl{*this, endo}, denl{*this, exo}, dsnl{*this, dnl, denl}, ss_tol{sstol}, varOrder{var_order},
  ll_Incidence{llincidence}, qz_criterium{criterium}, g1p{nullptr},
  g2p{nullptr}, g3p{nullptr}, dynamicModelFile{move(dynamicModelFile_arg)}
49 50 51 52 53 54 55 56
{
  ReorderDynareJacobianIndices();

  //	Initialise ModelDerivativeContainer(*this, this->md, nOrder);
  for (int iord = 1; iord <= nOrder; iord++)
    md.insert(new FSSparseTensor(iord, nY+nYs+nYss+nExog, nY));
}

57 58
KordpDynare::KordpDynare(const std::vector<std::string> &endo, int num_endo,
                         const std::vector<std::string> &exo, int nexog, int npar,
59 60 61
                         Vector &ysteady, TwoDMatrix &vcov, Vector &inParams, int nstat,
                         int npred, int nforw, int nboth, const int jcols, const Vector &nnzd,
                         const int nsteps, int norder,
62 63
                         Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg, double sstol,
                         const std::vector<int> &var_order, const TwoDMatrix &llincidence, double criterium,
64 65 66 67 68 69 70
                         TwoDMatrix *g1_arg, TwoDMatrix *g2_arg, TwoDMatrix *g3_arg) noexcept(false) :
  nStat{nstat}, nBoth{nboth}, nPred{npred}, nForw{nforw}, nExog{nexog}, nPar{npar},
  nYs{npred + nboth}, nYss{nboth + nforw}, nY{num_endo}, nJcols{jcols}, NNZD{nnzd}, nSteps{nsteps},
  nOrder{norder}, journal{jr}, ySteady{ysteady}, params{inParams}, vCov{vcov},
  md{1}, dnl{*this, endo}, denl{*this, exo}, dsnl{*this, dnl, denl}, ss_tol{sstol}, varOrder{var_order},
  ll_Incidence{llincidence}, qz_criterium{criterium},
  g1p{g1_arg}, g2p{g2_arg}, g3p{g3_arg}, dynamicModelFile{move(dynamicModelFile_arg)}
george's avatar
george committed
71
{
72
  ReorderDynareJacobianIndices();
73

sebastien's avatar
sebastien committed
74 75
  //	Initialise ModelDerivativeContainer(*this, this->md, nOrder);
  for (int iord = 1; iord <= nOrder; iord++)
76
    md.insert(new FSSparseTensor(iord, nY+nYs+nYss+nExog, nY));
george's avatar
george committed
77 78 79 80
}

KordpDynare::~KordpDynare()
{
81
  // No need to manually delete tensors in "md", they are deleted by the TensorContainer destructor
george's avatar
george committed
82 83
}

84
void
sebastien's avatar
sebastien committed
85
KordpDynare::solveDeterministicSteady()
george's avatar
george committed
86 87 88 89 90
{
  JournalRecordPair pa(journal);
  pa << "Non-linear solver for deterministic steady state By-passed " << endrec;
}

91
void
92
KordpDynare::evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) noexcept(false)
george's avatar
george committed
93
{
94 95
  // This method is only called when checking the residuals at steady state (Approximation::check), so return zero residuals
  out.zeros();
george's avatar
george committed
96 97
}

98
void
99 100
KordpDynare::evaluateSystem(Vector &out, const ConstVector &yym, const ConstVector &yy,
                            const ConstVector &yyp, const Vector &xx) noexcept(false)
george's avatar
george committed
101
{
102 103
  // This method is only called when checking the residuals at steady state (Approximation::check), so return zero residuals
  out.zeros();
george's avatar
george committed
104
}
sebastien's avatar
sebastien committed
105

george's avatar
george committed
106
/************************************************
107 108 109 110
 * this is main derivative calculation functin that indirectly calls dynamic.dll
 * which performs actual calculation and reorders
 ***************************************************/
void
111
KordpDynare::calcDerivativesAtSteady()
george's avatar
george committed
112
{
113
  if (g1p == nullptr)
114
    {
115 116
      g1p = new TwoDMatrix(nY, nJcols);
      g1p->zeros();
117

118
      if (nOrder > 1)
119 120 121 122 123 124
        {
          // allocate space for sparse Hessian
          g2p = new TwoDMatrix((int) NNZD[1], 3);
          g2p->zeros();
        }

125
      if (nOrder > 2)
126 127 128 129
        {
          g3p = new TwoDMatrix((int) NNZD[2], 3);
          g3p->zeros();
        }
130 131 132

      Vector xx(nexog());
      xx.zeros();
133

134 135 136 137
      Vector out(nY);
      out.zeros();
      Vector llxSteady(nJcols-nExog);
      LLxSteady(ySteady, llxSteady);
138

139
      dynamicModelFile->eval(llxSteady, xx, params, ySteady, out, g1p, g2p, g3p);
140

141
    }
142

143 144
  populateDerivativesContainer(*g1p, 1, JacobianIndices);
  delete g1p;
145

146
  if (nOrder > 1)
147
    {
148 149
      populateDerivativesContainer(*g2p, 2, JacobianIndices);
      delete g2p;
150
    }
151
  if (nOrder > 2)
152
    {
153 154
      populateDerivativesContainer(*g3p, 3, JacobianIndices);
      delete g3p;
155
    }
george's avatar
george committed
156 157 158
}

/*******************************************************************************
sebastien's avatar
sebastien committed
159 160
 * populateDerivatives to sparse Tensor and fit it in the Derivatives Container
 *******************************************************************************/
161
void
162
KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const std::vector<int> &vOrder)
george's avatar
george committed
163
{
164
  // model derivatives FSSparseTensor instance
165
  FSSparseTensor *mdTi = (new FSSparseTensor(ord, nJcols, nY));
166

167
  IntSequence s(ord, 0);
168

169 170
  if (ord == 1)
    {
171
      for (int i = 0; i < g.ncols(); i++)
sebastien's avatar
sebastien committed
172
        {
173
          for (int j = 0; j < g.nrows(); j++)
sebastien's avatar
sebastien committed
174 175 176
            {
              double x;
              if (s[0] < nJcols-nExog)
177
                x = g.get(j, vOrder[s[0]]);
sebastien's avatar
sebastien committed
178
              else
179
                x = g.get(j, s[0]);
sebastien's avatar
sebastien committed
180 181 182 183 184
              if (x != 0.0)
                mdTi->insert(s, j, x);
            }
          s[0]++;
        }
185
    }
sebastien's avatar
sebastien committed
186
  else if (ord == 2)
187
    {
sebastien's avatar
sebastien committed
188
      int nJcols1 = nJcols-nExog;
189
      std::vector<int> revOrder(nJcols1);
sebastien's avatar
sebastien committed
190
      for (int i = 0; i < nJcols1; i++)
191 192
        revOrder[vOrder[i]] = i;
      for (int i = 0; i < g.nrows(); i++)
sebastien's avatar
sebastien committed
193
        {
194 195
          int j = (int) g.get(i, 0)-1; // hessian indices start with 1
          int i1 = (int) g.get(i, 1) -1;
196 197
          if (j < 0 || i1 < 0)
            continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix())
198 199
          int s0 = i1 / nJcols;
          int s1 = i1 % nJcols;
sebastien's avatar
sebastien committed
200 201 202 203 204 205 206 207 208 209
          if (s0 < nJcols1)
            s[0] = revOrder[s0];
          else
            s[0] = s0;
          if (s1 < nJcols1)
            s[1] = revOrder[s1];
          else
            s[1] = s1;
          if (s[1] >= s[0])
            {
210
              double x = g.get(i, 2);
sebastien's avatar
sebastien committed
211 212 213
              mdTi->insert(s, j, x);
            }
        }
214
    }
sebastien's avatar
sebastien committed
215
  else if (ord == 3)
216 217 218
    {
      int nJcols1 = nJcols-nExog;
      int nJcols2 = nJcols*nJcols;
219
      std::vector<int> revOrder(nJcols1);
220
      for (int i = 0; i < nJcols1; i++)
221 222
        revOrder[vOrder[i]] = i;
      for (int i = 0; i < g.nrows(); i++)
sebastien's avatar
sebastien committed
223
        {
224 225
          int j = (int) g.get(i, 0)-1;
          int i1 = (int) g.get(i, 1) -1;
226 227
          if (j < 0 || i1 < 0)
            continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix())
228 229 230 231
          int s0 = i1 / nJcols2;
          int i2 = i1 % nJcols2;
          int s1 = i2 / nJcols;
          int s2 = i2 % nJcols;
sebastien's avatar
sebastien committed
232 233 234 235 236 237 238 239 240 241 242 243
          if (s0 < nJcols1)
            s[0] = revOrder[s0];
          else
            s[0] = s0;
          if (s1 < nJcols1)
            s[1] = revOrder[s1];
          else
            s[1] = s1;
          if (s2 < nJcols1)
            s[2] = revOrder[s2];
          else
            s[2] = s2;
244 245 246
          double x = g.get(i, 2);
          if (s.isSorted())
            mdTi->insert(s, j, x);
sebastien's avatar
sebastien committed
247
        }
248
    }
sebastien's avatar
sebastien committed
249

george's avatar
george committed
250
  // md container
251
  md.remove(Symmetry{ord});
252
  md.insert(mdTi);
253
  // No need to delete mdTi, it will be deleted by TensorContainer destructor
george's avatar
george committed
254 255 256
}

/*********************************************************
257 258 259 260
 * LLxSteady()
 * returns ySteady extended with leads and lags suitable for
 * passing to <model>_dynamic DLL
 *************************************************************/
261
void
262
KordpDynare::LLxSteady(const Vector &yS, Vector &llxSteady) noexcept(false)
263 264
{
  if ((nJcols-nExog) == yS.length())
sebastien's avatar
sebastien committed
265 266
    throw DynareException(__FILE__, __LINE__, "ySteady already of right size");

267
  // create temporary square 2D matrix size nEndo x nEndo (sparse)
george's avatar
george committed
268
  // for the lag, current and lead blocks of the jacobian
269 270 271 272
  if (llxSteady.length() != nJcols-nExog)
    throw DynareException(__FILE__, __LINE__, "llxSteady has wrong size");

  for (int ll_row = 0; ll_row < ll_Incidence.nrows(); ll_row++)
sebastien's avatar
sebastien committed
273 274 275
    {
      // populate (non-sparse) vector with ysteady values
      for (int i = 0; i < nY; i++)
276
        {
277 278
          if (ll_Incidence.get(ll_row, i))
            llxSteady[((int) ll_Incidence.get(ll_row, i))-1] = yS[i];
279
        }
sebastien's avatar
sebastien committed
280
    }
george's avatar
george committed
281 282 283
}

/************************************
sebastien's avatar
sebastien committed
284 285 286 287
 * Reorder DynareJacobianIndices of variables in a vector according to
 * given int * varOrder together with lead & lag incidence matrix and
 * any the extra columns for exogenous vars, and then,
 * reorders its blocks given by the varOrder and the Dynare++ expectations:
george's avatar
george committed
288

sebastien's avatar
sebastien committed
289 290
 * extra	nboth+ npred (t-1) lags
 * varOrder
291
                static:
george's avatar
george committed
292
    pred
293 294
    both
    forward
sebastien's avatar
sebastien committed
295 296
    * extra both + nforw (t+1) leads, and
    * extra exogen
297

sebastien's avatar
sebastien committed
298
    * so to match the jacobian organisation expected by the Appoximation class
george's avatar
george committed
299 300 301 302 303 304 305 306 307
      both + nforw (t+1) leads
      static
      pred
      both
      forward
      nboth+ npred  (t-1) lags
      exogen
************************************/

308
void
309
KordpDynare::ReorderDynareJacobianIndices() noexcept(false)
310 311
{
  // create temporary square 2D matrix size nEndo x nEndo (sparse)
george's avatar
george committed
312
  // for the lag, current and lead blocks of the jacobian
313
  JacobianIndices.resize(nJcols);
314
  std::vector<int> tmp(nY);
sebastien's avatar
sebastien committed
315 316
  int i, j, rjoff = nJcols-nExog-1;

317
  for (int ll_row = 0; ll_row < ll_Incidence.nrows(); ll_row++)
sebastien's avatar
sebastien committed
318 319 320 321
    {
      // reorder in orde-var order & populate temporary nEndo (sparse) vector with
      // the lag, current and lead blocks of the jacobian respectively
      for (i = 0; i < nY; i++)
322
        tmp[i] = ((int) ll_Incidence.get(ll_row, varOrder[i]-1));
sebastien's avatar
sebastien committed
323 324 325 326 327
      // write the reordered blocks back to the jacobian
      // in reverse order
      for (j = nY-1; j >= 0; j--)
        if (tmp[j])
          {
328
            JacobianIndices[rjoff] = tmp[j] -1;
sebastien's avatar
sebastien committed
329 330 331 332 333
            rjoff--;
            if (rjoff < 0)
              break;
          }
    }
sebastien's avatar
sebastien committed
334

george's avatar
george committed
335
  //add the indices for the nExog exogenous jacobians
336
  for (j = nJcols-nExog; j < nJcols; j++)
337
    JacobianIndices[j] = j;
george's avatar
george committed
338 339 340 341 342 343
}

/**************************************************************************************/
/*       DynareNameList class                                                         */
/**************************************************************************************/

344
DynareNameList::DynareNameList(const KordpDynare &dynare, const std::vector<std::string> &names_arg) : names(names_arg)
george's avatar
george committed
345 346 347
{
}

348
DynareStateNameList::DynareStateNameList(const KordpDynare &dynare, const DynareNameList &dnl,
349
                                         const DynareNameList &denl)
george's avatar
george committed
350
{
351
  for (int i = 0; i < dynare.nys(); i++)
352
    names.push_back(std::string{dnl.getName(i+dynare.nstat())});
353
  for (int i = 0; i < dynare.nexog(); i++)
354
    names.push_back(std::string{denl.getName(i)});
george's avatar
george committed
355
}