/*
 * Common functions for "Nash or Ramsey Inflation" by Sargent, Williams, and Zha (SWZ).
 *
 * Copywright (c) by Sargent, Williams, and Zha, June 2004.
 * Addition, 03/17/05 (SWZ Notes: pp.42-43);
*/
#include "swz_comfuns.h"
#include "numgradgen.h"   //Only used by CreateNumericalHessian().
#define N1STSET  4     //Number of parameters for the 1st set of parameters (vstar, theta0, theta1, tau1).
#define N2NDSET  1     //Number of parameters for the 2nd set of parameters (zeta1).
#define N3RDSET  1     //Number of parameters for the 3rd set of parameters (zeta2).
#define AFORMAT  " %.16e "      //A: accurate (most accurate).


static void OneDrawOf1stSet(struct TSconstgibbs_tag *constgibbs_ps, struct TSgovprob_tag *govprob_ps, struct TSconstmodpars_tag *constmodpars_ps);
static void OneDrawOf2nd3rdSet(struct TSconstgibbs_tag *constgibbs_ps, struct TSgovprob_tag *govprob_ps, struct TSconstmodpars_tag *constmodpars_ps);
static void OneDrawOf4a5a6thSets(struct TSconstgibbs_tag *constgibbs_ps, struct TSgovprob_tag *govprob_ps, struct TSconstmodpars_tag *constmodpars_ps);
static void ftd_UpdatelogMLH_mhm(struct TSconstgibbs_tag *constgibbs_ps);
//--- Added 03/17/05.
static void ftd_computeqvals(struct TSconstmodpars_tag *constmodpars_ps, struct TSgovprob_tag *govprob_ps, struct TSconstgibbs_tag *constgibbs_ps);
//
static struct TSforecasts_tag *CreateTSforecasts(int kx, int tv);
static struct TSforecasts_tag *DestroyTSforecasts(struct TSforecasts_tag *forecasts_ps);
static void InitializeKalman1datapoint(TSkalcvfurw *kalcvfurw1datapoint_ps, int locti, struct TSgovprob_tag *govprob_ps, void *constmod_ps, int pointflag);
//
static void ftd_histresids(void *constmod_ps, int pointflag, struct TSgovprob_tag *govprob_ps);

//************************************************************
// I. For both estimation and posterior draws.
//************************************************************
TSinput *CreateTSinput(FILE *fptr_input, FILE *fptr_output, const char *filename_input, const char *filenametag)
{
   //Outputs:
   //  input_ps: Input arguments.
   //Inputs:
   //  fptr_input:  Pointer to first input file.
   //  filenametag:  File name tag read from the command line, such as bmpce.

   int nData;  //Number of months or quarters for the raw data imported.
   int nrows, ncols;
   int nvec;
   //===
   TSinput *input_ps=tzMalloc(1, TSinput);


   //=== Memory initialization for pointers.
   input_ps->Dataraw_dm = (TSdmatrix *)NULL;

   //=== Points to a File.
   input_ps->fptr_input = fptr_input;
   input_ps->fptr_output = fptr_output;

   //--------- Setting up the model from the definition in swz_comfuns.h. ---------

   //--------- Reads data from fptr_input.  Order of the following lines may matter. ---------
   if ( !fn_SetFilePosition(fptr_input, "//== indxWhichModel ==//") || fscanf(fptr_input, " %d ", &input_ps->indxWhichModel) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== indxWhichModel ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== indxFixz10P10 ==//") || fscanf(fptr_input, " %d ", &input_ps->indxFixz10P10) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== indxFixz10P10 ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== tc_P10 ==//") || fscanf(fptr_input, " %lf ", &input_ps->tc_P10) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== tc_P10 ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== indxEstFinder ==//") || fscanf(fptr_input, " %d ", &input_ps->indxEstFinder) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== indxEstFinder ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== indxFcsts ==//") || fscanf(fptr_input, " %d ", &input_ps->indxFcsts) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== indxFcsts ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== locti ==//") || fscanf(fptr_input, " %d ", &input_ps->locti) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== locti ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== nsteps ==//") || fscanf(fptr_input, " %d ", &input_ps->nsteps) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== nsteps ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== indxCnterFact ==//") || fscanf(fptr_input, " %d ", &input_ps->indxCnterFact) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== indxCnterFact ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== indxNoUpdating ==//") || fscanf(fptr_input, " %d ", &input_ps->indxNoUpdating) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== indxNoUpdating ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if (fn_SetFilePosition(fptr_input, "//== sclhistshocks_dv ==//")) {
      if ( fscanf(fptr_input, " %d ", &nvec) != 1) {
         printf("\n\nCheck the first integer in the first row below the line //== sclhistshocks_dv ==// in the input data file %s!\n", filename_input);
         exit(EXIT_FAILURE);
      }
      input_ps->sclhistshocks_dv = CreateVector_lf(nvec);
      if ( !ReadVector_lf(fptr_input, input_ps->sclhistshocks_dv) )  {
         printf("\n\nCheck the data matrix or vector after the first row below the line //== sclhistshocks_dv ==// in the input data file %s!\n", filename_input);
         exit(EXIT_FAILURE);
      }
      input_ps->sclhistshocks_dv->flag = V_DEF;
   }
   if ( !fn_SetFilePosition(fptr_input, "//== store_sims ==//") || fscanf(fptr_input, " %d ", &input_ps->store_sims) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== store_sims ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== nsims ==//") || fscanf(fptr_input, " %d ", &input_ps->nsims) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== nsims ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== scl_epsilon ==//") || fscanf(fptr_input, " %lf ", &input_ps->scl_epsilon) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== scl_epsilon ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== randomseed ==//") || fscanf(fptr_input, " %d ", &input_ps->randomseed) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== randomseed ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== qm ==//") || fscanf(fptr_input, " %d ", &input_ps->qm) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== qm ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== yrbin ==//") || fscanf(fptr_input, " %d ", &input_ps->yrbin) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== yrbin ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== qmbin ==//") || fscanf(fptr_input, " %d ", &input_ps->qmbin) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== qmbin ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== yrfin ==//") || fscanf(fptr_input, " %d ", &input_ps->yrfin) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== yrfin ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== qmfin ==//") || fscanf(fptr_input, " %d ", &input_ps->qmfin) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== qmfin ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== yrstart ==//") || fscanf(fptr_input, " %d ", &input_ps->yrstart) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== yrstart ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== qmstart ==//") || fscanf(fptr_input, " %d ", &input_ps->qmstart) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== qmstart ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
//   if ( !fn_SetFilePosition(fptr_input, "//== kx ==//") || fscanf(fptr_input, " %d ", &input_ps->kx) != 1 ) {
//      printf("\n\nFatal Error: check the below the line //== kx ==// in the input data file %s!\n", filename_input);
//      exit(EXIT_FAILURE);
//   }
//   if ( !fn_SetFilePosition(fptr_input, "//== nlagstrue ==//") || fscanf(fptr_input, " %d ", &input_ps->nlagstrue) != 1 ) {
//      printf("\n\nFatal Error: check the below the line //== nlagstrue ==// in the input data file %s!\n", filename_input);
//      exit(EXIT_FAILURE);
//   }
//   if ( !fn_SetFilePosition(fptr_input, "//== nlagsinflation ==//") || fscanf(fptr_input, " %d ", &input_ps->nlagsinflation) != 1 ) {
//      printf("\n\nFatal Error: check the below the line //== nlagsinflation ==// in the input data file %s!\n", filename_input);
//      exit(EXIT_FAILURE);
//   }
//   if ( !fn_SetFilePosition(fptr_input, "//== govmodel ==//") || fscanf(fptr_input, " %d ", &input_ps->govmodel) != 1 ) {
//      printf("\n\nFatal Error: check the below the line //== govmodel ==// in the input data file %s!\n", filename_input);
//      exit(EXIT_FAILURE);
//   }
//   if ( !fn_SetFilePosition(fptr_input, "//== nlagsgov ==//") || fscanf(fptr_input, " %d ", &input_ps->nlagsgov) != 1 ) {
//      printf("\n\nFatal Error: check the below the line //== nlagsgov ==// in the input data file %s!\n", filename_input);
//      exit(EXIT_FAILURE);
//   }
   if ( !fn_SetFilePosition(fptr_input, "//== govpitarget ==//") || fscanf(fptr_input, " %lf ", &input_ps->govpitarget) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== govpitarget ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== govutarget ==//") || fscanf(fptr_input, " %lf ", &input_ps->govutarget) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== govutarget ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== govbeta ==//") || fscanf(fptr_input, " %lf ", &input_ps->govbeta) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== govbeta ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== govlambda1 ==//") || fscanf(fptr_input, " %lf ", &input_ps->govlambda1) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== govlambda1 ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== govlambda2 ==//") || fscanf(fptr_input, " %lf ", &input_ps->govlambda2) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== govlambda2 ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== _div_gov ==//") || fscanf(fptr_input, " %lf ", &input_ps->_div_gov) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== _div_gov ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== tc1 ==//") || fscanf(fptr_input, " %lf ", &input_ps->tc1) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== tc1 ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== tc2 ==//") || fscanf(fptr_input, " %lf ", &input_ps->tc2) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== tc2 ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== tc3 ==//") || fscanf(fptr_input, " %lf ", &input_ps->tc3) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== tc3 ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== tc4 ==//") || fscanf(fptr_input, " %lf ", &input_ps->tc4) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== tc4 ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== tc5 ==//") || fscanf(fptr_input, " %lf ", &input_ps->tc5) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== tc5 ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   //+ Posterior draws.
   if ( !fn_SetFilePosition(fptr_input, "//== ndraws1 ==//") || fscanf(fptr_input, " %d ", &input_ps->ndraws1) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== ndraws1 ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== ndrawsq ==//") || fscanf(fptr_input, " %d ", &input_ps->ndrawsq) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== ndrawsq ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== napartsave ==//") || fscanf(fptr_input, " %d ", &input_ps->napartsave) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== napartsave ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== timesk2b ==//") || fscanf(fptr_input, " %d ", &input_ps->timesk2b) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== timesk2b ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== sf_mtpl ==//") || fscanf(fptr_input, " %lf ", &input_ps->sf_mtpl) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== sf_mtpl ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== sf_mhm ==//") || fscanf(fptr_input, " %lf ", &input_ps->sf_mhm) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== sf_mhm ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== indx_readgibbs ==//") || fscanf(fptr_input, " %d ", &input_ps->indx_readgibbs) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== indx_readgibbs ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== indx_lrsacrats ==//") || fscanf(fptr_input, " %d ", &input_ps->indx_lrsacrats) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== indx_lrsacrats ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== indx_mhm ==//") || fscanf(fptr_input, " %d ", &input_ps->indx_mhm) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== indx_mhm ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== indxRecord ==//") || fscanf(fptr_input, " %d ", &input_ps->indxRecord) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== indxRecord ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== indxMHMPeak ==//") || fscanf(fptr_input, " %d ", &input_ps->indxMHMPeak) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== indxMHMPeak ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if ( !fn_SetFilePosition(fptr_input, "//== cutval_logpost ==//") || fscanf(fptr_input, " %lf ", &input_ps->cutval_logpost) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== cutval_logpost ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
//   else if ((input_ps->wlowlogpost>1.0) || (input_ps->wlowlogpost<0.0))
//   {
//      printf("\n\nFatal Error: number must be between 0.0 and 1.0 below the line //== wlowlogpost ==// in the input data file %s!\n", filename_input);
//      exit(EXIT_FAILURE);
//   }
   if (input_ps->indx_mhm && fn_SetFilePosition(fptr_input, "//== pcuts_dv ==//")) {
      if ( fscanf(fptr_input, " %d ", &nvec) != 1) {
         printf("\n\nCheck the first integer in the first row below the line //== pcuts_dv ==// in the input data file %s!\n", filename_input);
         exit(EXIT_FAILURE);
      }
      input_ps->pcuts_dv = CreateVector_lf(nvec);
      if ( !ReadVector_lf(fptr_input, input_ps->pcuts_dv) )  {
         printf("\n\nCheck the data matrix or vector after the first row below the line //== pcuts_dv ==// in the input data file %s!\n", filename_input);
         exit(EXIT_FAILURE);
      }
      input_ps->pcuts_dv->flag = V_DEF;
   }
   else  input_ps->pcuts_dv = NULL;
   //+ Added 03/17/05.
   if (input_ps->indx_mhm && fn_SetFilePosition(fptr_input, "//== zeta1cuts_dv ==//")) {
      if ( fscanf(fptr_input, " %d ", &nvec) != 1) {
         printf("\n\nCheck the first integer in the first row below the line //== zeta1cuts_dv ==// in the input data file %s!\n", filename_input);
         exit(EXIT_FAILURE);
      }
      input_ps->zeta1cuts_dv = CreateVector_lf(nvec);
      if ( !ReadVector_lf(fptr_input, input_ps->zeta1cuts_dv) )  {
         printf("\n\nCheck the data matrix or vector after the first row below the line //== zeta1cuts_dv ==// in the input data file %s!\n", filename_input);
         exit(EXIT_FAILURE);
      }
      input_ps->zeta1cuts_dv->flag = V_DEF;
   }
   else  input_ps->zeta1cuts_dv = NULL;
   if (input_ps->indx_mhm && fn_SetFilePosition(fptr_input, "//== zeta2cuts_dv ==//")) {
      if ( fscanf(fptr_input, " %d ", &nvec) != 1) {
         printf("\n\nCheck the first integer in the first row below the line //== zeta2cuts_dv ==// in the input data file %s!\n", filename_input);
         exit(EXIT_FAILURE);
      }
      input_ps->zeta2cuts_dv = CreateVector_lf(nvec);
      if ( !ReadVector_lf(fptr_input, input_ps->zeta2cuts_dv) )  {
         printf("\n\nCheck the data matrix or vector after the first row below the line //== zeta2cuts_dv ==// in the input data file %s!\n", filename_input);
         exit(EXIT_FAILURE);
      }
      input_ps->zeta2cuts_dv->flag = V_DEF;
   }
   else  input_ps->zeta2cuts_dv = NULL;
   //---
   if ( !fn_SetFilePosition(fptr_input, "//== locinflation ==//") || fscanf(fptr_input, " %d ", &input_ps->locinflation) != 1 ) {
      printf("\n\nFatal Error: check the below the line //== locinflation ==// in the input data file %s!\n", filename_input);
      exit(EXIT_FAILURE);
   }
   if (fn_SetFilePosition(fptr_input, "//== Dataraw_dm ==//")) {
      if ( fscanf(fptr_input, " %d ", &nrows) != 1) {
         printf("\n\nCheck the first integer in the first row below the line //== Dataraw_dm ==// in the input data file %s!\n", filename_input);
         exit(EXIT_FAILURE);
      }
      if ( fscanf(fptr_input, " %d ", &ncols) != 1) {
         printf("\n\nCheck the second integer in the first row below the line //== Dataraw_dm ==// in the input data file %s!\n", filename_input);
         exit(EXIT_FAILURE);
      }
      input_ps->Dataraw_dm = CreateMatrix_lf(nrows, ncols);
      if ( !ReadMatrix_lf(fptr_input, input_ps->Dataraw_dm) )  {
         printf("\n\nCheck the data matrix or vector after the first row below the line //== Dataraw_dm ==// in the input data file %s!\n", filename_input);
         exit(EXIT_FAILURE);
      }
      input_ps->Dataraw_dm->flag = M_GE;
   }
   else {
      printf("\n\n.../swz_comfuns.c/CreateTSinput(): the line with //== Dataraw_dm ==// in the input data file %s does not exist!\n", filename_input);
      exit(EXIT_FAILURE);
   }


   //--- Model specifications.
   if (input_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0 | NGO_CGP2TP1TI0)) {
      input_ps->kx = 6;   //Number of right-hand regressors in the government's Phillips curve.
                          //Classical right-hand variables:  {pi_t, pi_{t-1}, u_{t-1}, pi_{t-2}, u_{t-2}, 1}.
                          //Keynesian right-hand variables:  {u_t, u_{t-1}, pi_{t-1}, u_{t-2}, pi_{t-2}, 1}.
      input_ps->nlagsgov = 2;   //Number of lags or lag length for the government's misspecified Phillips curve.
      input_ps->nlagstrue = 1;  //Number of lags or lag length for the true Phillips curve.
      input_ps->nlagsinflation = 0;   //Number of lags or lag length for the true inflation process.

   }
   //+
   if (input_ps->nlagsinflation > input_ps->nlagsgov)  fn_DisplayError(".../swz_comfuns.c/CreateTSinput(): It must be nlagsinflation <= nlagsgov, i.e., lag length for the true inflation process must be less than that of the government regression");
   if (input_ps->nlagstrue > input_ps->nlagsgov)  fn_DisplayError(".../swz_comfuns.c/CreateTSinput(): It must be nlagstrue <= nlagsgov, i.e., lag length for the true Phillips curve must be less than that of the government regression");


//   //--------- Reads data from fptr_case.  Order of the following lines may matter. ---------
//   if ( !fn_SetFilePosition(fptr_case, "//== locinflation ==//") || fscanf(fptr_case, " %d ", &input_ps->locinflation) != 1 ) {
//      printf("\n\nFatal Error: check the below the line //== locinflation ==// in the input data file datainp_%s.prn!\n", filenametag);
//      exit(EXIT_FAILURE);
//   }


   //------- Checks the dates against the dimension of the data matrix. -------
   nData = fn_locofyearqm(input_ps->qm, input_ps->yrbin, input_ps->qmbin, input_ps->yrfin, input_ps->qmfin) + 1;
   if (nData != input_ps->Dataraw_dm->nrows) {
      printf("\n\nFatal Error in .../swz_comfuns.c/CreateTSinput(): [yrbin qmbin] and [yrfin qmfin] give the range (%d) that differs from the row number (%d) of Dataraw_dm in the input data file %s!\n", nData, input_ps->Dataraw_dm->nrows, filename_input);
      exit(EXIT_FAILURE);
   }

   return (input_ps);
}
//---
TSinput *DestroyTSinput(TSinput *input_ps)
{
   if (input_ps) {
      DestroyVector_lf(input_ps->sclhistshocks_dv);
      DestroyVector_lf(input_ps->pcuts_dv);
      DestroyVector_lf(input_ps->zeta1cuts_dv);  //Added 03/17/05.
      DestroyVector_lf(input_ps->zeta2cuts_dv);  //Added 03/17/05.
      DestroyMatrix_lf(input_ps->Dataraw_dm);

      //===
      free(input_ps);
      return ((TSinput *)NULL);
   }
   else  return (input_ps);
}


//-----------------------
// Some initializations of the outside-government parameters.
//-----------------------
TSconstmodpars *CreateTSconstmodpars(FILE *fptr_input, TSinput *input_ps, TFconstmodest *est_func, TFconstmodprob *prob_func)   //fptr_input is added 03/17/05.
{
   int worki, nvec;
   TSconstmodpars *constmodpars_ps = tzMalloc(1, TSconstmodpars);

   //------- Added, 03/17/05. -------
   //--- Reads integer.
   if ( !fn_SetFilePosition(fptr_input, "//== indxNoGovOpt ==//") || fscanf(fptr_input, " %d ", &constmodpars_ps->indxNoGovOpt) != 1 )
   {
      printf("\n\nFatal Error: check the below the line //== indxNoGovOpt ==// in the input data file!\n");
      exit(EXIT_FAILURE);
   }
   //--- Reading a double vector.
   if (fn_SetFilePosition(fptr_input, "//== ngo_x_dv ==//"))
   {
      if ( fscanf(fptr_input, " %d ", &nvec) != 1)
      {
         printf("\n\nCheck the first integer in the first row below the line //== ngo_x_dv ==// in the input data file!\n");
         exit(EXIT_FAILURE);
      }
      constmodpars_ps->ngo_x_dv = CreateVector_lf(nvec);
      if ( !ReadVector_lf(fptr_input, constmodpars_ps->ngo_x_dv) )
      {
         printf("\n\nCheck the data matrix or vector after the first row below the line //== ngo_x_dv ==// in the input data file!\n");
         exit(EXIT_FAILURE);
      }
      constmodpars_ps->ngo_x_dv->flag = V_DEF;
   }
   else
   {
      printf("\n\nCheck the data matrix or vector after the first row below the line //== ngo_x_dv ==// in the input data file!\n");
      exit(EXIT_FAILURE);
   }
   //--- Normalization, this number is arbitrary. See 42-43.
   constmodpars_ps->ngo_sigmasq = 1.0/3.5653810749685604e+001;  //Fixed at the estimate with the government's optimization problem.




//   constmodpars_ps->indxEstFinder = input_ps->indxEstFinder;
   constmodpars_ps->indxWhichModel = input_ps->indxWhichModel;
   constmodpars_ps->indxFixz10P10 = input_ps->indxFixz10P10;
   constmodpars_ps->indxFcsts = input_ps->indxFcsts;
   constmodpars_ps->indxCnterFact = input_ps->indxCnterFact;
   constmodpars_ps->indxNoUpdating = input_ps->indxNoUpdating;
   constmodpars_ps->sclhistshocks_dv = CreateVector_lf(input_ps->sclhistshocks_dv->n);
   CopyVector0(constmodpars_ps->sclhistshocks_dv, input_ps->sclhistshocks_dv);
   constmodpars_ps->tc_P10 = input_ps->tc_P10;
   constmodpars_ps->scl_epsilon = input_ps->scl_epsilon;
   constmodpars_ps->randomseed = input_ps->randomseed;

   //=== For the MLE or posterior estimate.
//   constmodpars_ps->govmodel = input_ps->govmodel;    //== govmodel ==//   //Government's model (0: Classical direction or regresssion; 1: Keynesian direction).
   constmodpars_ps->kx = input_ps->kx;  //== kx ==//   //Number of right-hand regressors in the government's Phillips curve.
   constmodpars_ps->nlagsinflation = input_ps->nlagsinflation;   //== nlagsinflation ==//   //Number of lags or lag length for the infla
   constmodpars_ps->nlagstrue = input_ps->nlagstrue;  //== nlagstrue ==//   //Number of lags or lag length for the true Phillips curve.
   worki = input_ps->nlagstrue > input_ps->nlagsinflation ? input_ps->nlagstrue : input_ps->nlagsinflation;
   constmodpars_ps->maxnlags = worki > input_ps->nlagsgov ? worki : input_ps->nlagsgov;  //max of nlagstrue, nlagsinflation, and nlagsgov.
   constmodpars_ps->fss = input_ps->Dataraw_dm->nrows - constmodpars_ps->maxnlags;  //Effective sample size (i.e., excluding max lag length in the whole model).
   if (!input_ps->nlagsinflation && constmodpars_ps->nlagstrue==1) {
      //No lag and thus no rho1 and rho2 in the true inflation process to be estimated.
      //=== 6+kx+(kx+1)*kx variables in the order of: //vstar, theta0, theta1, tau1, xi1 (sqrt(zeta1)), xi2 (sqrt(zeta2)), z10, P10, and V.
      if (constmodpars_ps->indxWhichModel & NGO_CGP2TP1TI0)  constmodpars_ps->totnpars = 6;   //Added 03/17/05.
      else if (constmodpars_ps->indxFixz10P10 == 1)  constmodpars_ps->totnpars = 6 + input_ps->kx*(input_ps->kx+1)/2;
      else if (constmodpars_ps->indxFixz10P10 == 2)  constmodpars_ps->totnpars = 6 + input_ps->kx*(input_ps->kx+1);
      else  constmodpars_ps->totnpars = 6 + input_ps->kx + input_ps->kx*(input_ps->kx+1);  //Total number of free parameters.  If rho1=rho2=0, equals to 6+kx+(kx+1)*kx.
   }
   else if (input_ps->nlagstrue != 1)  fn_DisplayError(".../CreateTSconstmodpars(): Have not got time to work out the case for nlagstrue != 1");
   else  fn_DisplayError(".../CreateTSconstmodpars(): Have not got time to code up the case for nlagsinflation > 0");
   constmodpars_ps->govpitarget = input_ps->govpitarget;    //Government's targeted inflation.

   //=== Memory allocation.
   constmodpars_ps->x_dv = CreateVector_lf(constmodpars_ps->totnpars);   //totnpars-by-1 vector of free parameters.
   constmodpars_ps->z10_dv = CreateVector_lf(input_ps->kx);
   constmodpars_ps->k_fcstinf_z10_dv = CreateVector_lf(input_ps->kx);
   constmodpars_ps->P10_dm = CreateMatrix_lf(input_ps->kx,input_ps->kx);
   constmodpars_ps->V_dm = CreateMatrix_lf(input_ps->kx,input_ps->kx);
   //
   constmodpars_ps->upred_dv = CreateVector_lf(constmodpars_ps->fss);   //T-by-1 vector of predicted unemployment rate.
   constmodpars_ps->StrResidstran_dm = CreateMatrix_lf(2, constmodpars_ps->fss);
   //
   constmodpars_ps->sampledates_dv = CreateDatesVector_lf(input_ps->qm, input_ps->yrstart, input_ps->qmstart, input_ps->yrfin, input_ps->qmfin);


   //=== Initialization.
   constmodpars_ps->ustar = 5.0;  //Natural rate of unemployment.
   constmodpars_ps->gtheta0 = -1.0;
   constmodpars_ps->gtheta1 = 0.0;
   constmodpars_ps->gtau1 = 0.0;
   constmodpars_ps->gzeta1 = 1.0;  //Inverse of residual variance in the true Phillips curve.
   //+
   if (!input_ps->nlagsinflation) {
      //No lag and thus no rho1 and rho2 in the true inflation process to be estimated.
      constmodpars_ps->grho1 = 0.0;
      constmodpars_ps->grho2 = 0.0;
   }
   else if (input_ps->nlagsinflation == 1) {
      //Lag 1 and thus estimating rho1 and rho2 in the true inflation process
      constmodpars_ps->grho1 = 0.0;  //0.98;
      constmodpars_ps->grho2 = 0.0;  //30;
   }
   else  fn_DisplayError(".../CreateTSconstmodpars(): So far the lag length, input_ps->nlagsinflation, for the true inflation process is only allowed to be 0 or 1");
   constmodpars_ps->gzeta2 = 1.0;
   //+
   constmodpars_ps->gzeta = constmodpars_ps->gzeta1;  //Normalization.
   //--- Wild guess.
   InitializeConstantVector_lf(constmodpars_ps->z10_dv, 1.0);     //The initial coefficients on the government's false updated Phillips-curve regression.
   InitializeDiagonalMatrix_lf(constmodpars_ps->P10_dm, 1.0);
   InitializeDiagonalMatrix_lf(constmodpars_ps->V_dm, 1.0);
   if (constmodpars_ps->indxWhichModel & (CGP2TP1TI0 | NGO_CGP2TP1TI0))    //NGO_CGP2TP1TI0 is added 03/17/04.
   {
      //--- Use the estimated coefficients from ...\DataArrange\govregress48_59.m (from 1948:01 to 1959:12).
      if (constmodpars_ps->indxFixz10P10 == 1) {
         constmodpars_ps->z10_dv->v[0] = -0.13236567947738;
         constmodpars_ps->z10_dv->v[1] =  0.14191629076828;
         constmodpars_ps->z10_dv->v[2] =  1.09275370936974;
         constmodpars_ps->z10_dv->v[3] = -0.02157077022996;
         constmodpars_ps->z10_dv->v[4] = -0.13378436709391;
         constmodpars_ps->z10_dv->v[5] =  0.21897028305873;

         constmodpars_ps->P10_dm->M[mos(0,0,6)] = 0.00267876931048*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(0,1,6)] = -0.00369045856092*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(0,2,6)] = 0.00013255563323*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(0,3,6)] = 0.00116570612453*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(0,4,6)] = -0.00007178711007*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(0,5,6)] = -0.00057084483436*constmodpars_ps->tc_P10;
         //
         constmodpars_ps->P10_dm->M[mos(1,1,6)] = 0.00734636595888*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(1,2,6)] = 0.00073958708312*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(1,3,6)] = -0.00374815005491*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(1,4,6)] = -0.00075481949242*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(1,5,6)] = 0.00027514643633*constmodpars_ps->tc_P10;
         //
         constmodpars_ps->P10_dm->M[mos(2,2,6)] = 0.00689027588484*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(2,3,6)] = -0.00081411533828*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(2,4,6)] = -0.00666515364393*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(2,5,6)] = -0.00123611890375*constmodpars_ps->tc_P10;
         //
         constmodpars_ps->P10_dm->M[mos(3,3,6)] = 0.00269200416549*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(3,4,6)] = 0.00090517462673*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(3,5,6)] = -0.00070599769660*constmodpars_ps->tc_P10;
         //
         constmodpars_ps->P10_dm->M[mos(4,4,6)] = 0.00694047472377*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->M[mos(4,5,6)] = -0.00135429545568*constmodpars_ps->tc_P10;
         //
         constmodpars_ps->P10_dm->M[mos(5,5,6)] = 0.01471395068969*constmodpars_ps->tc_P10;
         constmodpars_ps->P10_dm->flag = M_SU;
      }
      else if (constmodpars_ps->indxFixz10P10 == 2) {
         if (0) {//(constmodpars_ps->indxFcsts && constmodpars_ps->indxCnterFact) {
            constmodpars_ps->z10_dv->v[0] =  4.020000000000000e-002;
            constmodpars_ps->z10_dv->v[1] =  8.280000000000000e-002;
            constmodpars_ps->z10_dv->v[2] = -3.438000000000000e-001;
            constmodpars_ps->z10_dv->v[3] =  9.130000000000001e-002;
            constmodpars_ps->z10_dv->v[4] =  1.614000000000000e-001;
            constmodpars_ps->z10_dv->v[5] =  8.246200000000000e+000;
            //=== A limiting result with the initial condition at locti = 10.
            //  4.020000000000000e-002
            //  8.280000000000000e-002
            // -3.438000000000000e-001
            //  9.130000000000001e-002
            //  1.614000000000000e-001
            //  8.246200000000000e+000
         }
         else {
            constmodpars_ps->z10_dv->v[0] = -0.13236567947738;
            constmodpars_ps->z10_dv->v[1] =  0.14191629076828;
            constmodpars_ps->z10_dv->v[2] =  1.09275370936974;
            constmodpars_ps->z10_dv->v[3] = -0.02157077022996;
            constmodpars_ps->z10_dv->v[4] = -0.13378436709391;
            constmodpars_ps->z10_dv->v[5] =  0.21897028305873;
            //=== the estimated coefficients from ...\DataArrange\govregress48_59.m (from 1948:01 to 1959:12).
            // -0.13236567947738
            //  0.14191629076828
            //  1.09275370936974
            // -0.02157077022996
            // -0.13378436709391
            //  0.21897028305873
            //=== the estimated coefficients from ...\DataArrange\govregress48_91.m (from 1948:01 to 1991:12).
            //   -0.116749201518924
            //    0.132519395920216
            //     1.06582288792372
            // -0.00237964508526484
            //  -0.0862657390263775
            //   0.0691841583278944
         }
      }
   }
   else if (constmodpars_ps->indxWhichModel & KGP2TP1TI0) {
      //--- Use the estimated coefficients from ...\DataArrange\govregress48_59.m (from 1948:01 to 1959:12).
      if (constmodpars_ps->indxFixz10P10 == 1) {
         constmodpars_ps->z10_dv->v[0] = -0.13236567947738;
         constmodpars_ps->z10_dv->v[1] =  0.14191629076828;
         constmodpars_ps->z10_dv->v[2] =  1.09275370936974;
         constmodpars_ps->z10_dv->v[3] = -0.02157077022996;
         constmodpars_ps->z10_dv->v[4] = -0.13378436709391;
         constmodpars_ps->z10_dv->v[5] =  0.21897028305873;

         fn_DisplayError(".../swz_comfuns.c/CreateTSconstmodpars(): Have not got time to deal with Keynesian model with P10 fixed");
      }
      else if (constmodpars_ps->indxFixz10P10 == 2) {
         //--- WARNING: No conversion yet.  See (37.3) on p.37.  Converting the Keynesian direction (b0, b1, b2, b3, b4, b5) to the classical direction (a0, a1, a2, a3, a4, a5).
         constmodpars_ps->z10_dv->v[0] = -0.32823623917777;
         constmodpars_ps->z10_dv->v[1] =  0.31134751951751;
         constmodpars_ps->z10_dv->v[2] =  1.36439559257540;
         constmodpars_ps->z10_dv->v[3] = -0.01827865811705;
         constmodpars_ps->z10_dv->v[4] = -0.42333841556037;
         constmodpars_ps->z10_dv->v[5] =  0.27571505773370;
      }
   }




   //--- Almost converged results.
   // constmodpars_ps->z10_dv->v[0] = 14.81195;
   // constmodpars_ps->z10_dv->v[1] = -4.05202;
   // constmodpars_ps->z10_dv->v[2] = -4.22886;
   // constmodpars_ps->z10_dv->v[3] = -2.20691;
   // constmodpars_ps->z10_dv->v[4] = 0.95348 ;
   // constmodpars_ps->z10_dv->v[5] = -0.76895;
   //--- Almost converged results.
   // constmodpars_ps->P10_dm->M[mos(0,0,6)] = 106.32918;
   // constmodpars_ps->P10_dm->M[mos(0,1,6)] = -27.32784;
   // constmodpars_ps->P10_dm->M[mos(0,2,6)] = 2.79046;
   // constmodpars_ps->P10_dm->M[mos(0,3,6)] = 33.78822;
   // constmodpars_ps->P10_dm->M[mos(0,4,6)] = -4.92765;
   // constmodpars_ps->P10_dm->M[mos(0,5,6)] = -27.88119;
   // //
   // constmodpars_ps->P10_dm->M[mos(1,1,6)] = 114.09872;
   // constmodpars_ps->P10_dm->M[mos(1,2,6)] = 34.71816;
   // constmodpars_ps->P10_dm->M[mos(1,3,6)] = 25.18305;
   // constmodpars_ps->P10_dm->M[mos(1,4,6)] = -0.84149;
   // constmodpars_ps->P10_dm->M[mos(1,5,6)] = 72.03649;
   // //
   // constmodpars_ps->P10_dm->M[mos(2,2,6)] = 34.44133;
   // constmodpars_ps->P10_dm->M[mos(2,3,6)] = 47.35676;
   // constmodpars_ps->P10_dm->M[mos(2,4,6)] = 5.10799;
   // constmodpars_ps->P10_dm->M[mos(2,5,6)] = 13.44840;
   // //
   // constmodpars_ps->P10_dm->M[mos(3,3,6)] = 92.12431;
   // constmodpars_ps->P10_dm->M[mos(3,4,6)] = 8.44715;
   // constmodpars_ps->P10_dm->M[mos(3,5,6)] = -23.52167;
   // //
   // constmodpars_ps->P10_dm->M[mos(4,4,6)] = 6.10698;
   // constmodpars_ps->P10_dm->M[mos(4,5,6)] = -13.72212;
   // //
   // constmodpars_ps->P10_dm->M[mos(5,5,6)] = 216.33734;
   // constmodpars_ps->P10_dm->flag = M_SU;



   //=== Prior settings.
   constmodpars_ps->tc1 = input_ps->tc1;  //(0, 1].  Tightness control for vstar, theta, theta1, and tau1.
   constmodpars_ps->tc2 = input_ps->tc2;  //>0.  Tightness control for zeta1 and zeta2.
   constmodpars_ps->tc3 = input_ps->tc3;  //>0.  Tightness control for the elements of the upper triangular decomposition of P_{1|0}.
   constmodpars_ps->tc4 = input_ps->tc4;  //>0.  Tightness control for the elements of the upper triangular decomposition of V.
   constmodpars_ps->tc5 = input_ps->tc5;
   //+
   constmodpars_ps->mean1stset_dv = NULL;
   constmodpars_ps->Var1stSet_dm = NULL;
   constmodpars_ps->cholVar1stSet_dm = NULL;
   constmodpars_ps->mean4thset_dv = NULL;
   constmodpars_ps->Var4thSet_dm = NULL;
   constmodpars_ps->cholVar4thSet_dm = NULL;
   constmodpars_ps->mean5thset_dv = NULL;
   constmodpars_ps->Var5thSet_dm = NULL;
   constmodpars_ps->cholVar5thSet_dm = NULL;
   constmodpars_ps->mean6thset_dv = NULL;
   constmodpars_ps->Var6thSet_dm = NULL;
   constmodpars_ps->cholVar6thSet_dm = NULL;


   constmodpars_ps->modest = est_func;
   constmodpars_ps->modprob = prob_func;

   return (constmodpars_ps);
}
//---
TSconstmodpars *DestroyTSconstmodpars(TSconstmodpars *constmodpars_ps)
{
   if (constmodpars_ps)
   {
      //--- Added 03/17/05.
      DestroyVector_lf(constmodpars_ps->ngo_x_dv);


      DestroyVector_lf(constmodpars_ps->x_dv);
      DestroyMatrix_lf(constmodpars_ps->V_dm);
      DestroyVector_lf(constmodpars_ps->z10_dv);
      DestroyVector_lf(constmodpars_ps->k_fcstinf_z10_dv);
      DestroyMatrix_lf(constmodpars_ps->P10_dm);
      //
      DestroyVector_lf(constmodpars_ps->upred_dv);
      DestroyMatrix_lf(constmodpars_ps->StrResidstran_dm);
      DestroyVector_lf(constmodpars_ps->sclhistshocks_dv);
      //
      DestroyDatesVector_lf(constmodpars_ps->sampledates_dv);


      DestroyVector_lf(constmodpars_ps->mean1stset_dv);
      DestroyMatrix_lf(constmodpars_ps->Var1stSet_dm);
      DestroyMatrix_lf(constmodpars_ps->cholVar1stSet_dm);
      DestroyVector_lf(constmodpars_ps->mean4thset_dv);
      DestroyMatrix_lf(constmodpars_ps->Var4thSet_dm);
      DestroyMatrix_lf(constmodpars_ps->cholVar4thSet_dm);
      DestroyVector_lf(constmodpars_ps->mean5thset_dv);
      DestroyMatrix_lf(constmodpars_ps->Var5thSet_dm);
      DestroyMatrix_lf(constmodpars_ps->cholVar5thSet_dm);
      DestroyVector_lf(constmodpars_ps->mean6thset_dv);
      DestroyMatrix_lf(constmodpars_ps->Var6thSet_dm);
      DestroyMatrix_lf(constmodpars_ps->cholVar6thSet_dm);


      //===
      free(constmodpars_ps);
      return ((TSconstmodpars *)NULL);
   }
   else  return (constmodpars_ps);
}


