nlopr.cc

一个统计模块的C++源程序！应该能用到的！

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

nlmdl: nlopr.cc

nlmdl is a C++ implementation of the statistical methods in A. Ronald 
Gallant, "Nonlinear Statistical Models," New York: John Wiley and Sons, 
1987, ISBN 0-471-80260-3, using a matrix class realmat that is distributed 
with it.  The header files nlmdl.h and realmat.h describe the use of the 
program and matrix class, respectively.  

Copyright (C) 1990.

A. Ronald Gallant
P.O. Box 5513 
Raleigh NC 27650-5513 
USA   

Permission to use, copy, modify, and distribute this software and its 
documentation for any purpose and without fee is hereby granted, provided 
that the above copyright notice appear in all copies and that both that 
copyright notice and this permission notice appear in supporting 
documentation.  

This software is provided "as is" without any expressed or implied warranty.

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

This is the collection of operators used by nlmdl.  They act on s, of class 
status, using m, of class model; s must have been filled in and the function 
initialize(), of class model, must have been called before using any of these 
operators.  

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

int line_search(char* msg)
{
  REAL    norm_D, norm_theta;
  REAL    step_length;
  REAL    save_obj = s.obj;
  realmat save_theta = s.theta;
  realmat ss(1,1);
  
  ss = T(s.D) * s.D;
  norm_D = sqrt(ss[1]);

  ss = T(s.theta) * s.theta;
  norm_theta = sqrt(ss[1]);

  if ( norm_D < s.tol*(norm_theta + s.eps) ) {
    strcpy(msg,"Tolerence check passed.\n");
    return 1;
  }

  for (step_length=1; step_length >= 0.5; step_length -= 0.1) {

    s.theta = save_theta + step_length * s.D;

    (*opr_obj)();

    if (s.obj < save_obj) {
      sprintf(msg,"Step length = %g \n",step_length);
      return 0 ;
    }
  }

  for (step_length=0.5; step_length > 1.0e-13; step_length *= 0.5) {

    s.theta = save_theta + step_length * s.D;

    (*opr_obj)();

    if (s.obj < save_obj) {
      sprintf(msg,"Step length = %g \n",step_length);
      return 0;
    }
  }

  strcpy(msg,"A line search did not improve this estimate.\n");
  s.obj = save_obj;
  s.theta = save_theta;
  return 1;
}


void SUR_obj()
{
  INTEGER t; 
  realmat sse(1,1,(REAL)0 );
  realmat varinv = invpsd(s.var,s.eps);

  for (t=1; t<=s.n; t++) 
    sse = sse + T(m.e(t)) * varinv * m.e(t);

  s.obj = sse[1];
}


void SUR_mgn()
{
  INTEGER t; 
  realmat sse;
  realmat varinv = invpsd(s.var,s.eps);

  sse.resize(1,1,(REAL)0 );
  s.V.resize(s.p, s.p, (REAL)0 );
  s.D.resize(s.p, 1, (REAL)0 );

  for(t=1; t<=s.n; t++) {
    sse = sse + T(m.e(t))    * varinv * m.e(t);
    s.V = s.V + T(m.dele(t)) * varinv * m.dele(t);
    s.D = s.D + T(m.dele(t)) * varinv * m.e(t);
  }

  s.obj = sse[1];
  s.V = invpsd(s.V,s.eps);
  s.D = - s.V * s.D;
}


void SUR_var(int var_loop)
{
  INTEGER i,j,t;
  REAL fn;
  
  if (rows(s.var)==0 || cols(s.var)==0) {
    s.var.resize(s.M,s.M,(REAL)0 );
    for (i=1; i<=s.M; i++) s.var.elem(i,i)=1.0;
  }

  if (var_loop == 0) return;

  for (i=1; i<=s.M; i++) 
  for (j=1; j<=s.M; j++) 
    s.var.elem(i,j) = (REAL)0 ;

  for (t=1; t<=s.n; t++) 
    s.var = s.var + m.e(t) * T(m.e(t));

  if (s.M == 1 && s.n > s.p) {
    fn = s.n - s.p;
    strcpy(s.df,"corrected");
  }
  else {
    fn = s.n;
    strcpy(s.df,"uncorrected");
  }

  for (i=1; i<=s.M; i++) 
  for (j=1; j<=s.M; j++) 
    s.var.elem(i,j) = s.var.elem(i,j)/fn;
}


