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

// GP, based on work by O.Kamenik

#include <vector>
#include "first_order.h"
24 25
#include "k_ord_dynare.hh"
#include "dynamic_dll.hh"
george's avatar
george committed
26

27
#include <cmath>
sebastien's avatar
sebastien committed
28
#include <sstream>
george's avatar
george committed
29 30 31 32 33 34 35

#include "memory_file.h"

/**************************************************************************************/
/*       Dynare DynamicModel class                                                                 */
/**************************************************************************************/

36 37
KordpDynare::KordpDynare(const vector<string> &endo, int num_endo,
                         const vector<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,
sebastien's avatar
sebastien committed
41
                         Journal &jr, DynamicModelDLL &dynamicDLL, double sstol,
42
                         const vector<int> &var_order, const TwoDMatrix &llincidence, double criterium) throw (TLException) :
sebastien's avatar
sebastien committed
43 44
  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),
45
  nOrder(norder), journal(jr), ySteady(ysteady), params(inParams), vCov(vcov),
46
  md(1), dnl(*this, endo), denl(*this, exo), dsnl(*this, dnl, denl), ss_tol(sstol), varOrder(var_order),
47
  ll_Incidence(llincidence), qz_criterium(criterium), dynamicDLL(dynamicDLL)
george's avatar
george committed
48
{
49
  ReorderDynareJacobianIndices();
50

sebastien's avatar
sebastien committed
51 52
  //	Initialise ModelDerivativeContainer(*this, this->md, nOrder);
  for (int iord = 1; iord <= nOrder; iord++)
53
    md.insert(new FSSparseTensor(iord, nY+nYs+nYss+nExog, nY));
george's avatar
george committed
54 55 56 57
}

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

61
void
sebastien's avatar
sebastien committed
62
KordpDynare::solveDeterministicSteady()
george's avatar
george committed
63 64 65 66 67
{
  JournalRecordPair pa(journal);
  pa << "Non-linear solver for deterministic steady state By-passed " << endrec;
}

68
void
sebastien's avatar
sebastien committed
69
KordpDynare::evaluateSystem(Vector &out, const Vector &yy, const Vector &xx) throw (DynareException)
george's avatar
george committed
70
{
71 72
  // 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
73 74
}

75 76
void
KordpDynare::evaluateSystem(Vector &out, const Vector &yym, const Vector &yy,
sebastien's avatar
sebastien committed
77
                            const Vector &yyp, const Vector &xx) throw (DynareException)
george's avatar
george committed
78
{
79 80
  // 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
81
}
sebastien's avatar
sebastien committed
82

george's avatar
george committed
83
/************************************************
84 85 86 87
 * this is main derivative calculation functin that indirectly calls dynamic.dll
 * which performs actual calculation and reorders
 ***************************************************/
void
88
KordpDynare::calcDerivativesAtSteady() throw (DynareException)
george's avatar
george committed
89
{
90 91
  TwoDMatrix g1(nY, nJcols);
  g1.zeros();
92

93
  TwoDMatrix *g2p = NULL, *g3p = NULL;
94 95 96

  if (nOrder > 1)
    {
97 98 99
      // allocate space for sparse Hessian
      g2p = new TwoDMatrix((int) NNZD[1], 3);
      g2p->zeros();
100
    }
101

102 103
  if (nOrder > 2)
    {
104 105
      g3p = new TwoDMatrix((int) NNZD[2], 3);
      g3p->zeros();
106
    }
107 108 109 110

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

sebastien's avatar
sebastien committed
111
  Vector out(nY);
george's avatar
george committed
112
  out.zeros();
113 114
  Vector llxSteady(nJcols-nExog);
  LLxSteady(ySteady, llxSteady);
sebastien's avatar
sebastien committed
115

116
  dynamicDLL.eval(llxSteady, xx, params, ySteady, out, &g1, g2p, g3p);
117

118
  populateDerivativesContainer(g1, 1, JacobianIndices);
119

120
  if (nOrder > 1)
george's avatar
george committed
121
    {
122 123
      populateDerivativesContainer(*g2p, 2, JacobianIndices);
      delete g2p;
george's avatar
george committed
124
    }
125
  if (nOrder > 2)
george's avatar
george committed
126
    {
127 128
      populateDerivativesContainer(*g3p, 3, JacobianIndices);
      delete g3p;
george's avatar
george committed
129
    }
george's avatar
george committed
130 131 132
}