void ftd_freepars2modpars(void *constmod_ps, int pointflag)
{
   //Copying values of the vectorized free (minimization) parameters to regular model parameters.
   struct TSconstmodpars_tag *constmodpars_ps;
   struct TSconstgibbs_tag *constgibbs_ps;
   int _i, _j;
   int cnt;
   int kx;
   double workd;
   TSdvector *x_dv;
   //===
   TSdmatrix *cholMkbyk_dm = NULL;

   if (pointflag & POINTEST) {
      constmodpars_ps = (struct TSconstmodpars_tag *)constmod_ps;
      kx = constmodpars_ps->kx;
      x_dv = constmodpars_ps->x_dv;
      cholMkbyk_dm = CreateConstantMatrix_lf(kx,kx,0.0);
      cholMkbyk_dm->flag = M_GE | M_UT;  //Will become upper triangular and general as well.

      if (!x_dv->flag)  fn_DisplayError(".../ftd_freepars2modpars(): x_dv must have legal values");

//      if (!constmodpars_ps->nlagsinflation && (constmodpars_ps->nlagstrue==1) && !constmodpars_ps->govmodel) {
      if (constmodpars_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0 | NGO_CGP2TP1TI0))  //NGO_CGP2TP1TI0 is added 03/17/05.
      {
         //No lag and thus no rho1 and rho2 in the true inflation process to be estimated.  See p.20a.
         //6+kx+(kx+1)*kx variables in the order of: //vstar, theta0, theta1, tau1, xi1 (sqrt(zeta1)), xi2 (sqrt(zeta2)), z10, P10, and V.
         constmodpars_ps->gtau1 = x_dv->v[3];
         if (fabs(workd=1.0-constmodpars_ps->gtau1)>MACHINEZERO)  constmodpars_ps->ustar = x_dv->v[0]/workd;
         else {
            //Sets ustar to any arbitrary number and x_dv->v[0] to zero to make it compatible.
            constmodpars_ps->ustar = 5.0;
            x_dv->v[0] = 0.0;
         }
         // else if (fabs(x_dv->v[0])<MACHINEZERO)  constmodpars_ps->ustar = 5.0;   //Sets ustar to any arbitrary number.
         // else  fn_DisplayError(".../ftd_freepars2modpars(): 1.0-gtau1 is zero but vstar is not, thus no solution for ustar (see p.20a)");
                                                //Done with all work* arrays.
         constmodpars_ps->gtheta0 = x_dv->v[1];
         constmodpars_ps->gtheta1 = x_dv->v[2];
         constmodpars_ps->gzeta1 = square(x_dv->v[4]);
         constmodpars_ps->gzeta2 = square(x_dv->v[5]);
         //+
         if (constmodpars_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0))  //Added 03/17/05.
         {
            if (!constmodpars_ps->indxFixz10P10) {
               for (_i=0; _i<kx; _i++)  constmodpars_ps->z10_dv->v[_i] = x_dv->v[6+_i];
               constmodpars_ps->z10_dv->flag = V_DEF;
               for (cnt=0, _j=0; _j<kx; _j++)
                  for (_i=0; _i<=_j; _i++)
                     cholMkbyk_dm->M[mos(_i,_j,kx)] = x_dv->v[6+kx+cnt++];
               MatrixTimesSelf(constmodpars_ps->P10_dm, 'U', cholMkbyk_dm, 'T', 1.0, 0.0);  //P = C'*C;
               for (_j=0; _j<kx; _j++)
                  for (_i=0; _i<=_j; _i++)
                     cholMkbyk_dm->M[mos(_i,_j,kx)] = x_dv->v[6+kx+cnt++];
               MatrixTimesSelf(constmodpars_ps->V_dm, 'U', cholMkbyk_dm, 'T', 1.0, 0.0);  //V = C'*C;
            }
            else if (constmodpars_ps->indxFixz10P10 == 2) {
               for (cnt=0, _j=0; _j<kx; _j++)
                  for (_i=0; _i<=_j; _i++)
                     cholMkbyk_dm->M[mos(_i,_j,kx)] = x_dv->v[6+cnt++];
               MatrixTimesSelf(constmodpars_ps->P10_dm, 'U', cholMkbyk_dm, 'T', 1.0, 0.0);  //P = C'*C;
               for (_j=0; _j<kx; _j++)
                  for (_i=0; _i<=_j; _i++)
                     cholMkbyk_dm->M[mos(_i,_j,kx)] = x_dv->v[6+cnt++];
               MatrixTimesSelf(constmodpars_ps->V_dm, 'U', cholMkbyk_dm, 'T', 1.0, 0.0);  //V = C'*C;
            }
            else if (constmodpars_ps->indxFixz10P10 == 1) {
               for (cnt=0, _j=0; _j<kx; _j++)
                  for (_i=0; _i<=_j; _i++)
                     cholMkbyk_dm->M[mos(_i,_j,kx)] = x_dv->v[6+cnt++];
               MatrixTimesSelf(constmodpars_ps->V_dm, 'U', cholMkbyk_dm, 'T', 1.0, 0.0);  //V = C'*C;
            }
            else  fn_DisplayError(".../swz_comfuns.c/ftd_freepars2modpars(): indxFixz10P10 must be 0, 1, or 2");
         }
      }
      else {
         printf(".../swz_comfuns.c/ftd_freepars2modpars(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d!", CGP2TP1TI0, KGP2TP1TI0);
         exit(EXIT_FAILURE);
      }
//      else  fn_DisplayError(".../swz_comfuns.c/ftd_freepars2modpars(): Have not got time to do the case for indxWhichModel be other than %dnlagstrue != 1 or nlagsinflation !=0 or govmodel != 0");
   }
   else {
      //For GIBBSDRAWS (Gibbs draws)
      constgibbs_ps = (struct TSconstgibbs_tag *)constmod_ps;
      kx = constgibbs_ps->kx;
      x_dv = constgibbs_ps->x_dv;
      cholMkbyk_dm = CreateConstantMatrix_lf(kx,kx,0.0);
      cholMkbyk_dm->flag = M_GE | M_UT;  //Will become upper triangular and general as well.

      if (!x_dv->flag)  fn_DisplayError(".../ftd_freepars2modpars(): x_dv must have legal values");

//      if (!constgibbs_ps->nlagsinflation && (constgibbs_ps->nlagstrue==1) && !constgibbs_ps->govmodel) {
      if (constgibbs_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0)) {
         //No lag and thus no rho1 and rho2 in the true inflation process to be estimated.  See p.20a.
         //6+kx+(kx+1)*kx variables in the order of: //vstar, theta0, theta1, tau1, xi1 (sqrt(zeta1)), xi2 (sqrt(zeta2)), z10, P10, and V.
         constgibbs_ps->gtau1 = x_dv->v[3];
         if (fabs(workd=1.0-constgibbs_ps->gtau1)>MACHINEZERO)  constgibbs_ps->ustar = x_dv->v[0]/workd;
         else {
            //Sets ustar to any arbitrary number and x_dv->v[0] to zero to make it compatible.
            constgibbs_ps->ustar = 5.0;
            x_dv->v[0] = 0.0;
         }
                                                //Done with all work* arrays.
         constgibbs_ps->gtheta0 = x_dv->v[1];
         constgibbs_ps->gtheta1 = x_dv->v[2];
         constgibbs_ps->gzeta1 = square(x_dv->v[4]);
         constgibbs_ps->gzeta2 = square(x_dv->v[5]);
         //+
         if (!constgibbs_ps->indxFixz10P10) {
            for (_i=0; _i<kx; _i++)  constgibbs_ps->z10_dv->v[_i] = x_dv->v[6+_i];
            constgibbs_ps->z10_dv->flag = V_DEF;
            for (cnt=0, _j=0; _j<kx; _j++)
               for (_i=0; _i<=_j; _i++)
                  cholMkbyk_dm->M[mos(_i,_j,kx)] = x_dv->v[6+kx+cnt++];
            MatrixTimesSelf(constgibbs_ps->P10_dm, 'U', cholMkbyk_dm, 'T', 1.0, 0.0);  //P = C'*C;
            for (_j=0; _j<kx; _j++)
               for (_i=0; _i<=_j; _i++)
                  cholMkbyk_dm->M[mos(_i,_j,kx)] = x_dv->v[6+kx+cnt++];
            MatrixTimesSelf(constgibbs_ps->V_dm, 'U', cholMkbyk_dm, 'T', 1.0, 0.0);  //V = C'*C;
         }
         else if (constgibbs_ps->indxFixz10P10 == 2) {
            for (cnt=0, _j=0; _j<kx; _j++)
               for (_i=0; _i<=_j; _i++)
                  cholMkbyk_dm->M[mos(_i,_j,kx)] = x_dv->v[6+cnt++];
            MatrixTimesSelf(constgibbs_ps->P10_dm, 'U', cholMkbyk_dm, 'T', 1.0, 0.0);  //P = C'*C;
            for (_j=0; _j<kx; _j++)
               for (_i=0; _i<=_j; _i++)
                  cholMkbyk_dm->M[mos(_i,_j,kx)] = x_dv->v[6+cnt++];
            MatrixTimesSelf(constgibbs_ps->V_dm, 'U', cholMkbyk_dm, 'T', 1.0, 0.0);  //V = C'*C;
         }
         else if (constgibbs_ps->indxFixz10P10 == 1) {
            for (cnt=0, _j=0; _j<kx; _j++)
               for (_i=0; _i<=_j; _i++)
                  cholMkbyk_dm->M[mos(_i,_j,kx)] = x_dv->v[6+cnt++];
            MatrixTimesSelf(constgibbs_ps->V_dm, 'U', cholMkbyk_dm, 'T', 1.0, 0.0);  //V = C'*C;
         }
         else  fn_DisplayError(".../swz_comfuns.c/ftd_freepars2modpars(): indxFixz10P10 must be 0, 1, or 2");
      }
      else {
         printf(".../swz_comfuns.c/ftd_freepars2modpars(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d!", CGP2TP1TI0, KGP2TP1TI0);
         exit(EXIT_FAILURE);
      }
//      else  fn_DisplayError(".../swz_comfuns.c/ftd_freepars2modpars(): Have not got time to do the case for nlagstrue != 1 or nlagsinflation !=0 or govmodel != 0");
   }


   //===
   DestroyMatrix_lf(cholMkbyk_dm);
}
//---
void ftd_modpars2freepars(void *constmod_ps, int pointflag)
{
   //Copying values of regular model parameters to the vectorized free (minimization) parameters.
   struct TSconstmodpars_tag *constmodpars_ps;
   struct TSconstgibbs_tag *constgibbs_ps;
   int _i, _j;
   int cnt;
   int kx;
   TSdvector *x_dv;
   //===
   TSdmatrix *cholMkbyk_dm = NULL;

   if (pointflag & POINTEST) {
      constmodpars_ps = (TSconstmodpars *)constmod_ps;
      kx = constmodpars_ps->kx;
      x_dv = constmodpars_ps->x_dv;
      cholMkbyk_dm = CreateMatrix_lf(kx,kx);   //Will become upper triangular.

//      if (!constmodpars_ps->nlagsinflation && (constmodpars_ps->nlagstrue==1) && !constmodpars_ps->govmodel) {
      if (constmodpars_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0 | NGO_CGP2TP1TI0))  //NGO_CGP2TP1TI0 is added 03/17/05.
      {
         //No lag and thus no rho1 and rho2 in the true inflation process to be estimated.  See p.20a.
         //6+kx+(kx+1)*kx variables in the order of: //vstar, theta0, theta1, tau1, xi1 (sqrt(zeta1)), xi2 (sqrt(zeta2)), z10, P10, and V.
         x_dv->v[0] = constmodpars_ps->ustar*(1.0-constmodpars_ps->gtau1);   //vstar in the notation of p.20a.
         x_dv->v[1] = constmodpars_ps->gtheta0;
         x_dv->v[2] = constmodpars_ps->gtheta1;
         x_dv->v[3] = constmodpars_ps->gtau1;
         x_dv->v[4] = sqrt(constmodpars_ps->gzeta1);
         x_dv->v[5] = sqrt(constmodpars_ps->gzeta2);
         //+
         if (constmodpars_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0))  //Added 03/17/05.
         {
            if (!constmodpars_ps->indxFixz10P10) {
               for (_i=0; _i<kx; _i++)  x_dv->v[6+_i] = constmodpars_ps->z10_dv->v[_i];
               if (chol(cholMkbyk_dm, constmodpars_ps->P10_dm, 'U'))  fn_DisplayError(".../swz_comfuns.c/modpars2freepars(): Choleski decomposition of P10_dm fails when calling chol()");
               for (cnt=0, _j=0; _j<kx; _j++)
                  for (_i=0; _i<=_j; _i++)
                     x_dv->v[6+kx+cnt++] = cholMkbyk_dm->M[mos(_i,_j,kx)];
               if (chol(cholMkbyk_dm, constmodpars_ps->V_dm, 'U'))  fn_DisplayError(".../swz_comfuns.c/modpars2freepars(): Choleski decomposition of V_dm fails when calling chol()");
               for (_j=0; _j<kx; _j++)
                  for (_i=0; _i<=_j; _i++)
                     x_dv->v[6+kx+cnt++] = cholMkbyk_dm->M[mos(_i,_j,kx)];
            }
            else if (constmodpars_ps->indxFixz10P10 == 2) {
               if (chol(cholMkbyk_dm, constmodpars_ps->P10_dm, 'U')<0) {
                  //fprintf(fptr_debug, "cholP10_dm with error %d:\n", tmpi);
                  //WriteMatrix(fptr_debug, cholMkbyk_dm, " %.16e ");
                  fn_DisplayError(".../swz_comfuns.c/modpars2freepars(): Choleski decomposition of P10_dm fails when calling chol()");
               }
               for (cnt=0, _j=0; _j<kx; _j++)
                  for (_i=0; _i<=_j; _i++)
                     x_dv->v[6+cnt++] = cholMkbyk_dm->M[mos(_i,_j,kx)];
               if (chol(cholMkbyk_dm, constmodpars_ps->V_dm, 'U'))  fn_DisplayError(".../swz_comfuns.c/modpars2freepars(): Choleski decomposition of V_dm fails when calling chol()");
               for (_j=0; _j<kx; _j++)
                  for (_i=0; _i<=_j; _i++)
                     x_dv->v[6+cnt++] = cholMkbyk_dm->M[mos(_i,_j,kx)];
            }
            else if (constmodpars_ps->indxFixz10P10 == 1) {
               if (chol(cholMkbyk_dm, constmodpars_ps->V_dm, 'U'))  fn_DisplayError(".../swz_comfuns.c/modpars2freepars(): Choleski decomposition of V_dm fails when calling chol()");
               for (cnt=0, _j=0; _j<kx; _j++)
                  for (_i=0; _i<=_j; _i++)
                     x_dv->v[6+cnt++] = cholMkbyk_dm->M[mos(_i,_j,kx)];
            }
            else  fn_DisplayError(".../swz_comfuns.c/ftd_modpars2freepars(): indxFixz10P10 must be 0, 1, or 2");
         }
         x_dv->flag = V_DEF;
      }
      else {
         printf(".../swz_comfuns.c/ftd_modpars2freepars(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d!", CGP2TP1TI0, KGP2TP1TI0);
         exit(EXIT_FAILURE);
      }
//      else  fn_DisplayError(".../swz_comfuns.c/ftd_modpars2freepars(): Have not got time to do the case for nlagstrue != 1 or nlagsinflation !=0 or govmodel != 0");
   }
   else {
      //For GIBBSDRAWS (Gibbs draws).
      constgibbs_ps = (struct TSconstgibbs_tag *)constmod_ps;
      kx = constgibbs_ps->kx;
      x_dv = constgibbs_ps->x_dv;
      cholMkbyk_dm = CreateMatrix_lf(kx,kx);   //Will become upper triangular.

//      if (!constgibbs_ps->nlagsinflation && (constgibbs_ps->nlagstrue==1) && !constgibbs_ps->govmodel) {
      if (!constgibbs_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0)) {
         //No lag and thus no rho1 and rho2 in the true inflation process to be estimated.  See p.20a.
         //6+kx+(kx+1)*kx variables in the order of: //vstar, theta0, theta1, tau1, xi1 (sqrt(zeta1)), xi2 (sqrt(zeta2)), z10, P10, and V.
         x_dv->v[0] = constgibbs_ps->ustar*(1.0-constgibbs_ps->gtau1);   //vstar in the notation of p.20a.
         x_dv->v[1] = constgibbs_ps->gtheta0;
         x_dv->v[2] = constgibbs_ps->gtheta1;
         x_dv->v[3] = constgibbs_ps->gtau1;
         x_dv->v[4] = sqrt(constgibbs_ps->gzeta1);
         x_dv->v[5] = sqrt(constgibbs_ps->gzeta2);
         //+
         if (!constgibbs_ps->indxFixz10P10) {
            for (_i=0; _i<kx; _i++)  x_dv->v[6+_i] = constgibbs_ps->z10_dv->v[_i];
            if (chol(cholMkbyk_dm, constgibbs_ps->P10_dm, 'U'))  fn_DisplayError(".../modpars2freepars(): Choleski decomposition of P10_dm fails when calling chol()");
            for (cnt=0, _j=0; _j<kx; _j++)
               for (_i=0; _i<=_j; _i++)
                  x_dv->v[6+kx+cnt++] = cholMkbyk_dm->M[mos(_i,_j,kx)];
            if (chol(cholMkbyk_dm, constgibbs_ps->V_dm, 'U'))  fn_DisplayError(".../modpars2freepars(): Choleski decomposition of V_dm fails when calling chol()");
            for (_j=0; _j<kx; _j++)
               for (_i=0; _i<=_j; _i++)
                  x_dv->v[6+kx+cnt++] = cholMkbyk_dm->M[mos(_i,_j,kx)];
         }
         else if (constgibbs_ps->indxFixz10P10 == 2) {
            if (chol(cholMkbyk_dm, constgibbs_ps->P10_dm, 'U'))  fn_DisplayError(".../modpars2freepars(): Choleski decomposition of P10_dm fails when calling chol()");
            for (cnt=0, _j=0; _j<kx; _j++)
               for (_i=0; _i<=_j; _i++)
                  x_dv->v[6+cnt++] = cholMkbyk_dm->M[mos(_i,_j,kx)];
            if (chol(cholMkbyk_dm, constgibbs_ps->V_dm, 'U'))  fn_DisplayError(".../modpars2freepars(): Choleski decomposition of V_dm fails when calling chol()");
            for (_j=0; _j<kx; _j++)
               for (_i=0; _i<=_j; _i++)
                  x_dv->v[6+cnt++] = cholMkbyk_dm->M[mos(_i,_j,kx)];
         }
         else if (constgibbs_ps->indxFixz10P10 == 1) {
            if (chol(cholMkbyk_dm, constgibbs_ps->V_dm, 'U'))  fn_DisplayError(".../modpars2freepars(): Choleski decomposition of V_dm fails when calling chol()");
            for (cnt=0, _j=0; _j<kx; _j++)
               for (_i=0; _i<=_j; _i++)
                  x_dv->v[6+cnt++] = cholMkbyk_dm->M[mos(_i,_j,kx)];
         }
         else  fn_DisplayError(".../swz_comfuns.c/ftd_modpars2freepars(): indxFixz10P10 must be 0, 1, or 2");
         x_dv->flag = V_DEF;
      }
      else {
         printf(".../swz_comfuns.c/ftd_modpars2freepars(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d!", CGP2TP1TI0, KGP2TP1TI0);
         exit(EXIT_FAILURE);
      }
//      else  fn_DisplayError(".../swz_comfuns.c/ftd_modpars2freepars(): Have not got time to do the case for nlagstrue != 1 or nlagsinflation !=0 or govmodel != 0");
   }


   //===
   DestroyMatrix_lf(cholMkbyk_dm);
}


//======= 03/17/05.  No government optimization problem but estimating chol(V), chol(P10), and sigma directly. =======
void ngo_freepars2modpars(void *constmod_ps, int pointflag)
{
   //Copying values of the vectorized free (minimization) parameters to regular model parameters.
   struct TSconstmodpars_tag *constmodpars_ps;
   int _i, _j;
   int cnt;
   int kx;
   TSdvector *x_dv;
   //===
   TSdmatrix *cholMkbyk_dm = NULL;

   if (pointflag & POINTEST)
   {
      constmodpars_ps = (struct TSconstmodpars_tag *)constmod_ps;
      kx = constmodpars_ps->kx;
      x_dv = constmodpars_ps->ngo_x_dv;
      cholMkbyk_dm = CreateConstantMatrix_lf(kx,kx,0.0);
      cholMkbyk_dm->flag = M_GE | M_UT;  //Will become upper triangular and general as well.

      if (!x_dv->flag)  fn_DisplayError(".../ngo_freepars2modpars(): ngo_x_dv must have legal values");

      //(kx+1)*kx variables in the order of: 6*7/3 free elements for chol(P10) and 6*7/3 free elements for chol(V).
      for (cnt=0, _j=0; _j<kx; _j++)
         for (_i=0; _i<=_j; _i++)
            cholMkbyk_dm->M[mos(_i,_j,kx)] = x_dv->v[cnt++];
      MatrixTimesSelf(constmodpars_ps->P10_dm, 'U', cholMkbyk_dm, 'T', 1.0, 0.0);  //P = C'*C;
      for (_j=0; _j<kx; _j++)
         for (_i=0; _i<=_j; _i++)
            cholMkbyk_dm->M[mos(_i,_j,kx)] = x_dv->v[cnt++];
      MatrixTimesSelf(constmodpars_ps->V_dm, 'U', cholMkbyk_dm, 'T', 1.0, 0.0);  //V = C'*C;
   }
   else  fn_DisplayError(".../swz_comfuns.c/ngo_freepars2modpars(): no Gibbs draws for this revision");


   //===
   DestroyMatrix_lf(cholMkbyk_dm);
}
//---
void ngo_modpars2freepars(void *constmod_ps, int pointflag)
{
   //Copying values of regular model parameters to the vectorized free (minimization) parameters.
   struct TSconstmodpars_tag *constmodpars_ps;
   int _i, _j;
   int cnt;
   int kx;
   TSdvector *x_dv;
   //===
   TSdmatrix *cholMkbyk_dm = NULL;

   if (pointflag & POINTEST) {
      constmodpars_ps = (TSconstmodpars *)constmod_ps;
      kx = constmodpars_ps->kx;
      x_dv = constmodpars_ps->ngo_x_dv;
      cholMkbyk_dm = CreateMatrix_lf(kx,kx);   //Will become upper triangular.

      //(kx+1)*kx variables in the order of: sigma, 6*7/3 free elements for chol(P10), and 6*7/3 free elements for chol(V).
      if (chol(cholMkbyk_dm, constmodpars_ps->P10_dm, 'U'))  fn_DisplayError(".../swz_comfuns.c/ngo_modpars2freepars(): Choleski decomposition of P10_dm fails when calling chol()");
      for (cnt=0, _j=0; _j<kx; _j++)
         for (_i=0; _i<=_j; _i++)
            x_dv->v[cnt++] = cholMkbyk_dm->M[mos(_i,_j,kx)];
      if (chol(cholMkbyk_dm, constmodpars_ps->V_dm, 'U'))  fn_DisplayError(".../swz_comfuns.c/ngo_modpars2freepars(): Choleski decomposition of V_dm fails when calling chol()");
      for (_j=0; _j<kx; _j++)
         for (_i=0; _i<=_j; _i++)
            x_dv->v[cnt++] = cholMkbyk_dm->M[mos(_i,_j,kx)];
   }
   else  fn_DisplayError(".../swz_comfuns.c/ngo_modpars2freepars(): no Gibbs draws for this revision");


   //===
   DestroyMatrix_lf(cholMkbyk_dm);
}



//----------------------------
// Updating the model's free parameters from the subset, i.e. the government's free parameters.
//----------------------------
void updatemodfreepars5govfreepars(void *modpars_ps, int pointflag, TSdvector *xgov_dv)
{
   int locgovpars = N1STSET + N2NDSET + N3RDSET;
   struct TSconstmodpars_tag *constmodpars_ps;
   struct TSconstgibbs_tag *constgibbs_ps;


   if (pointflag & POINTEST) {
      constmodpars_ps = (struct TSconstmodpars_tag *)modpars_ps;
      if ((locgovpars+xgov_dv->n) > constmodpars_ps->x_dv->n)  fn_DisplayError(".../swz_comfuns.c/updatemodfreepars5govfreepars(): dimension of xgov_dv exceeds the length of constmodpars_ps->x_dv when copying");
      CopySubvector(constmodpars_ps->x_dv, locgovpars, xgov_dv, 0, xgov_dv->n);
   }
   else {
      constgibbs_ps = (struct TSconstgibbs_tag *)modpars_ps;
      if ((locgovpars+xgov_dv->n) > constgibbs_ps->x_dv->n)  fn_DisplayError(".../swz_comfuns.c/updatemodfreepars5govfreepars(): dimension of xgov_dv exceeds the length of constgibbs_ps->x_dv when copying");
      CopySubvector(constgibbs_ps->x_dv, locgovpars, xgov_dv, 0, xgov_dv->n);
   }
}
//----------------------------
// Updating the government's free parameters from the model's free parameters.
//----------------------------
void updategovfreepars5modfreepars(void *modpars_ps, int pointflag, TSdvector *xgov_dv)
{
   int locgovpars = N1STSET + N2NDSET + N3RDSET;
   struct TSconstmodpars_tag *constmodpars_ps;
   struct TSconstgibbs_tag *constgibbs_ps;


   if (pointflag & POINTEST) {
      constmodpars_ps = (struct TSconstmodpars_tag *)modpars_ps;
      if ((locgovpars+xgov_dv->n) > constmodpars_ps->x_dv->n)  fn_DisplayError(".../swz_comfuns.c/updategovfreepars5modfreepars(): dimension of xgov_dv exceeds the length of constmodpars_ps->x_dv when copying");
      CopySubvector(xgov_dv, 0, constmodpars_ps->x_dv, locgovpars, xgov_dv->n);
   }
   else {
      constgibbs_ps = (struct TSconstgibbs_tag *)modpars_ps;
      if ((locgovpars+xgov_dv->n) > constgibbs_ps->x_dv->n)  fn_DisplayError(".../swz_comfuns.c/updategovfreepars5modfreepars(): dimension of xgov_dv exceeds the length of constgibbs_ps->x_dv when copying");
      CopySubvector(xgov_dv, 0, constgibbs_ps->x_dv, locgovpars, xgov_dv->n);
   }
}