void SUR_V() 
{ 
  INTEGER t; 
  INTEGER tau; 
  INTEGER i; 
  REAL    x,weight;
  realmat varinv = invpsd(s.var,s.eps);
  
  if (strcmp(s.vartype,"heteroskedastic") == 0 || s.MA > 0) {

    realmat I(s.p,s.p,(REAL)0 );
    realmat I_tau;

    for (tau=0; tau<=s.MA; tau++) {

      I_tau.resize(s.p,s.p,(REAL)0 );

      for (t=tau+1; t<=s.n; t++)
        I_tau = I_tau 
                + T(m.dele(t)) * varinv * m.e(t)  
                  * T(m.e(t-tau)) * varinv * m.dele(t-tau);

      if (strcmp(s.weights,"Parzen") == 0 && tau > 0) {
        x = tau/(REAL)s.MA;

        if ( x < 0.5 ) 
          weight = 1.0 - 6.0*pow(x,2) + 6.0*pow(x,3);
        else
          weight = 2.0*pow((1.0 - x),3);
      }
      else {
        weight = 1.0;
      }

      I = I + weight*I_tau;

      if (tau > 0)
        I = I + weight*T(I_tau);
    }

    s.V = s.V * I * s.V;
  }

  s.rank = 0;
  for (i=1; i<=s.p; i++)
    if (s.V.elem(i,i) > (REAL)0  ) s.rank++;

  return;
}


realmat qZ;

void opr_qZ(INTEGER t)

{
  INTEGER alpha, i, ii;

  realmat q_tmp(s.M,1);
  realmat Z_tmp(s.K,1);

  qZ.resize(s.M*s.K,1);

  q_tmp = m.e(t);
  Z_tmp = m.Z(t);

  for (alpha=1; alpha<=s.M; alpha++) {
  for (i=1; i<=s.K; i++)             {

    ii = s.K*(alpha-1) + i;

    qZ[ii] = q_tmp[alpha] * Z_tmp[i];
  }
  }
}


realmat QZ;

void opr_QZ(INTEGER t)
{
  INTEGER alpha, i, ii, k;

  realmat Q_tmp(s.M,s.p);
  realmat Z_tmp(s.K,1);

  QZ.resize(s.M*s.K,s.p);

  Q_tmp = m.dele(t);
  Z_tmp = m.Z(t);

  for (k=1; k<=s.p; k++)             {
  for (alpha=1; alpha<=s.M; alpha++) {
  for (i=1; i<=s.K; i++)             {

    ii = s.K*(alpha-1) + i;

    QZ.elem(ii,k) = Q_tmp.elem(alpha,k) * Z_tmp[i];
  }
  }
  }
}


realmat qZ_mts;

void opr_qZ_mts()
{ 

  INTEGER t;
  
  qZ_mts.resize(s.M*s.K, 1, (REAL)0 );

  for (t=1; t<=s.n; t++) {

    opr_qZ(t);

    qZ_mts = qZ_mts + qZ;

  }
}


realmat QZ_mts;

void opr_QZ_mts()
{ 
  INTEGER t;
  
  QZ_mts.resize(s.M*s.K, s.p, (REAL)0 );

  for (t=1; t<=s.n; t++) {

    opr_QZ(t);

    QZ_mts = QZ_mts + QZ;

  }
}


realmat ZZ;

void opr_ZZ()
{
  INTEGER t;

  ZZ.resize(s.K, s.K, (REAL)0 );

  for (t=1; t<=s.n; t++)
    ZZ = ZZ + m.Z(t) * T(m.Z(t));
}


void TSLS_obj() 
{

  INTEGER alpha, i, ii;
  INTEGER beta,  j, jj;

  realmat varinv = invpsd(s.var,s.eps);

  opr_qZ_mts();
  opr_ZZ();

  realmat ZZinv = invpsd(ZZ,s.eps);
  s.obj = (REAL)0 ;

  for (alpha=1; alpha<=s.M; alpha++) {
  for (i=1; i<=s.K; i++)             {
  for (beta=1; beta<=s.M; beta++)    {
  for (j=1; j<=s.K; j++)             {

    ii = s.K*(alpha-1) + i;
    jj = s.K*(beta-1)  + j;

    s.obj += qZ_mts[ii] * qZ_mts[jj] 
             * varinv.elem(alpha,beta) * ZZinv.elem(i,j);
  }
  }
  }
  }
}


