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

Copyright (C) 2018

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      varcoef - compute VAR estimates 

Syntax        #include "libscl.h"
              #include "kronprd.h"
              INTEGER varcoef(const realmat& data, const intvec& idx, 
                INTEGER lag, realmat& b0, realmat& B, realmat& S,
                realmat& C)

Prototype in  libscl.h

Description   Computes least squares estimates for a VAR defined by the 
              inputs data, idx, and lag. The elements of the intvec idx 
              specify the columns of data to be used and the order in 
              which they are to be used. The dimension M of the left hand 
              side y_t of the VAR, the sample size n, and the definition 
              of the VAR are inferred from the inputs as follows:
              M = idx.size()
              n = data.nrow()
              y_t = (data(t,idx[1], ..., data(t,idx[M])    a column vector
              x_{t-1} = (y_{t-1}, y_{t-2), ..., y_{t-lag}) a column vector
              e_t ~ N(0,S) of dimension M
              y_t = b0 + B*x_{t-1} + e_t,  t = lag+1, ..., n

Return value  Estimated rank of the design matrix; full rank is 1 + M*lag.

Remark        S is unbiased, i.e., S = SSE/(n-lag-rank). The mle estimate 
              is of S is obtined by the usage
              realmat S_mle = REAL(n-lag-rank)/REAL(n-lag)*S. 
              The estimated variance of vec[b0|B] is obtained with the usage 
              realmat V = kronprd(C,S).

Functions     Library: (none)
called        libscl:  realmat

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


#include "libscl.h"
using namespace std;

namespace scl {

  INTEGER varcoef(const realmat& data, const intvec& idx, INTEGER lag,
    realmat& b0, realmat& B, realmat& S, realmat& C)
  {
    INTEGER M = idx.size();
    INTEGER K = 1+M*lag;
    INTEGER n = data.nrow();
    INTEGER df = n - lag - K;
  
    if (lag < 1) error("Error, varcoef, lag < 1");
    if (M > data.ncol()) error("Error, varcoef, size(idx) > data.ncol()");
    if (df <= 0.0) error("Error, varcoef, degress freedom not positive");
    for (INTEGER i=1; i<=M; ++i) {
      if (idx[i] < 1 || data.ncol() < idx[i]) {
        error("Error, varcoef, idx[i] out of bounds for i =" + fmt('d',2,i)());
      }
    }
  
    realmat yX(M,K,0.0);
    realmat XX(K,K,0.0);
  
    realmat y(M,1);
    realmat x(K,1);
    x[1] = 1.0;
    for (INTEGER t=lag+1; t<=n; ++t) {
      for (INTEGER j=1; j<=M; ++j) y[j] = data(t,idx[j]);
      for (INTEGER k=1; k<=lag; ++k) {
        for (INTEGER j=1; j<=M; ++j) x[1+M*(k-1)+j] = data(t-k,idx[j]);
      }
      for (INTEGER j=1; j<=K; ++j) {
        for (INTEGER i=1; i<=M; ++i) yX(i,j) += y[i]*x[j];
        for (INTEGER i=1; i<=K; ++i) XX(i,j) += x[i]*x[j];
      }
    }
  
    INTEGER rank;
    C = invpsd(XX,rank);
  
    realmat A = yX*C;
  
    S.resize(M,M,0.0);
  
    for (INTEGER t=lag+1; t<=n; ++t) {
      for (INTEGER j=1; j<=M; ++j) y[j] = data(t,idx[j]);
      for (INTEGER k=1; k<=lag; ++k) {
        for (INTEGER j=1; j<=M; ++j) x[1+M*(k-1)+j] = data(t-k,idx[j]);
      }
  
      realmat e = y - A*x;
      
      for (INTEGER j=1; j<=M; ++j) {
        for (INTEGER i=1; i<=M; ++i) S(i,j) += e[i]*e[j];
      }
    }
  
    //S = S/df;
  
    S = S/(n-lag-rank);
  
    b0.resize(M,1);
    B.resize(M,K-1);
  
    for (INTEGER i=1; i<=M; ++i) b0[i] = A[i];
    for (INTEGER i=1; i<=M*K-M; ++i) B[i] = A[M+i];
  
    return rank;
  }
}