//----------------------------
// Computes the log value of prior pdf (not just kernal).  See p.20b.
//----------------------------
double value_logpriorpdf(struct TSconstmodpars_tag *constmodpars_ps, struct TSconstgibbs_tag *constgibbs_ps, int pointflag)
{
   //Must be comformable to ftd_freepars2modpars().
   int cntloc, _m;
   double valuelogpdf;
   //---
   double gzeta1p2;
   double gzeta1p2_prioralpha;
   double gzeta1p2_priorbetainv;
   //---
   TSdvector *x_dv;
   TSdvector xmveset_sdv;
   TSdvector *mean1stset_dv;
   TSdmatrix *cholVar1stSet_dm;
   TSdvector *mean4thset_dv;
   TSdmatrix *cholVar4thSet_dm;
   TSdvector *mean5thset_dv;
   TSdmatrix *cholVar5thSet_dm;
   TSdvector *mean6thset_dv;
   TSdmatrix *cholVar6thSet_dm;
   //+
   TSdvector *work_dv = NULL;

   if (pointflag & POINTEST)  x_dv = constmodpars_ps->x_dv;
   else  x_dv = constgibbs_ps->x_dv;

//   if (!constmodpars_ps->nlagsinflation && constmodpars_ps->nlagstrue==1) {
   if (constmodpars_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0 | NGO_CGP2TP1TI0))  //NGO_CGP2TP1TI0 is added 03/17/05.
   {
      //No lag and thus no rho1 and rho2 in the true inflation process to be estimated.
      //6+kx+(kx+1)*kx variables in the order of: //vstar, theta0, theta1, tau1, xi1 (sqrt(zeta1)), xi2 (sqrt(zeta2)), z10, P10, and V.

      //=== Prior settings.
      mean1stset_dv = constmodpars_ps->mean1stset_dv;
      cholVar1stSet_dm = constmodpars_ps->cholVar1stSet_dm;
      mean4thset_dv = constmodpars_ps->mean4thset_dv;
      cholVar4thSet_dm = constmodpars_ps->cholVar4thSet_dm;
      mean5thset_dv = constmodpars_ps->mean5thset_dv;
      cholVar5thSet_dm = constmodpars_ps->cholVar5thSet_dm;
      mean6thset_dv = constmodpars_ps->mean6thset_dv;
      cholVar6thSet_dm = constmodpars_ps->cholVar6thSet_dm;

      //--------------------------
      // The order of the following lines matters!
      //--------------------------
      valuelogpdf = 0.0;   //Initialization.
      xmveset_sdv.flag = V_DEF;
      //------- logpriorpdf for 1st set of parameters (vstar, theta0, theta1, and tau1). See p.20b. -------
      cntloc = 0;
      xmveset_sdv.n = N1STSET;
      xmveset_sdv.v = x_dv->v;
      work_dv = CreateVector_lf(N1STSET);
      //+
      VectorMinusVector(work_dv, &xmveset_sdv, mean1stset_dv);
      SolveTriSysVector(work_dv, cholVar1stSet_dm, work_dv, 'T', 'N');
      valuelogpdf += constmodpars_ps->term1st1 - 0.5*VectorDotVector(work_dv, work_dv);
      //+
      work_dv = DestroyVector_lf(work_dv);

      //------- logpriorpdf for 2nd and 3rd sets of parameters (i.e., zeta1 and zeta2). See p.20b. -------
      gzeta1p2_prioralpha = constmodpars_ps->gzeta1p2_prioralpha;
      gzeta1p2_priorbetainv = 1.0/constmodpars_ps->gzeta1p2_priorbeta;
      //
      gzeta1p2 = x_dv->v[cntloc += N1STSET];
      gzeta1p2 = square(gzeta1p2);
      valuelogpdf += constmodpars_ps->term2nd1 + (gzeta1p2_prioralpha-1.0)*log(gzeta1p2) - (gzeta1p2*gzeta1p2_priorbetainv);
      //+
      gzeta1p2 = x_dv->v[cntloc += N2NDSET];
      gzeta1p2 = square(gzeta1p2);
      valuelogpdf += constmodpars_ps->term3rd1 + (gzeta1p2_prioralpha-1.0)*log(gzeta1p2) - (gzeta1p2*gzeta1p2_priorbetainv);


      if (constmodpars_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0))   //Added 03/17/05.
      {
         if (!constmodpars_ps->indxFixz10P10) {
             //------- logpriorpdf for 4th set of parameters (i.e., z_{1|0}, coefficients on pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and Constant). See p.20b. -------
            xmveset_sdv.n = constmodpars_ps->kx;
            xmveset_sdv.v = x_dv->v + (cntloc += N3RDSET);
            work_dv = CreateVector_lf(constmodpars_ps->kx);
            //+
            VectorMinusVector(work_dv, &xmveset_sdv, mean4thset_dv);
            SolveTriSysVector(work_dv, cholVar4thSet_dm, work_dv, 'T', 'N');
            valuelogpdf += constmodpars_ps->term4th1 - 0.5*VectorDotVector(work_dv, work_dv);
            //+
            work_dv = DestroyVector_lf(work_dv);

            //------- logpriorpdf for 5th set of parameters (i.e., only upper triangular part of the Choleski decompostion of P_{1|0}). See p.20b. -------
            xmveset_sdv.n = _m = constmodpars_ps->kx*(constmodpars_ps->kx+1)/2;
            xmveset_sdv.v = x_dv->v + (cntloc += constmodpars_ps->kx);
            work_dv = CreateVector_lf(_m);
            //+
            VectorMinusVector(work_dv, &xmveset_sdv, mean5thset_dv);
            SolveTriSysVector(work_dv, cholVar5thSet_dm, work_dv, 'T', 'N');
            valuelogpdf += constmodpars_ps->term5th1 - 0.5*VectorDotVector(work_dv, work_dv);
            //+
            work_dv = DestroyVector_lf(work_dv);
         }
         else if (constmodpars_ps->indxFixz10P10 == 2) {
            //------- logpriorpdf for 5th set of parameters (i.e., only upper triangular part of the Choleski decompostion of P_{1|0}). See p.20b. -------
            xmveset_sdv.n = _m = constmodpars_ps->kx*(constmodpars_ps->kx+1)/2;
            xmveset_sdv.v = x_dv->v + (cntloc += N3RDSET);
            work_dv = CreateVector_lf(_m);
            //+
            VectorMinusVector(work_dv, &xmveset_sdv, mean5thset_dv);
            SolveTriSysVector(work_dv, cholVar5thSet_dm, work_dv, 'T', 'N');
            valuelogpdf += constmodpars_ps->term5th1 - 0.5*VectorDotVector(work_dv, work_dv);
            //+
            work_dv = DestroyVector_lf(work_dv);
         }

         //------- logpriorpdf for 6th set of parameters (i.e., only upper triangular part of the Choleski decompostion of V). See p.20b. -------
         xmveset_sdv.n = _m = constmodpars_ps->kx*(constmodpars_ps->kx+1)/2;
         xmveset_sdv.v = x_dv->v + (cntloc += N3RDSET);
         work_dv = CreateVector_lf(_m);
         //+
         VectorMinusVector(work_dv, &xmveset_sdv, mean6thset_dv);
         SolveTriSysVector(work_dv, cholVar6thSet_dm, work_dv, 'T', 'N');
         valuelogpdf +=  constmodpars_ps->term6th1 - 0.5*VectorDotVector(work_dv, work_dv);
         //+
         work_dv = DestroyVector_lf(work_dv);
      }
   }
   else if (constmodpars_ps->nlagstrue != 1)  fn_DisplayError(".../ftd_logprior(): Have not got time to work out the case for nlagstrue != 1");
   else if (constmodpars_ps->nlagsinflation > 0)  fn_DisplayError(".../ftd_logprior(): Have not got time to code up the case for nlagsinflation > 0");
   else {
      printf(".../swz_comfuns.c/value_logpriorpdf(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d or %d!\n", CGP2TP1TI0, KGP2TP1TI0, NGO_CGP2TP1TI0);
      exit(EXIT_FAILURE);
   }

   //===
   DestroyVector_lf(work_dv);
   //
   return (valuelogpdf);
}


//-------------------------------
// (1) Copies xchange_dv to constmodpars_ps->x_dv or constgibbs_ps->x_dv, depending on the flag "pointflag."
// (2) Refreshes all the parameters in constmodpars_ps or constgibbs_ps, using ftd_freepars2modpars().
// (3) Update the government's controled inflation and expected inflation, using govprob_ps->govestimator().
// (4) Refreshes one-step predictions of unemployment, upred_dv, for the case of POINTEST.
// (5) Computes the log value of LHpdf * Priorpdf.  See p.20a.
//-------------------------------
#define NEGHALFLOG2PI  (-0.91893853320467)   //-0.5*log(2*pi)
double value_logpostkernal_refresh(struct TSconstmodpars_tag *constmodpars_ps, struct TSconstgibbs_tag *constgibbs_ps, int pointflag, TSgovprob *govprob_ps, TSdvector *xchange_dv)
{
   int ti;
   int fss;  //Effective sample size.
   double workdterm, z1t, z2t;
   double logval;

   //--- Initialization. The order matters!
   TSkalcvfurw *kalcvfurw_ps = govprob_ps->kalcvfurw_ps;
   TSdvector *ylhtran_dv = kalcvfurw_ps->ylhtran_dv;
   TSdmatrix *Xrhtran_dm = kalcvfurw_ps->Xrhtran_dm;  //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.
               //For the benchmark case, 6 variables for Xrhtran_dm in the order of pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and constant terms.
   int nrhvars = Xrhtran_dm->nrows;
   TSdvector *expinf_dv = govprob_ps->expinf_dv;  //Will be updated by govprob_ps->govestimator() below.
   TSdvector *upred_dv;


   if (pointflag & POINTEST) {
      upred_dv = constmodpars_ps->upred_dv;   //Will be updated below.

      //------- Getting log LH or log Posterior first, NOT its negative value.  See (20a.4) on p.20a.
//      if (!constmodpars_ps->nlagsinflation && (constmodpars_ps->nlagstrue==1) && !constmodpars_ps->govmodel) {
      if (constmodpars_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0 | NGO_CGP2TP1TI0)) {
         //indxWhichModel: government's model (0: Classical direction or regresssion; 1: Keynesian direction).
         //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.
         //If indxWhichModel == 1 (classical), 6 variables for Xrhtran_dm are in the order of pi_t, pi_{t-1}, u_{t-1}, pi_{t-2}, u_{t-2}, and constant terms.
         //If indxWhichModel == 2 (Keynesian), 6 variables for Xrhtran_dm are in the order of u_t, u_{t-1}, pi_{t-1}, u_{t-2}, pi_{t-2}, and constant terms.
         fss = constmodpars_ps->fss;
         //--- Converting free paramaters to model (meaningful) parameters.
         CopyVector0(constmodpars_ps->x_dv, xchange_dv);
         ftd_freepars2modpars((void *)constmodpars_ps, POINTEST);

         //--- $$$WARNING$$$: should be always commented out.  Scaling down V to get a self-confirming equilibrium.  ONLY used to see if there is a self-confirming equilibrium.
         // if (constmodpars_ps->scl_epsilon != 1.0) {
         //    ScalarTimesMatrix(constmodpars_ps->V_dm, constmodpars_ps->scl_epsilon, constmodpars_ps->V_dm, 0.0);
         //    ftd_modpars2freepars((void *)constmodpars_ps, POINTEST);   //Converting back to x_dv to be consistent with the changed V.
         //    CopyVector0(xchange_dv, constmodpars_ps->x_dv);
         // }
         //--- Normalization.
         ftd_normalizationforgov((void *)constmodpars_ps, POINTEST);
         //--- Solving the government's problem for expected inflation and government's optimization decision for x_t (which controls at least part of the inflation process).
         govprob_ps->govestimator(govprob_ps, (void *)constmodpars_ps, POINTEST);

         //=== Getting logLH or logPosterior.
         logval = fss*(NEGHALFLOG2PI + 0.5*(log(constmodpars_ps->gzeta1) + log(constmodpars_ps->gzeta2)));
         if (constmodpars_ps->indxWhichModel & (CGP2TP1TI0 | NGO_CGP2TP1TI0))  //NGO_CGP2TP1TI0 is added 03/17/05.
         {
            for (workdterm=0.0, ti=0; ti<fss; ti++) {
               //--- Classical regression: ylhtran_dv is u_t; 6 variables for Xrhtran_dm in the order of pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and Constant.
               //--- Note that this order of lagged variables IS different from that in the government's gensys problem.
               z1t = ylhtran_dv->v[ti] - (upred_dv->v[ti] =  xchange_dv->v[0] + constmodpars_ps->gtau1*Xrhtran_dm->M[mos(2,ti,nrhvars)]
                     + constmodpars_ps->gtheta0*(z2t=(Xrhtran_dm->M[mos(0,ti,nrhvars)] - expinf_dv->v[ti]))
                     + constmodpars_ps->gtheta1*(Xrhtran_dm->M[mos(1,ti,nrhvars)] - ((ti<1) ? Xrhtran_dm->M[mos(1,ti,nrhvars)] : expinf_dv->v[ti-1])));
                                        //xchange_dv->v[0] is vstar on p.20a.
               workdterm += constmodpars_ps->gzeta1*square(z1t) + constmodpars_ps->gzeta2*square(z2t);
            }
         }
         else if (constmodpars_ps->indxWhichModel & KGP2TP1TI0)
         {
            for (workdterm=0.0, ti=0; ti<fss; ti++) {
               //--- Keynesian regression: ylhtran_dv is pi_t; 6 variables for Xrhtran_dm in the order of u_t, u_{t-1}, pi_{t-1}, u_{t-2}, pi_{t-2}, and Constant.
               //--- Note that this order of lagged variables IS different from that in the government's gensys problem.
               z1t = Xrhtran_dm->M[mos(0,ti,nrhvars)] - (upred_dv->v[ti] =  xchange_dv->v[0] + constmodpars_ps->gtau1*Xrhtran_dm->M[mos(1,ti,nrhvars)]
                     + constmodpars_ps->gtheta0*(z2t=(ylhtran_dv->v[ti] - expinf_dv->v[ti]))
                     + constmodpars_ps->gtheta1*(Xrhtran_dm->M[mos(2,ti,nrhvars)] - ((ti<1) ? Xrhtran_dm->M[mos(2,ti,nrhvars)] : expinf_dv->v[ti-1])));
                                        //xchange_dv->v[0] is vstar on p.20a.
               workdterm += constmodpars_ps->gzeta1*square(z1t) + constmodpars_ps->gzeta2*square(z2t);
            }
         }
         else
         {
            printf(".../swz_comfuns.c/value_logpostkernal_refresh(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d or %d!\n", CGP2TP1TI0, KGP2TP1TI0, NGO_CGP2TP1TI0);
            exit(EXIT_FAILURE);
         }
         upred_dv->flag = V_DEF;
         logval -= 0.5*workdterm;
                                               //Done with all work* arrays.
      }
//      else if (constmodpars_ps->indxWhichModel & NGO_CGP2TP1TI0)
//      {
//      }
      else {
         printf(".../swz_comfuns.c/value_logpostkernal_refresh(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d or %d!\n", CGP2TP1TI0, KGP2TP1TI0, NGO_CGP2TP1TI0);  //NGO_CGP2TP1TI0 is added 03/17/05.
         exit(EXIT_FAILURE);
      }
//      else  fn_DisplayError(".../value_logpostkernal_refresh():  Have not got time to handle the case for nlagstrue != 1 or nlagsinflation !=0 or Keyesian model (govmodel != 0)");

      if (!constmodpars_ps->indxMLE) {
         //Computes log posterior kernal with an informative prior, NOT its negative value yet.
         constmodpars_ps->LHvalatpostpeak = logval;
         logval += value_logpriorpdf(constmodpars_ps, (struct TSconstgibbs_tag *)NULL, POINTEST);
      }
   }
   else {
      upred_dv = constgibbs_ps->upred_dv;        //Will be updated below.

      //------- Getting log LH or log Posterior first, NOT its negative value.  See (20a.4) on p.20a.
//      if (!constgibbs_ps->nlagsinflation && (constgibbs_ps->nlagstrue==1) && !constgibbs_ps->govmodel) {
      if (constgibbs_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0)) {
         //indxWhichModel: government's model (0: Classical direction or regresssion; 1: Keynesian direction).
         //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.
         //6 variables for Xrhtran_dm in the order of pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and constant terms.
         fss = constgibbs_ps->fss;
         //--- Converting free paramaters to model (meaningful) parameters.
         CopyVector0(constgibbs_ps->x_dv, xchange_dv);
         ftd_freepars2modpars((void *)constgibbs_ps, GIBBSDRAWS);
         //--- $$$WARNING$$$: should be always commented out.  Scaling down V to get a self-confirming equilibrium.  ONLY used to see if there is a self-confirming equilibrium.
         // ScalarTimesMatrix(constgibbs_ps->V_dm, constgibbs_ps->scl_epsilon, constgibbs_ps->V_dm, 0.0);
         // ftd_modpars2freepars((void *)constgibbs_ps, GIBBSDRAWS);   //Converting back to x_dv to be consistent with the changed V.
         // CopyVector0(xchange_dv, constgibbs_ps->x_dv);
         //--- Normalization.
         ftd_normalizationforgov((void *)constgibbs_ps, GIBBSDRAWS);
         //--- Solving the government's problem for expected inflation and government's optimization decision for x_t (which controls at least part of the inflation process).
         govprob_ps->govestimator(govprob_ps, (void *)constgibbs_ps, GIBBSDRAWS);

         //=== Getting logLH or logPosterior.
         logval = fss*(NEGHALFLOG2PI + 0.5*(log(constgibbs_ps->gzeta1) + log(constgibbs_ps->gzeta2)));
         if (constgibbs_ps->indxWhichModel & CGP2TP1TI0) {
            for (workdterm=0.0, ti=0; ti<fss; ti++) {
               //--- Classical regression: ylhtran_dv is u_t; 6 variables for Xrhtran_dm in the order of pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and Constant.
               //--- Note that this order of lagged variables IS different from that in the government's gensys problem.
               z1t = ylhtran_dv->v[ti] - (upred_dv->v[ti] = xchange_dv->v[0] + constgibbs_ps->gtau1*Xrhtran_dm->M[mos(2,ti,nrhvars)]
                     + constgibbs_ps->gtheta0*(z2t=(Xrhtran_dm->M[mos(0,ti,nrhvars)] - expinf_dv->v[ti]))
                     + constgibbs_ps->gtheta1*(Xrhtran_dm->M[mos(1,ti,nrhvars)] - ((ti<1) ? Xrhtran_dm->M[mos(1,ti,nrhvars)] : expinf_dv->v[ti-1])));
                                        //xchange_dv->v[0] is vstar on p.20a.
               workdterm += constgibbs_ps->gzeta1*square(z1t) + constgibbs_ps->gzeta2*square(z2t);
            }
         }
         else if (constgibbs_ps->indxWhichModel & KGP2TP1TI0) {
            for (workdterm=0.0, ti=0; ti<fss; ti++) {
               //--- Keynesian regression: ylhtran_dv is pi_t; 6 variables for Xrhtran_dm in the order of u_t, u_{t-1}, pi_{t-1}, u_{t-2}, pi_{t-2}, and Constant.
               //--- Note that this order of lagged variables IS different from that in the government's gensys problem.
               z1t = Xrhtran_dm->M[mos(0,ti,nrhvars)] - (upred_dv->v[ti] =  xchange_dv->v[0] + constgibbs_ps->gtau1*Xrhtran_dm->M[mos(1,ti,nrhvars)]
                     + constgibbs_ps->gtheta0*(z2t=(ylhtran_dv->v[ti] - expinf_dv->v[ti]))
                     + constgibbs_ps->gtheta1*(Xrhtran_dm->M[mos(2,ti,nrhvars)] - ((ti<1) ? Xrhtran_dm->M[mos(2,ti,nrhvars)] : expinf_dv->v[ti-1])));
                                        //xchange_dv->v[0] is vstar on p.20a.
               workdterm += constgibbs_ps->gzeta1*square(z1t) + constgibbs_ps->gzeta2*square(z2t);
            }
         }
         upred_dv->flag = V_DEF;
         logval -= 0.5*workdterm;
                                               //Done with all work* arrays.
      }
      else {
         printf(".../swz_comfuns.c/value_logpostkernal_refresh(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d!", CGP2TP1TI0, KGP2TP1TI0);
         exit(EXIT_FAILURE);
      }
//      else  fn_DisplayError(".../value_logpostkernal_refresh():  Have not got time to handle the case for nlagstrue != 1 or nlagsinflation !=0 or Keyesian model (govmodel != 0)");

      //Computes log posterior kernal with an informative prior, NOT its negative value yet.
      logval += value_logpriorpdf(constmodpars_ps, constgibbs_ps, GIBBSDRAWS);
   }

   return (logval);
}
#undef NEGHALFLOG2PI


//------------------------------- 03/17/05 Revision.
// (1) Copies xchange_dv to constmodpars_ps->ngo_x_dv.
// (2) Refreshes all the parameters in constmodpars_ps, using ngo_freepars2modpars().
// (3) Computes the log value of LHvalue.  See pp.42-43.
//-------------------------------
#define NEGHALFLOG2PI  (-0.91893853320467)   //-0.5*log(2*pi)
double ngo_value_loglh_refresh(struct TSconstmodpars_tag *constmodpars_ps, TSgovprob *govprob_ps, TSdvector *xchange_dv)
{
   int ti;
   int fss;  //Effective sample size.
   double logval, acumnegterm1, acumterm2;
   //--- Initialization. The order matters!
   TSkalcvfurw *kalcvfurw_ps = govprob_ps->kalcvfurw_ps;
   TSdvector *ylhtran_dv = kalcvfurw_ps->ylhtran_dv;
   TSdmatrix *Xrhtran_dm = kalcvfurw_ps->Xrhtran_dm;  //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.
               //For the benchmark case, 6 variables for Xrhtran_dm in the order of pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and constant terms.
   int nrhvars = Xrhtran_dm->nrows;
   TSdvector *expinf_dv = govprob_ps->expinf_dv;  //Will be updated by govprob_ps->govestimator() below.
   TSdvector *upred_dv;
   //--- For Kalman filter.
   int kx = kalcvfurw_ps->kx;
   double workd, workdenominv;
   double sigmasq_fix = kalcvfurw_ps->sigmasq;
   TSdmatrix *V_dm = kalcvfurw_ps->V_dm;
   TSdmatrix *Zpredtran_dm = kalcvfurw_ps->Zpredtran_dm;
   TSdcell *Ppred_dc = kalcvfurw_ps->Ppred_dc;
   //=== For Kalman filter as well.
   TSdvector *workkxby1_dv = NULL;  //kx-by-1.
   TSdmatrix *workkxbykx_dm = NULL;  //kx-by-kx symmetric and positive positive.
   //******* WARNING: Some dangerous pointer movement to gain efficiency *******
   TSdvector xt_sdv;
   TSdvector zbefore_sdv;
   TSdvector zafter_sdv;

   //=== Memory allocation.
   workkxby1_dv = CreateVector_lf(kx);
   workkxbykx_dm = CreateMatrix_lf(kx, kx);
   //---
   xt_sdv.n = kx;
   xt_sdv.flag = V_DEF;
   zbefore_sdv.n = kx;
   zbefore_sdv.flag = V_DEF;
   zafter_sdv.n = kx;
   zafter_sdv.flag = V_DEF;


   upred_dv = constmodpars_ps->upred_dv;   //Will be updated below.
   logval = (fss = constmodpars_ps->fss)*NEGHALFLOG2PI;  //-0.5*fss*log(2*pi);
   acumnegterm1 = 0.0;
   acumterm2 = 0.0;

   //========== Getting log LHfirst, NOT its negative value.  See pp.42-43. ==========
   //With indxWhichModel == 1 (classical), 6 variables for Xrhtran_dm are in the order of pi_t, pi_{t-1}, u_{t-1}, pi_{t-2}, u_{t-2}, and constant terms.
   //--- Converting free paramaters to model (meaningful) parameters.
   CopyVector0(constmodpars_ps->ngo_x_dv, xchange_dv);
   ngo_freepars2modpars((void *)constmodpars_ps, POINTEST);
   sigmasq_fix = constmodpars_ps->ngo_sigmasq;
   //--- Refreshing kalcvfurw_ps->P10_dm and kalcvfurw_ps->V_dm.
   CopyMatrix0(kalcvfurw_ps->P10_dm, constmodpars_ps->P10_dm);
   CopyMatrix0(V_dm, constmodpars_ps->V_dm);


   //------- The first period (ti=0). -------
   zbefore_sdv.v = kalcvfurw_ps->z10_dv->v;
   zafter_sdv.v = Zpredtran_dm->M;
   xt_sdv.v = Xrhtran_dm->M;
   //--- Getting the current time-varying coefficient and its co-variance matrix.
   workd = ylhtran_dv->v[0] - VectorDotVector(&xt_sdv, &zbefore_sdv);   //y_t - x_t'*z_{t-1}.
   SymmatrixTimesVector(workkxby1_dv, kalcvfurw_ps->P10_dm, &xt_sdv, 1.0, 0.0);   //P_{t|t-1} x_t;
   workdenominv = 1.0/(sigmasq_fix + VectorDotVector(&xt_sdv, workkxby1_dv));   //1/[sigma^2 + x_t' P_{t|t-1} x_t]
   //--- Forming the LH.
   acumnegterm1 += log(workdenominv);  //(sum of) -log[sigma^2 + x_t' P_{t|t-1} x_t];
   acumterm2 += square(workd)*workdenominv;  //(sum of) (y_t - x_t'*z_{t-1})^2 / [sigma^2 + x_t' P_{t|t-1} x_t];
   //--- Updating z_{t+1|t}.
   CopyVector0(&zafter_sdv, &zbefore_sdv);
   VectorPlusMinusVectorUpdate(&zafter_sdv, workkxby1_dv, workd*workdenominv);  //z_{t+1|t} = z_{t|t-1} + P_{t|t-1} x_t [y_t - x_t'*z_{t-1}] / [sigma^2 + x_t' P_{t|t-1} x_t];
   //--- Updating P_{t+1|t}.
   CopyMatrix0(workkxbykx_dm, V_dm);
   VectorTimesSelf(workkxbykx_dm, workkxby1_dv, -workdenominv, 1.0, (V_dm->flag & M_SU) ? 'U' : 'L');
                                     // - P_{t|t-1}*x_t * xt'*P_{t|t-1} / [sigma^2 + x_t' P_{t|t-1} x_t] + V;
   MatrixPlusMatrix(Ppred_dc->C[0], kalcvfurw_ps->P10_dm, workkxbykx_dm);
                                     //P_{t|t-1} - P_{t|t-1}*x_t * xt'*P_{t|t-1} / [sigma^2 + x_t' P_{t|t-1} x_t] + V;
   Ppred_dc->C[0]->flag = M_GE | M_SU | M_SL;   //This is necessary because if P10_dm is initialized as diagonal, it will have M_GE | M_SU | M_SL | M_UT | M_LT,
                                          //  which is no longer true for workkxbykx_dm and therefore gives Ppred_dc->C[0] with M_GE only as a result of MatrixPlusMatrix().
                            //Done with all work* arrays.

   //------- The rest of the periods (ti=1:T-1). -------
   for (ti=1; ti<fss; ti++) {
      //NOTE: ti=0 has been taken care of outside of this loop.
      zbefore_sdv.v = Zpredtran_dm->M + (ti-1)*kx;
      zafter_sdv.v = Zpredtran_dm->M + ti*kx;
      xt_sdv.v = Xrhtran_dm->M + ti*kx;
      //--- Getting the current time-varying coefficient and its co-variance matrix.
      workd = ylhtran_dv->v[ti] - VectorDotVector(&xt_sdv, &zbefore_sdv);   //y_t - x_t'*z_{t-1}.
      SymmatrixTimesVector(workkxby1_dv, Ppred_dc->C[ti-1], &xt_sdv, 1.0, 0.0);   //P_{t|t-1} x_t;
      workdenominv = 1.0/(sigmasq_fix + VectorDotVector(&xt_sdv, workkxby1_dv));   //1/[sigma^2 + x_t' P_{t|t-1} x_t]
      //--- Forming the LH.
      acumnegterm1 += log(workdenominv);  //(sum of) -log[sigma^2 + x_t' P_{t|t-1} x_t];
      acumterm2 += square(workd)*workdenominv;  //(sum of) (y_t - x_t'*z_{t-1})^2 / [sigma^2 + x_t' P_{t|t-1} x_t];
      //--- Updating z_{t+1|t}.
      CopyVector0(&zafter_sdv, &zbefore_sdv);
      VectorPlusMinusVectorUpdate(&zafter_sdv, workkxby1_dv, workd*workdenominv);  //z_{t+1|t} = z_{t|t-1} + P_{t|t-1} x_t [y_t - x_t'*z_{t-1}] / [sigma^2 + x_t' P_{t|t-1} x_t];
      //--- Updating P_{t+1|t}.
      CopyMatrix0(workkxbykx_dm, V_dm);
      VectorTimesSelf(workkxbykx_dm, workkxby1_dv, -workdenominv, 1.0, (V_dm->flag & M_SU) ? 'U' : 'L');
                                        // - P_{t|t-1}*x_t * xt'*P_{t|t-1} / [sigma^2 + x_t' P_{t|t-1} x_t] + V;
      MatrixPlusMatrix(Ppred_dc->C[ti], Ppred_dc->C[ti-1], workkxbykx_dm);
                                        //P_{t|t-1} - P_{t|t-1}*x_t * xt'*P_{t|t-1} / [sigma^2 + x_t' P_{t|t-1} x_t] + V;
      Ppred_dc->C[ti]->flag = M_GE | M_SU | M_SL;   //This is necessary because if P10_dm is initialized as diagonal, it will have M_GE | M_SU | M_SL | M_UT | M_LT,
                                             //  which is no longer true for workkxbykx_dm and therefore gives Ppred_dc->C[0] with M_GE only as a result of MatrixPlusMatrix().
                               //Done with all work* arrays.
   }
   CopyVector0(kalcvfurw_ps->zupdate_dv, &zafter_sdv);
   Zpredtran_dm->flag = M_GE;

   //=== Finishing off with the logLH.  See p.43.
   logval += 0.5*(acumnegterm1 - acumterm2);

   DestroyVector_lf(workkxby1_dv);
   DestroyMatrix_lf(workkxbykx_dm);

   return (logval);
}
#undef NEGHALFLOG2PI



//-----------------------
// Normalization.
//-----------------------
void ftd_normalizationforgov(void *constmod_ps, int pointflag)
{
   struct TSconstmodpars_tag *constmodpars_ps;
   struct TSconstgibbs_tag *constgibbs_ps;


   //------- Sargent-William normalization. -------
   if (pointflag & POINTEST) {
      constmodpars_ps = (struct TSconstmodpars_tag *)constmod_ps;
      constmodpars_ps->gzeta = constmodpars_ps->gzeta1;  //Normalization.
//      if (constmodpars_ps->indxWhichModel & CGP2TP1TI0)  constmodpars_ps->gzeta = constmodpars_ps->gzeta1;  //Normalization.
//      else if (constmodpars_ps->indxWhichModel & KGP2TP1TI0)  constmodpars_ps->gzeta = constmodpars_ps->gzeta1/square(constmodpars_ps->gtheta0);  //Normalization.
//      else {
//         printf(".../swz_comfuns.c/ftd_normalizationforgov(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d!", CGP2TP1TI0, KGP2TP1TI0);
//         exit(EXIT_FAILURE);
//      }
   }
   else {
      //For GIBBSDRAWS (Gibbs draws).
      constgibbs_ps = (struct TSconstgibbs_tag *)constmod_ps;
      constgibbs_ps->gzeta = constgibbs_ps->gzeta1;  //Normalization.
//      if (constgibbs_ps->indxWhichModel & CGP2TP1TI0)  constgibbs_ps->gzeta = constgibbs_ps->gzeta1;  //Normalization.
//      else if (constgibbs_ps->indxWhichModel & KGP2TP1TI0)  constgibbs_ps->gzeta = constgibbs_ps->gzeta1/square(constgibbs_ps->gtheta0);  //Normalization.
//      else {
//         printf(".../swz_comfuns.c/ftd_normalizationforgov(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d!", CGP2TP1TI0, KGP2TP1TI0);
//         exit(EXIT_FAILURE);
//      }
   }
}