void TSLS_mgn() 
{ 
  INTEGER alpha, i, ii;
  INTEGER beta,  j, jj;
  INTEGER k,l;

  realmat varinv = invpsd(s.var,s.eps);

  opr_qZ_mts();
  opr_QZ_mts();
  opr_ZZ();

  realmat ZZinv = invpsd(ZZ,s.eps);

  s.obj = (REAL)0 ;
  s.V.resize(s.p, s.p, (REAL)0 );
  s.D.resize(s.p, 1, (REAL)0 );

  for (alpha=1; alpha<=s.M; alpha++) {
  for (i=1; i<=s.K; i++)             {
  for (beta=1; beta<=s.M; beta++)    {
  for (j=1; j<=s.K; j++)             {

    ii = s.K*(alpha-1) + i;
    jj = s.K*(beta-1)  + j;

    s.obj += qZ_mts[ii] * qZ_mts[jj] 
             * varinv.elem(alpha,beta) * ZZinv.elem(i,j);

  for (k=1; k<=s.p; k++)             {

    s.D[k] += QZ_mts.elem(ii,k) * qZ_mts[jj] 
              * varinv.elem(alpha,beta) * ZZinv.elem(i,j);

  for (l=1; l<=s.p; l++)             {

    s.V.elem(k,l) += QZ_mts.elem(ii,k) * QZ_mts.elem(jj,l) 
                     * varinv.elem(alpha,beta) * ZZinv.elem(i,j);

  }
  }
  }
  }
  }
  }

  s.V = invpsd(s.V,s.eps);

  s.D = - s.V * s.D;

}


void TSLS_V() 
{ 
  INTEGER i;

  s.rank = 0;
  for (i=1; i<=s.p; i++)
    if (s.V.elem(i,i) > (REAL)0  ) s.rank++;

  return;
};


void GMM_obj() 
{ 

  realmat sse(1,1);
  realmat varinv = invpsd(s.var,s.eps);

  opr_qZ_mts();

  sse = T(qZ_mts) * varinv * qZ_mts ;

  s.obj = sse[1];
}


void GMM_mgn() { 

  realmat sse(1,1);
  realmat varinv = invpsd(s.var,s.eps);

  opr_qZ_mts();
  opr_QZ_mts();

  sse = T(qZ_mts) * varinv * qZ_mts ;

  s.obj = sse[1];

  s.D = T(QZ_mts) * varinv * qZ_mts ;

  s.V = T(QZ_mts) * varinv * QZ_mts ;

  s.V = invpsd(s.V,s.eps);

  s.D = - s.V * s.D;

}


void GMM_var(int var_loop) 
{
  INTEGER l = s.M * s.K;
  INTEGER t; 
  INTEGER tau; 
  INTEGER alpha, i, j, ii, jj; 
  REAL    x,weight;

  strcpy(s.df,"uncorrected");

  if (rows(s.var)==0 || cols(s.var)==0) {

    s.var.resize(l,l,(REAL)0 );

    opr_ZZ();

    for (alpha=1; alpha<=s.M; alpha++) {
    for (i=1; i<=s.K; i++)             {
    for (j=1; j<=s.K; j++)             {

      ii = s.K*(alpha-1) + i;
      jj = s.K*(alpha-1) + j;

      s.var.elem(ii,jj) = ZZ.elem(i,j);

    }
    }
    }
  }
  if (var_loop == 0) return;

  realmat I(l,l);
  realmat I_tau(l,l);
  realmat qZ_lag(l,1);

  for (tau=0; tau<=s.MA; tau++) {

    for (i=1; i<=l; i++)
    for (j=1; j<=l; j++)
      I_tau.elem(i,j)=(REAL)0 ;

    for (t=tau+1; t<=s.n; t++) {

      opr_qZ(t-tau);
      qZ_lag = qZ;
      opr_qZ(t);
      
      I_tau = I_tau + qZ * T(qZ_lag);
    }

    if (strcmp(s.weights,"Parzen") == 0 && tau > 0) {
      x = tau/(REAL)s.MA;
      if ( x < 0.5 ) 
        weight = 1.0 - 6.0*pow(x,2) + 6.0*pow(x,3);
      else
        weight = 2.0*pow((1.0 - x),3);
        
    }
    else {
      weight = 1.0;
    }

    I = I + weight*I_tau;
    if (tau > 0) {
      I = I + weight*T(I_tau);
    }
  }

  s.var = I;
 
  return;
}


void GMM_V() 
{ 
  INTEGER i;

  s.rank = 0;
  for (i=1; i<=s.p; i++)
    if (s.V.elem(i,i) > (REAL)0  ) s.rank++;

  return;
}