dynare/dynare++/kord/global_check.cc

// Copyright 2005, Ondra Kamenik

#include "SymSchurDecomp.hh"

#include "global_check.hh"
#include "seed_generator.hh"

#include "smolyak.hh"
#include "product.hh"
#include "quasi_mcarlo.hh"

#include <utility>
#include <cmath>

/* Here we just set a reference to the approximation, and create a new
   |DynamicModel|. */

ResidFunction::ResidFunction(const Approximation &app)
  : VectorFunction(app.getModel().nexog(), app.getModel().numeq()), approx(app),
    model(app.getModel().clone())
{
}

ResidFunction::ResidFunction(const ResidFunction &rf)
  : VectorFunction(rf), approx(rf.approx), model(rf.model->clone())
{
  if (rf.yplus)
    yplus = std::make_unique<Vector>(*(rf.yplus));
  if (rf.ystar)
    ystar = std::make_unique<Vector>(*(rf.ystar));
  if (rf.u)
    u = std::make_unique<Vector>(*(rf.u));
  if (rf.hss)
    hss = std::make_unique<FTensorPolynomial>(*(rf.hss));
}

/* This sets $y^*$ and $u$. We have to create |ystar|, |u|, |yplus| and
   |hss|. */

void
ResidFunction::setYU(const ConstVector &ys, const ConstVector &xx)
{
  ystar = std::make_unique<Vector>(ys);
  u = std::make_unique<Vector>(xx);
  yplus = std::make_unique<Vector>(model->numeq());
  approx.getFoldDecisionRule().evaluate(DecisionRule::emethod::horner,
                                        *yplus, *ystar, *u);

  // make a tensor polynomial of in-place subtensors from decision rule
  /* Note that the non-|const| polynomial will be used for a
     construction of |hss| and will be used in |const| context. So this
     const_cast is safe.

     Note, that there is always a folded decision rule in |Approximation|. */
  const FoldDecisionRule &dr = approx.getFoldDecisionRule();
  FTensorPolynomial dr_ss(model->nstat()+model->npred(), model->nboth()+model->nforw(),
                          const_cast<FoldDecisionRule &>(dr));

  // make |ytmp_star| be a difference of |yplus| from steady
  Vector ytmp_star(ConstVector(*yplus, model->nstat(), model->npred()+model->nboth()));
  ConstVector ysteady_star(dr.getSteady(), model->nstat(),
                           model->npred()+model->nboth());
  ytmp_star.add(-1.0, ysteady_star);

  // make |hss| and add steady to it
  /* Here is the |const| context of |dr_ss|. */
  hss = std::make_unique<FTensorPolynomial>(dr_ss, ytmp_star);
  ConstVector ysteady_ss(dr.getSteady(), model->nstat()+model->npred(),
                         model->nboth()+model->nforw());
  if (hss->check(Symmetry{0}))
    hss->get(Symmetry{0}).getData().add(1.0, ysteady_ss);
  else
    {
      auto ten = std::make_unique<FFSTensor>(hss->nrows(), hss->nvars(), 0);
      ten->getData() = ysteady_ss;
      hss->insert(std::move(ten));
    }
}

/* Here we evaluate the residual $F(y^*,u,u')$. We have to evaluate |hss|
   for $u'=$|point| and then we evaluate the system $f$. */

void
ResidFunction::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
{
  KORD_RAISE_IF(point.length() != hss->nvars(),
                "Wrong dimension of input vector in ResidFunction::eval");
  KORD_RAISE_IF(out.length() != model->numeq(),
                "Wrong dimension of output vector in ResidFunction::eval");
  Vector yss(hss->nrows());
  hss->evalHorner(yss, point);
  model->evaluateSystem(out, *ystar, *yplus, yss, *u);
}

/* This checks the $E[F(y^*,u,u')]$ for a given $y^*$ and $u$ by
   integrating with a given quadrature. Note that the input |ys| is $y^*$
   not whole $y$. */

void
GlobalChecker::check(const Quadrature &quad, int level,
                     const ConstVector &ys, const ConstVector &x, Vector &out)
{
  for (int ifunc = 0; ifunc < vfs.getNum(); ifunc++)
    dynamic_cast<GResidFunction &>(vfs.getFunc(ifunc)).setYU(ys, x);
  quad.integrate(vfs, level, out);
}

/* This method is a bulk version of |@<|GlobalChecker::check| vector
   code@>|. It decides between Smolyak and product quadrature according
   to |max_evals| constraint.

   Note that |y| can be either full (all endogenous variables including
   static and forward looking), or just $y^*$ (state variables). The
   method is able to recognize it. */