//************************************************************
// II. For the minimization (estimation) problem only.
//************************************************************
/**
void ftd_freepars2modpars(TSconstmodpars *constmodpars_ps)
{
   //Copying values of the vectorized free (minimization) parameters to regular model parameters.
   int _i, _j;
   int cnt;
   int kx = constmodpars_ps->kx;
   double workd;
   TSdvector *x_dv = constmodpars_ps->x_dv;
   //===
   TSdmatrix *cholMkbyk_dm = CreateConstantMatrix_lf(kx,kx,0.0);
   cholMkbyk_dm->flag = M_GE | M_UT;  //Will become upper triangular and general as well.

   if (!x_dv->flag)  fn_DisplayError(".../freepars2modpars(): x_dv must have legal values");

   if (!constmodpars_ps->nlagsinflation && (constmodpars_ps->nlagstrue==1) && !constmodpars_ps->govmodel) {
      //No lag and thus no rho1 and rho2 in the true inflation process to be estimated.  See p.20a.
      //6+kx+(kx+1)*kx variables in the order of: //vstar, theta0, theta1, tau1, xi1 (sqrt(zeta1)), xi2 (sqrt(zeta2)), z10, P10, and V.
      constmodpars_ps->gtau1 = x_dv->v[3];
      if (fabs(workd=1.0-constmodpars_ps->gtau1)>MACHINEZERO)  constmodpars_ps->ustar = x_dv->v[0]/workd;
      else {
         //Sets ustar to any arbitrary number and x_dv->v[0] to zero to make it compatible.
         constmodpars_ps->ustar = 5.0;
         x_dv->v[0] = 0.0;
      }
      // else if (fabs(x_dv->v[0])<MACHINEZERO)  constmodpars_ps->ustar = 5.0;   //Sets ustar to any arbitrary number.
      // else  fn_DisplayError(".../ftd_freepars2modpars(): 1.0-gtau1 is zero but vstar is not, thus no solution for ustar (see p.20a)");
                                             //Done with all work* arrays.
      constmodpars_ps->gtheta0 = x_dv->v[1];
      constmodpars_ps->gtheta1 = x_dv->v[2];
      constmodpars_ps->gzeta1 = square(x_dv->v[4]);
      constmodpars_ps->gzeta2 = square(x_dv->v[5]);
      //+
      for (_i=0; _i<kx; _i++)  constmodpars_ps->z10_dv->v[_i] = x_dv->v[6+_i];
      for (cnt=0, _j=0; _j<kx; _j++)
         for (_i=0; _i<=_j; _i++)
            cholMkbyk_dm->M[mos(_i,_j,kx)] = x_dv->v[6+kx+cnt++];
      MatrixTimesSelf(constmodpars_ps->P10_dm, 'U', cholMkbyk_dm, 'T', 1.0, 0.0);  //P = C'*C;
      for (_j=0; _j<kx; _j++)
         for (_i=0; _i<=_j; _i++)
            cholMkbyk_dm->M[mos(_i,_j,kx)] = x_dv->v[6+kx+cnt++];
      MatrixTimesSelf(constmodpars_ps->V_dm, 'U', cholMkbyk_dm, 'T', 1.0, 0.0);  //V = C'*C;
   }
   else  fn_DisplayError(".../ftd_freepars2modpars(): Have not got time to do the case for nlagstrue != 1 or nlagsinflation !=0 or govmodel != 0");

   //===
   DestroyMatrix_lf(cholMkbyk_dm);
}
//---
void ftd_modpars2freepars(TSconstmodpars *constmodpars_ps)
{
   //Copying values of regular model parameters to the vectorized free (minimization) parameters.
   int _i, _j;
   int cnt;
   int kx = constmodpars_ps->kx;
   TSdvector *x_dv = constmodpars_ps->x_dv;
   //===
   TSdmatrix *cholMkbyk_dm = CreateMatrix_lf(kx,kx);   //Will become upper triangular.


   if (!constmodpars_ps->nlagsinflation && (constmodpars_ps->nlagstrue==1) && !constmodpars_ps->govmodel) {
      //No lag and thus no rho1 and rho2 in the true inflation process to be estimated.  See p.20a.
      //6+kx+(kx+1)*kx variables in the order of: //vstar, theta0, theta1, tau1, xi1 (sqrt(zeta1)), xi2 (sqrt(zeta2)), z10, P10, and V.
      x_dv->v[0] = constmodpars_ps->ustar*(1.0-constmodpars_ps->gtau1);   //vstar in the notation of p.20a.
      x_dv->v[1] = constmodpars_ps->gtheta0;
      x_dv->v[2] = constmodpars_ps->gtheta1;
      x_dv->v[3] = constmodpars_ps->gtau1;
      x_dv->v[4] = sqrt(constmodpars_ps->gzeta1);
      x_dv->v[5] = sqrt(constmodpars_ps->gzeta2);
      //+
      for (_i=0; _i<kx; _i++)  x_dv->v[6+_i] = constmodpars_ps->z10_dv->v[_i];
      if (chol(cholMkbyk_dm, constmodpars_ps->P10_dm, 'U'))  fn_DisplayError(".../modpars2freepars(): Choleski decomposition of P10_dm fails when calling chol()");
      for (cnt=0, _j=0; _j<kx; _j++)
         for (_i=0; _i<=_j; _i++)
            x_dv->v[6+kx+cnt++] = cholMkbyk_dm->M[mos(_i,_j,kx)];
      if (chol(cholMkbyk_dm, constmodpars_ps->V_dm, 'U'))  fn_DisplayError(".../modpars2freepars(): Choleski decomposition of V_dm fails when calling chol()");
      for (_j=0; _j<kx; _j++)
         for (_i=0; _i<=_j; _i++)
            x_dv->v[6+kx+cnt++] = cholMkbyk_dm->M[mos(_i,_j,kx)];
      x_dv->flag = V_DEF;
   }
   else  fn_DisplayError(".../ftd_modpars2freepars(): Have not got time to do the case for nlagstrue != 1 or nlagsinflation !=0 or govmodel != 0");

   //===
   DestroyMatrix_lf(cholMkbyk_dm);
}


//-----------------------
// Normalization.
//-----------------------
ftd_normalizationforgov(TSconstmodpars *constmodpars_ps)
{
   //--- Sargent-William normalization.
   constmodpars_ps->gzeta = constmodpars_ps->gzeta1;  //Normalization.
}
/**/


//-----------------------
// Sets (manually) the values of prior hyperparameters.
//-----------------------
#define NEGHALFLOG2PI (-0.91893853320467)   //-0.5*log(2*pi).
void ftd_priorsettings(TSconstmodpars *constmodpars_ps)    //Added 03/17/05.
{
   int _m, _i, cnter, kx;


   //Hyperparameter settings for the prior.  See pp. 20c-d.
//   if (!constmodpars_ps->nlagsinflation && constmodpars_ps->nlagstrue==1) {
   if (constmodpars_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0 | NGO_CGP2TP1TI0))  //NGO_CGP2TP1TI0 is added 03/17/05.
   {
      //No lag and thus no rho1 and rho2 in the true inflation process to be estimated.
      //6+kx+(kx+1)*kx variables in the order of: //vstar, theta0, theta1, tau1, xi1 (sqrt(zeta1)), xi2 (sqrt(zeta2)), z10, P10, and V.
      constmodpars_ps->vstar_priormean = 0.12; //0.6;  //0.12;
      constmodpars_ps->vstar_priorvar =  constmodpars_ps->tc1*square(0.06); //(0.3);  //(0.06);
      constmodpars_ps->gtheta0_priormean = -0.20;  //-1.0;  //-0.20;
      constmodpars_ps->gtheta0_priorvar = constmodpars_ps->tc1*square(0.10);  //(0.5);  //(0.10);
      constmodpars_ps->gtheta1_priormean = -0.16; //-1.0; //-0.16;
      constmodpars_ps->gtheta1_priorvar = constmodpars_ps->tc1*square(0.08); //(0.5);  //(0.08);
      constmodpars_ps->gtau1_priormean = 0.98; //0.9;  //0.98;
      constmodpars_ps->gtau1_priorvar = constmodpars_ps->tc1*square(0.01);  //(0.05);  //(0.01);
      //+
      constmodpars_ps->corr_theta0_theta1 = 0.8;  //Correlation between theta0 and theta1.
      constmodpars_ps->corr_vstar_theta0 = 0.0;  //Correlation between vstar and theta0.
      constmodpars_ps->corr_vstar_theta1 = 0.0;  //Correlation between vstar and theta1.
      constmodpars_ps->corr_vstar_tau1 = -0.9;  //Correlation between vstar and tau1.
      constmodpars_ps->corr_theta0_tau1 = 0.0;  //-0.5 will make the matrix negative definite;  //Correlation between theta0 and tau1.
      constmodpars_ps->corr_theta1_tau1 = 0.0;  //0.5 may make the matrix negative definte;  //Correlation between theta1 and tau1.
      //+
      constmodpars_ps->gzeta1p2_prioralpha = 4.0;   //The parameter alpha in a Gamma prior.
      constmodpars_ps->gzeta1p2_priorbeta = constmodpars_ps->tc2*12.5;    //The parameter beta in a Gamma prior.
      //+
      if (constmodpars_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0))  //Added 03/17/05.
      {
         constmodpars_ps->z10_priormean = 0.0;  //Mean for all elements in z_{1|0} or alpha^{\hat}_{1|0}.
         constmodpars_ps->z10_priorvar = constmodpars_ps->tc5*square(100.0);  //Variance for all elements in z_{1|0} or alpha^{\hat}_{1|0}.
         constmodpars_ps->Cpdiag_priormean = 0.0; //20.0;  //Mean for the diagonal of the upper Choleski decomposition of P_{1|0} = Cp'*Cp where Cp is upper triangular.
         constmodpars_ps->Cpoff_priormean = 0.0;  //Mean for the off-diagonal elements of the upper Choleski decomposition of P_{1|0} = Cp'*Cp where Cp is upper triangular.
         constmodpars_ps->Cpdiag_priorvar = constmodpars_ps->tc3*square(5.0); //(10.0);  //Variance for the diagonal of the upper Choleski decomposition of P_{1|0} = Cp'*Cp where Cp is upper triangular.
         constmodpars_ps->Cpoff_priorvar = constmodpars_ps->tc3*square(2.5);  //(10.0);  //Variance for the off-diagonal elements of the upper Choleski decomposition of P_{1|0} = Cp'*Cp where Cp is upper triangular.
         constmodpars_ps->Cvdiag_priormean = 0.0;   //10.0;  //Mean for the diagonal of the upper Choleski decomposition of V = Cv'*Cv where Cv is upper triangular.
         constmodpars_ps->Cvoff_priormean = 0.0;  //Mean for the off-diagonal elements of the upper Choleski decomposition of V = Cv'*Cv where Cv is upper triangular.
         constmodpars_ps->Cvdiag_priorvar = constmodpars_ps->tc4*square(5.0); //(5.0) //Variance for the diagonal of the upper Choleski decomposition of V = Cv'*Cv where Cv is upper triangular.
         constmodpars_ps->Cvoff_priorvar = constmodpars_ps->tc4*square(2.5);  //(5.0) //Variance for the off-diagonal elements of the upper Choleski decomposition of V = Cv'*Cv where Cv is upper triangular.
      }

      //------- Setting mean and covariance for the 1st set of parameters.
      if (!constmodpars_ps->Var1stSet_dm)  constmodpars_ps->Var1stSet_dm = CreateMatrix_lf(N1STSET,N1STSET);  //4: vstar, theta0, theta1, and tau1.
      else  fn_DisplayError(".../ftd_priorsettings():  Must sure that constmodpars_ps->Var1stSet_dm is created with correct dimension");
      constmodpars_ps->Var1stSet_dm->M[mos(0,0,N1STSET)] = constmodpars_ps->vstar_priorvar;
      constmodpars_ps->Var1stSet_dm->M[mos(1,1,N1STSET)] = constmodpars_ps->gtheta0_priorvar;
      constmodpars_ps->Var1stSet_dm->M[mos(2,2,N1STSET)] = constmodpars_ps->gtheta1_priorvar;
      constmodpars_ps->Var1stSet_dm->M[mos(3,3,N1STSET)] = constmodpars_ps->gtau1_priorvar;
      //+
      constmodpars_ps->Var1stSet_dm->M[mos(0,1,N1STSET)] = constmodpars_ps->corr_vstar_theta0*sqrt(constmodpars_ps->vstar_priorvar*constmodpars_ps->gtheta0_priorvar);
      constmodpars_ps->Var1stSet_dm->M[mos(0,2,N1STSET)] = constmodpars_ps->corr_vstar_theta1*sqrt(constmodpars_ps->vstar_priorvar*constmodpars_ps->gtheta1_priorvar);
      constmodpars_ps->Var1stSet_dm->M[mos(0,3,N1STSET)] = constmodpars_ps->corr_vstar_tau1*sqrt(constmodpars_ps->vstar_priorvar*constmodpars_ps->gtau1_priorvar);
      constmodpars_ps->Var1stSet_dm->M[mos(1,2,N1STSET)] = constmodpars_ps->corr_theta0_theta1*sqrt(constmodpars_ps->gtheta0_priorvar*constmodpars_ps->gtheta1_priorvar);
      constmodpars_ps->Var1stSet_dm->M[mos(1,3,N1STSET)] = constmodpars_ps->corr_theta0_tau1*sqrt(constmodpars_ps->gtheta0_priorvar*constmodpars_ps->gtau1_priorvar);
      constmodpars_ps->Var1stSet_dm->M[mos(2,3,N1STSET)] = constmodpars_ps->corr_theta1_tau1*sqrt(constmodpars_ps->gtheta1_priorvar*constmodpars_ps->gtau1_priorvar);
      constmodpars_ps->Var1stSet_dm->flag = M_SU;
      constmodpars_ps->cholVar1stSet_dm = CreateConstantMatrix_lf(N1STSET,N1STSET, 0.0);
      if (chol(constmodpars_ps->cholVar1stSet_dm, constmodpars_ps->Var1stSet_dm, 'U'))  fn_DisplayError(".../ftd_priorsettings(): Make sure Var1stSet_dm is positive definite");
      constmodpars_ps->term1st1 = NEGHALFLOG2PI*N1STSET - tracelogfabs(constmodpars_ps->cholVar1stSet_dm);
      //+
      if (!constmodpars_ps->mean1stset_dv) constmodpars_ps->mean1stset_dv = CreateVector_lf(N1STSET);
      else  fn_DisplayError(".../ftd_priorsettings():  Must sure that constmodpars_ps->mean1stset_dv is created with correct dimension");
      constmodpars_ps->mean1stset_dv->v[0] = constmodpars_ps->vstar_priormean;
      constmodpars_ps->mean1stset_dv->v[1] = constmodpars_ps->gtheta0_priormean;
      constmodpars_ps->mean1stset_dv->v[2] = constmodpars_ps->gtheta1_priormean;
      constmodpars_ps->mean1stset_dv->v[3] = constmodpars_ps->gtau1_priormean;
      constmodpars_ps->mean1stset_dv->flag = V_DEF;


      //------- Some terms for the 2nd and 3rd sets of parameters.   -------
      constmodpars_ps->term2nd1 = constmodpars_ps->term3rd1 = -fn_gammalog(constmodpars_ps->gzeta1p2_prioralpha) - constmodpars_ps->gzeta1p2_prioralpha*log(constmodpars_ps->gzeta1p2_priorbeta);



      if (constmodpars_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0))  //Added 03/17/05.
      {
         if (!constmodpars_ps->indxFixz10P10) {
            //------- Setting mean and covariance for the 4th set of parameters (coefficients on pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and Constant).
            if (!constmodpars_ps->mean4thset_dv) constmodpars_ps->mean4thset_dv = CreateConstantVector_lf(constmodpars_ps->kx, constmodpars_ps->z10_priormean);
            else  fn_DisplayError(".../ftd_priorsettings():  Must sure that constmodpars_ps->mean4thset_dv is created with correct dimension");
            constmodpars_ps->mean4thset_dv->flag = V_DEF;
            //+
            if (!constmodpars_ps->Var4thSet_dm) {
               constmodpars_ps->Var4thSet_dm = CreateMatrix_lf(constmodpars_ps->kx,constmodpars_ps->kx);
               InitializeDiagonalMatrix_lf(constmodpars_ps->Var4thSet_dm, constmodpars_ps->z10_priorvar);
            }
            else  fn_DisplayError(".../ftd_priorsettings():  Must sure that constmodpars_ps->Var4thSet_dm is created with correct dimension");
            constmodpars_ps->Var4thSet_dm->flag = M_GE | M_SU;
            constmodpars_ps->cholVar4thSet_dm = CreateConstantMatrix_lf(constmodpars_ps->kx,constmodpars_ps->kx, 0.0);
            if (chol(constmodpars_ps->cholVar4thSet_dm, constmodpars_ps->Var4thSet_dm, 'U'))  fn_DisplayError(".../ftd_priorsettings(): Make sure Var4thSet_dm is positive definite");
            constmodpars_ps->term4th1 = NEGHALFLOG2PI*constmodpars_ps->kx - tracelogfabs(constmodpars_ps->cholVar4thSet_dm);


            //------- Setting mean and covariance for the 5th set of parameters (only upper triangular part of the Choleski decompostion of P_{1|0}).
            _m = (kx=constmodpars_ps->kx)*(constmodpars_ps->kx+1)/2;
            if (!constmodpars_ps->mean5thset_dv) {
               constmodpars_ps->mean5thset_dv = CreateConstantVector_lf(_m, constmodpars_ps->Cpoff_priormean);
               for (cnter=-1, _i=0; _i<kx; _i++) {
                  //Diagonal elements.
                  constmodpars_ps->mean5thset_dv->v[cnter += _i+1] = constmodpars_ps->Cpdiag_priormean;
               }
            }
            else  fn_DisplayError(".../ftd_priorsettings():  Must sure that constmodpars_ps->mean5thset_dv is created with correct dimension");
            constmodpars_ps->mean5thset_dv->flag = V_DEF;
            //+
            if (!constmodpars_ps->Var5thSet_dm) {
               constmodpars_ps->Var5thSet_dm = CreateMatrix_lf(_m, _m);
               InitializeDiagonalMatrix_lf(constmodpars_ps->Var5thSet_dm, constmodpars_ps->tc3*constmodpars_ps->Cpoff_priorvar);
               for (cnter=-1, _i=0; _i<kx; _i++) {
                  //Diagonal elements for the Choleski of P_{1|0}, NOT for constmodpars_ps->Var5thSet_dm itself.
                  cnter += _i+1;
                  constmodpars_ps->Var5thSet_dm->M[mos(cnter,cnter,_m)] = constmodpars_ps->tc3*constmodpars_ps->Cpdiag_priorvar;
               }
            }
            else  fn_DisplayError(".../ftd_priorsettings():  Must sure that constmodpars_ps->Var5thSet_dm is created with correct dimension");
            constmodpars_ps->Var5thSet_dm->flag = M_GE | M_SU;
            constmodpars_ps->cholVar5thSet_dm = CreateConstantMatrix_lf(_m, _m, 0.0);
            if (chol(constmodpars_ps->cholVar5thSet_dm, constmodpars_ps->Var5thSet_dm, 'U'))  fn_DisplayError(".../ftd_priorsettings(): Make sure Var5thSet_dm is positive definite");
            constmodpars_ps->term5th1 = NEGHALFLOG2PI*_m - tracelogfabs(constmodpars_ps->cholVar5thSet_dm);
         }
         if (constmodpars_ps->indxFixz10P10 == 2) {
            //------- Setting mean and covariance for the 5th set of parameters (only upper triangular part of the Choleski decompostion of P_{1|0}).
            _m = (kx=constmodpars_ps->kx)*(constmodpars_ps->kx+1)/2;
            if (!constmodpars_ps->mean5thset_dv) {
               constmodpars_ps->mean5thset_dv = CreateConstantVector_lf(_m, constmodpars_ps->Cpoff_priormean);
               for (cnter=-1, _i=0; _i<kx; _i++) {
                  //Diagonal elements.
                  constmodpars_ps->mean5thset_dv->v[cnter += _i+1] = constmodpars_ps->Cpdiag_priormean;
               }
            }
            else  fn_DisplayError(".../ftd_priorsettings():  Must sure that constmodpars_ps->mean5thset_dv is created with correct dimension");
            constmodpars_ps->mean5thset_dv->flag = V_DEF;
            //+
            if (!constmodpars_ps->Var5thSet_dm) {
               constmodpars_ps->Var5thSet_dm = CreateMatrix_lf(_m, _m);
               InitializeDiagonalMatrix_lf(constmodpars_ps->Var5thSet_dm, constmodpars_ps->tc3*constmodpars_ps->Cpoff_priorvar);
               for (cnter=-1, _i=0; _i<kx; _i++) {
                  //Diagonal elements for the Choleski of P_{1|0}, NOT for constmodpars_ps->Var5thSet_dm itself.
                  cnter += _i+1;
                  constmodpars_ps->Var5thSet_dm->M[mos(cnter,cnter,_m)] = constmodpars_ps->tc3*constmodpars_ps->Cpdiag_priorvar;
               }
            }
            else  fn_DisplayError(".../ftd_priorsettings():  Must sure that constmodpars_ps->Var5thSet_dm is created with correct dimension");
            constmodpars_ps->Var5thSet_dm->flag = M_GE | M_SU;
            constmodpars_ps->cholVar5thSet_dm = CreateConstantMatrix_lf(_m, _m, 0.0);
            if (chol(constmodpars_ps->cholVar5thSet_dm, constmodpars_ps->Var5thSet_dm, 'U'))  fn_DisplayError(".../ftd_priorsettings(): Make sure Var5thSet_dm is positive definite");
            constmodpars_ps->term5th1 = NEGHALFLOG2PI*_m - tracelogfabs(constmodpars_ps->cholVar5thSet_dm);
         }

         //------- Setting mean and covariance for the 6th set of parameters (only upper triangular part of the Choleski decompostion of P_{1|0}).
         _m = (kx=constmodpars_ps->kx)*(constmodpars_ps->kx+1)/2;
         if (!constmodpars_ps->mean6thset_dv) {
            constmodpars_ps->mean6thset_dv = CreateConstantVector_lf(_m, constmodpars_ps->Cvoff_priormean);
            for (cnter=-1, _i=0; _i<kx; _i++) {
               //Diagonal elements.
               constmodpars_ps->mean6thset_dv->v[cnter += _i+1] = constmodpars_ps->Cvdiag_priormean;
            }
         }
         else  fn_DisplayError(".../ftd_priorsettings():  Must sure that constmodpars_ps->mean6thset_dv is created with correct dimension");
         constmodpars_ps->mean6thset_dv->flag = V_DEF;
         //+
         if (!constmodpars_ps->Var6thSet_dm) {
            constmodpars_ps->Var6thSet_dm = CreateMatrix_lf(_m, _m);
            InitializeDiagonalMatrix_lf(constmodpars_ps->Var6thSet_dm, constmodpars_ps->tc4*constmodpars_ps->Cvoff_priorvar);
            for (cnter=-1, _i=0; _i<kx; _i++) {
               //Diagonal elements for the Choleski of V, NOT for constmodpars_ps->Var6thSet_dm itself.
               cnter += _i+1;
               constmodpars_ps->Var6thSet_dm->M[mos(cnter,cnter,_m)] = constmodpars_ps->tc4*constmodpars_ps->Cvdiag_priorvar;
            }
         }
         else  fn_DisplayError(".../ftd_priorsettings():  Must sure that constmodpars_ps->Var6thSet_dm is created with correct dimension");
         constmodpars_ps->Var6thSet_dm->flag = M_GE | M_SU;
         constmodpars_ps->cholVar6thSet_dm = CreateConstantMatrix_lf(_m, _m, 0.0);
         if (chol(constmodpars_ps->cholVar6thSet_dm, constmodpars_ps->Var6thSet_dm, 'U'))  fn_DisplayError(".../ftd_priorsettings(): Make sure Var6thSet_dm is positive definite");
         constmodpars_ps->term6th1 = NEGHALFLOG2PI*_m - tracelogfabs(constmodpars_ps->cholVar6thSet_dm);
      }
   }
   else if (constmodpars_ps->nlagstrue != 1)  fn_DisplayError(".../swz_comfuns.c/ftd_priorsettings(): Have not got time to work out the case for nlagstrue != 1");
   else if (constmodpars_ps->nlagsinflation > 0)  fn_DisplayError(".../swz_comfuns.c/ftd_priorsettings(): Have not got time to code up the case for nlagsinflation > 0");
   else {
      printf(".../swz_comfuns.c/ftd_priorsettings(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d!", CGP2TP1TI0, KGP2TP1TI0);
      exit(EXIT_FAILURE);
   }
}
#undef NEGHALFLOG2PI
/**
//---
double value_logpriorpdf(TSconstmodpars *constmodpars_ps)
{
   //Log prior pdf (not just kernal).  See p.20b.
   //Must be comformable to ftd_freepars2modpars().
   int cntloc, _m;
   double valuelogpdf;
   //---
   double gzeta1p2;
   double gzeta1p2_prioralpha;
   double gzeta1p2_priorbetainv;
   //---
   TSdvector *x_dv = constmodpars_ps->x_dv;
   TSdvector xmveset_sdv;
   TSdvector *mean1stset_dv;
   TSdmatrix *cholVar1stSet_dm;
   TSdvector *mean4thset_dv;
   TSdmatrix *cholVar4thSet_dm;
   TSdvector *mean5thset_dv;
   TSdmatrix *cholVar5thSet_dm;
   TSdvector *mean6thset_dv;
   TSdmatrix *cholVar6thSet_dm;
   //+
   TSdvector *work_dv = NULL;

   if (!constmodpars_ps->nlagsinflation && constmodpars_ps->nlagstrue==1) {
      //No lag and thus no rho1 and rho2 in the true inflation process to be estimated.
      //6+kx+(kx+1)*kx variables in the order of: //vstar, theta0, theta1, tau1, xi1 (sqrt(zeta1)), xi2 (sqrt(zeta2)), z10, P10, and V.

      //=== Prior settings.
      mean1stset_dv = constmodpars_ps->mean1stset_dv;
      cholVar1stSet_dm = constmodpars_ps->cholVar1stSet_dm;
      mean4thset_dv = constmodpars_ps->mean4thset_dv;
      cholVar4thSet_dm = constmodpars_ps->cholVar4thSet_dm;
      mean5thset_dv = constmodpars_ps->mean5thset_dv;
      cholVar5thSet_dm = constmodpars_ps->cholVar5thSet_dm;
      mean6thset_dv = constmodpars_ps->mean6thset_dv;
      cholVar6thSet_dm = constmodpars_ps->cholVar6thSet_dm;

      //--------------------------
      // The order of the following lines matters!
      //--------------------------
      valuelogpdf = 0.0;   //Initialization.
      xmveset_sdv.flag = V_DEF;
      //------- logpriorpdf for 1st set of parameters (vstar, theta0, theta1, and tau1). See p.20b. -------
      cntloc = 0;
      xmveset_sdv.n = N1STSET;
      xmveset_sdv.v = x_dv->v;
      work_dv = CreateVector_lf(N1STSET);
      //+
      VectorMinusVector(work_dv, &xmveset_sdv, mean1stset_dv);
      SolveTriSysVector(work_dv, cholVar1stSet_dm, work_dv, 'T', 'N');
      valuelogpdf += constmodpars_ps->term1st1 - 0.5*VectorDotVector(work_dv, work_dv);
      //+
      work_dv = DestroyVector_lf(work_dv);

      //------- logpriorpdf for 2nd and 3rd sets of parameters (i.e., zeta1 and zeta2). See p.20b. -------
      gzeta1p2_prioralpha = constmodpars_ps->gzeta1p2_prioralpha;
      gzeta1p2_priorbetainv = 1.0/constmodpars_ps->gzeta1p2_priorbeta;
      //
      gzeta1p2 = x_dv->v[cntloc += N1STSET];
      gzeta1p2 = square(gzeta1p2);
      valuelogpdf += constmodpars_ps->term2nd1 + (gzeta1p2_prioralpha-1.0)*log(gzeta1p2) - (gzeta1p2*gzeta1p2_priorbetainv);
      //+
      gzeta1p2 = x_dv->v[cntloc += N2NDSET];
      gzeta1p2 = square(gzeta1p2);
      valuelogpdf += constmodpars_ps->term3rd1 + (gzeta1p2_prioralpha-1.0)*log(gzeta1p2) - (gzeta1p2*gzeta1p2_priorbetainv);

      //------- logpriorpdf for 4th set of parameters (i.e., z_{1|0}, coefficients on pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and Constant). See p.20b. -------
      xmveset_sdv.n = constmodpars_ps->kx;
      xmveset_sdv.v = x_dv->v + (cntloc += N3RDSET);
      work_dv = CreateVector_lf(constmodpars_ps->kx);
      //+
      VectorMinusVector(work_dv, &xmveset_sdv, mean4thset_dv);
      SolveTriSysVector(work_dv, cholVar4thSet_dm, work_dv, 'T', 'N');
      valuelogpdf += constmodpars_ps->term4th1 - 0.5*VectorDotVector(work_dv, work_dv);
      //+
      work_dv = DestroyVector_lf(work_dv);

      //------- logpriorpdf for 5th set of parameters (i.e., only upper triangular part of the Choleski decompostion of P_{1|0}). See p.20b. -------
      xmveset_sdv.n = _m = constmodpars_ps->kx*(constmodpars_ps->kx+1)/2;
      xmveset_sdv.v = x_dv->v + (cntloc += constmodpars_ps->kx);
      work_dv = CreateVector_lf(_m);
      //+
      VectorMinusVector(work_dv, &xmveset_sdv, mean5thset_dv);
      SolveTriSysVector(work_dv, cholVar5thSet_dm, work_dv, 'T', 'N');
      valuelogpdf += constmodpars_ps->term5th1 - 0.5*VectorDotVector(work_dv, work_dv);
      //+
      work_dv = DestroyVector_lf(work_dv);

      //------- logpriorpdf for 6th set of parameters (i.e., only upper triangular part of the Choleski decompostion of V). See p.20b. -------
      xmveset_sdv.n = _m;
      xmveset_sdv.v = x_dv->v + (cntloc += _m);
      work_dv = CreateVector_lf(_m);
      //+
      VectorMinusVector(work_dv, &xmveset_sdv, mean6thset_dv);
      SolveTriSysVector(work_dv, cholVar6thSet_dm, work_dv, 'T', 'N');
      valuelogpdf +=  constmodpars_ps->term6th1 - 0.5*VectorDotVector(work_dv, work_dv);
      //+
      work_dv = DestroyVector_lf(work_dv);
   }
   else if (constmodpars_ps->nlagstrue != 1)  fn_DisplayError(".../ftd_logprior(): Have not got time to work out the case for nlagstrue != 1");
   else  fn_DisplayError(".../ftd_logprior(): Have not got time to code up the case for nlagsinflation > 0");

   //===
   DestroyVector_lf(work_dv);
   //
   return (valuelogpdf);
}
/**/

