#include "libscl.h"
using namespace scl;
using namespace std;

namespace {
  realmat corr(const realmat& var){
    INTEGER d = var.nrow();
    realmat c(d,d);
    for (INTEGER j=1; j<=d; ++j) {
      c(j,j) = sqrt(var(j,j));
      for (INTEGER i=1; i<j; ++i) {
        c(i,j) = var(i,j)/sqrt(var(i,i)*var(j,j));
        c(j,i) = c(i,j);
      }
    } 
    return c;
  }  
}


int main(int argc, char** argp, char** envp)
{
  INTEGER d = 3;
  INTEGER n = 10000;
  REAL rho = 0.9;

  INT_32BIT seed = 8943;
  
  realmat data(d,n);

  for (INTEGER i=1; i<=d; ++i) {
    data(i,1) = sqrt(1.0/(1.0-pow(rho,2)))*unsk(seed);
  }

  for (INTEGER t=2; t<=n; ++t) {
    for  (INTEGER i=1; i<=d; ++i) {
      data(i,t) = unsk(seed) + rho*data(i,t-1);
    }
  }

  for (INTEGER i=1; i<=d; ++i) {
    data(i,1) = data(i,n);
  }

  for (INTEGER t=2; t<=n; ++t) {
    for  (INTEGER i=1; i<=d; ++i) {
      data(i,t) = unsk(seed) + rho*data(i,t-1);
    }
  }


  realmat mean;
  realmat var;
  INTEGER lag=0;
  lag = hac(T(data),mean,var,lag);
  cout << mean << var << corr(var) << lag << '\n';
  lag=1;
  lag = hac(T(data),mean,var,lag);
  cout << mean << var << corr(var) << lag << '\n';
  lag = hac(T(data),mean,var);
  cout << mean << var << corr(var) << lag << '\n';


  realmat tdata;
  tdata = T(data);

  lag = hac(tdata,mean,var);
  cout << mean << var << corr(var) << lag << '\n';

  data = tdata;
  realmat X;
  X = T(data);

  lag = hac(data,mean,var);
  cout << mean << var << corr(var) << lag << '\n';
  lag = hac(T(X),mean,var);
  cout << mean << var << corr(var) << lag << '\n';


  return 0;
}