void
GlobalChecker::check(int max_evals, const ConstTwoDMatrix &y,
                     const ConstTwoDMatrix &x, TwoDMatrix &out)
{
  JournalRecordPair pa(journal);
  pa << "Checking approximation error for " << y.ncols()
     << " states with at most " << max_evals << " evaluations" << endrec;

  // decide about type of quadrature
  GaussHermite gh;

  SmolyakQuadrature dummy_sq(model.nexog(), 1, gh);
  int smol_evals;
  int smol_level;
  dummy_sq.designLevelForEvals(max_evals, smol_level, smol_evals);

  ProductQuadrature dummy_pq(model.nexog(), gh);
  int prod_evals;
  int prod_level;
  dummy_pq.designLevelForEvals(max_evals, prod_level, prod_evals);

  bool take_smolyak = (smol_evals < prod_evals) && (smol_level >= prod_level-1);

  std::unique_ptr<Quadrature> quad;
  int lev;

  // create the quadrature and report the decision
  if (take_smolyak)
    {
      quad = std::make_unique<SmolyakQuadrature>(model.nexog(), smol_level, gh);
      lev = smol_level;
      JournalRecord rec(journal);
      rec << "Selected Smolyak (level,evals)=(" << smol_level << ","
          << smol_evals << ") over product (" << prod_level << ","
          << prod_evals << ")" << endrec;
    }
  else
    {
      quad = std::make_unique<ProductQuadrature>(model.nexog(), gh);
      lev = prod_level;
      JournalRecord rec(journal);
      rec << "Selected product (level,evals)=(" << prod_level << ","
          << prod_evals << ") over Smolyak (" << smol_level << ","
          << smol_evals << ")" << endrec;
    }

  // check all column of |y| and |x|
  int first_row = (y.nrows() == model.numeq()) ? model.nstat() : 0;
  ConstTwoDMatrix ysmat(y, first_row, 0, model.npred()+model.nboth(), y.ncols());
  for (int j = 0; j < y.ncols(); j++)
    {
      ConstVector yj{ysmat.getCol(j)};
      ConstVector xj{x.getCol(j)};
      Vector outj{out.getCol(j)};
      check(*quad, lev, yj, xj, outj);
    }
}

/* This method checks an error of the approximation by evaluating
   residual $E[F(y^*,u,u')\vert y^*,u]$ for $y^*$ being the steady state, and
   changing $u$. We go through all elements of $u$ and vary them from
   $-mult\cdot\sigma$ to $mult\cdot\sigma$ in |m| steps. */

void
GlobalChecker::checkAlongShocksAndSave(mat_t *fd, const std::string &prefix,
                                       int m, double mult, int max_evals)
{
  JournalRecordPair pa(journal);
  pa << "Calculating errors along shocks +/- "
     << mult << " std errors, granularity " << m << endrec;

  // setup |y_mat| of steady states for checking
  TwoDMatrix y_mat(model.numeq(), 2*m*model.nexog()+1);
  for (int j = 0; j < 2*m*model.nexog()+1; j++)
    {
      Vector yj{y_mat.getCol(j)};
      yj = model.getSteady();
    }

  // setup |exo_mat| for checking
  TwoDMatrix exo_mat(model.nexog(), 2*m*model.nexog()+1);
  exo_mat.zeros();
  for (int ishock = 0; ishock < model.nexog(); ishock++)
    {
      double max_sigma = sqrt(model.getVcov().get(ishock, ishock));
      for (int j = 0; j < 2*m; j++)
        {
          int jmult = (j < m) ? j-m : j-m+1;
          exo_mat.get(ishock, 1+2*m*ishock+j) = mult*jmult*max_sigma/m;
        }
    }

  TwoDMatrix errors(model.numeq(), 2*m*model.nexog()+1);
  check(max_evals, y_mat, exo_mat, errors);

  // report errors along shock and save them
  TwoDMatrix res(model.nexog(), 2*m+1);
  JournalRecord rec(journal);
  rec << "Shock    value         error" << endrec;
  ConstVector err0{errors.getCol(0)};
  for (int ishock = 0; ishock < model.nexog(); ishock++)
    {
      TwoDMatrix err_out(model.numeq(), 2*m+1);
      for (int j = 0; j < 2*m+1; j++)
        {
          int jj;
          Vector error{err_out.getCol(j)};
          if (j != m)
            {
              if (j < m)
                jj = 1 + 2*m*ishock+j;
              else
                jj = 1 + 2*m*ishock+j-1;
              ConstVector coljj{errors.getCol(jj)};
              error = coljj;
            }
          else
            {
              jj = 0;
              error = err0;
            }
          JournalRecord rec1(journal);
          std::string shockname{model.getExogNames().getName(ishock)};
          shockname.resize(8, ' ');
          rec1 << shockname << ' ' << exo_mat.get(ishock, jj)
               << "\t" << error.getMax() << endrec;
        }
      err_out.writeMat(fd, prefix + "_shock_" + model.getExogNames().getName(ishock) + "_errors");
    }
}

/* This method checks errors on ellipse of endogenous states
   (predetermined variables). The ellipse is shaped according to
   covariance matrix of endogenous variables based on the first order
   approximation and scaled by |mult|. The points on the
   ellipse are chosen as polar images of the low discrepancy grid in a
   cube.

   The method works as follows: First we calculate symmetric Schur factor of
   covariance matrix of the states. Second we generate low discrepancy
   points on the unit sphere. Third we transform the sphere with the
   variance-covariance matrix factor and multiplier |mult| and initialize
   matrix of $u_t$ to zeros. Fourth we run the |check| method and save
   the results. */