/**
double value_logpostkernal_refresh(TSconstmodpars *constmodpars_ps, TSgovprob *govprob_ps, TSdvector *xchange_dv)
{
   int ti;
   int fss;  //Effective sample size.
   double workdterm, z1t, z2t;
   double logval;
   //--- Initialization. The order matters!
   TSkalcvfurw *kalcvfurw_ps = govprob_ps->kalcvfurw_ps;
   TSdvector *ylhtran_dv = kalcvfurw_ps->ylhtran_dv;
   TSdmatrix *Xrhtran_dm = kalcvfurw_ps->Xrhtran_dm;  //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->govmodel.
   int nrhvars = Xrhtran_dm->nrows;
   TSdvector *expinf_dv = govprob_ps->expinf_dv;
   TSdvector *upred_dv = constmodpars_ps->upred_dv;



   //------- Getting log LH or log Posterior first, NOT its negative value.  See (20a.4) on p.20a.
   if (!constmodpars_ps->nlagsinflation && (constmodpars_ps->nlagstrue==1) && !constmodpars_ps->govmodel) {
      //govmodel: government's model (0: Classical direction or regresssion; 1: Keynesian direction).
      //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->govmodel.
      fss = constmodpars_ps->fss;
      //--- Converting free paramaters to model (meaningful) parameters.
      CopyVector0(constmodpars_ps->x_dv, xchange_dv);
      ftd_freepars2modpars((void *)constmodpars_ps, POINTEST);
      ftd_normalizationforgov((void *)constmodpars_ps, POINTEST);  //Normalization.
      //--- Solving the government's problem for expected inflation and government's optimization decision for x_t (which controls at least part of the inflation process).
      govprob_ps->govestimator(govprob_ps, (void *)constmodpars_ps, POINTEST);

      //=== Getting logLH or logPosterior.
      logval = fss*(NEGHALFLOG2PI + 0.5*(log(constmodpars_ps->gzeta1) + log(constmodpars_ps->gzeta2)));
      for (workdterm=0.0, ti=0; ti<fss; ti++) {
         //--- 6 variables for Xrhtran_dm in the order of pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and Constant.
         //--- Note that this order of lagged variables IS different from that in the government's gensys problem.
         z1t = ylhtran_dv->v[ti] - (upred_dv->v[ti] =  xchange_dv->v[0] + constmodpars_ps->gtau1*Xrhtran_dm->M[mos(2,ti,nrhvars)]
               + constmodpars_ps->gtheta0*(z2t=(Xrhtran_dm->M[mos(0,ti,nrhvars)] - expinf_dv->v[ti]))
               + constmodpars_ps->gtheta1*(Xrhtran_dm->M[mos(1,ti,nrhvars)] - ((ti<1) ? Xrhtran_dm->M[mos(1,ti,nrhvars)] : expinf_dv->v[ti-1])));
                                  //xchange_dv->v[0] is vstar on p.20a.
         workdterm += constmodpars_ps->gzeta1*square(z1t) + constmodpars_ps->gzeta2*square(z2t);
      }
      upred_dv->flag = V_DEF;
      logval -= 0.5*workdterm;
                                            //Done with all work* arrays.
   }
   else  fn_DisplayError(".../minneglogpost():  Have not got time to handle the case for nlagstrue != 1 or nlagsinflation !=0 or Keyesian model (govmodel != 0)");

   if (!constmodpars_ps->indxMLE) {
      //Computes log posterior kernal with an informative prior, NOT its negative value yet.
      logval += value_logpriorpdf(constmodpars_ps, (struct TSconstgibbs_tag *)NULL, POINTEST);
   }

   return (logval);
}
/**/




//************************************************************
// III. For the posterior probability (inferences).
//************************************************************
struct TSconstgibbs_tag *CreateTSconstgibbs(struct TSinput_tag *input_ps, struct TSconstmodpars_tag *constmodpars_ps)
{
   int fourdim;
   struct TSconstgibbs_tag *constgibbs_ps = tzMalloc(1, struct TSconstgibbs_tag);

   constgibbs_ps->indxWhichModel = input_ps->indxWhichModel;
   constgibbs_ps->indxFixz10P10 = input_ps->indxFixz10P10;
   constgibbs_ps->indxCnterFact = input_ps->indxCnterFact;
   constgibbs_ps->indxNoUpdating = input_ps->indxNoUpdating;
   constgibbs_ps->sclhistshocks_dv = CreateVector_lf(input_ps->sclhistshocks_dv->n);
   CopyVector0(constgibbs_ps->sclhistshocks_dv, input_ps->sclhistshocks_dv);
   constgibbs_ps->scl_epsilon = input_ps->scl_epsilon;

   //=== Fixed values used for evaluating the logLH or logPosterior that is formed according to nlagstrue.
//   constgibbs_ps->indxWhichModel = constmodpars_ps->indxWhichModel;
   constgibbs_ps->kx = constmodpars_ps->kx;
   constgibbs_ps->nlagsinflation = constmodpars_ps->nlagsinflation;
   constgibbs_ps->nlagstrue = constmodpars_ps->nlagstrue;
   constgibbs_ps->fss = constmodpars_ps->fss;
   constgibbs_ps->govpitarget = input_ps->govpitarget;    //Government's targeted inflation.

   //=== Draws of free parameters.
   constgibbs_ps->ndraws1 = input_ps->ndraws1;
   constgibbs_ps->ndrawsq = input_ps->ndrawsq;
   constgibbs_ps->ndraws2 = input_ps->timesk2b*constgibbs_ps->ndraws1;
   constgibbs_ps->napartsave = input_ps->napartsave;
   constgibbs_ps->jumpingratio = 0.0;   //Cumulative.
   constgibbs_ps->sf_mtpl = input_ps->sf_mtpl;
   constgibbs_ps->sf_mhm = input_ps->sf_mhm;
   constgibbs_ps->indx_readgibbs = input_ps->indx_readgibbs;
   constgibbs_ps->indx_lrsacrats = input_ps->indx_lrsacrats;
   constgibbs_ps->indx_mhm = input_ps->indx_mhm;
   constgibbs_ps->indxRecord = input_ps->indxRecord;
   constgibbs_ps->indxMHMPeak = input_ps->indxMHMPeak;
   constgibbs_ps->cutval_logpost = input_ps->cutval_logpost;
   constgibbs_ps->maxgibbs_logpostkernel = -MACHINEINFINITY;
   constgibbs_ps->mingibbs_logpostkernel = MACHINEINFINITY;
   //
   constgibbs_ps->totnpars = constmodpars_ps->totnpars;
   constgibbs_ps->x_dv = CreateVector_lf(constgibbs_ps->totnpars);   //totnpars-by-1 vector of free parameters. For indxWhichModel == 1, we have 6+kx+(kx+1)*kx variables in the order of
          //vstar, theta0, theta1, tau1, zeta1, zeta2, z10 (k-by-1), c_p1, ..., cpf (upper triangular of P10), c_v1, ..., c_vf (upper triangular of V).
   constgibbs_ps->xmax_dv = CreateVector_lf(constgibbs_ps->totnpars);
   constgibbs_ps->xmin_dv = CreateVector_lf(constgibbs_ps->totnpars);
   //--- Other crucial (e.g. government) parameters.
   constgibbs_ps->z10_dv = CreateVector_lf(input_ps->kx);
   constgibbs_ps->P10_dm = CreateMatrix_lf(input_ps->kx,input_ps->kx);
   if (constgibbs_ps->indxFixz10P10 == 1) {
      CopyVector0(constgibbs_ps->z10_dv, constmodpars_ps->z10_dv);
      CopyMatrix0(constgibbs_ps->P10_dm, constmodpars_ps->P10_dm);
   }
   if (constgibbs_ps->indxFixz10P10 == 2)  CopyVector0(constgibbs_ps->z10_dv, constmodpars_ps->z10_dv);


   constgibbs_ps->V_dm = CreateMatrix_lf(input_ps->kx,input_ps->kx);
   if (constgibbs_ps->indx_mhm) {
      constgibbs_ps->xzetapostest_dv = CreateVector_lf(constgibbs_ps->totnpars);
      constgibbs_ps->xzeta_dv = CreateVector_lf(constgibbs_ps->totnpars);
      constgibbs_ps->xmean_dv = CreateVector_lf(constgibbs_ps->totnpars);
      constgibbs_ps->cholXcov_dm = CreateMatrix_lf(constgibbs_ps->totnpars, constgibbs_ps->totnpars);
      //+
      constgibbs_ps->pcuts_dv = CreateVector_lf(input_ps->pcuts_dv->n);
      CopyVector0(constgibbs_ps->pcuts_dv, input_ps->pcuts_dv);
      //Commented out 03/17/05.  tzDestroy(input_ps->pcuts_dv);  //$$$$$$$ Make sure no memory leak.  8/9/04.
      //+ Added 03/17/05.
      constgibbs_ps->zeta1cuts_dv = CreateVector_lf(input_ps->zeta1cuts_dv->n);
      constgibbs_ps->zeta2cuts_dv = CreateVector_lf(input_ps->zeta2cuts_dv->n);
      CopyVector0(constgibbs_ps->zeta1cuts_dv, input_ps->zeta1cuts_dv);
      CopyVector0(constgibbs_ps->zeta2cuts_dv, input_ps->zeta2cuts_dv);
      //Commented out 03/17/05.  tzDestroy(input_ps->zetacuts_dv);  //$$$$$$$ Make sure no memory leak.  8/9/04.
      //+
      constgibbs_ps->logqvals_dv = CreateVector_lf(input_ps->pcuts_dv->n);
      constgibbs_ps->logpcuts_dv = CreateVector_lf(constgibbs_ps->pcuts_dv->n);
      constgibbs_ps->nzeros_pdf_iv = CreateConstantVector_int(constgibbs_ps->pcuts_dv->n, 0);
      constgibbs_ps->veclogMLH_ps = CreateVeclogsum(constgibbs_ps->pcuts_dv->n);
   }
   else {
      constgibbs_ps->xzetapostest_dv = NULL;
      constgibbs_ps->xzeta_dv = NULL;
      constgibbs_ps->xmean_dv = NULL;
      constgibbs_ps->cholXcov_dm =  NULL;
      constgibbs_ps->pcuts_dv = NULL;
      constgibbs_ps->zeta1cuts_dv = NULL;
      constgibbs_ps->zeta2cuts_dv = NULL;
      constgibbs_ps->logqvals_dv = NULL;
      constgibbs_ps->logpcuts_dv = NULL;
      constgibbs_ps->nzeros_pdf_iv = NULL;
      constgibbs_ps->veclogMLH_ps = NULL;
   }
   constgibbs_ps->upred_dv = CreateVector_lf(constgibbs_ps->fss);   //T-by-1 vector of predicted unemployment rate.
   constgibbs_ps->lreu_dv = CreateVector_lf(constmodpars_ps->fss);   //T-by-1 vector of predicted unemployment rate.
   constgibbs_ps->StrResidstran_dm = CreateMatrix_lf(2, constgibbs_ps->fss);

   //--- For the 1st set of parameters: x_dv->v[0:3].  vstar, theta0, theta1, and tau1.
   constgibbs_ps->postmean1stset_dv = CreateVector_lf(N1STSET);    //Posterior mean.
   constgibbs_ps->postCov1stSetinv_dm = CreateMatrix_lf(N1STSET, N1STSET);   //inverse of posterior covaraince matrix.
   constgibbs_ps->cholpostCov1stSetinv_dm = CreateConstantMatrix_lf(N1STSET, N1STSET, 0.0);   //choleski of inverse of posterior covaraince matrix.
   //--- For the 4th set of parameters: x_dv->v[6:end-1]: z_{1|0} and only upper triangular parts of the Choleski decompostions of P_{1|0} and V,  See p.20d.
   if (constgibbs_ps->indxFixz10P10 == 0)
      fourdim = constgibbs_ps->kx+constgibbs_ps->kx*(constgibbs_ps->kx+1);
   else if (constgibbs_ps->indxFixz10P10 == 1)
      fourdim = constgibbs_ps->kx*(constgibbs_ps->kx+1)/2;
   else if (constgibbs_ps->indxFixz10P10 == 2)
      fourdim = constgibbs_ps->kx*(constgibbs_ps->kx+1);
   else  fn_DisplayError(".../swz_comfuns.c/CreateTSconstgibbs(): indxFixz10P10 must be either 0, 1. or 2");

//   constgibbs_ps->postmean4a5a6thsets_dv = CreateVector_lf(fourdim);    //Posterior mean.
//   constgibbs_ps->postCov4a5a6thSetsinv_dm = CreateMatrix_lf(fourdim, fourdim);   //Inverse of posterior covaraince matrix.
//   constgibbs_ps->cholpostCov4a5a6thSetsinv_dm = CreateConstantMatrix_lf(fourdim, fourdim, 0.0);   //Choleski of inverse of posterior covaraince matrix.
   constgibbs_ps->CovMplsinv_dm = NULL;  //Will be created when govprob_ps and constmodpars_ps have legal values.
   constgibbs_ps->cholCovMplsinv_dm = CreateConstantMatrix_lf(fourdim, fourdim, 0.0);

   return (constgibbs_ps);
}
//---
struct TSconstgibbs_tag *DestroyTSconstgibbs(struct TSconstgibbs_tag *constgibbs_ps)
{
   if (constgibbs_ps) {
      DestroyVector_lf(constgibbs_ps->x_dv);
      DestroyVector_lf(constgibbs_ps->xmax_dv);
      DestroyVector_lf(constgibbs_ps->xmin_dv);
      //+
      DestroyVector_lf(constgibbs_ps->z10_dv);
      DestroyMatrix_lf(constgibbs_ps->P10_dm);
      DestroyMatrix_lf(constgibbs_ps->V_dm);
      //
      DestroyVector_lf(constgibbs_ps->xzetapostest_dv);
      DestroyVector_lf(constgibbs_ps->xzeta_dv);
      DestroyVector_lf(constgibbs_ps->xmean_dv);
      DestroyMatrix_lf(constgibbs_ps->cholXcov_dm);
      DestroyVector_lf(constgibbs_ps->pcuts_dv);
      DestroyVector_lf(constgibbs_ps->zeta1cuts_dv);  //Added 03/17/05.
      DestroyVector_lf(constgibbs_ps->zeta2cuts_dv);  //Added 03/17/05.
      DestroyVector_lf(constgibbs_ps->logqvals_dv);      //Added 03/17/05.
      DestroyVector_lf(constgibbs_ps->logpcuts_dv);
      DestroyVector_int(constgibbs_ps->nzeros_pdf_iv);
      DestroyVeclogsum(constgibbs_ps->veclogMLH_ps);
      //
      DestroyVector_lf(constgibbs_ps->upred_dv);
      DestroyVector_lf(constgibbs_ps->lreu_dv);
      DestroyMatrix_lf(constgibbs_ps->StrResidstran_dm);
      DestroyVector_lf(constgibbs_ps->sclhistshocks_dv);
      //+
      DestroyVector_lf(constgibbs_ps->postmean1stset_dv);
      DestroyMatrix_lf(constgibbs_ps->postCov1stSetinv_dm);
      DestroyMatrix_lf(constgibbs_ps->cholpostCov1stSetinv_dm);
      //+
//      DestroyVector_lf(constgibbs_ps->postmean4a5a6thsets_dv);
//      DestroyMatrix_lf(constgibbs_ps->postCov4a5a6thSetsinv_dm);
//      DestroyMatrix_lf(constgibbs_ps->cholpostCov4a5a6thSetsinv_dm);
      DestroyMatrix_lf(constgibbs_ps->CovMplsinv_dm);
      DestroyMatrix_lf(constgibbs_ps->cholCovMplsinv_dm);

      //===
      free(constgibbs_ps);
      return ((struct TSconstgibbs_tag *)NULL);
   }
   else  return (constgibbs_ps);
}


//----------------------------
// Read the posterior estimate from the starting-point file.
//----------------------------
void RefreshWPosteriorEstimates(struct TSconstmodpars_tag *constmodpars_ps, TSinput *input_ps, TSgovprob *govprob_ps, FILE *fptr_sp)
{
   //Refreshes constmodpars->x_dv with the posterior estimate.
   //Note fptr_sp must point to the starting-point file that contains the posterior estimates.
   int _n, _i;
   TSdvector *x_dv = constmodpars_ps->x_dv;


   rewind(fptr_sp);   //Puts the pointer at the beginning of the file.

   for (_n=x_dv->n, _i=0; _i<_n; _i++) {
      if (fscanf(fptr_sp, " %lf ", x_dv->v+_i) != 1) {
         printf("Fatal Error: RefreshWPosteriorEstimates() -- cannot read the number from the file fptr_sp");
         exit(EXIT_FAILURE);
      }
   }
   x_dv->flag = V_DEF;
   //--- Giving the overall posterior peak value.
   constmodpars_ps->peaklogpost = value_logpostkernal_refresh(constmodpars_ps, (struct TSconstgibbs_tag *)NULL, POINTEST, govprob_ps, x_dv);
}




//----------------------------
// Gibbs sampling of the 1st set of parameters (vstar, theta0, theta1, and tau1).
//----------------------------
static void OneDrawOf1stSet(struct TSconstgibbs_tag *constgibbs_ps, struct TSgovprob_tag *govprob_ps, struct TSconstmodpars_tag *constmodpars_ps)
{
   //Note that (1) constmodpars_ps is used only for getting **prior** mean and variance;
   //          (2) govprob_ps is used only for solving the government's controlled inflation x_{t-1}.
   //See pp.20e-20f.
   int ti;
   int fss = constgibbs_ps->fss;
   //--- Initialization. The order matters!
   TSkalcvfurw *kalcvfurw_ps = govprob_ps->kalcvfurw_ps;
   TSdvector *ylhtran_dv = kalcvfurw_ps->ylhtran_dv;
   TSdmatrix *Xrhtran_dm = kalcvfurw_ps->Xrhtran_dm;  //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.
   int nrhvars = Xrhtran_dm->nrows;
   TSdvector *expinf_dv = govprob_ps->expinf_dv;   //Will be updated with a Gibbs draw with govprob_ps->govestimator().
   //
   TSdmatrix *priorCov1stSetinv_dm = constmodpars_ps->Var1stSet_dm;   //Will be an inverse later.
   TSdmatrix *postCov1stSetinv_dm = constgibbs_ps->postCov1stSetinv_dm;
   TSdmatrix *cholpostCov1stSetinv_dm = constgibbs_ps->cholpostCov1stSetinv_dm;
   TSdvector *priormean1stset_dv = constmodpars_ps->mean1stset_dv;
   TSdvector *postmean1stset_dv = constgibbs_ps->postmean1stset_dv;
   TSdvector randnN1_sdv;
   //=== Memory allocated here.
   TSdvector *yt_dv = CreateConstantVector_lf(N1STSET, 1.0);     //1st element is 1.0.
   TSdvector *sumytut_dv = CreateConstantVector_lf(N1STSET, 0.0);   //Cumulative.


   //--- Initialization.
   InitializeConstantMatrix_lf(postCov1stSetinv_dm, 0.0);  //Cumulative.
   postCov1stSetinv_dm->flag = M_SU;   //Upper part symmetric.
   //+
   randnN1_sdv.v = constgibbs_ps->x_dv->v;
   randnN1_sdv.n = N1STSET;
   randnN1_sdv.flag = V_DEF;


//   if (!constgibbs_ps->nlagsinflation && (constgibbs_ps->nlagstrue==1) && !constgibbs_ps->govmodel) {
   if (constgibbs_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0)) {
      //indxWhichModel: government's model (0: Classical direction or regresssion; 1: Keynesian direction).
      //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.
      fss = constgibbs_ps->fss;
      ftd_freepars2modpars((void *)constgibbs_ps, GIBBSDRAWS);
      ftd_normalizationforgov((void *)constgibbs_ps, GIBBSDRAWS);  //Normalization.

      //--- Conditioning on last draw of the government's parameters, solves the government's problem for expected inflation and government's optimization decision for x_t (which controls at least part of the inflation process).
      govprob_ps->govestimator(govprob_ps, (void *)constgibbs_ps, GIBBSDRAWS);

      //------- Forming posterior mean and variance. -------
      if (N1STSET != 4)  fn_DisplayError(".../OneDrawOf1stSet(): Make sure N1STSET matches the specification of yt_dv");
      if (constgibbs_ps->indxWhichModel & CGP2TP1TI0) {
         for (ti=0; ti<fss; ti++) {
            //--- Classical regression: ylhtran_dv is ut; 6 variables for Xrhtran_dm in the order of pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and Constant.
            //--- Note that this order of lagged variables *is* different from that in the government's gensys problem.
            sumytut_dv->v[0] += ylhtran_dv->v[ti] * yt_dv->v[0];  //u_t*1.0.
            sumytut_dv->v[1] += ylhtran_dv->v[ti] * (yt_dv->v[1] = Xrhtran_dm->M[mos(0,ti,nrhvars)] - expinf_dv->v[ti]); //u_t*z_{2t}, p.20e.
            sumytut_dv->v[2] += ylhtran_dv->v[ti] * (yt_dv->v[2] = Xrhtran_dm->M[mos(1,ti,nrhvars)] - ((ti<1) ? Xrhtran_dm->M[mos(1,ti,nrhvars)] : expinf_dv->v[ti-1])); //u_t*z_{2t-1}, p.20e.
            sumytut_dv->v[3] += ylhtran_dv->v[ti] * (yt_dv->v[3] = Xrhtran_dm->M[mos(2,ti,nrhvars)]);  //u_t*u_{t-1}.

            VectorTimesSelf(postCov1stSetinv_dm, yt_dv, 1.0, 1.0, 'U');  //Cumulative.  sum(y_t*y_t').
         }
      }
      else if (constgibbs_ps->indxWhichModel & KGP2TP1TI0) {
         for (ti=0; ti<fss; ti++) {
            //--- Keynesian regression: ylhtran_dv is pi_t; 6 variables for Xrhtran_dm in the order of u_t, u_{t-1}, pi_{t-1}, u_{t-2}, pi_{t-2}, and Constant.
            //--- Note that this order of lagged variables *is* different from that in the government's gensys problem.
            sumytut_dv->v[0] += Xrhtran_dm->M[mos(0,ti,nrhvars)] * yt_dv->v[0];  //u_t*1.0.
            sumytut_dv->v[1] += Xrhtran_dm->M[mos(0,ti,nrhvars)] * (yt_dv->v[1] = ylhtran_dv->v[ti] - expinf_dv->v[ti]); //u_t*z_{2t}, p.20e.
            sumytut_dv->v[2] += Xrhtran_dm->M[mos(0,ti,nrhvars)] * (yt_dv->v[2] = Xrhtran_dm->M[mos(2,ti,nrhvars)] - ((ti<1) ? Xrhtran_dm->M[mos(2,ti,nrhvars)] : expinf_dv->v[ti-1])); //u_t*z_{2t-1}, p.20e.
            sumytut_dv->v[3] += Xrhtran_dm->M[mos(0,ti,nrhvars)] * (yt_dv->v[3] = Xrhtran_dm->M[mos(1,ti,nrhvars)]);  //u_t*u_{t-1}.

            VectorTimesSelf(postCov1stSetinv_dm, yt_dv, 1.0, 1.0, 'U');  //Cumulative.  sum(y_t*y_t').
         }
      }
      //+
      if (invspd(priorCov1stSetinv_dm, priorCov1stSetinv_dm, (priorCov1stSetinv_dm->flag & M_SU) ? 'U' : 'L'))
         fn_DisplayError(".../OneDrawOf1stSet(): constmodpars_ps->Var1stSet_dm is not invertible when calling invspd()");
      ScalarTimesMatrix(postCov1stSetinv_dm, 1.0, priorCov1stSetinv_dm, constgibbs_ps->gzeta1);
      //+ The order matters!
      ScalarTimesVector(postmean1stset_dv, constgibbs_ps->gzeta1, sumytut_dv, 0.0);  //zeta1*sum(yt*ut).
      SymmatrixTimesVector(postmean1stset_dv, priorCov1stSetinv_dm, priormean1stset_dv, 1.0, 1.0);  //inv(Sigmag^bar)*theta^bar + zeta1*sum(yt*ut).  See p.20f.
      if (chol(cholpostCov1stSetinv_dm, postCov1stSetinv_dm, (postCov1stSetinv_dm->flag & M_SU) ? 'U' : 'L'))
         fn_DisplayError(".../OneDrawOf1stSet(): constgibbs_ps->postCov1stSet_dm must be spd when calling chol()");
      //--- postmean1stset_dv = Sigma^tld * (inv(Sigmag^bar)*theta^bar + zeta1*sum(yt*ut)).  See p.20f.
      SolveTriSysVector(postmean1stset_dv, cholpostCov1stSetinv_dm, postmean1stset_dv, 'T', 'N');
      SolveTriSysVector(postmean1stset_dv, cholpostCov1stSetinv_dm, postmean1stset_dv, 'N', 'N');


      //------- Normal posterior draw. --------
      StandardNormalVector(&randnN1_sdv);
      SolveTriSysVector(&randnN1_sdv, cholpostCov1stSetinv_dm, &randnN1_sdv, 'N', 'N');
      VectorPlusMinusVectorUpdate(&randnN1_sdv, postmean1stset_dv, 1.0);
   }
   else {
      printf(".../swz_comfuns.c/OneDrawOf1stSet(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d!", CGP2TP1TI0, KGP2TP1TI0);
      exit(EXIT_FAILURE);
   }
//   else  fn_DisplayError(".../OneDrawOf1stSet():  Have not got time to handle the case for nlagstrue != 1 or nlagsinflation !=0 or Keyesian model (govmodel != 0)");


   //===
   DestroyVector_lf(yt_dv);
   DestroyVector_lf(sumytut_dv);
}
//----------------------------
// Gibbs sampling of the 2nd and 3rd parameters (zeta1 and zeta2).
//----------------------------
static void OneDrawOf2nd3rdSet(struct TSconstgibbs_tag *constgibbs_ps, struct TSgovprob_tag *govprob_ps, struct TSconstmodpars_tag *constmodpars_ps)
{
   //Note that (1) constmodpars_ps is used only for getting *prior* settings.
   //          (2) govprob_ps is used only for solving the government's controlled inflation x_{t-1}.
   //See pp.20e-20f.
   int ti;
   int fss = constgibbs_ps->fss;
   TSdvector *x_dv = constgibbs_ps->x_dv;
   //--- Initialization. The order matters!
   TSkalcvfurw *kalcvfurw_ps = govprob_ps->kalcvfurw_ps;
   TSdvector *ylhtran_dv = kalcvfurw_ps->ylhtran_dv;
   TSdmatrix *Xrhtran_dm = kalcvfurw_ps->Xrhtran_dm;  //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.
   int nrhvars = Xrhtran_dm->nrows;
   TSdvector *expinf_dv = govprob_ps->expinf_dv;    //Will be updated with a Gibbs draw with govprob_ps->govestimator().
   //---
   double gzeta1p2_postalpha;
   double gzeta1_postbeta, gzeta2_postbeta;
   double z1t, z2t, sum_z1tsq, sum_z2tsq;   //Cumulative.


   //------- Getting gzeta1p2_postbeta. -------
//   if (!constgibbs_ps->nlagsinflation && (constgibbs_ps->nlagstrue==1) && !constgibbs_ps->govmodel) {
   if (constgibbs_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0)) {
      //indxWhichModel: government's model (0: Classical direction or regresssion; 1: Keynesian direction).
      //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.
      fss = constgibbs_ps->fss;
      ftd_freepars2modpars((void *)constgibbs_ps, GIBBSDRAWS);
      ftd_normalizationforgov((void *)constgibbs_ps, GIBBSDRAWS);  //Normalization.

      //--- Conditioning on last draw of the government's parameters, solves the government's problem for expected inflation and government's optimization decision for x_t (which controls at least part of the inflation process).
      govprob_ps->govestimator(govprob_ps, (void *)constgibbs_ps, GIBBSDRAWS);

      //------- Forming posterior alpha and beta in Gamma distribution. -------
      gzeta1p2_postalpha = 0.5*constgibbs_ps->fss + constmodpars_ps->gzeta1p2_prioralpha;
      sum_z1tsq = sum_z2tsq = 0.0;   //Cumulative.
      if (constgibbs_ps->indxWhichModel & CGP2TP1TI0) {
         for (ti=0; ti<fss; ti++) {
            //--- Classical regression: ylhtran_dv is u_t; 6 variables for Xrhtran_dm in the order of pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and Constant.
            //--- Note that this order of lagged variables *is* different from that in the government's gensys problem.
            z1t = ylhtran_dv->v[ti] - (x_dv->v[0] + constgibbs_ps->gtau1*Xrhtran_dm->M[mos(2,ti,nrhvars)]
                  + constgibbs_ps->gtheta0*(z2t=(Xrhtran_dm->M[mos(0,ti,nrhvars)] - expinf_dv->v[ti]))
                  + constgibbs_ps->gtheta1*(Xrhtran_dm->M[mos(1,ti,nrhvars)] - ((ti<1) ? Xrhtran_dm->M[mos(1,ti,nrhvars)] : expinf_dv->v[ti-1])));
                                     //x_dv->v[0] is vstar on p.20a.
            sum_z1tsq += square(z1t); //sum(z_{1t}^2).
            sum_z2tsq += square(z2t); //sum(z_{2t}^2), p.20e.
         }
      }
      else if (constgibbs_ps->indxWhichModel & KGP2TP1TI0) {
         for (ti=0; ti<fss; ti++) {
            //--- Keynesian regression: ylhtran_dv is pi_t; 6 variables for Xrhtran_dm in the order of u_t, u_{t-1}, pi_{t-1}, u_{t-2}, pi_{t-2}, and Constant.
            //--- Note that this order of lagged variables *is* different from that in the government's gensys problem.
            z1t = Xrhtran_dm->M[mos(0,ti,nrhvars)] - (x_dv->v[0] + constgibbs_ps->gtau1*Xrhtran_dm->M[mos(1,ti,nrhvars)]
                  + constgibbs_ps->gtheta0*(z2t=(ylhtran_dv->v[ti] - expinf_dv->v[ti]))
                  + constgibbs_ps->gtheta1*(Xrhtran_dm->M[mos(2,ti,nrhvars)] - ((ti<1) ? Xrhtran_dm->M[mos(2,ti,nrhvars)] : expinf_dv->v[ti-1])));
                                     //x_dv->v[0] is vstar on p.20a.
            sum_z1tsq += square(z1t); //sum(z_{1t}^2).
            sum_z2tsq += square(z2t); //sum(z_{2t}^2), p.20e.
         }
      }
      gzeta1_postbeta = 1.0/(0.5*sum_z1tsq + 1.0/constmodpars_ps->gzeta1p2_priorbeta);
      gzeta2_postbeta = 1.0/(0.5*sum_z2tsq + 1.0/constmodpars_ps->gzeta1p2_priorbeta);

      //=== Gamma posterior draw.
      *(x_dv->v + N1STSET) = sqrt(GammaDouble(gzeta1p2_postalpha, gzeta1_postbeta));
      *(x_dv->v + N1STSET + N2NDSET) = sqrt(GammaDouble(gzeta1p2_postalpha, gzeta2_postbeta));
   }
   else {
      printf(".../swz_comfuns.c/OneDrawOf2nd3rdSet(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d!", CGP2TP1TI0, KGP2TP1TI0);
      exit(EXIT_FAILURE);
   }