/*******************************************************************************
sebastien's avatar
sebastien committed
133 134
 * populateDerivatives to sparse Tensor and fit it in the Derivatives Container
 *******************************************************************************/
135
void
136
KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const vector<int> &vOrder)
george's avatar
george committed
137
{
138
  // model derivatives FSSparseTensor instance
139
  FSSparseTensor *mdTi = (new FSSparseTensor(ord, nJcols, nY));
140 141 142

  IntSequence s(ord, 0);

143 144
  if (ord == 1)
    {
145
      for (int i = 0; i < g.ncols(); i++)
sebastien's avatar
sebastien committed
146
        {
147
          for (int j = 0; j < g.nrows(); j++)
sebastien's avatar
sebastien committed
148 149 150
            {
              double x;
              if (s[0] < nJcols-nExog)
151
                x = g.get(j, vOrder[s[0]]);
sebastien's avatar
sebastien committed
152
              else
153
                x = g.get(j, s[0]);
sebastien's avatar
sebastien committed
154 155 156 157 158
              if (x != 0.0)
                mdTi->insert(s, j, x);
            }
          s[0]++;
        }
159
    }
sebastien's avatar
sebastien committed
160
  else if (ord == 2)
161
    {
sebastien's avatar
sebastien committed
162 163 164
      int nJcols1 = nJcols-nExog;
      vector<int> revOrder(nJcols1);
      for (int i = 0; i < nJcols1; i++)
165 166
        revOrder[vOrder[i]] = i;
      for (int i = 0; i < g.nrows(); i++)
sebastien's avatar
sebastien committed
167
        {
168 169
          int j = (int) g.get(i, 0)-1; // hessian indices start with 1
          int i1 = (int) g.get(i, 1) -1;
sebastien's avatar
sebastien committed
170 171 172 173 174 175 176 177 178 179 180 181
          int s0 = (int) floor(((double) i1)/((double) nJcols));
          int s1 = i1- (nJcols*s0);
          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])
            {
182
              double x = g.get(i, 2);
sebastien's avatar
sebastien committed
183 184 185
              mdTi->insert(s, j, x);
            }
        }
186
    }
sebastien's avatar
sebastien committed
187
  else if (ord == 3)
188 189 190 191 192
    {
      int nJcols1 = nJcols-nExog;
      int nJcols2 = nJcols*nJcols;
      vector<int> revOrder(nJcols1);
      for (int i = 0; i < nJcols1; i++)
193 194
        revOrder[vOrder[i]] = i;
      for (int i = 0; i < g.nrows(); i++)
sebastien's avatar
sebastien committed
195
        {
196 197
          int j = (int) g.get(i, 0)-1;
          int i1 = (int) g.get(i, 1) -1;
sebastien's avatar
sebastien committed
198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213
          int s0 = (int) floor(((double) i1)/((double) nJcols2));
          int i2 = i1 - nJcols2*s0;
          int s1 = (int) floor(((double) i2)/((double) nJcols));
          int s2 = i2 - nJcols*s1;
          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;
214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
	  // reordering symmetry indices in increasing order
	  if (s[0] > s[1])
	    {
	      int temp = s[0];
	      s[0] = s[1];
	      s[1] = temp;
	    }
	  if (s[1] > s[2])
	    {
	      int temp = s[1];
	      s[1] = s[2];
	      s[2] = temp;
	    }
	  if (s[0] > s[1])
	    {
	      int temp = s[0];
	      s[0] = s[1];
	      s[1] = temp;
	    }
	  double x = g.get(i, 2);
	  mdTi->insert(s, j, x);
sebastien's avatar
sebastien committed
235
        }
236
    }
sebastien's avatar
sebastien committed
237

george's avatar
george committed
238 239
  // md container
  md.remove(Symmetry(ord));
