/*-----------------------------------------------------------------------------

Copyright (C) 2009.

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

-----------------------------------------------------------------------------*/

#include "libscl.h"

using namespace std;
using namespace scl;

int main(int argc, char** argp, char** envp)
{

  ifstream fin("annual4.dat");
  ofstream fout("bgt.dat");

  realmat raw;
  vecread(fin,raw,11,73);

  INTEGER M = 4;
  INTEGER n = raw.ncol() - 1;

  realmat bgt(M,n);

  REAL p1, d1, c1;
  p1 = d1 = c1 = 0.0;
  for (INTEGER t=1; t<=raw.ncol(); t++) {
    /*
    REAL yyyy   = raw( 1,t);   Not used below.
    */
    REAL pr     = raw( 2,t);
    REAL div    = raw( 3,t);
    REAL ret    = raw( 4,t);
    /*
    REAL qvar   = raw( 5,t);   Not used below.
    REAL dur    = raw( 6,t);   Not used below.
    */
    REAL ndur   = raw( 7,t);
    REAL serv   = raw( 8,t);
    REAL pop    = raw( 9,t);
    /*
    REAL expinf = raw(10,t);   Not used below.
    REAL irate  = raw(11,t);   Not used below.
    */
    REAL p = log(pr/pop);
    REAL d = log(div/pop);
    REAL c = log(ndur+serv);
    if (t > 1) {
      REAL dc = d-c;
      REAL delc = c-c1;
      REAL pd = p-d;
      bgt(1,t-1) = dc;
      bgt(2,t-1) = delc;
      bgt(3,t-1) = pd;
      bgt(4,t-1) = ret;
      fout << fmt('e',27,16,dc) << fmt('e',27,16,delc) 
           << fmt('e',27,16,pd) << fmt('e',27,16,ret) << '\n';
    }
    p1 = p;
    d1 = d;
    c1 = c;
  }

  realmat mean = bgt*realmat(n,1,1.0)/n;

  realmat sdev(M,1,0.0);

  for (INTEGER t=1; t<=n; ++t) {
    for (INTEGER i=1; i<=M; ++i) {
      sdev[i] += pow(bgt(i,t)-mean[i],2);
    }
  }

  for (INTEGER i=1; i<=M; ++i) {
    sdev[i] = sqrt(sdev[i]/(n-1));
  }

  cout << "\nCompare to Table 3 of Bansal, Gallant, and Tauchen (2007)\n";

  cout << mean << sdev << '\n';

  return 0;
}