//   else  fn_DisplayError(".../OneDrawOf2nd3rdSet():  Have not got time to handle the case for nlagstrue != 1 or nlagsinflation !=0 or Keyesian model (govmodel != 0)");
}
//----------------------------
// Gibbs sampling of the 4th, 5th, and 6th sets of parameters: z_{1|0} and only upper triangular parts of the Choleski decompostions of P_{1|0} and V,  See p.20d.
//----------------------------
static void OneDrawOf4a5a6thSets(struct TSconstgibbs_tag *constgibbs_ps, struct TSgovprob_tag *govprob_ps, struct TSconstmodpars_tag *constmodpars_ps)
{
   //Note that (1) constmodpars_ps is used only for getting *prior* settings and posterior estimate gphihat.
   //          (2) govprob_ps is used only for solving the government's controlled inflation x_{t-1}.
   //See pp.20h, 20e-20f.
   int fss = constgibbs_ps->fss;
   int dim4a5a6thset;
   double proposallogval, lastlogval;
   //--- Initialization. The order matters!
   TSkalcvfurw *kalcvfurw_ps = govprob_ps->kalcvfurw_ps;
   TSdvector *ylhtran_dv = kalcvfurw_ps->ylhtran_dv;
   TSdmatrix *Xrhtran_dm = kalcvfurw_ps->Xrhtran_dm;  //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.
   int nrhvars = Xrhtran_dm->nrows;
   TSdvector *expinf_dv = govprob_ps->expinf_dv;    //Will be updated with a Gibbs draw with govprob_ps->govestimator().
   //===
   TSdvector *gphidrawold_dv = NULL;
   TSdvector *gphidrawstar_dv = NULL;

   if (constgibbs_ps->indxFixz10P10 == 0)
      dim4a5a6thset = constgibbs_ps->kx+constgibbs_ps->kx*(constgibbs_ps->kx+1);
   else if (constgibbs_ps->indxFixz10P10 == 1)
      dim4a5a6thset = constgibbs_ps->kx*(constgibbs_ps->kx+1)/2;
   else if (constgibbs_ps->indxFixz10P10 == 2)
      dim4a5a6thset = constgibbs_ps->kx*(constgibbs_ps->kx+1);
   else  fn_DisplayError(".../swz_comfuns.c/OneDrawOf4a5a6thSets(): indxFixz10P10 must be either 0, 1. or 2");

   gphidrawold_dv = CreateVector_lf(dim4a5a6thset);  //Last (previous) Metropolis draw, to be decided whether to be kept.
   gphidrawstar_dv = CreateVector_lf(dim4a5a6thset);  //Metropolis proposal draw for gphidraw_sdv.


//   if (!constgibbs_ps->nlagsinflation && (constgibbs_ps->nlagstrue==1) && !constgibbs_ps->govmodel) {
   if (constgibbs_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0)) {
      //indxWhichModel: government's model (0: Classical direction or regresssion; 1: Keynesian direction).
      //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.

      //----------------
      // The order matters!
      //----------------
      //=== Previous Metropolis draw.
      fss = constgibbs_ps->fss;
      ftd_freepars2modpars((void *)constgibbs_ps, GIBBSDRAWS);
      ftd_normalizationforgov((void *)constgibbs_ps, GIBBSDRAWS);  //Normalization.
      //+
      lastlogval = value_logpostkernal_refresh(constmodpars_ps, constgibbs_ps, GIBBSDRAWS, govprob_ps, constgibbs_ps->x_dv);
      if (constgibbs_ps->indxRecord && (lastlogval>constgibbs_ps->maxgibbs_logpostkernel)) {
         CopyVector0(constgibbs_ps->xmax_dv, constgibbs_ps->x_dv);
         constgibbs_ps->maxgibbs_logpostkernel = lastlogval;
      }
      if (constgibbs_ps->indxRecord && (lastlogval<constgibbs_ps->mingibbs_logpostkernel)) {
         CopyVector0(constgibbs_ps->xmin_dv, constgibbs_ps->x_dv);
         constgibbs_ps->mingibbs_logpostkernel = lastlogval;
      }
      updategovfreepars5modfreepars((void *)constgibbs_ps, GIBBSDRAWS, gphidrawold_dv);

      //=== Proposl normal draw.
      StandardNormalVector(gphidrawstar_dv);
      SolveTriSysVector(gphidrawstar_dv, constgibbs_ps->cholCovMplsinv_dm, gphidrawstar_dv, 'N', 'N');
      ScalarTimesVector(gphidrawstar_dv, 1.0, gphidrawold_dv, constgibbs_ps->sf_mtpl);  //Proposal Metropolis draw.
      // VectorPlusMinusVectorUpdate(gphidrawstar_dv, gphidrawold_dv, 1.0);       //Proposal Metropolis draw.

      //=== Testing for acceptance or jumping -- constgibbs_ps will be affected in the middle.
      updatemodfreepars5govfreepars((void *)constgibbs_ps, GIBBSDRAWS, gphidrawstar_dv);
      proposallogval = value_logpostkernal_refresh(constmodpars_ps, constgibbs_ps, GIBBSDRAWS, govprob_ps, constgibbs_ps->x_dv);
      if (log(UniformDouble()) <= (proposallogval - lastlogval) ) {
         //Jumping or accepting.
         constgibbs_ps->jumpingratio += 1.0;
         constgibbs_ps->logLHtimespriorpdf_lastval = proposallogval;  //For MHM method of computing Bayes factors.
      }
      else {
         //Sets the current draw back to be the previous draw,
         updatemodfreepars5govfreepars((void *)constgibbs_ps, GIBBSDRAWS, gphidrawold_dv);
         constgibbs_ps->logLHtimespriorpdf_lastval = lastlogval;      //For MHM method of computing Bayes factors.
      }
   }
   else {
      printf(".../swz_comfuns.c/OneDrawOf4a5a6thSets(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d!", CGP2TP1TI0, KGP2TP1TI0);
      exit(EXIT_FAILURE);
   }


   //===
   DestroyVector_lf(gphidrawold_dv);
   DestroyVector_lf(gphidrawstar_dv);
}
//----------------------------
// Updating towards getting the log marginal likelihood using the MHM method.
//----------------------------
static void ftd_UpdatelogMLH_mhm(struct TSconstgibbs_tag *constgibbs_ps)
{
   int _i, _j;
   int loc_printoutp = 0;  //$$$$$$$$Hard coded.$$$$$$$$$$
   double quadterm, logh;
   double mhmterm;
   //---
   TSdvector *xzeta_dv = constgibbs_ps->xzeta_dv;
   int totpars = xzeta_dv->n;
   TSdvector *xmean_dv = constgibbs_ps->xmean_dv;
   TSdmatrix *cholXcov_dm = constgibbs_ps->cholXcov_dm;
   TSdvector *pcuts_dv = constgibbs_ps->pcuts_dv;
   int npcuts = pcuts_dv->n;
   TSdvector *logpcuts_dv = constgibbs_ps->logpcuts_dv;
   double neghalfklog2pi_minus_dethalfcholXcov = constgibbs_ps->neghalfklog2pi_minus_dethalfcholXcov;
                                       //= -(k/2)*log(2pi) - log(det(cholXcov)).
   double logLHtimespriorpdf_lastval = constgibbs_ps->logLHtimespriorpdf_lastval;
   TSivector *nzeros_pdf_iv = constgibbs_ps->nzeros_pdf_iv;
   struct TSveclogsum_tag *veclogMLH_ps = constgibbs_ps->veclogMLH_ps;
   TSivector *N_iv = veclogMLH_ps->N_iv;
   TSdvector *logsum_dv = veclogMLH_ps->logsum_dv;
   TSdvector *logmax_dv = veclogMLH_ps->logmax_dv;
   TSdvector *logqvals_dv = constgibbs_ps->logqvals_dv;
   //===
   TSdvector *gepsilon_dv = CreateVector_lf(xzeta_dv->n);



   //=== DDDDDDDDDDebug.
   // fprintf(fptr_debug, "Draw xzeta_dv:\n");
   // WriteVector(fptr_debug, xzeta_dv, " %.16e ");
   // fprintf(fptr_debug, "Choleski U transpose:\n");
   // WriteMatrix(fptr_debug, cholXcov_dm, " %.16e ");
   // fprintf(fptr_debug, "quadterm: %d\n", quadterm);
   // fflush(fptr_debug);
   //--- Printing out log importance pdf and log posterior kernal.
   // logh = neghalfklog2pi_minus_dethalfcholXcov - 0.5*quadterm;
   // fprintf(fptr_debug, "logpdf_importance: %10.6f;   logpostkernal: %10.6f\n", logh ,logLHtimespriorpdf_lastval);


   if ( (logLHtimespriorpdf_lastval<constgibbs_ps->cutval_logpost) ||
        (xzeta_dv->v[4]<constgibbs_ps->zeta1cuts_dv->v[0]) || (xzeta_dv->v[5]<constgibbs_ps->zeta2cuts_dv->v[0]) ||
        (xzeta_dv->v[4]>constgibbs_ps->zeta1cuts_dv->v[1]) || (xzeta_dv->v[5]>constgibbs_ps->zeta2cuts_dv->v[1]) )   //Added 03/17/05.
   {
//      printf("\n*******ZEROS*********\n");  //DDDDDDDebug.
//      fflush(stdout);

      //Early exit by setting the importance pdf h() = 0 for all _j's.
      for (_j=npcuts-1; _j>=0; _j--) {
         mhmterm = -MACHINEINFINITY;
         //--- DDDDDebug.
         if (_j==loc_printoutp)
         {
            fprintf(FPTR_DEBUG, " %.16e \n", mhmterm);
            fflush(FPTR_DEBUG);
         }

         ++nzeros_pdf_iv->v[_j];
         ++N_iv->v[_j];  //Yet, the log value of the sum remains the same because we add 0.0 to the sum (or logh() - log(LH*prior) = -infinity).
      }
   }
   else
   {
      if (constgibbs_ps->indxMHMPeak)  VectorMinusVector(gepsilon_dv, xzeta_dv, constgibbs_ps->xzetapostest_dv);  //gepsilon = xdraw - xhat;
      else  VectorMinusVector(gepsilon_dv, xzeta_dv, xmean_dv);  //gepsilon = xdraw - xmean;
      SolveTriSysVector(gepsilon_dv, cholXcov_dm, gepsilon_dv, 'T', 'N');  //Now epsilon = inv(U')*epsilon, becoming a standard normal.
      quadterm = VectorDotVector(gepsilon_dv, gepsilon_dv);

//      for ( _i=npcuts-1; _i>=0; _i--) {
//         if (quadterm < fn_chi2inv(pcuts_dv->v[_i], (double)totpars)) {

            _i=0;
            //--- Computing the log importance density h().
            //logh = -logqvals_dv->v[_i] - logpcuts_dv->v[_i] + neghalfklog2pi_minus_dethalfcholXcov - 0.5*quadterm;  //log h(). -logqvals_dv->v[_i] is added 03/17/05.
            logh = -logqvals_dv->v[_i] + neghalfklog2pi_minus_dethalfcholXcov - 0.5*quadterm;  //log h(). -logqvals_dv->v[_i] is added 03/17/05.
            mhmterm = logh-logLHtimespriorpdf_lastval;
            //--- Accumulating log of sum of h()/(LH*prior).
            N_iv->v[_i] = fn_update_logofsum(N_iv->v[_i], mhmterm, logsum_dv->v+_i, logmax_dv->v+_i);

            //--- DDDDDebug.
            if (_i==loc_printoutp)
            {
               fprintf(FPTR_DEBUG, " %.16e \n", mhmterm);
               fflush(FPTR_DEBUG);
            }
//         }
//         else {
//            //Early exit by setting the importance pdf h() = 0 for the rest of _i's.
//            for (_j=_i; _j>=0; _j--) {
//               mhmterm = -MACHINEINFINITY;
//               //--- DDDDDebug.
//               if (_j==loc_printoutp)
//               {
//                  fprintf(FPTR_DEBUG, " %.16e \n", mhmterm);
//                  fflush(FPTR_DEBUG);
//               }

//               ++nzeros_pdf_iv->v[_j];
//               ++N_iv->v[_j];  //Yet, the log value of the sum remains the same because we add 0.0 to the sum (or logh() - log(LH*prior) = -infinity).
//            }
//            break;   //Early exit.
//         }
//      }
   }



   //===
   DestroyVector_lf(gepsilon_dv);
}



#define NHLOG2PI (-0.91893853320467)       // Negative half of log(2*pi) or -0.5*log(2*pi).  Note pi = 3.1415926535897931.
//----------------------------
// Gibbs draws of all model parameters.
//----------------------------
void GibbsSimulations(struct TSconstgibbs_tag *constgibbs_ps, struct TSgovprob_tag *govprob_ps, struct TSconstmodpars_tag *constmodpars_ps, struct TSinput_tag *input_ps, char *filename_sp, FILE *fptr_inpfor_matlab, char *filename_hess)
{
   int _i;
   int drawi;
   int ndraws1 = constgibbs_ps->ndraws1,
       ndraws2 = constgibbs_ps->ndraws2;
   int _n, errflag;
   int cntbuffer = 0;  //Cumulative.  Must be set to 0.
   int cntsave = 0,  //Cumulative.  Must be set to 0.
       napartsave = constgibbs_ps->napartsave;
   double workd;
   double ndraws1inv = 1.0/((double)ndraws1);
   time_t sim_endtime;
   FILE *fptr_sp = NULL;
   FILE *fptr_sp_hess = NULL;


   if (govprob_ps->indxFixz10P10 == 0)
      _n = govprob_ps->kx+govprob_ps->kx*(govprob_ps->kx+1);
   else if (govprob_ps->indxFixz10P10 == 1)
      _n = govprob_ps->kx*(govprob_ps->kx+1)/2;
   else if (govprob_ps->indxFixz10P10 == 2)
      _n = govprob_ps->kx*(govprob_ps->kx+1);
   else  fn_DisplayError(".../swz_comfuns.c/GibbsSimulations(): indxFixz10P10 must be either 0, 1. or 2");


   //------- (1) The posterior estimates is copied to constmodpars_ps->x_dv, which leads to constgibbs_ps->xzetapostest_dv. -------
   //------- (2) The numerical Hessian at the posterior estimates is created or read. -------
   if (constgibbs_ps->indxHess) {
      //=== Reads the numerical Hessian for the Metropolis algorithm.
      fptr_sp_hess = tzFopen(filename_hess,"r");
      rewind(fptr_sp_hess);   //Must put the pointer at the beginning of the file.
      if (!ReadMatrix_lf(fptr_sp_hess, constgibbs_ps->CovMplsinv_dm = CreateConstantMatrix_lf(_n, _n, 0.0))) {
         printf(".../GibbsSimulations(): cannot read the spd matrix from the file %s", filename_hess);
         exit(EXIT_FAILURE);
      }
      constgibbs_ps->CovMplsinv_dm->flag = M_SU;
      tzFclose(fptr_sp_hess);

      //=== Reads the posterior estimates into constmodpars_ps->x_dv.
      fptr_sp = tzFopen(filename_sp,"r");
      rewind(fptr_sp);   //Must put the pointer at the beginning of the file.
      RefreshWPosteriorEstimates(constmodpars_ps, input_ps, govprob_ps, fptr_sp);
      tzFclose(fptr_sp);
   }
   else {
      fptr_sp_hess = tzFopen(filename_hess,"w");
      //+
      fptr_sp = tzFopen(filename_sp,"r");
      rewind(fptr_sp);   //Must put the pointer at the beginning of the file.
      constgibbs_ps->CovMplsinv_dm = CreateNumericalHessian(govprob_ps, constmodpars_ps, input_ps, fptr_sp);
      WriteMatrix(fptr_sp_hess, constgibbs_ps->CovMplsinv_dm, " %.16e ");
      //
      tzFclose(fptr_sp_hess);
      tzFclose(fptr_sp);
   }
   chol(constgibbs_ps->cholCovMplsinv_dm, constgibbs_ps->CovMplsinv_dm, 'U');


   //--- Initializes the Gibbs draw with the posterior estimates constmodpars_ps->x_dv.
   value_logpostkernal_refresh(constmodpars_ps, constgibbs_ps, GIBBSDRAWS, govprob_ps, constmodpars_ps->x_dv);


   //=== 1st set of Gibbs draws to be tossed away.
   constgibbs_ps->indxRecord = 0;  //Does not record, say, the max of logPostKernel.
   constgibbs_ps->jumpingratio = 0.0;  //Cumulative.
   for (drawi=0; drawi<ndraws1; drawi++) {
      OneDrawOf1stSet(constgibbs_ps, govprob_ps, constmodpars_ps);
      OneDrawOf2nd3rdSet(constgibbs_ps, govprob_ps, constmodpars_ps);
      OneDrawOf4a5a6thSets(constgibbs_ps, govprob_ps, constmodpars_ps);
   }
   constgibbs_ps->jumpingratio /= (double)ndraws1;
   constgibbs_ps->jumpingratio *= 100.0;
   printf("\nThe jumping ratio for the 1st tossed-up loop is %5.2f\n", constgibbs_ps->jumpingratio);
   fflush( stdout );


   //=== A set of Gibbs draws to get the mean and covariance of the importance Gaussian pdf using the MHM method.
   if (constgibbs_ps->indx_mhm) {
      CopyVector0(constgibbs_ps->xzetapostest_dv, constmodpars_ps->x_dv);
      constgibbs_ps->xzetapostest_dv->v[4] = square(constmodpars_ps->x_dv->v[4]);
      constgibbs_ps->xzetapostest_dv->v[5] = square(constmodpars_ps->x_dv->v[5]);
      constgibbs_ps->indxRecord = 1;  //Begins to record, say, the max of logPostKernel.  Used by OneDrawOf4a5a6thSets().
      for (drawi=0; drawi<ndraws1; drawi++) {
         OneDrawOf1stSet(constgibbs_ps, govprob_ps, constmodpars_ps);
         OneDrawOf2nd3rdSet(constgibbs_ps, govprob_ps, constmodpars_ps);
         OneDrawOf4a5a6thSets(constgibbs_ps, govprob_ps, constmodpars_ps);

         //--- Updating 1st and 2nd moments for the importance Gaussian pdf.
         CopyVector0(constgibbs_ps->xzeta_dv, constgibbs_ps->x_dv);
         constgibbs_ps->xzeta_dv->v[4] = square(constgibbs_ps->x_dv->v[4]);
         constgibbs_ps->xzeta_dv->v[5] = square(constgibbs_ps->x_dv->v[5]);
         UpdateSumFor1st2ndMoments(constgibbs_ps->xmean_dv, constgibbs_ps->cholXcov_dm, constgibbs_ps->xzeta_dv);
      }
      ScalarTimesVector(constgibbs_ps->xmean_dv, ndraws1inv, constgibbs_ps->xmean_dv, 0.0);  //Converting it to mean.

      //=== Converting 2nd moment to covariance matrix.
      if (constgibbs_ps->indxMHMPeak) {
         //Centers the importance pdf at the posterior peak.
         VectorTimesSelf(constgibbs_ps->cholXcov_dm, constgibbs_ps->xzeta_dv, 1.0, ndraws1inv, 'U'); //X'X/n + xhat*xhat';
         VectorTimesVector(constgibbs_ps->cholXcov_dm, constgibbs_ps->xzetapostest_dv, constgibbs_ps->xmean_dv, -1.0, 1.0);  //X'X/n + xhat*xhat' - xhat*xmean';
         VectorTimesVector(constgibbs_ps->cholXcov_dm, constgibbs_ps->xmean_dv, constgibbs_ps->xzetapostest_dv, -1.0, 1.0);  //X'X/n + xhat*xhat' - xhat*xmean' - xmean*xhat';
         ScalarTimesMatrix(constgibbs_ps->cholXcov_dm, constgibbs_ps->sf_mhm, constgibbs_ps->cholXcov_dm, 0.0);
         constgibbs_ps->cholXcov_dm->flag = M_SU;  //Getting rid of the flag M_GE to make sure the upper symmetric part is used later.
      }
      else {
         //Centers the importance pdf at the posterior mean -- a traditional method.
         VectorTimesSelf(constgibbs_ps->cholXcov_dm, constgibbs_ps->xmean_dv, -1.0, ndraws1inv, 'U'); //Converting it to covariance matrix.
         ScalarTimesMatrix(constgibbs_ps->cholXcov_dm, constgibbs_ps->sf_mhm, constgibbs_ps->cholXcov_dm, 0.0);
         constgibbs_ps->cholXcov_dm->flag = M_SU;  //Getting rid of the flag M_GE to make sure the upper symmetric part is used later.
      }


      //=== Getting the upper Choleski of the covariance matrix.
      if (errflag=chol(constgibbs_ps->cholXcov_dm, constgibbs_ps->cholXcov_dm, 'U')) {
         //cholXcov_dm is now upper triangular Choleski so that chol'chol = Xcov.
         printf("\n.../swz_comfuns.c/GibbsSimulations(): No chol decomp for cholXcov_dm with error flag %d\n", errflag);
         exit(EXIT_FAILURE);
      }
      //--- Taking log of pcuts_dv and computing a normal pdf term outside the second loop.
      for (_i=constgibbs_ps->pcuts_dv->n-1; _i>=0; _i--) {
         constgibbs_ps->logpcuts_dv->v[_i] = log(constgibbs_ps->pcuts_dv->v[_i]);
         constgibbs_ps->logpcuts_dv->flag = V_DEF;
      }
      constgibbs_ps->neghalfklog2pi_minus_dethalfcholXcov = NHLOG2PI*(double)constgibbs_ps->x_dv->n - tracelogfabs(constgibbs_ps->cholXcov_dm);
                                   //-(k/2)*log(2pi) - log(det(cholXcov)).
   }

   //=== Printing out the matters.
   time(&sim_endtime);  //End time for this first loop.
   printf("The minutes of completing the 1st tossed-up loop (and setting up the importance pdf under MHM) is %5.2f\n", difftime(sim_endtime, constgibbs_ps->prog_begtime)/60.0);
   fflush( stdout );



   //=== 2nd set of Gibbs draws to be kept and written to the file pointed by fptr_inpfor_matlab.
   constgibbs_ps->indxRecord = 1;  //Begins to record, say, the max of logPostKernel.  Used by OneDrawOf4a5a6thSets().
   constgibbs_ps->jumpingratio = 0.0;  //Cumulative.
   rewind(fptr_inpfor_matlab);   //Puts the pointer at the beginning of the file.
   for (drawi=0; drawi<ndraws2; drawi++) {
      OneDrawOf1stSet(constgibbs_ps, govprob_ps, constmodpars_ps);
      OneDrawOf2nd3rdSet(constgibbs_ps, govprob_ps, constmodpars_ps);
      OneDrawOf4a5a6thSets(constgibbs_ps, govprob_ps, constmodpars_ps);
      //--- Check the condition (e.g., u*>0) before writing the draw to a file.
      if (fabs(workd=1.0-constgibbs_ps->x_dv->v[3])>MACHINEZERO)
         if (constgibbs_ps->x_dv->v[0]/workd < 0.0) {
            drawi--;
            continue;
         }



      //--- Updating towards getting the log marginal likelihood using the MHM method.
      if (constgibbs_ps->indx_mhm) {
         CopyVector0(constgibbs_ps->xzeta_dv, constgibbs_ps->x_dv);
         constgibbs_ps->xzeta_dv->v[4] = square(constgibbs_ps->x_dv->v[4]);
         constgibbs_ps->xzeta_dv->v[5] = square(constgibbs_ps->x_dv->v[5]);
         ftd_UpdatelogMLH_mhm(constgibbs_ps);
      }


      //--- Writes the draw to a file.
      if (++cntsave == napartsave) {
         WriteVector(fptr_inpfor_matlab, constgibbs_ps->x_dv, AFORMAT);
         fflush(fptr_inpfor_matlab);
         cntsave = 0;  //Cumulative.  Must be reset to 0.
      }

      if (++cntbuffer == ndraws1) {
         printf("\n\n==============================================");
         printf("\n (A) Total number of kept draws for the 2nd loop: %d\n (B) Number of draws so far: %d\n (C) The peak log posterior kernal at this point is: %g\n (D) The lowest log posterior kernal at this point is: %g",
                  ndraws2, drawi, constgibbs_ps->maxgibbs_logpostkernel, constgibbs_ps->mingibbs_logpostkernel);
         fflush(stdout);          //Flush the buffer to get out this message without delay.
         cntbuffer = 0;  //Cumulative.  Must be reset to 0.
      }

   }
   constgibbs_ps->jumpingratio /= (double)ndraws2;
   constgibbs_ps->jumpingratio *= 100.0;
}
//----------------------------
// This function is used to
//   (1) read the posterior draws xdraw_dv to search the Gibbs max and export xmax_dv and to compute 1st and 2nd moments of xzetadraw_dv; //Added 03/17/05.
//   (2) compute the normalizing constant q for each cut-off p value of Gaussain weight function.
//   (3) read the posterior draws xdraw_dv to compute the marginal data density;
//   (4) read the posterior draws xdraw_dv to draw perceived long-run excess unemployment under Ramsey policy.
//----------------------------
void ReadGibbsDraws_Compute(struct TSconstgibbs_tag *constgibbs_ps, struct TSgovprob_tag *govprob_ps, struct TSconstmodpars_tag *constmodpars_ps, struct TSinput_tag *input_ps, char *filename_sp, char *filename_readgibbs, FILE *fptr_lreutradeoffs, FILE *fptr_inpfor_maxGibbs)     //fptr_inpfor_maxGibbs is added 03/17/05.
{
   int _i, _n, drawi;
   int ndraws2, cntdraws;
   int nbuffer = 1000;
   int cntbuffer = 0;  //Cumulative.  Must be set to 0.
   int indic_end = 0;
   int errflag;
   int indxGibbsmax = 0;  //$$$$$Hard coded.$$$$$$$ Added 03/17/05.  1: using the max in the Gibbs draws; 0: using the posterior estimate.
   double ndraws2inv;
   //---
   TSdvector *xdraw_dv = constgibbs_ps->x_dv;
   TSdvector *xzetadraw_dv = constgibbs_ps->xzeta_dv;
   //
   TSdvector *lreudraw_dv = constgibbs_ps->lreu_dv;
   TSdmatrix *Zpredtran_dm = govprob_ps->kalcvfurw_ps->Zpredtran_dm;  //6-by-fss(T).
   TSdvector *z10_dv = govprob_ps->kalcvfurw_ps->z10_dv;   //fss(T)-by-1.
   double sumofcoeff_pi, sumofcoeff_u;
   int ti;
   int fss = constgibbs_ps->fss;
   int kx = constgibbs_ps->kx;
//   int ndraws1 = constgibbs_ps->ndraws1,
//       ndraws2 = constgibbs_ps->ndraws2;
//   int _n, errflag;
//   int  cntbuffer = 0;  //Cumulative.  Must be set to 0.
//   int  cntsave = 0,  //Cumulative.  Must be set to 0.
//        napartsave = constgibbs_ps->napartsave;
//   double workd;
//   double ndraws1inv = 1.0/((double)ndraws1);
//   time_t sim_endtime;
   //===
   FILE *fptr_sp = NULL;
   FILE *fptr_readgibbs = NULL;


   //--- Reads the posterior estimates into constmodpars_ps->x_dv.
   fptr_sp = tzFopen(filename_sp,"r");
   rewind(fptr_sp);   //Must put the pointer at the beginning of the file.
   RefreshWPosteriorEstimates(constmodpars_ps, input_ps, govprob_ps, fptr_sp);
   tzFclose(fptr_sp);

   //--- Opens the file that contains Gibbs draws.
   fptr_readgibbs = tzFopen(filename_readgibbs,"r");


   _n = xdraw_dv->n;
   xdraw_dv->flag = V_DEF;  //Because the data will be read (at least one would hope).
   if (constgibbs_ps->indx_mhm)
   {
      //======= (1) Added 03/17/05.  Reading the posterior draws, computing cutval_logpost (which is set at 150.00 below Gibbmax of logPostKernel, hard coded in ReadGibbsDraws_Compute() in sze_comfuns.c) for MHM, and recording and exporting xmax_dv for the max logPosteriorKernal. =======
      if (constgibbs_ps->indxRecord)   //1: computing cutval_logpost and recording the max and min logPosteriorKernal; 0: no such recording.
      {
         constgibbs_ps->maxgibbs_logpostkernel = -MACHINEINFINITY;
         constgibbs_ps->mingibbs_logpostkernel = MACHINEINFINITY;
         //
         ndraws2 = 0;
         cntdraws = -1;
         rewind(fptr_readgibbs);   //Must put the pointer at the beginning of the file.
         for (; ;) {
            //--- Reads x_dv one row at a time from filename_readgibbs.  It's the user's responsiblity to match _n to the number of columns in filename_readgibbs.
            for (_i=0; _i<_n; _i++)
               if ( fscanf(fptr_readgibbs, " %lf ", xdraw_dv->v+_i) !=1 ) {
                  indic_end = 1;
                  break;
               }
            if (indic_end)  break;

            //=== Recording and exporting xdraw_dv for the max and min logPosteriorKernal.
            constgibbs_ps->logLHtimespriorpdf_lastval = value_logpostkernal_refresh(constmodpars_ps, constgibbs_ps, GIBBSDRAWS, govprob_ps, xdraw_dv);
            if (constgibbs_ps->logLHtimespriorpdf_lastval>constgibbs_ps->maxgibbs_logpostkernel)
            {
               CopyVector0(constgibbs_ps->xmax_dv, constgibbs_ps->x_dv);
               constgibbs_ps->maxgibbs_logpostkernel = constgibbs_ps->logLHtimespriorpdf_lastval;
            }
            if (constgibbs_ps->logLHtimespriorpdf_lastval<constgibbs_ps->mingibbs_logpostkernel)
            {
               CopyVector0(constgibbs_ps->xmin_dv, constgibbs_ps->x_dv);
               constgibbs_ps->mingibbs_logpostkernel = constgibbs_ps->logLHtimespriorpdf_lastval;
            }

            //=== Updating 1st and 2nd moments for the importance Gaussian pdf.
            CopyVector0(xzetadraw_dv, xdraw_dv);
            xzetadraw_dv->v[4] = square(xdraw_dv->v[4]);
            xzetadraw_dv->v[5] = square(xdraw_dv->v[5]);
            UpdateSumFor1st2ndMoments(constgibbs_ps->xmean_dv, constgibbs_ps->cholXcov_dm, xzetadraw_dv);

            //--- Printing out stuff for monitoring.
            ++cntdraws;
            if (++cntbuffer == nbuffer) {
               printf("\n\n============= Reading the draws from a file: 1st loop to (1) compute cutval_logpost, xmax_dv, and xmin_dv and (2) get 1st and 2nd moments for importance pdf under MHM: =================================");
               printf("\n Number of draws so far: %d\n", cntdraws);
               fflush(stdout);          //Flush the buffer to get out this message without delay.
               cntbuffer = 0;  //Cumulative.  Must be reset to 0.
            }

            //---
            ++ndraws2;
         }
         WriteVector(fptr_inpfor_maxGibbs, constgibbs_ps->xmax_dv, " %.16e ");
         if (indxGibbsmax)   //Gibbs max.
         {
            CopyVector0(constgibbs_ps->xzetapostest_dv, constgibbs_ps->xmax_dv);
            constgibbs_ps->xzetapostest_dv->v[4] = square(constmodpars_ps->x_dv->v[4]);
            constgibbs_ps->xzetapostest_dv->v[5] = square(constmodpars_ps->x_dv->v[5]);
         }
         else     //Posterior peak.
         {
            CopyVector0(constgibbs_ps->xzetapostest_dv, constmodpars_ps->x_dv);
            constgibbs_ps->xzetapostest_dv->v[4] = square(constmodpars_ps->x_dv->v[4]);
            constgibbs_ps->xzetapostest_dv->v[5] = square(constmodpars_ps->x_dv->v[5]);
         }
         //--- For MHM, records the cut-off values for logPostKernel to toss away the posterior draws with low values of logPostKernel.
         constgibbs_ps->cutval_logpost = constgibbs_ps->maxgibbs_logpostkernel - 150.00;
      }
      else
      {
         //======= (1) Reading the posterior draws of constgibbs_ps->x_dv and computing 1st and 2nd moments of xzetadraw_dv for weight Gaussian pdf using the MHM method. =======
         //--- Getting constgibbs_ps->xzetapostest_dv from constmodpars_ps->x_dv.  Commented out 03/17/05.
         // CopyVector0(constgibbs_ps->xzetapostest_dv, constmodpars_ps->x_dv);
         // constgibbs_ps->xzetapostest_dv->v[4] = square(constmodpars_ps->x_dv->v[4]);
         // constgibbs_ps->xzetapostest_dv->v[5] = square(constmodpars_ps->x_dv->v[5]);
         //--- Getting constgibbs_ps->xzetapostest_dv from constgibbs_ps->xmax_dv.  Added 03/17/05.
         if (indxGibbsmax)    //Gibbs max.
         {
            if ( !ReadVector_lf(fptr_inpfor_maxGibbs, constgibbs_ps->xmax_dv) )
            {
               printf("\n\nCheck the data matrix or vector in the input data file outdatainp_for_xmaxGibbs_filenametag.prn!\n");
               exit(EXIT_FAILURE);
            }
            constgibbs_ps->xmax_dv->flag = V_DEF;
            constgibbs_ps->maxgibbs_logpostkernel = value_logpostkernal_refresh(constmodpars_ps, constgibbs_ps, GIBBSDRAWS, govprob_ps, constgibbs_ps->xmax_dv);
            //+
            CopyVector0(constgibbs_ps->xzetapostest_dv, constgibbs_ps->xmax_dv);
            constgibbs_ps->xzetapostest_dv->v[4] = square(constmodpars_ps->x_dv->v[4]);
            constgibbs_ps->xzetapostest_dv->v[5] = square(constmodpars_ps->x_dv->v[5]);
         }
         else   //Posterior peak.
         {
            CopyVector0(constgibbs_ps->xzetapostest_dv, constmodpars_ps->x_dv);
            constgibbs_ps->xzetapostest_dv->v[4] = square(constmodpars_ps->x_dv->v[4]);
            constgibbs_ps->xzetapostest_dv->v[5] = square(constmodpars_ps->x_dv->v[5]);
         }
         //+
         ndraws2 = 0;
         cntdraws = -1;
         rewind(fptr_readgibbs);   //Must put the pointer at the beginning of the file.
         for (; ;) {
            //--- Reads x_dv one row at a time from filename_readgibbs.  It's the user's responsiblity to match _n to the number of columns in filename_readgibbs.
            for (_i=0; _i<_n; _i++)
               if ( fscanf(fptr_readgibbs, " %lf ", xdraw_dv->v+_i) !=1 ) {
                  indic_end = 1;
                  break;
               }
            if (indic_end)  break;

            //--- Updating 1st and 2nd moments for the importance Gaussian pdf.
            CopyVector0(xzetadraw_dv, xdraw_dv);
            xzetadraw_dv->v[4] = square(xdraw_dv->v[4]);
            xzetadraw_dv->v[5] = square(xdraw_dv->v[5]);
            UpdateSumFor1st2ndMoments(constgibbs_ps->xmean_dv, constgibbs_ps->cholXcov_dm, xzetadraw_dv);

            //--- Printing out stuff for monitoring.
            ++cntdraws;
            if (++cntbuffer == nbuffer) {
               printf("\n\n============= Reading the draws from a file: 1st loop to read xmax_dv and get 1st and 2nd moments for importance pdf under MHM: =================================");
               printf("\n Number of draws so far: %d\n", cntdraws);
               fflush(stdout);          //Flush the buffer to get out this message without delay.
               cntbuffer = 0;  //Cumulative.  Must be reset to 0.
            }


            //---
            ++ndraws2;
         }
      }
      //End of the short loop.
      constgibbs_ps->ndraws2 = ndraws2;
      if (!ndraws2) {
         printf(".../swz_comfuns.c/ReadGibbsDraws_Compute(): the data file %s is empty!\n", filename_readgibbs);
         exit(EXIT_FAILURE);
      }
      // fprintf(fptr_debug, "Draw x_dv:\n");
      // WriteVector(fptr_debug, constgibbs_ps->x_dv, " %.16e ");
      ndraws2inv = 1.0/(double)ndraws2;
      //--- Converting 1st moment to mean.
      ScalarTimesVector(constgibbs_ps->xmean_dv, ndraws2inv, constgibbs_ps->xmean_dv, 0.0);
      //--- Converting 2nd moment to covariance matrix.
      if (constgibbs_ps->indxMHMPeak) {
         //Centers the importance pdf at the posterior peak.
         VectorTimesSelf(constgibbs_ps->cholXcov_dm, constgibbs_ps->xzetapostest_dv, 1.0, ndraws2inv, 'U'); //X'X/n + xhat*xhat';
         VectorTimesVector(constgibbs_ps->cholXcov_dm, constgibbs_ps->xzetapostest_dv, constgibbs_ps->xmean_dv, -1.0, 1.0);  //X'X/n + xhat*xhat' - xhat*xmean';
         VectorTimesVector(constgibbs_ps->cholXcov_dm, constgibbs_ps->xmean_dv, constgibbs_ps->xzetapostest_dv, -1.0, 1.0);  //X'X/n + xhat*xhat' - xhat*xmean' - xmean*xhat';
         //--- Adjusts the covariance matrix.
         ScalarTimesMatrix(constgibbs_ps->cholXcov_dm, constgibbs_ps->sf_mhm, constgibbs_ps->cholXcov_dm, 0.0);
         constgibbs_ps->cholXcov_dm->flag = M_SU;  //Getting rid of the flag M_GE to make sure the upper symmetric part is used later.
      }
      else {
         //Centers the importance pdf at the posterior mean -- a traditional method.
         VectorTimesSelf(constgibbs_ps->cholXcov_dm, constgibbs_ps->xmean_dv, -1.0, ndraws2inv, 'U'); //Converting it to covariance matrix.
         //--- Adjusts the covariance matrix.
         ScalarTimesMatrix(constgibbs_ps->cholXcov_dm, constgibbs_ps->sf_mhm, constgibbs_ps->cholXcov_dm, 0.0);
         constgibbs_ps->cholXcov_dm->flag = M_SU;  //Getting rid of the flag M_GE to make sure the upper symmetric part is used later.
      }
      //--- Getting the upper Choleski of the covariance matrix.
      if (errflag=chol(constgibbs_ps->cholXcov_dm, constgibbs_ps->cholXcov_dm, 'U')) {
         //cholXcov_dm is now upper triangular Choleski so that chol'chol = Xcov.
         printf("\n.../swz_comfuns.c/GibbsSimulations(): No chol decomp for cholXcov_dm with error flag %d\n", errflag);
         exit(EXIT_FAILURE);
      }
      //--- Taking log of pcuts_dv and computing a normal pdf term outside the second loop.
      for (_i=constgibbs_ps->pcuts_dv->n-1; _i>=0; _i--) {
         constgibbs_ps->logpcuts_dv->v[_i] = log(constgibbs_ps->pcuts_dv->v[_i]);
         constgibbs_ps->logpcuts_dv->flag = V_DEF;
      }
      constgibbs_ps->neghalfklog2pi_minus_dethalfcholXcov = NHLOG2PI*(double)constgibbs_ps->x_dv->n - tracelogfabs(constgibbs_ps->cholXcov_dm);
                                   //-(k/2)*log(2pi) - log(det(cholXcov)).


      //======= (2) Computing the normalizing constant q for each cut-off p value of Gaussain weight function.
      //InitializeConstantVector_lf(constgibbs_ps->logqvals_dv, 0.0);   //Cumulative.
      if (1)  ftd_computeqvals(constmodpars_ps, govprob_ps, constgibbs_ps);  //$$$$$Hard coded.$$$$$$$


      //======= (3) Reading the posterior draws of constgibbs_ps->x_dv and computing objects like marginal data density. =======
      constgibbs_ps->maxgibbs_logpostkernel = -MACHINEINFINITY;
      constgibbs_ps->mingibbs_logpostkernel = MACHINEINFINITY;
      cntbuffer = 0;
      indic_end = 0;
      rewind(fptr_readgibbs);   //Must put the pointer at the beginning of the file.
      for (drawi=0; drawi<ndraws2; drawi++)
      {
         //--- Reads x_dv one row at a time from filename_readgibbs.  It's the user's responsiblity to match _n to the number of columns in filename_readgibbs.
         for (_i=0; _i<_n; _i++)
            if ( fscanf(fptr_readgibbs, " %lf ", xdraw_dv->v+_i) !=1 ) {
               indic_end = 1;
               break;
            }
         if (indic_end) {
            printf(".../swz_comfuns.c/ReadGibbsDraws_Compute(): at 2nd loop and draw %d, we cannot propertly read the data from %s", drawi, filename_readgibbs);
            exit(EXIT_FAILURE);
         }

         //--- Updating towards getting the log marginal likelihood using the MHM method.
         constgibbs_ps->logLHtimespriorpdf_lastval = value_logpostkernal_refresh(constmodpars_ps, constgibbs_ps, GIBBSDRAWS, govprob_ps, xdraw_dv);
         CopyVector0(xzetadraw_dv, xdraw_dv);
         xzetadraw_dv->v[4] = square(xdraw_dv->v[4]);
         xzetadraw_dv->v[5] = square(xdraw_dv->v[5]);
         ftd_UpdatelogMLH_mhm(constgibbs_ps);
         //
         if (!constgibbs_ps->indxRecord && (constgibbs_ps->logLHtimespriorpdf_lastval>constgibbs_ps->maxgibbs_logpostkernel)) {
            CopyVector0(constgibbs_ps->xmax_dv, constgibbs_ps->x_dv);
            constgibbs_ps->maxgibbs_logpostkernel = constgibbs_ps->logLHtimespriorpdf_lastval;
         }
         if (!constgibbs_ps->indxRecord && (constgibbs_ps->logLHtimespriorpdf_lastval<constgibbs_ps->mingibbs_logpostkernel)) {
            CopyVector0(constgibbs_ps->xmin_dv, constgibbs_ps->x_dv);
            constgibbs_ps->mingibbs_logpostkernel = constgibbs_ps->logLHtimespriorpdf_lastval;
         }
         //--- Printing out stuff for monitoring.
         if (++cntbuffer == nbuffer) {
            printf("\n\n============= Reading the draws from a file: 2nd loop to finish computing marginal data density for MHM: =================================");
            printf("\n (A) Total number of kept draws: %d\n (B) Number of draws so far: %d\n",
                     ndraws2, drawi);
            fflush(stdout);          //Flush the buffer to get out this message without delay.
            cntbuffer = 0;  //Cumulative.  Must be reset to 0.
         }
      }
   }


   //======= (4) Drawing perceived long-run excess unemployment under Ramsey policy. =======
   if (constgibbs_ps->indx_lrsacrats) {
      rewind(fptr_lreutradeoffs);    //Must put the pointer at the beginning of the file.

      //=== First draw is at the posterior estimate.
      //--- Refreshing the posterior estimate and Zpredtran_dm.
      value_logpostkernal_refresh(constmodpars_ps, (struct TSconstgibbs_tag *)NULL, POINTEST, govprob_ps, constmodpars_ps->x_dv);
      //--- At time 0.
      sumofcoeff_pi = z10_dv->v[0] + z10_dv->v[1] + z10_dv->v[3];
      sumofcoeff_u = 1.0 - z10_dv->v[2] - z10_dv->v[4];
      lreudraw_dv->v[0] = (constmodpars_ps->govpitarget*sumofcoeff_pi + z10_dv->v[5]) / sumofcoeff_u - constmodpars_ps->ustar;
      //--- From 1 to T-1.
      for (ti=0; ti<fss-1; ti++) {
         sumofcoeff_pi = Zpredtran_dm->M[mos(0,ti,kx)] + Zpredtran_dm->M[mos(1,ti,kx)] + Zpredtran_dm->M[mos(3,ti,kx)];
         sumofcoeff_u = 1.0 - Zpredtran_dm->M[mos(2,ti,kx)] - Zpredtran_dm->M[mos(4,ti,kx)];
         lreudraw_dv->v[ti+1] = (constmodpars_ps->govpitarget*sumofcoeff_pi + Zpredtran_dm->M[mos(5,ti,kx)]) / sumofcoeff_u - constmodpars_ps->ustar;
      }
      WriteVector(fptr_lreutradeoffs, lreudraw_dv, " %.16e ");

      //=== The rest of draws are from Gibbs.
      ndraws2 = 0;
      cntdraws = -1;
      rewind(fptr_readgibbs);   //Must put the pointer at the beginning of the file.
      for (; ;) {
         //--- Reads x_dv one row at a time from filename_readgibbs.  It's the user's responsiblity to match _n to the number of columns in filename_readgibbs.
         for (_i=0; _i<_n; _i++)
            if ( fscanf(fptr_readgibbs, " %lf ", xdraw_dv->v+_i) !=1 ) {
               indic_end = 1;
               break;
            }
         if (indic_end)  break;

         //--- Refreshing the posterior estimate and Zpredtran_dm.
         value_logpostkernal_refresh(constmodpars_ps, constgibbs_ps, GIBBSDRAWS, govprob_ps, xdraw_dv);
         //--- At time 0.
         sumofcoeff_pi = z10_dv->v[0] + z10_dv->v[1] + z10_dv->v[3];
         sumofcoeff_u = 1.0 - z10_dv->v[2] - z10_dv->v[4];
         lreudraw_dv->v[0] = (constgibbs_ps->govpitarget*sumofcoeff_pi + z10_dv->v[5]) / sumofcoeff_u - constgibbs_ps->ustar;
         //--- From 1 to T-1.
         for (ti=0; ti<fss-1; ti++) {
            sumofcoeff_pi = Zpredtran_dm->M[mos(0,ti,kx)] + Zpredtran_dm->M[mos(1,ti,kx)] + Zpredtran_dm->M[mos(3,ti,kx)];
            sumofcoeff_u = 1.0 - Zpredtran_dm->M[mos(2,ti,kx)] - Zpredtran_dm->M[mos(4,ti,kx)];
            lreudraw_dv->v[ti+1] = (constgibbs_ps->govpitarget*sumofcoeff_pi + Zpredtran_dm->M[mos(5,ti,kx)]) / sumofcoeff_u - constgibbs_ps->ustar;
         }
         WriteVector(fptr_lreutradeoffs, lreudraw_dv, " %.16e ");

         //--- Printing out stuff for monitoring.
         ++cntdraws;
         if (++cntbuffer == nbuffer) {
            printf("\n\n============= Reading the draws from a file: getting the long-run excess unemployment tradeoffs under Ramsey: =================================");
            printf("\n Number of draws so far: %d\n", cntdraws);
            fflush(stdout);          //Flush the buffer to get out this message without delay.
            cntbuffer = 0;  //Cumulative.  Must be reset to 0.
         }
         //---
         ++ndraws2;
      }
      //End of the short loop.
      constgibbs_ps->ndraws2 = ndraws2;
      if (!ndraws2) {
         printf(".../swz_comfuns.c/ReadGibbsDraws_Compute(): the data file %s is empty!\n", filename_readgibbs);
         exit(EXIT_FAILURE);
      }
   }


   //===
   tzFclose(fptr_readgibbs);
}