240
  md.insert(mdTi);
241
  // No need to delete mdTi, it will be deleted by TensorContainer destructor
george's avatar
george committed
242 243 244
}

/*********************************************************
245 246 247 248
 * LLxSteady()
 * returns ySteady extended with leads and lags suitable for
 * passing to <model>_dynamic DLL
 *************************************************************/
249 250
void
KordpDynare::LLxSteady(const Vector &yS, Vector &llxSteady) throw (DynareException, TLException)
251 252
{
  if ((nJcols-nExog) == yS.length())
sebastien's avatar
sebastien committed
253 254
    throw DynareException(__FILE__, __LINE__, "ySteady already of right size");

255
  // create temporary square 2D matrix size nEndo x nEndo (sparse)
george's avatar
george committed
256
  // for the lag, current and lead blocks of the jacobian
257 258 259 260
  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
261 262 263
    {
      // populate (non-sparse) vector with ysteady values
      for (int i = 0; i < nY; i++)
264
        {
265 266
          if (ll_Incidence.get(ll_row, i))
            llxSteady[((int) ll_Incidence.get(ll_row, i))-1] = yS[i];
267
        }
sebastien's avatar
sebastien committed
268
    }
george's avatar
george committed
269 270 271
}

/************************************
sebastien's avatar
sebastien committed
272 273 274 275
 * 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
276

sebastien's avatar
sebastien committed
277 278
 * extra	nboth+ npred (t-1) lags
 * varOrder
279
                static:
george's avatar
george committed
280
    pred
281 282
    both
    forward
sebastien's avatar
sebastien committed
283 284
    * extra both + nforw (t+1) leads, and
    * extra exogen
285

sebastien's avatar
sebastien committed
286
    * so to match the jacobian organisation expected by the Appoximation class
george's avatar
george committed
287 288 289 290 291 292 293 294 295
      both + nforw (t+1) leads
      static
      pred
      both
      forward
      nboth+ npred  (t-1) lags
      exogen
************************************/

296 297
void
KordpDynare::ReorderDynareJacobianIndices() throw (TLException)
298 299
{
  // create temporary square 2D matrix size nEndo x nEndo (sparse)
george's avatar
george committed
300
  // for the lag, current and lead blocks of the jacobian
301
  JacobianIndices.resize(nJcols);
302
  vector <int> tmp(nY);
sebastien's avatar
sebastien committed
303 304
  int i, j, rjoff = nJcols-nExog-1;

305
  for (int ll_row = 0; ll_row < ll_Incidence.nrows(); ll_row++)
sebastien's avatar
sebastien committed
306 307 308 309
    {
      // 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++)
310
        tmp[i] = ((int) ll_Incidence.get(ll_row, varOrder[i]-1));
sebastien's avatar
sebastien committed
311 312 313 314 315
      // write the reordered blocks back to the jacobian
      // in reverse order
      for (j = nY-1; j >= 0; j--)
        if (tmp[j])
          {
316
            JacobianIndices[rjoff] = tmp[j] -1;
sebastien's avatar
sebastien committed
317 318 319 320 321
            rjoff--;
            if (rjoff < 0)
              break;
          }
    }
sebastien's avatar
sebastien committed
322

george's avatar
george committed
323
  //add the indices for the nExog exogenous jacobians
324
  for (j = nJcols-nExog; j < nJcols; j++)
325
    JacobianIndices[j] = j;
george's avatar
george committed
326 327 328 329 330 331
}

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

332
DynareNameList::DynareNameList(const KordpDynare &dynare, const vector<string> &names_arg) : names(names_arg)
george's avatar
george committed
333 334 335
{
}

336
DynareStateNameList::DynareStateNameList(const KordpDynare &dynare, const DynareNameList &dnl,
337
                                         const DynareNameList &denl)
george's avatar
george committed
338
{
339
  for (int i = 0; i < dynare.nys(); i++)
340
    names.push_back(string(dnl.getName(i+dynare.nstat())));
341
  for (int i = 0; i < dynare.nexog(); i++)
342
    names.push_back(string(denl.getName(i)));
george's avatar
george committed
343
}