void
GlobalChecker::checkOnEllipseAndSave(mat_t *fd, const std::string &prefix,
                                     int m, double mult, int max_evals)
{
  JournalRecordPair pa(journal);
  pa << "Calculating errors at " << m
     << " ellipse points scaled by " << mult << endrec;

  // make factor of covariance of variables
  /* Here we set |ycovfac| to the symmetric Schur decomposition factor of
     a submatrix of covariances of all endogenous variables. The submatrix
     corresponds to state variables (predetermined plus both). */
  TwoDMatrix ycov{approx.calcYCov()};
  TwoDMatrix ycovpred(const_cast<const TwoDMatrix &>(ycov), model.nstat(), model.nstat(),
                      model.npred()+model.nboth(), model.npred()+model.nboth());
  SymSchurDecomp ssd(ycovpred);
  ssd.correctDefinitness(1.e-05);
  TwoDMatrix ycovfac(ycovpred.nrows(), ycovpred.ncols());
  KORD_RAISE_IF(!ssd.isPositiveSemidefinite(),
                "Covariance matrix of the states not positive \
				  semidefinite in GlobalChecker::checkOnEllipseAndSave");
  ssd.getFactor(ycovfac);

  // put low discrepancy sphere points to |ymat|
  /* Here we first calculate dimension |d| of the sphere, which is a
     number of state variables minus one. We go through the |d|-dimensional
     cube $\langle 0,1\rangle^d$ by |QMCarloCubeQuadrature| and make a
     polar transformation to the sphere. The polar transformation $f^i$ can
     be written recursively wrt. the dimension $i$ as:
     $$\eqalign{
     f^0() &= \left[1\right]\cr
     f^i(x_1,\ldots,x_i) &=
     \left[\matrix{cos(2\pi x_i)\cdot f^{i-1}(x_1,\ldots,x_{i-1})\cr sin(2\pi x_i)}\right]
     }$$ */
  int d = model.npred()+model.nboth()-1;
  TwoDMatrix ymat(model.npred()+model.nboth(), (d == 0) ? 2 : m);
  if (d == 0)
    {
      ymat.get(0, 0) = 1;
      ymat.get(0, 1) = -1;
    }
  else
    {
      int icol = 0;
      ReversePerScheme ps;
      QMCarloCubeQuadrature qmc(d, m, ps);
      qmcpit beg = qmc.start(m);
      qmcpit end = qmc.end(m);
      for (qmcpit run = beg; run != end; ++run, icol++)
        {
          Vector ycol{ymat.getCol(icol)};
          Vector x(run.point());
          x.mult(2*M_PI);
          ycol[0] = 1;
          for (int i = 0; i < d; i++)
            {
              Vector subsphere(ycol, 0, i+1);
              subsphere.mult(cos(x[i]));
              ycol[i+1] = sin(x[i]);
            }
        }
    }

  // transform sphere |ymat| and prepare |umat| for checking
  /* Here we multiply the sphere points in |ymat| with the Cholesky
     factor to obtain the ellipse, scale the ellipse by the given |mult|,
     and initialize matrix of shocks |umat| to zero. */
  TwoDMatrix umat(model.nexog(), ymat.ncols());
  umat.zeros();
  ymat.mult(mult);
  ymat.multLeft(ycovfac);
  ConstVector ys(model.getSteady(), model.nstat(),
                 model.npred()+model.nboth());
  for (int icol = 0; icol < ymat.ncols(); icol++)
    {
      Vector ycol{ymat.getCol(icol)};
      ycol.add(1.0, ys);
    }

  // check on ellipse and save
  /* Here we check the points and save the results to MAT-4 file. */
  TwoDMatrix out(model.numeq(), ymat.ncols());
  check(max_evals, ymat, umat, out);

  ymat.writeMat(fd, prefix + "_ellipse_points");
  out.writeMat(fd, prefix + "_ellipse_errors");
}

/* Here we check the errors along a simulation. We simulate, then set
   |x| to zeros, check and save results. */

void
GlobalChecker::checkAlongSimulationAndSave(mat_t *fd, const std::string &prefix,
                                           int m, int max_evals)
{
  JournalRecordPair pa(journal);
  pa << "Calculating errors at " << m
     << " simulated points" << endrec;
  RandomShockRealization sr(model.getVcov(), seed_generator::get_new_seed());
  TwoDMatrix y{approx.getFoldDecisionRule().simulate(DecisionRule::emethod::horner,
                                                     m, model.getSteady(), sr)};
  TwoDMatrix x(model.nexog(), m);
  x.zeros();
  TwoDMatrix out(model.numeq(), m);
  check(max_evals, y, x, out);

  y.writeMat(fd, prefix + "_simul_points");
  out.writeMat(fd, prefix + "_simul_errors");
}