#undef NHLOG2PI
//+
#undef N1STSET
#undef N2NDSET
#undef N3RDSET
#undef AFORMAT



//************************************************************
// IV. For forecasts of inflation and unemployment only.
//************************************************************
typedef struct TSforecasts_tag {
//           TSdmatrix *DynForesuinfs_dm;  //nsteps-by-2 where 2 concerns inflation and unemployment.
        int nsteps;   //Number of forecast steps in simulating inflation and unemployment.
        TSdvector *foreuinf_k_dv;   //4-by-1.  k-step forecast of expected inflation, government inflation (= expected inflation in this case), inflation, and unemployment.

        //--- For government's updated inflation given the forecast data (1 data point at a time).
        TSkalcvfurw *kalcvfurw1datapoint_ps;  //For kalcvfurw1datapoint_ps->Xrhtran_dm have variable in the order of, say, pi_{t}, pi_{t-1}, u_{t-1}, pi_{t-2}, u_{t-2}, and 1.0 for indxWhitchModel == CGP2TP1TI0.

        //--- Probability of nonconvergence.

} TSforecasts;

static struct TSforecasts_tag *CreateTSforecasts(int kx, int tv)
{
   //kth-step forecasts of gov_inf, exp_inf (=gov_inf), inf_t, and u_t.
   struct TSforecasts_tag *forecasts_ps = tzMalloc(1, struct TSforecasts_tag);

//   forecasts_ps->DynForesuinfs_dm = CreateMatrix_lf(nsteps, 2);
   forecasts_ps->foreuinf_k_dv = CreateVector_lf(4);   //gov_inf_t, exp_inf_t (=gov_inf_t in most cases), inf_t, u_t.
   forecasts_ps->kalcvfurw1datapoint_ps = CreateTSkalcvfurw(kalcvf_urw, 1, kx, tv);

   return (forecasts_ps);
};
//---
static struct TSforecasts_tag *DestroyTSforecasts(struct TSforecasts_tag *forecasts_ps)
{
   if (forecasts_ps) {
      DestroyVector_lf(forecasts_ps->foreuinf_k_dv);
      DestroyTSkalcvfurw(forecasts_ps->kalcvfurw1datapoint_ps);

      //===
      free(forecasts_ps);
      return ((struct TSforecasts_tag *)NULL);
   }
   else  return (forecasts_ps);

}

static void InitializeKalman1datapoint(TSkalcvfurw *kalcvfurw1datapoint_ps, int locti, struct TSgovprob_tag *govprob_ps, void *constmod_ps, int pointflag)
{
   //locti: location of the forecast date locti+ki for ki=0 (which is one step after the sample end date).
   struct TSconstmodpars_tag *constmodpars_ps;
   struct TSconstgibbs_tag *constgibbs_ps;
   //
   TSkalcvfurw *kalcvfurw_ps = govprob_ps->kalcvfurw_ps;
   TSdmatrix *Xrhtran_dm = kalcvfurw_ps->Xrhtran_dm;  //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.
               //For the benchmark case, 6 variables for Xrhtran_dm in the order of pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and constant terms.
   int nrhvars = Xrhtran_dm->nrows;


   //=== Initialization.
   if (pointflag & POINTEST) {
      constmodpars_ps = (struct TSconstmodpars_tag *)constmod_ps;
      kalcvfurw1datapoint_ps->sigmasq = 1.0/constmodpars_ps->gzeta;
      if (locti==0) {
         CopyVector0(kalcvfurw1datapoint_ps->z10_dv, constmodpars_ps->z10_dv);
         CopyMatrix0(kalcvfurw1datapoint_ps->P10_dm, constmodpars_ps->P10_dm);
      }
      else {
         CopySubmatrix2vector(kalcvfurw1datapoint_ps->z10_dv, 0, kalcvfurw_ps->Zpredtran_dm, 0, locti-1, nrhvars);
         CopyMatrix0(kalcvfurw1datapoint_ps->P10_dm, kalcvfurw_ps->Ppred_dc->C[locti-1]);
      }
      CopyMatrix0(kalcvfurw1datapoint_ps->V_dm, constmodpars_ps->V_dm);
   }
   else {
      //For GIBBSDRAWS (Gibbs draws)
      constgibbs_ps = (struct TSconstgibbs_tag *)constmod_ps;
      kalcvfurw1datapoint_ps->sigmasq = 1.0/constgibbs_ps->gzeta;
      if (locti==0) {
         CopyVector0(kalcvfurw1datapoint_ps->z10_dv, constgibbs_ps->z10_dv);
         CopyMatrix0(kalcvfurw1datapoint_ps->P10_dm, constgibbs_ps->P10_dm);
      }
      else {
         CopySubmatrix2vector(kalcvfurw1datapoint_ps->z10_dv, 0, kalcvfurw_ps->Zpredtran_dm, 0, locti-1, nrhvars);
         CopyMatrix0(kalcvfurw1datapoint_ps->P10_dm, kalcvfurw_ps->Ppred_dc->C[locti-1]);
      }
      CopyMatrix0(kalcvfurw1datapoint_ps->V_dm, constgibbs_ps->V_dm);
   }
   //=== Order matters!  The following stores lagged values that will be used in a loop by ftd_dynfores_uinfs()!
   if (govprob_ps->indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0)) {
      //Classical regression: Xrhtran_dm in the order of pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and constant terms.
      //Keynesian regression: Xrhtran_dm in the order of u_t, u_{t-1}, pi_{t-1}, u_{t-2}, pi_{t-2}, and constant terms.
      kalcvfurw1datapoint_ps->ylhtran_dv->v[0] = Xrhtran_dm->M[mos(2,locti,nrhvars)];   //u_{locti-1} for classical; piu_{locti-1} for Keynesian.
      kalcvfurw1datapoint_ps->ylhtran_dv->flag = V_DEF;
      //+
      kalcvfurw1datapoint_ps->Xrhtran_dm->M[0] = Xrhtran_dm->M[mos(1,locti,nrhvars)];  //pi_{locti-1} for classical; u_{locti-1} for Keynesian.
      kalcvfurw1datapoint_ps->Xrhtran_dm->M[1] = Xrhtran_dm->M[mos(3,locti,nrhvars)];  //pi_{locti-2} for classical; u_{locti-2} for Keynesian.
      kalcvfurw1datapoint_ps->Xrhtran_dm->M[2] = Xrhtran_dm->M[mos(4,locti,nrhvars)];  //u_{locti-2} for classical;  pi_{locti-2} for Keynesian.
      kalcvfurw1datapoint_ps->Xrhtran_dm->M[5] = 1.0;
      kalcvfurw1datapoint_ps->Xrhtran_dm->flag = M_GE;
              //Initial values for M[3] and M[4] are not needed but will be filled in the next loop in ftd_dynfores_uinfs().
   }
   else  {
      printf(".../swz_comfuns.c/InitializeKalman1datapoint): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d!", CGP2TP1TI0, KGP2TP1TI0);
      exit(EXIT_FAILURE);
   }
}


