dynare/dynare++/kord/global_check.cc

// Copyright 2005, Ondra Kamenik

#include "SymSchurDecomp.hh"

#include "global_check.hh"

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

#ifdef __MINGW32__
# define __CROSS_COMPILATION__
#endif

#ifdef __MINGW64__
# define __CROSS_COMPILATION__
#endif

#ifdef __CROSS_COMPILATION__
# define M_PI 3.14159265358979323846
#endif

/* 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()),
    yplus(NULL), ystar(NULL), u(NULL), hss(NULL)
{
}

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

ResidFunction::~ResidFunction()
{
  delete model;
  // delete |y| and |u| dependent data
  if (yplus)
    delete yplus;
  if (ystar)
    delete ystar;
  if (u)
    delete u;
  if (hss)
    delete hss;
}

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

void
ResidFunction::setYU(const Vector &ys, const Vector &xx)
{
  // delete |y| and |u| dependent data
  /* NB: code shared with the destructor */
  if (yplus)
    delete yplus;
  if (ystar)
    delete ystar;
  if (u)
    delete u;
  if (hss)
    delete hss;

  ystar = new Vector(ys);
  u = new Vector(xx);
  yplus = new Vector(model->numeq());
  approx.getFoldDecisionRule().evaluate(DecisionRule::horner,
                                        *yplus, *ystar, *u);

  // make a tensor polynomial of in-place subtensors from decision rule
  /* Here we use a dirty tricky of converting |const| to non-|const| to
     obtain a polynomial of subtensor corresponding to non-predetermined
     variables. However, this new non-|const| polynomial will be used for a
     construction of |hss| and will be used in |const| context. So this
     dirty thing is safe.

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

  // 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.c->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 = new FTensorPolynomial(dr_ss, ytmp_star);
  ConstVector ysteady_ss(dr.c->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
    {
      FFSTensor *ten = new FFSTensor(hss->nrows(), hss->nvars(), 0);
      ten->getData() = ysteady_ss;
      hss->insert(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++)
    ((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);

  Quadrature *quad;
  int lev;

  // create the quadrature and report the decision
  if (take_smolyak)
    {
      quad = new 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 = new 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, j);
      ConstVector xj(x, j);
      Vector outj(out, j);
      check(*quad, lev, yj, xj, outj);
    }

  delete quad;
}

/* 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 char *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, j);
      yj = (const Vector &) 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, 0);
  char shock[9];
  char erbuf[17];
  for (int ishock = 0; ishock < model.nexog(); ishock++)
    {
      TwoDMatrix err_out(model.numeq(), 2*m+1);
      sprintf(shock, "%-8s", model.getExogNames().getName(ishock));
      for (int j = 0; j < 2*m+1; j++)
        {
          int jj;
          Vector error(err_out, j);
          if (j != m)
            {
              if (j < m)
                jj = 1 + 2*m*ishock+j;
              else
                jj = 1 + 2*m*ishock+j-1;
              ConstVector coljj(errors, jj);
              error = coljj;
            }
          else
            {
              jj = 0;
              error = err0;
            }
          JournalRecord rec1(journal);
          sprintf(erbuf, "%12.7g    ", error.getMax());
          rec1 << shock << " " << exo_mat.get(ishock, jj)
               << "\t" << erbuf << endrec;
        }
      char tmp[100];
      sprintf(tmp, "%s_shock_%s_errors", prefix, model.getExogNames().getName(ishock));
      err_out.writeMat(fd, tmp);
    }
}

/* 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 char *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 TwoDMatrix &)*ycov, model.nstat(), model.nstat(),
                      model.npred()+model.nboth(), model.npred()+model.nboth());
  delete ycov;
  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, 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, 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);

  char tmp[100];
  sprintf(tmp, "%s_ellipse_points", prefix);
  ymat.writeMat(fd, tmp);
  sprintf(tmp, "%s_ellipse_errors", prefix);
  out.writeMat(fd, tmp);
}

/* 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 char *prefix,
                                           int m, int max_evals)
{
  JournalRecordPair pa(journal);
  pa << "Calculating errors at " << m
     << " simulated points" << endrec;
  RandomShockRealization sr(model.getVcov(), system_random_generator.int_uniform());
  TwoDMatrix *y = approx.getFoldDecisionRule().simulate(DecisionRule::horner,
                                                        m, model.getSteady(), sr);
  TwoDMatrix x(model.nexog(), m);
  x.zeros();
  TwoDMatrix out(model.numeq(), m);
  check(max_evals, *y, x, out);

  char tmp[100];
  sprintf(tmp, "%s_simul_points", prefix);
  y->writeMat(fd, tmp);
  sprintf(tmp, "%s_simul_errors", prefix);
  out.writeMat(fd, tmp);

  delete y;
}