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

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

Copyright (C) 2020

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.

-------------------------------------------------------------------------------

Function      hac - computes a HAC matrix for data stored columnwise

Syntax        #include "libscl.h"
              INTEGER hac(const realmat& data, realmat& mean, realmat& var)
              INTEGER hac
                (const realmat& data, realmat& mean, realmat& var, INTEGER lag)
              INTEGER hac(const trealmat& data, realmat& mean, realmat& var)
              INTEGER hac
                (const realmat& data, trealmat& mean, realmat& var, INTEGER lag)

Prototype in  libscl.h 

Description   Each row of data is an observation on a multivariat time series. 
              If lag is not specified it is computed as lag=pow(n,1/5).

Remark        The data matrix has n rows and d columns; n is the number of
              observations. 
              The computed matrix var estimates the variance of sqrt(n)*mean. 
              The estimated variance of the mean would therefore be var/n.

Return value  lag

Sample code   realmat data(n,d); // n is number of observations, d is dimension
              //fill data with something
              INTEGER lag = hac(data,mean,var);
              cout << mean << var << corr(var) << lag << '\n';
              realmat X;
              X = T(data); // X = [x_1, ..., x_n] where x_i is d by 1
              lag = hac(T(X),mean,var);
              cout << mean << var << corr(var) << lag << '\n';
              //both usages should produce the same output.

Reference     Gallant, A. Ronald (1987), Nonlinear Statistical Models,
              Wiley, New York, pp. 447, 533

Functions     Library: pow, fabs
called        libscl : realmat

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

namespace {

  REAL parzen(REAL x)   // Note that function is zero for x >= 1
  {
    REAL z = fabs(x);
    if (z <= 0.5) {
      return 1.0 - 6.0*pow(z,2) + 6.0*pow(z,3);
    }
    if (z <= 1.0) {
      return 2.0*pow((1.0 - z),3);
    }
    return 0.0;
  }

}

namespace scl {

  INTEGER hac(const realmat& data, realmat& mean, realmat& var)
  {
    INTEGER n = data.nrow();
    INTEGER lag = pow(n,1.0/5.0);
    return hac(data,mean,var,lag);
  }
  
  INTEGER hac(const realmat& data, realmat& mean, realmat& var, INTEGER lag)
  {
    INTEGER d = data.ncol();
    INTEGER n = data.nrow();
  
    if (lag < 0) error("Error, hac, lag is negative");
    if (lag > n-1) error("Error, hac, lag > n-1");
    
    realmat sum(d,1,0.0);
  
    for (INTEGER j=1; j<=d; ++j) {
      for (INTEGER t=1; t<=n; ++t) {
        sum[j] += data(t,j);
      }
    }
  
    mean = sum/n;
  
    realmat R0(d,d,0.0);
  
    for (INTEGER j=1; j<=d; ++j) {
      for (INTEGER i=1; i<=d; ++i) {
        for (INTEGER t=1; t<=n; t++) {
          R0(i,j) += (data(t,i)-mean[i])*(data(t,j)-mean[j]);
        }
      }
    }
  
    R0 = R0/n;
  
    var = R0;
  
    realmat Rl(d,d);
    
    for (INTEGER l=1; l<=lag; l++) {
     
      fill(Rl);
  
      for (INTEGER j=1; j<=d; ++j) {
        for (INTEGER i=1; i<=d; ++i) {
          for (INTEGER t=1+lag; t<=n; t++) {
            Rl(i,j) += (data(t,i)-mean[i])*(data(t-l,j)-mean[j]);
          }
        }
      }
      
      Rl = Rl/n; 
  
      var += parzen(REAL(l)/REAL(lag+1))*(Rl + T(Rl));
    }
  
    return lag;
  }

  INTEGER hac(const trealmat& data, realmat& mean, realmat& var)
  {
    INTEGER n = data.nrow();
    INTEGER lag = pow(n,1.0/5.0);
    return hac(data,mean,var,lag);
  }
  
  INTEGER hac(const trealmat& data, realmat& mean, realmat& var, INTEGER lag)
  {
    INTEGER d = data.ncol();
    INTEGER n = data.nrow();
  
    if (lag < 0) error("Error, hac, lag is negative");
    if (lag > n-1) error("Error, hac, lag > n-1");
    
    realmat sum(d,1,0.0);
  
    for (INTEGER j=1; j<=d; ++j) {
      for (INTEGER t=1; t<=n; ++t) {
        sum[j] += data(t,j);
      }
    }
  
    mean = sum/n;
  
    realmat R0(d,d,0.0);
  
    for (INTEGER j=1; j<=d; ++j) {
      for (INTEGER i=1; i<=d; ++i) {
        for (INTEGER t=1; t<=n; t++) {
          R0(i,j) += (data(t,i)-mean[i])*(data(t,j)-mean[j]);
        }
      }
    }
  
    R0 = R0/n;
  
    var = R0;
  
    realmat Rl(d,d);
    
    for (INTEGER l=1; l<=lag; l++) {
     
      fill(Rl);
  
      for (INTEGER j=1; j<=d; ++j) {
        for (INTEGER i=1; i<=d; ++i) {
          for (INTEGER t=1+lag; t<=n; t++) {
            Rl(i,j) += (data(t,i)-mean[i])*(data(t-l,j)-mean[j]);
          }
        }
      }
      
      Rl = Rl/n; 
  
      var += parzen(REAL(l)/REAL(lag+1))*(Rl + T(Rl));
    }
  
    return lag;
  }

}