//----------------------------
// Dynamic forecasts of inflation and unemployment.
//----------------------------
void ftd_dynfores_uinfs(FILE *fptr_fcsts, int nsteps, int locti, struct TSconstmodpars_tag *constmodpars_ps, struct TSgovprob_tag *govprob_ps, struct TSconstgibbs_tag *constgibbs_ps, int pointflag, int indxWhichModel, int store_sims)
{
   //locti: location of the forecast date locti+ki for ki=0 (which is one step after the sample end date).
   //Structures govprob_ps and constmodpars_ps must be inputs.  Depending on pointflag, constgibbs_ps may or may not be NULL.
   //TSdvector *inivars_dv: initial conditions or variables, inivars_dv, contain variables in the order of x_{t-1}, x_{t-2}, pi_{t-1}, u_{t-1}, pi_{t-2}, and u_{t-2}.
   int ki;
   int chosenvar;   //This chosen variable is manually keyed in below.
   int  cntsave = 0,  //Cumulative.  Must be set to 0.
        napartsave;
   double x_1, x_2, pi_1, u_1, pi_2, u_2; //Initial conditions: x_{t-1}, x_{t-2}, pi_{t-1}, u_{t-1}, pi_{t-2}, and u_{t-2}.
   double pi_0, u_0;
   //--- Initialization. The order matters!
   TSkalcvfurw *kalcvfurw_ps = govprob_ps->kalcvfurw_ps;
   TSdvector *ylhtran_dv = kalcvfurw_ps->ylhtran_dv;
   TSdmatrix *Xrhtran_dm = kalcvfurw_ps->Xrhtran_dm;  //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.
               //For classical regression, the order of 6 variables is pi_t, pi_{t-1}, u_{t-1}, pi_{t-2}, u_{t-2}, and 1.
               //For Keynesian regression, the order of 6 variables is u_t, u_{t-1}, pi_{t-1}, u_{t-2}, pi_{t-2}, and 1.
   int nrhvars = Xrhtran_dm->nrows;
   TSdvector *expinf_dv = govprob_ps->expinf_dv;  //Will be updated by govprob_ps->govestimator() below.
   //
   TSdvector *x_dv;   //Will be pointed to constmodpars_ps->x_dv or constgibbs_ps->x_dv, depending on pointflag.
   TSkalcvfurw *kalcvfurw1datapoint_ps;   //Will be pointed to forecasts_ps->kalcvfurw1datapoint_ps.
   //
   TSdvector *sclhistshocks_dv;
   TSdmatrix *StrResidstran_dm;
   //===
   struct TSforecasts_tag *forecasts_ps = CreateTSforecasts(nrhvars, govprob_ps->indx_tvsigmasq);
   TSdvector *gov_expinf_t_dv = CreateVector_lf(2); //2-by-1.  Will be copied to the first 2 elements of forecasts_ps->foreuinf_k_dv.
   TSdvector *onepath_dv = NULL;  //Allocated when store_sims is on.
   TSdvector *z10old_dv = CreateVector_lf(nrhvars);
   TSdmatrix *P10old_dm = CreateConstantMatrix_lf(nrhvars, nrhvars, 0.0);
   TSdmatrix *StrResids_dm = NULL;
   // TSdvector *XrhtranOld_dv = NULL;

   if (store_sims) {
      onepath_dv = CreateVector_lf(nsteps);
      onepath_dv->flag = V_DEF;
      //=== 0: gov inflation; 1: expected inflation; 2: inflation; 3: unemployment.
      if ((chosenvar = 2) >= forecasts_ps->foreuinf_k_dv->n)  fn_DisplayError(".../swz_comfuns.c/ftd_dynfores_uinfs(): the location of the chosen variable, chosenvar, must be within the length of forecasts_ps->foreuinf_k_dv");
   }


   if (locti<0 || locti >= constmodpars_ps->fss)  fn_DisplayError(".../swz_comfuns.c/ftd_dynfores_uinfs(): the forecast date, locti, must be between 0 and fss-1 inclusive");
   // XrhtranOld_dv = CreateVector_lf(nrhvars);
   kalcvfurw1datapoint_ps = forecasts_ps->kalcvfurw1datapoint_ps;


   if (indxWhichModel & (CGP2TP1TI0 | KGP2TP1TI0)) {
      if (pointflag & POINTEST) {
         //Updating government's decision rules through t=1:T.
         //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.
         //For classical regression, the order of 6 variables is pi_t, pi_{t-1}, u_{t-1}, pi_{t-2}, u_{t-2}, and 1.
         //For Keynesian regression, the order of 6 variables is u_t, u_{t-1}, pi_{t-1}, u_{t-2}, pi_{t-2}, and 1.
         x_dv = constmodpars_ps->x_dv;   //No updating, but will be used below.

         //--- Converting free paramaters to model (meaningful) parameters.
         ftd_freepars2modpars((void *)constmodpars_ps, POINTEST);
         //--- Scaling down V to get a self-confirming equilibrium.
         if (constmodpars_ps->scl_epsilon != 1.0) {
            ScalarTimesMatrix(constmodpars_ps->V_dm, constmodpars_ps->scl_epsilon, constmodpars_ps->V_dm, 0.0);
            ftd_modpars2freepars((void *)constmodpars_ps, POINTEST);   //Converting back to x_dv to be consistent with the changed V.
         }
         //--- Normalization.
         ftd_normalizationforgov((void *)constmodpars_ps, POINTEST);
         //--- Refreshed by solving the government's problem for (1) expected inflation and (2) government's inflation x_t (which controls at least part of the inflation process).
         govprob_ps->govestimator(govprob_ps, (void *)constmodpars_ps, POINTEST);
         //--- Initial conditions before locti.
         if (indxWhichModel & CGP2TP1TI0) {
            //--- Classical regression: ylhtran_dv is u_t; 6 variables for Xrhtran_dm in the order of pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and Constant.
            x_1 = expinf_dv->v[locti];
            x_2 = (locti<1) ? Xrhtran_dm->M[mos(1,locti,nrhvars)] : expinf_dv->v[locti-1];
            pi_1 = Xrhtran_dm->M[mos(1,locti,nrhvars)];
            u_1 = Xrhtran_dm->M[mos(2,locti,nrhvars)];
            pi_2 = Xrhtran_dm->M[mos(3,locti,nrhvars)];
            u_2 = Xrhtran_dm->M[mos(4,locti,nrhvars)];
         }
         else if (indxWhichModel & KGP2TP1TI0) {
            //--- Keynesian regression: ylhtran_dv is pi_t; 6 variables for Xrhtran_dm in the order of u_t, u_{t-1}, pi_{t-1}, u_{t-2}, pi_{t-2}, and Constant.
            x_1 = expinf_dv->v[locti];
            x_2 = (locti<1) ? Xrhtran_dm->M[mos(2,locti,nrhvars)] : expinf_dv->v[locti-1];
            pi_1 = Xrhtran_dm->M[mos(2,locti,nrhvars)];
            u_1 = Xrhtran_dm->M[mos(1,locti,nrhvars)];
            pi_2 = Xrhtran_dm->M[mos(4,locti,nrhvars)];
            u_2 = Xrhtran_dm->M[mos(3,locti,nrhvars)];
         }
         //+
         // CopyVector0(z10old_dv, constmodpars_ps->z10_dv);
         // CopyMatrix0(P10old_dm, constmodpars_ps->P10_dm);

         //DDDDDDDDDDDDDDDDebug.
         // fprintf(fptr_debug, "govprob_ps->expinf_dv:\n");
         // WriteVector(fptr_debug, expinf_dv, " %.16e ");
         // fprintf(fptr_debug, "gov_expinf_t_dv at each time:\n");

         #if defined (DEGUG_SWZ)
            if (nsteps+locti > constmodpars_ps->fss) {
               printf(".../swz_comfuns.c/ftd_dynfores_uinfs(): In a debug mode DEGUG_SWZ, nsteps+locti, %d, cannot exceed the effective sample size %d", nsteps+locti, constmodpars_ps->fss);
               exit(EXIT_FAILURE);
            }
         #endif

         if (!constmodpars_ps->indxCnterFact) {     //No counterfactual exercises.
            //======= Forecasts for ki=locti+0:locti+nsteps-1. ========
            InitializeKalman1datapoint(kalcvfurw1datapoint_ps, locti, govprob_ps, (void *)constmodpars_ps, POINTEST);  //At ki=-1.
            //--- z_{t|t-1} at 73:12 (163).
            // kalcvfurw1datapoint_ps->z10_dv->v[0] = -8.2576116831132182e+000;
            // kalcvfurw1datapoint_ps->z10_dv->v[1] = -5.9170702577384295e+000;
            // kalcvfurw1datapoint_ps->z10_dv->v[2] = -4.0719350360103029e-001;
            // kalcvfurw1datapoint_ps->z10_dv->v[3] =  1.3267510383654432e+001;
            // kalcvfurw1datapoint_ps->z10_dv->v[4] =  5.5790782283342677e-001;
            // kalcvfurw1datapoint_ps->z10_dv->v[5] =  2.1165972013644286e+001;
            //
            // fprintf(fptr_debug, "\n**** Sum of pi coefficients:\n");
            for (ki=0; ki<nsteps; ki++) {
               // fprintf(fptr_debug, " %g\n", kalcvfurw1datapoint_ps->z10_dv->v[0]+kalcvfurw1datapoint_ps->z10_dv->v[1]+kalcvfurw1datapoint_ps->z10_dv->v[3]);
               //--- Government's inflation and expected inflation.
               forecasts_ps->foreuinf_k_dv->v[0] = forecasts_ps->foreuinf_k_dv->v[1] = x_1;
               //--- Inflation forecasts at ki.
               #if defined(DEGUG_SWZ)
                  if (indxWhichModel & CGP2TP1TI0)         //Classical direction of fit.
                     forecasts_ps->foreuinf_k_dv->v[2] = pi_0 = x_1 + Xrhtran_dm->M[mos(0, locti+ki, nrhvars)] - expinf_dv->v[locti+ki];
                  else if (indxWhichModel & KGP2TP1TI0)    //Keynesian direction of fit.
                     forecasts_ps->foreuinf_k_dv->v[2] = pi_0 = x_1 + ylhtran_dv->v[locti+ki] - expinf_dv->v[locti+ki];
               #else
                  forecasts_ps->foreuinf_k_dv->v[2] = pi_0 = constmodpars_ps->draw_k_fcstinf = x_1 + sqrt(1.0/constmodpars_ps->gzeta2)*StandardNormalDouble();
                                  //constmodpars_ps->draw_k_fcstinf is added 03/17/05.
               #endif
               //--- Unemployment forecasts at ki.
               #if defined (DEGUG_SWZ)
                  if (indxWhichModel & CGP2TP1TI0)      //Classical direction of fit.
                     forecasts_ps->foreuinf_k_dv->v[3] = u_0 = x_dv->v[0] + constmodpars_ps->gtau1*u_1
                        + constmodpars_ps->gtheta0*(pi_0 - x_1) + constmodpars_ps->gtheta1*(pi_1 - x_2) + ylhtran_dv->v[locti+ki] - constmodpars_ps->upred_dv->v[locti+ki];
                  else if (indxWhichModel & KGP2TP1TI0)  //Keynesian direction of fit.
                     forecasts_ps->foreuinf_k_dv->v[3] = u_0 = x_dv->v[0] + constmodpars_ps->gtau1*u_1
                        + constmodpars_ps->gtheta0*(pi_0 - x_1) + constmodpars_ps->gtheta1*(pi_1 - x_2) + Xrhtran_dm->M[mos(0, locti+ki, nrhvars)] - constmodpars_ps->upred_dv->v[locti+ki];
               #else
                  forecasts_ps->foreuinf_k_dv->v[3] = u_0 = x_dv->v[0] + constmodpars_ps->gtau1*u_1
                     + constmodpars_ps->gtheta0*(pi_0 - x_1) + constmodpars_ps->gtheta1*(pi_1 - x_2) + sqrt(1.0/constmodpars_ps->gzeta1)*StandardNormalDouble();
               #endif

               //------- For the next period or loop. --------
               //=== Updates x_t and expinf_t for the government's inflation and expected inflation, which are the same if govprob_ps->grho1=govprob_ps->grho2=0.0.
               //=== The following order matters!
               // CopySubmatrix2vector(XrhtranOld_dv, 0, kalcvfurw1datapoint_ps->Xrhtran_dm, 0, 0, nrhvars);
               //--- Classical regression: ylhtran_dv is u_t; 6 variables for Xrhtran_dm in the order of pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and Constant.
               //--- Keynesian regression: ylhtran_dv is pi_t; 6 variables for Xrhtran_dm in the order of u_t, u_{t-1}, pi_{t-1}, u_{t-2}, pi_{t-2}, and Constant.
               kalcvfurw1datapoint_ps->Xrhtran_dm->M[4] = kalcvfurw1datapoint_ps->Xrhtran_dm->M[2];  //[4] is now u_{ki-2} for classical.   [4] is now pi_{ki-2} for Keynesian.
               kalcvfurw1datapoint_ps->Xrhtran_dm->M[3] = kalcvfurw1datapoint_ps->Xrhtran_dm->M[1];  //[3] is now pi_{ki-2} for classical.  [3] is now u_{ki-2} for Keynesian.
               kalcvfurw1datapoint_ps->Xrhtran_dm->M[2] = kalcvfurw1datapoint_ps->ylhtran_dv->v[0];  //[2] is now u_{ki-1} for classical.   [2] is now pi_{ki-1} for Keynesian.
               kalcvfurw1datapoint_ps->Xrhtran_dm->M[1] = kalcvfurw1datapoint_ps->Xrhtran_dm->M[0];  //[1] is now pi_{ki-1} for classical.  [1] is now u_{ki-1} for Keynesian.
               if (indxWhichModel & CGP2TP1TI0) {
                  kalcvfurw1datapoint_ps->Xrhtran_dm->M[0] = pi_0;                                      //[0] is now pi_{ki}.
                  //+ Must be after all the above.
                  kalcvfurw1datapoint_ps->ylhtran_dv->v[0] = u_0;
               }
               else if (indxWhichModel & KGP2TP1TI0) {
                  kalcvfurw1datapoint_ps->Xrhtran_dm->M[0] = u_0;                                      //[0] is now u_{ki}.
                  //+ Must be after all the above.
                  kalcvfurw1datapoint_ps->ylhtran_dv->v[0] = pi_0;
               }
               //--- Getting government's inflation and expected inflation for the next ki (or ki++), stored in gov_expinf_t_dv.
               govexpinf_1update(gov_expinf_t_dv, kalcvfurw1datapoint_ps, govprob_ps, indxWhichModel);
               // if (govprob_ps->gensys_ps->eu_dv->v[0] !=1 || govprob_ps->gensys_ps->eu_dv->v[1] !=1) {
               //    //=== If non-unique or non-existent, force no Kalman updating.
               //    CopyVector0(kalcvfurw1datapoint_ps->z10_dv, z10old_dv);
               //    CopyMatrix0(kalcvfurw1datapoint_ps->P10_dm, P10old_dm);
               //    govexpinf_1update(gov_expinf_t_dv, kalcvfurw1datapoint_ps, govprob_ps, indxWhichModel);
               // }
               //--- Added 03/17/05.
               if (ki==nsteps-1)  CopyVector0(constmodpars_ps->k_fcstinf_z10_dv, kalcvfurw1datapoint_ps->zupdate_dv);


               //DDDDDDDDDDDDDDDDebug.
               // if (govprob_ps->gensys_ps->eu_dv->v[0] !=1 || govprob_ps->gensys_ps->eu_dv->v[1] !=1) {
               //    fprintf(fptr_debug, "\nAt forecast date ki = %d.\n", ki);
               //    fprintf(fptr_debug, "govprob_ps->gensys_ps->eu_dv: [%g  %g]\n", govprob_ps->gensys_ps->eu_dv->v[0], govprob_ps->gensys_ps->eu_dv->v[1]);
               //    fprintf(fptr_debug, "kalcvfurw1datapoint_ps->ylhtran_dv:\n");
               //    WriteVector(fptr_debug, kalcvfurw1datapoint_ps->ylhtran_dv, " %10.4f ");
               //    fprintf(fptr_debug, "Xrhtran_dm:\n");
               //    WriteMatrix(fptr_debug, kalcvfurw1datapoint_ps->Xrhtran_dm, " %.16e ");
               //    fprintf(fptr_debug, "z10_dv:\n");
               //    WriteVector(fptr_debug, kalcvfurw1datapoint_ps->z10_dv, " %.16e ");
               //    fprintf(fptr_debug, "zupdate_dv:\n");
               //    WriteVector(fptr_debug, kalcvfurw1datapoint_ps->zupdate_dv, " %10.4f ");
               //    fprintf(fptr_debug, "P10_dm:\n");
               //    WriteMatrix(fptr_debug, kalcvfurw1datapoint_ps->P10_dm, " %10.4f ");
               //    fprintf(fptr_debug, "V_dm:\n");
               //    WriteMatrix(fptr_debug, constmodpars_ps->V_dm, " %10.4f ");
               // }

               //=== A priori constraints on the bounds of inflation.
               // if (govprob_ps->gensys_ps->eu_dv->v[0] !=1 || govprob_ps->gensys_ps->eu_dv->v[1] !=1) {
               // if (gov_expinf_t_dv->v[0] < -10.0)  gov_expinf_t_dv->v[0] = gov_expinf_t_dv->v[1] = -10.0;
               // if (gov_expinf_t_dv->v[0] > 10.0)  gov_expinf_t_dv->v[0] = gov_expinf_t_dv->v[1] = 10.0;
               // }
               /**
               if (govprob_ps->gensys_ps->eu_dv->v[0] !=1 || govprob_ps->gensys_ps->eu_dv->v[1] !=1) {
                  //Tossing away non-existent or non-unique gensys solution.
                  //printf("\ngovprob_ps->gensys_ps->eu_dv: [%g  %g]\n", govprob_ps->gensys_ps->eu_dv->v[0], govprob_ps->gensys_ps->eu_dv->v[1]);
                  //printf("\nCONTINUE\n");

                  //               //?????????? Debug.
                  //               fprintf(fptr_debug, "*********** Debug ***************:\n");
                  //               fprintf(fptr_debug, "Xrhtran_dm:\n");
                  //               WriteMatrix(fptr_debug, kalcvfurw1datapoint_ps->Xrhtran_dm, " %10.4f ");
                  //               fprintf(fptr_debug, "kalcvfurw1datapoint_ps->ylhtran_dv:\n");
                  //               WriteVector(fptr_debug, kalcvfurw1datapoint_ps->ylhtran_dv, " %10.4f ");
                  //               fprintf(fptr_debug, "z10_dv:\n");
                  //               WriteVector(fptr_debug, kalcvfurw1datapoint_ps->z10_dv, " %10.4f ");
                  //               fprintf(fptr_debug, "zupdate_dv:\n");
                  //               WriteVector(fptr_debug, kalcvfurw1datapoint_ps->zupdate_dv, " %10.4f ");
                  //               fprintf(fptr_debug, "P10_dm:\n");
                  //               WriteMatrix(fptr_debug, kalcvfurw1datapoint_ps->P10_dm, " %10.4f ");
                  //               fprintf(fptr_debug, "V_dm:\n");
                  //               WriteMatrix(fptr_debug, constmodpars_ps->V_dm, " %10.4f ");
                  //               fflush(fptr_debug);

                  //
                  kalcvfurw1datapoint_ps->ylhtran_dv->v[0] = kalcvfurw1datapoint_ps->Xrhtran_dm->M[2];  //[0] is now u_{ki-1}.
                  CopySubvector2matrix(kalcvfurw1datapoint_ps->Xrhtran_dm, 0, 0, XrhtranOld_dv, 0, nrhvars);
                  ki--;
                  continue;
               }
               else  WriteVector(fptr_fcsts, forecasts_ps->foreuinf_k_dv, " %.16e ");     //Prints the results to an output file.
               /**/
               if (!store_sims)  WriteVector(fptr_fcsts, forecasts_ps->foreuinf_k_dv, " %.16e ");     //Prints the results to an output file.
               else onepath_dv->v[ki] = forecasts_ps->foreuinf_k_dv->v[chosenvar];




               //--- Save thid debug for ever!  DDDDDDDDDDDebug.
               // WriteVector(fptr_debug, gov_expinf_t_dv, " %.16e ");
               // if ((gov_expinf_t_dv->v[1]<-15.0) || (gov_expinf_t_dv->v[1]>15.0)) {
               //    //printf("\nHits the boundary\n");
               //    printf("\ngov_expinf_t_dv->v[1] is %g\n", gov_expinf_t_dv->v[1]);
               //    kalcvfurw1datapoint_ps->ylhtran_dv->v[0] = kalcvfurw1datapoint_ps->Xrhtran_dm->M[2];  //[0] is now u_{ki-1}.
               //    CopySubvector2matrix(kalcvfurw1datapoint_ps->Xrhtran_dm, 0, 0, XrhtranOld_dv, 0, nrhvars);
               //    ki--;
               //    continue;
               // }
               // if (gov_expinf_t_dv->v[0] < -3.0)  gov_expinf_t_dv->v[0] = gov_expinf_t_dv->v[1] = -0.0;
               // if (gov_expinf_t_dv->v[0] > 10.0)  gov_expinf_t_dv->v[0] = gov_expinf_t_dv->v[1] =10.0;
               // if ( (ki>370) && (ki<400) ) {
               //    //WriteVector(fptr_debug, gov_expinf_t_dv, " %.16e ");
               //    WriteVector(fptr_debug, kalcvfurw1datapoint_ps->z10_dv, " %.16e ");
               // }

               u_2 = u_1;
               pi_2 = pi_1;
               u_1 = u_0;
               pi_1 = pi_0;
               x_2 = x_1;
               x_1 = gov_expinf_t_dv->v[1];

               //--- Kalman-updating government's regression coefficients.
               kalcvf_urw(kalcvfurw1datapoint_ps, (void *)NULL);
               //--- In case we force no Kalman updating.
               // CopyVector0(z10old_dv, kalcvfurw1datapoint_ps->z10_dv);
               // CopyMatrix0(P10old_dm, kalcvfurw1datapoint_ps->P10_dm);
               //+
               CopyVector0(kalcvfurw1datapoint_ps->z10_dv, kalcvfurw1datapoint_ps->zupdate_dv);
               CopyMatrix0(kalcvfurw1datapoint_ps->P10_dm, kalcvfurw1datapoint_ps->Ppred_dc->C[0]);
            }
         }
         else {       //Counterfactual exercises.
            sclhistshocks_dv = constmodpars_ps->sclhistshocks_dv;
            StrResidstran_dm = constmodpars_ps->StrResidstran_dm;
            if (0) {  // (nsteps+locti > constmodpars_ps->fss) {
               printf(".../swz_comfuns.c/ftd_dynfores_uinfs(): In the counterfactual situation indxCnterFact==1, nsteps+locti, %d, cannot exceed the effective sample size %d", nsteps+locti, constmodpars_ps->fss);
               exit(EXIT_FAILURE);
            }
            //--- Getting all the historical residuals at x_dv.
            ftd_histresids((void *)constmodpars_ps, POINTEST, govprob_ps);
            // if (!store_sims) {
            //    fprintf(fptr_debug, "\n**** Historical Residuals:\n");
            //    if (!StrResids_dm) StrResids_dm = CreateMatrix_lf(constmodpars_ps->fss, 2);
            //    TransposeRegular(StrResids_dm, StrResidstran_dm);
            //    WriteMatrix(fptr_debug, StrResids_dm, " %.16e ");
            // }

            //======= Forecasts for ki=locti+0:locti+nsteps-1. ========
            InitializeKalman1datapoint(kalcvfurw1datapoint_ps, locti, govprob_ps, (void *)constmodpars_ps, POINTEST);  //At ki=-1.
            //--- z_{t|t-1} at 73:12 (163).
            kalcvfurw1datapoint_ps->z10_dv->v[0] = -8.2576116831132182e+000;
            kalcvfurw1datapoint_ps->z10_dv->v[1] = -5.9170702577384295e+000;
            kalcvfurw1datapoint_ps->z10_dv->v[2] = -4.0719350360103029e-001;
            kalcvfurw1datapoint_ps->z10_dv->v[3] =  1.3267510383654432e+001;
            kalcvfurw1datapoint_ps->z10_dv->v[4] =  5.5790782283342677e-001;
            kalcvfurw1datapoint_ps->z10_dv->v[5] =  2.1165972013644286e+001;
            //
            fprintf(fptr_debug, "\n**** Updated government coefficients Ztran:\n");
            for (ki=0; ki<nsteps; ki++) {
               WriteVector(fptr_debug, kalcvfurw1datapoint_ps->z10_dv, " %.16e ");

               //--- Government's inflation and expected inflation.
               forecasts_ps->foreuinf_k_dv->v[0] = forecasts_ps->foreuinf_k_dv->v[1] = x_1;
               if (locti+ki>=235 && locti+ki<=280) {    //118 (70:1); 202 (77:01), 249 (80:12), 273 (82:12).
                  //--- Inflation forecasts at ki.
                  forecasts_ps->foreuinf_k_dv->v[2] = pi_0 = x_1 + 0.0*StrResidstran_dm->M[mos(1, locti+ki, 2)];
                  //--- Unemployment forecasts at ki.
                  forecasts_ps->foreuinf_k_dv->v[3] = u_0 = x_dv->v[0] + constmodpars_ps->gtau1*u_1
                     + constmodpars_ps->gtheta0*(pi_0 - x_1) + constmodpars_ps->gtheta1*(pi_1 - x_2) + 0.0*StrResidstran_dm->M[mos(0, locti+ki, 2)];
               }
               else {
                  forecasts_ps->foreuinf_k_dv->v[2] = pi_0 = x_1 + sclhistshocks_dv->v[1]*StrResidstran_dm->M[mos(1, locti+ki, 2)];
                  //--- Unemployment forecasts at ki.
                  forecasts_ps->foreuinf_k_dv->v[3] = u_0 = x_dv->v[0] + constmodpars_ps->gtau1*u_1
                     + constmodpars_ps->gtheta0*(pi_0 - x_1) + constmodpars_ps->gtheta1*(pi_1 - x_2) + sclhistshocks_dv->v[0]*StrResidstran_dm->M[mos(0, locti+ki, 2)];
               }
               //------- For the next period or loop. --------
               //=== Updates x_t and expinf_t for the government's inflation and expected inflation, which are the same if govprob_ps->grho1=govprob_ps->grho2=0.0.
               //=== The following order matters!
               //--- Classical regression: ylhtran_dv is u_t; 6 variables for Xrhtran_dm in the order of pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and Constant.
               //--- Keynesian regression: ylhtran_dv is pi_t; 6 variables for Xrhtran_dm in the order of u_t, u_{t-1}, pi_{t-1}, u_{t-2}, pi_{t-2}, and Constant.
               kalcvfurw1datapoint_ps->Xrhtran_dm->M[4] = kalcvfurw1datapoint_ps->Xrhtran_dm->M[2];  //[4] is now u_{ki-2} for classical.   [4] is now pi_{ki-2} for Keynesian.
               kalcvfurw1datapoint_ps->Xrhtran_dm->M[3] = kalcvfurw1datapoint_ps->Xrhtran_dm->M[1];  //[3] is now pi_{ki-2} for classical.  [3] is now u_{ki-2} for Keynesian.
               kalcvfurw1datapoint_ps->Xrhtran_dm->M[2] = kalcvfurw1datapoint_ps->ylhtran_dv->v[0];  //[2] is now u_{ki-1} for classical.   [2] is now pi_{ki-1} for Keynesian.
               kalcvfurw1datapoint_ps->Xrhtran_dm->M[1] = kalcvfurw1datapoint_ps->Xrhtran_dm->M[0];  //[1] is now pi_{ki-1} for classical.  [1] is now u_{ki-1} for Keynesian.
               if (indxWhichModel & CGP2TP1TI0) {
                  kalcvfurw1datapoint_ps->Xrhtran_dm->M[0] = pi_0;                                      //[0] is now pi_{ki}.
                  //+ Must be after all the above.
                  kalcvfurw1datapoint_ps->ylhtran_dv->v[0] = u_0;
               }
               else if (indxWhichModel & KGP2TP1TI0) {
                  kalcvfurw1datapoint_ps->Xrhtran_dm->M[0] = u_0;                                      //[0] is now u_{ki}.
                  //+ Must be after all the above.
                  kalcvfurw1datapoint_ps->ylhtran_dv->v[0] = pi_0;
               }
               //--- Getting government's inflation and expected inflation for the next ki (or ki++), stored in gov_expinf_t_dv.
               govexpinf_1update(gov_expinf_t_dv, kalcvfurw1datapoint_ps, govprob_ps, indxWhichModel);

               //=== A priori constraints on the bounds of inflation.
               if (govprob_ps->gensys_ps->eu_dv->v[0] !=1 || govprob_ps->gensys_ps->eu_dv->v[1] !=1) {
                  if (gov_expinf_t_dv->v[0] < -10.0)  gov_expinf_t_dv->v[0] = gov_expinf_t_dv->v[1] = -10.0;
                  if (gov_expinf_t_dv->v[0] > 10.0)  gov_expinf_t_dv->v[0] = gov_expinf_t_dv->v[1] = 10.0;
               }
               if (!store_sims)  WriteVector(fptr_fcsts, forecasts_ps->foreuinf_k_dv, " %.16e ");     //Prints the results to an output file.
               else onepath_dv->v[ki] = forecasts_ps->foreuinf_k_dv->v[chosenvar];

               //--- Updating the data.
               u_2 = u_1;
               pi_2 = pi_1;
               u_1 = u_0;
               pi_1 = pi_0;
               x_2 = x_1;
               x_1 = gov_expinf_t_dv->v[1];
               //fprintf(fptr_debug, "%.16e  %.16e\n", expinf_dv->v[locti+ki+1], gov_expinf_t_dv->v[1]);



               if (!constmodpars_ps->indxNoUpdating) {
                  //--- Kalman-updating government's regression coefficients.
                  kalcvf_urw(kalcvfurw1datapoint_ps, (void *)NULL);
                  CopyVector0(kalcvfurw1datapoint_ps->z10_dv, kalcvfurw1datapoint_ps->zupdate_dv);
                  CopyMatrix0(kalcvfurw1datapoint_ps->P10_dm, kalcvfurw1datapoint_ps->Ppred_dc->C[0]);
               }
               //Else, we force no Kalman updating.
            }
         }
         //--- Writes the simulated path of inflation (or any chosen variable) to an output file.
         if (store_sims)  WriteVector(fptr_fcsts, onepath_dv, " %.4e ");  //Use %.4e to save space in the output file.
         else {
            //=== Writing the limiting result for understanding.
            fprintf(fptr_debug, "At forecast date ki = %d.\n", ki);
            fprintf(fptr_debug, "Xrhtran_dm:\n");
            WriteMatrix(fptr_debug, kalcvfurw1datapoint_ps->Xrhtran_dm, " %10.4f ");
            fprintf(fptr_debug, "z10_dv:\n");
            WriteVector(fptr_debug, kalcvfurw1datapoint_ps->z10_dv, " %10.4f ");
            fprintf(fptr_debug, "P10_dm:\n");
            WriteMatrix(fptr_debug, kalcvfurw1datapoint_ps->P10_dm, " %10.4f ");
            //=== $$$WARNING$$$: should be always commented out!  Uses it as a new initial condition to see how bad it fits to the data.
            // CopyVector0(constmodpars_ps->z10_dv,  kalcvfurw1datapoint_ps->z10_dv);
            // CopyMatrix0(constmodpars_ps->P10_dm,  kalcvfurw1datapoint_ps->P10_dm);
            // ftd_modpars2freepars((void *)constmodpars_ps, POINTEST);
            // fprintf(fptr_debug, "End-of-simulation x_dv:\n");
            // WriteVector(fptr_debug, constmodpars_ps->x_dv, " %.16e ");
         }
      }
      else {
         //Gibbs draws.
         napartsave = constgibbs_ps->napartsave;

         //Updating government's decision rules through t=1:T.
         //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.
         //6 variables for Xrhtran_dm in the order of pi_t, pi_{t-1}, U_{t-1}, pi_{t-2}, U_{t-2}, and constant terms.
         x_dv = constgibbs_ps->x_dv;   //No updating, but will be used below.

         fn_DisplayError(".../swz_comfuns.c/ftd_dynfores_uinfs():  Have not got time to deal with Gibbs draws yet");
      }
   }
   else {
      printf(".../swz_comfuns.c/ftd_dynfores_uinfs(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d!", CGP2TP1TI0, KGP2TP1TI0);
      exit(EXIT_FAILURE);
   }
//   else fn_DisplayError(".../ftd_dynfores_uinfs(): Check CGP2TP1TI0 in swz_comfuns.h. Have not got time to do the case for nlagstrue != 1 or nlagsinflation !=0 or govmodel != 0");


   //===
   DestroyTSforecasts(forecasts_ps);
   DestroyVector_lf(gov_expinf_t_dv);
   DestroyVector_lf(onepath_dv);
   // DestroyVector_lf(XrhtranOld_dv);
   DestroyVector_lf(z10old_dv);
   DestroyMatrix_lf(P10old_dm);

   DestroyMatrix_lf(StrResids_dm);
}


//----------------------------
// Backing out the historical shocks for the true Phillips equation and inflation process.
//----------------------------
static void ftd_histresids(void *constmod_ps, int pointflag, struct TSgovprob_tag *govprob_ps)
{
   //Output: StrResidstran_dm.  1st row: unemployment shocks; 2nd row: inflation shocks.
   int ti;
   TSkalcvfurw *kalcvfurw_ps = govprob_ps->kalcvfurw_ps;
   TSdvector *ylhtran_dv = kalcvfurw_ps->ylhtran_dv;
   TSdmatrix *Xrhtran_dm = kalcvfurw_ps->Xrhtran_dm;  //Order of the variables in Xrhtran_dm is dependent on the value of input_ps->indxWhichModel.
               //For classical regression, the order of 6 variables is pi_t, pi_{t-1}, u_{t-1}, pi_{t-2}, u_{t-2}, and 1.
               //For Keynesian regression, the order of 6 variables is u_t, u_{t-1}, pi_{t-1}, u_{t-2}, pi_{t-2}, and 1.
   int nrhvars = Xrhtran_dm->nrows;
   TSdvector *expinf_dv = govprob_ps->expinf_dv;
   //---
   struct TSconstmodpars_tag *constmodpars_ps;
   struct TSconstgibbs_tag *constgibbs_ps;
   TSdmatrix *StrResidstran_dm;

   if (pointflag & POINTEST) {
      constmodpars_ps = (struct TSconstmodpars_tag *)constmod_ps;
      StrResidstran_dm = constmodpars_ps->StrResidstran_dm;


      for (ti=constmodpars_ps->fss-1; ti>=0; ti--) {
         if (constmodpars_ps->indxWhichModel & CGP2TP1TI0) {        //Classical direction of fit.
            StrResidstran_dm->M[ti*2] = ylhtran_dv->v[ti] - constmodpars_ps->upred_dv->v[ti];  //U shocks -- 1st row of StrResidstran_dm.
            StrResidstran_dm->M[ti*2+1] = Xrhtran_dm->M[mos(0, ti, nrhvars)] - expinf_dv->v[ti];  //Inflation shocks -- 2nd row of StrResidstran_dm.
            StrResidstran_dm->flag = M_GE;
         }
         else if (constmodpars_ps->indxWhichModel & KGP2TP1TI0) {   //Keynesian direction of fit.
            StrResidstran_dm->M[ti*2] = Xrhtran_dm->M[mos(0, ti, nrhvars)] - constmodpars_ps->upred_dv->v[ti]; //U shocks -- 1st row of StrResidstran_dm.
            StrResidstran_dm->M[ti*2+1] = ylhtran_dv->v[ti] - expinf_dv->v[ti];  //Inflation shocks -- 2nd row of StrResidstran_dm.
            StrResidstran_dm->flag = M_GE;
         }
         else {
            printf(".../swz_comfuns.c/ftd_histresids(): Have not got time to do the case for indxWhichModel (defined in swz_comfuns.h and datainp_swz.prn) being other than %d or %d!", CGP2TP1TI0, KGP2TP1TI0);
            exit(EXIT_FAILURE);
         }
      }
   }
   else {
      constgibbs_ps = (struct TSconstgibbs_tag *)constmod_ps;
      StrResidstran_dm = constgibbs_ps->StrResidstran_dm;

   }
}




//      //===????????DDDDDDDDDDebug.
//      fprintf(fptr_debug, "\n********** Some matrices:\n\n");
//      fprintf(fptr_debug, "Total numbers of parameters (must equal): %d (correct) and %d (check)\n", constmodpars_ps->totnpars, cntloc);


//      fprintf(fptr_debug, "x_dv:\n");
//      WriteVector(fptr_debug, x_dv, " %10.4f ");

//      fprintf(fptr_debug, "mean1stset_dv:\n");
//      WriteVector(fptr_debug, mean1stset_dv, " %10.4f ");
//      fprintf(fptr_debug, "Var1stSet_dm:\n");
//      WriteMatrix(fptr_debug, Var1stSet_dm, " %10.4f ");

//      fprintf(fptr_debug, "mean4thset_dv:\n");
//      WriteVector(fptr_debug, mean4thset_dv, " %10.4f ");
//      fprintf(fptr_debug, "Var4thSet_dm:\n");
//      WriteMatrix(fptr_debug, Var4thSet_dm, " %10.4f ");

//      fprintf(fptr_debug, "mean5thset_dv:\n");
//      WriteVector(fptr_debug, mean5thset_dv, " %10.4f ");
//      fprintf(fptr_debug, "Var5thSet_dm:\n");
//      WriteMatrix(fptr_debug, Var5thSet_dm, " %10.4f ");

//      fprintf(fptr_debug, "mean6thset_dv:\n");
//      WriteVector(fptr_debug, mean6thset_dv, " %10.4f ");
//      fprintf(fptr_debug, "Var6thSet_dm:\n");
//      WriteMatrix(fptr_debug, Var6thSet_dm, " %10.4f ");



//----------------------------
// Added 03/17/05. Computing the normalizing constant q for each cut-off p value of Gaussain weight function.
//----------------------------
static void ftd_computeqvals(struct TSconstmodpars_tag *constmodpars_ps, struct TSgovprob_tag *govprob_ps, struct TSconstgibbs_tag *constgibbs_ps)
{
   int ndrawsq = constgibbs_ps->ndrawsq;
   int drawi, _i;
   int passonce_gaussiandraw, passonce_logpostkernel;
   double quadterm, logpostvalue;
   double cutval_logpost = constgibbs_ps->cutval_logpost;
   TSdvector *xmean_dv = constgibbs_ps->xmean_dv;
   int totpars = constgibbs_ps->xmean_dv->n;
   TSdmatrix *cholXcov_dm = constgibbs_ps->cholXcov_dm;
   TSdvector *pcuts_dv = constgibbs_ps->pcuts_dv;
   int npcuts = pcuts_dv->n;
   TSdvector *logqvals_dv = constgibbs_ps->logqvals_dv;
   //---
   TSdvector *teval_dv = NULL;   //n-by-1, the vector of the parameters that are treaed as a join normal.
   TSdvector *txzetadraw_dv = NULL;  //Tempory draw of xzeta.
   TSdvector *txdraw_dv = NULL;  //Tempory draw of x.
   TSivector *qtotn_iv = NULL;  //A vector of total number of draws from the right-truncated normal (pcuts_dv).  Will be <= ndrawsq.


   //=== Memory allocated for this function only.
   teval_dv = CreateVector_lf(totpars);
   txzetadraw_dv = CreateVector_lf(totpars);
   txdraw_dv = CreateVector_lf(totpars);
   qtotn_iv = CreateVector_int(npcuts);

   //=== Computing q values (unlogged at this stage!) for different p values.
   InitializeConstantVector_int(qtotn_iv, 0);   //Cumulative.
   InitializeConstantVector_lf(logqvals_dv, 0.0);   //Cumulative.
   for (drawi=ndrawsq-1; drawi>=0; drawi--)
   {
      passonce_gaussiandraw = 0;  //Reset to "no pass yet."
      passonce_logpostkernel = 0;  //Reset to "no pass yet."
      StandardNormalVector(teval_dv);   //Note: teval_dv = randn(k,1).
      quadterm = VectorDotVector(teval_dv, teval_dv);
      for ( _i=npcuts-1; _i>=0; _i--)
      {
         if ( quadterm < fn_chi2inv(pcuts_dv->v[_i], (double)totpars) )
         {
            qtotn_iv->v[_i]++;

            if (!passonce_gaussiandraw)
            {
               //--- Getting a draw from the trucated Gaussian.
               TrimatrixTimesVector(txzetadraw_dv, cholXcov_dm, teval_dv, 'T', 'N');   //'T' modified 03/17/05.
               if (constgibbs_ps->indxMHMPeak)  VectorPlusMinusVectorUpdate(txzetadraw_dv, constgibbs_ps->xzetapostest_dv, 1.0);  //x0draw = x0peak + L*randn(k,1).
               else  VectorPlusMinusVectorUpdate(txzetadraw_dv, xmean_dv, 1.0);  //x0draw = x0mean + L*randn(k,1).
               passonce_gaussiandraw = 1;
            }
            //--- Check if the restricted region is satisfied.
            if ( (txzetadraw_dv->v[4]>=constgibbs_ps->zeta1cuts_dv->v[0]) && (txzetadraw_dv->v[4]<=constgibbs_ps->zeta1cuts_dv->v[1]) &&
                 (txzetadraw_dv->v[5]>=constgibbs_ps->zeta2cuts_dv->v[0]) && (txzetadraw_dv->v[5]<=constgibbs_ps->zeta2cuts_dv->v[1]) )
            {
               if (!passonce_logpostkernel)
               {
                  //--- Taking sqrt to make it a draw.
                  CopyVector0(txdraw_dv, txzetadraw_dv);
                  txdraw_dv->v[4] = sqrt(txzetadraw_dv->v[4]);
                  txdraw_dv->v[5] = sqrt(txzetadraw_dv->v[5]);
                  logpostvalue = value_logpostkernal_refresh(constmodpars_ps, constgibbs_ps, GIBBSDRAWS, govprob_ps, txdraw_dv);
                  passonce_logpostkernel = 1;
               }
               if (logpostvalue>cutval_logpost)  logqvals_dv->v[_i] += 1.0;  //Recording.
            }

         }
         else  break;  //Early exit.
      }
   }

   for (_i=logqvals_dv->n-1; _i>=0; _i--)  {
      if ( logqvals_dv->v[_i] <= MACHINEZERO )
      {
         printf("\n********WARNING: Arbitraily set logqvals_dv->v[_i] to -infty because we do not get rectrictions satistifed by any draws from Gaussian weight functin with truncted %dth p values!***********\n", _i);
         logqvals_dv->v[_i] = -MACHINEINFINITY;
      }
      else
      {
         //--- Converting it to log values.
         logqvals_dv->v[_i] /= (double)qtotn_iv->v[_i];
         logqvals_dv->v[_i] = log(logqvals_dv->v[_i]);
      }
   }

   printf("\n--------Total numbers of truncated draws for different p's (always < %d, which is the total number of un-truncated draws):-------\n", ndrawsq);
   PrintVector_int(qtotn_iv);
   //
   printf("\n--------log q values for different p (the closer to -infty, the more inaccurate)-------\n");
   PrintVector(logqvals_dv, " %.16e \n");
   //
   fflush(stdout);

   //===
   DestroyVector_lf(teval_dv);
   DestroyVector_lf(txzetadraw_dv);
   DestroyVector_lf(txdraw_dv);
   DestroyVector_int(qtotn_iv);
}