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

Copyright (C) 2004, 2006.

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.

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

Class         rfunc - Implements the BEKK mulivariate GARCH for a mean 
                      equipped with a Jacobian.
                      
Syntax        #include "libsnp.h"
              #include "squash.h"

              class rfunc {
                //...
              public:
                rfunc(INTEGER ily, INTEGER ilx,
                      INTEGER squasher, REAL inflection, 
                      INTEGER ip, char iPtype, INTEGER iq, char iQtype,
                      INTEGER iv, char iVtype, INTEGER iw, char iWtype);
                rfunc(const rfunc& irf, INTEGER ily, INTEGER ilx,
                      INTEGER squasher, REAL inflection, 
                      INTEGER ip, char iPtype, INTEGER iq, char iQtype,
                      INTEGER iv, char iVtype, INTEGER iw, char iWtype);
                rfunc();
                rfunc(const rfunc& irf);
                ~rfunc();
                rfunc&  operator=(const rfunc& irf);
                INTEGER get_ly() const;
                INTEGER get_lx() const;
                INTEGER get_lR() const;
                const realmat& get_R0() const;
                INTEGER get_p() const;
                char get_Ptype() const;
                INTEGER get_rowsP() const;
                INTEGER get_colsP() const;
                const realmat& get_P() const;
                INTEGER get_q() const;
                char get_Qtype() const;
                INTEGER get_rowsQ() const;
                INTEGER get_colsQ() const;
                const realmat& get_Q() const;
                INTEGER get_v() const;
                char get_Vtype() const;
                INTEGER get_rowsV() const;
                INTEGER get_colsV() const;
                const realmat& get_V() const;
                INTEGER get_w() const;
                char get_Wtype() const;
                INTEGER get_rowsW() const;
                INTEGER get_colsW() const;
                INTEGER get_lRparms() const;
                const realmat& get_Rparms() const;
                const realmat& get_W() const:
                void initialize_state();
                void set_R0(const realmat& iR0);
                void set_P(const realmat& iP);
                void set_Q(const realmat& iQ);
                void set_V(const realmat& iV);
                void set_W(const realmat& iW);
                void set_Rparms(const realmat& iRparms);
                const realmat& operator() (
                  const realmat& iy, const realmat& iu, const realmat& ix,
                  const realmat& iduwb0, const kronprd& iduwB, 
                  realmat& dRwb0, realmat& dRwB, 
                  realmat& dRwR0, realmat& dRwP, realmat& dRwQ, 
                  realmat& dRwV,  realmat& dRwW);
                const realmat& operator() (
                  const realmat& iy, const realmat& iu, const realmat& ix,
                  const realmat& iduwb0, const kronprd& iduwB, 
                  realmat& dRwb0, realmat& dRwB, 
                  realmat& dRwRparms);
                const realmat& operator() (
                  const realmat& iy, const realmat& iu, const realmat& ix);
              };

Declared in   libsnp.h

Description   The BEKK parameterization here is, for p = q = v = w = 1, 
                
                H_t = (R0)(R0)' + P(e_t-1)(e_t-1)'P' + Q(H_t-l)Q'
                      + maxV(e_t-1)max'V(e_t-1) + W((x_t-1)(x_t-1)'W
     
              For more lags, the last four terms are sums.  P, Q, V, and 
              W can be scalar, diagonal, or full matrices; put Ptype, Qtype, 
              Vtype, Wtype to 's', 'd', or 'f', respectively.  R0 is upper 
              triangular, P is stored as [P_lag_1 | ... | P_lag_p] where
              each P_k is scalar for 's', ly by 1 for 'd' and ly by ly for
              'f'; similarly Q, V, W.  V is a leverage effect.  W is a
              level effect that depends on input x of length lx.  Most 
              methods are parameter input and output at a fine level of 
              detail.  The combined parameter vector Rparms is the more 
              convenient interface.  The main method is the application
              operator which computes rfunc and the derivatives of R with 
              respect to b0, B, Rparms.  These are anlytic derivatives.
              Either or both b0 and B can be null.  If squasher is 1 or 2
              then e_t_1 is replased by sqsh(e_t_1) in formulas above.

Remarks       The input parameter ily sets the dimension of y and ilx is
              the same for x.  The second contructor will impute coefficient 
              values from another instance of rfunc.  The only restriction 
              is that ly and lx be the same in both old and new.  This BEKK 
              is meant as an adjuct to SNP estimation.

Warning       The lags are tricky.  The variance for y_t is computed from 
              residual e_{t-1}.  This means that the input y is y_{t-1} and 
              the input u is u_{t-1} from an evaluation of ufunc at ufunc's
              x = y_{t-2}.  Because ufunc handles its own recursion, this 
              must be done by saving the result of a previous call to ufunc. 
              I.e., in a loop indexed by t, ufunc cannot be called with 
              x = y_{t-2} and again with x = y_{t-1}.  The input x is y_{t-1}.  

References    Engle, R. F, and K. F. Kroner (1995), "Multivariate
              Simultaneous Generalized ARCH," Econometric Theory 11, 122-150.
              Gallant, A. Ronald, and George Tauchen (1992), "A
              Nonparametric Approach to Nonlinear Time Series Analysis:
              Estimation and Simulation," in Brillinger, David, Peter
              Caines, John Geweke, Emanuel Parzen, Murray Rosenblatt, and
              Murad S. Taqqu eds. (1992),  New Directions in Time Series
              Analysis, Part II.  Springer--Verlag, New York, 71-92.

Sample        #include "libsnp.h"
program       using namespace scl; 
              using namespace std;
              using namespace libsnp;
              //...
              INTEGER Kz=4; INTEGER Iz=2; INTEGER ly=3;
              INTEGER maxKz=2; INTEGER maxIz=0;
              INTEGER Kx=1; INTEGER Ix=0; INTEGER lx=3;
              INTEGER p=1;  INTEGER q=1; INTEGER v=1; INTEGER w=1;
              INTEGER alag = 1; INTEGER ulag = 2;
              snpden f(Kz, Iz, ly);
              afunc af(f.get_alpha(),maxKz,maxIz,Kx,Ix,lx,alag)
              INTEGER lA = af.get_nrowA()*af.get_ncolA();
              INTEGER la0 = af.get_nrow0();  //Might be zero.
              ufunc uf(f.get_ly(),lx,ulag,true);
              INTEGER lB=uf.get_ly()*uf_get_lx()*uf_get_lag(); //Might be zero
              INTEGER lb0 = 0;
              if (uf.is_intercept()) lb0=uf.get_ly();
              rfunc rf(f.get_ly(),f.get_ly(),sqsh,infl,p,'d',q,'s',v,'s',w,'s');
              INTEGER lRparms = rf.get_lRparms;
              INTEGER lRtheta = lA+la0+lB+lb0+lRparms;
              theta[1]=1.0;  // Optimizer must hold theta[1] fixed.
              //... optimization loop
                INTEGER of=0;
                af.set_A(theta(seq(of+1,of+lA),1));
                of += lA;
                if (la0>0) af.set_a0(theta(seg(of+1,of+la0,1));
                of += la0;
                if (uf.is_intercept()) uf.set_b0(theta(seq(of+1,of+lb0),1));
                of += lb0;
                if (uf.is_regression()) uf.set_B(theta(seq(of+1,of+lB),1));
                of += lB;
                rf.set_Rparms(theta(seq(of+1,of+lR),1));
                of += lR;
                REAL log_likelihood = 0.0;
                realmat dlfwa, dlfwu, dlfwR;
                realmat dawa0,dawA;
                realmat duwb0, duwb0_lag;
                kronprd duwB, duwb_lag;
                realmat dRwb0,dRwB;
                realmat dRwRparms;
                realmat dllwa0;
                if (la0>0) dllwa0.resize(1,la0,0.0);
                realmat dllwA(1,lA,0.0);
                realmat dllwb0;
                if (uf.is_intercept()) dllwb0.resize(1,lb0,0.0);
                realmat dllwB;
                if (uf.is_regression()) dllwB.resize(1,lB,0.0);
                realmat dllwRparms(1,lRparms,0.0);
                realmat x = data("",1);
                realmat u = uf(x,duwb0,duwB)
                af.initialize_state();
                uf.initialize_state();
                rf.initialize_state();
                for (INTEGER t=3; t<=sample_size; ++t) {
                  realmat y = data("",t);
                  realmat y_lag = data("",t-1);
                  realmat x_lag = y_lag;
                  realmat u_lag = u;
                  duwb0_lag = duwb0;
                  duwB_lag = duwB;
                  realmat u = uf(x_lag,duwb0,duwB);
                  f.set_R(rf(y_lag,u_lag,x_lag,duwb0_lag,duwB_lag,
                      dRwb0,dRwB,dRwRparm));
                  f.set_a(af(x_lag,dawa0,dawA));
                  f.set_u(u);
                  log_likelihood += f.log_f(y,dlfwa,dlfwu,dlfwR);
                  dllwA += dlfwa*dawA;
                  if ((la0>0) dllwa0 += dlfwa*dawa0; 
                  if (uf.is_intercept()) {
                    dllwb0 += dlfwu*duwb0;
                    dllwb0 += dlfwR*dRwb0;
                  }
                  if (uf.is_regression()) {
                    dllwB += dlfwu*duwbB;
                    dllwB += dlfwR*dRwbB;
                  }
                  dllwRparms += dlfwR*dRwRparms;
                }
                realmat dllwtheta;
                if (la0>0) dllwtheta = cbind(dllwA,dllwa0);
                else dllwtheta = dllwA;
                if (uf.is_intercept()) cbind(dllwtheta,dllwb0);
                if (uf.is_regression()) cbind(dllwtheta,dllwB);
                dllwtheta = cbind(dllwtheta,dllwRparms);
              //... end optimization loop


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

#include "libsnp.h"
#include "squash.h"

using namespace scl;
using namespace std;

using namespace libsnp;

void libsnp::rfunc::initialize_state()
{
  realmat null;
  kronprd null_null;
  if (p > 0) {
    for (INTEGER k=0; k<p; ++k) {
      Pstate[k] = sE_k(null,null,null,null_null,true);
    }
  }
  if (q > 0) {
    for (INTEGER k=0; k<q; ++k) {
      Qstate[k] = sH_k(null,null,null,null,null,null,null,null,true);
    }
  }
  if (v > 0) {
    for (INTEGER k=0; k<v; ++k) {
      Vstate[k] = sE_k(null,null,null,null_null,true);
    }
  }
  if (w > 0) {
    for (INTEGER k=0; k<w; ++k) {
      Wstate[k] = sX_k(null,true);
    }
  }
}

void libsnp::rfunc::make_private()
{
  if ( q > 0 && p <= 0 ) {
    error("Error, rfunc, rfunc, q cannot be positive when p is zero");
  }
  if ( v > 0 && p <= 0 ) {
    error("Error, rfunc, rfunc, v cannot be positive when p is zero");
  }

  for (INTEGER i=1; i<=ly; ++i) {
    R[(i*i-i)/2+i]=1.0;
    R0[(i*i-i)/2+i]=1.0;
  }

  if (p > 0) {

    Pstate.resize(p);

    switch (Ptype) {
      case 's':
        rowsP = 1;
        colsP = p;
        break;
      case 'd':
        rowsP = ly;
        colsP = p;
        break;
      case 'f':
        rowsP = ly;
        colsP = ly*p;
        break;
      default:
        error("Error, rfunc, rfunc, Ptype must be 's', 'd', or 'f'");
        break;
    }
    P.resize(rowsP,colsP,0.0);
  }

  if (q > 0) {

    Qstate.resize(q);

    switch (Qtype) {
      case 's':
        rowsQ = 1;
        colsQ = q;
        break;
      case 'd':
        rowsQ = ly;
        colsQ = q;
        break;
      case 'f':
        rowsQ = ly;
        colsQ = ly*q;
        break;
      default:
        error("Error, rfunc, rfunc, Ptype must be 's', 'd', or 'f'");
        break;
    }
    Q.resize(rowsQ,colsQ,0.0);
  }

  if (v > 0) {

    Vstate.resize(v);

    switch (Vtype) {
      case 's':
        rowsV = 1;
        colsV = v;
        break;
      case 'd':
        rowsV = ly;
        colsV = v;
        break;
      case 'f':
        rowsV = ly;
        colsV = ly*v;
        break;
      default:
        error("Error, rfunc, rfunc, Vtype must be 's', 'd', or 'f'");
        break;
    }
    V.resize(rowsV,colsV,0.0);
  }

  if (w > 0) {

    Wstate.resize(w);

    switch (Wtype) {
      case 's':
        rowsW = 1;
        colsW = w;
        break;
      case 'd':
        rowsW = lx;
        colsW = w;
        break;
      case 'f':
        rowsW = lx;
        colsW = lx*w;
        break;
      default:
        error("Error, rfunc, rfunc, Wtype must be 's', 'd', or 'f'");
        break;
    }
    W.resize(rowsW,colsW,0.0);
  }

  lRparms = lR;
  if (p > 0) lRparms += rowsP*colsP;
  if (q > 0) lRparms += rowsQ*colsQ;
  if (v > 0) lRparms += rowsV*colsV;
  if (w > 0) lRparms += rowsW*colsW;
  Rparms.resize(lRparms,1);

  initialize_state();

}

libsnp::rfunc::rfunc(INTEGER ily, INTEGER ilx, 
    INTEGER sqsh, REAL infl,
    INTEGER ip, char iPtype, INTEGER iq, char iQtype,
    INTEGER iv, char iVtype, INTEGER iw, char iWtype)
  : ly(ily), lx(ilx), squasher(sqsh), inflection(infl), 
    lR((ly*ly+ly)/2), R(lR,1,0.0), R0(lR,1,0.0), 
    p(ip), Ptype(iPtype), rowsP(0), colsP(0), P(), 
    q(iq), Qtype(iQtype), rowsQ(0), colsQ(0), Q(),
    v(iv), Vtype(iVtype), rowsV(0), colsV(0), V(),
    w(iw), Wtype(iWtype), rowsW(0), colsW(0), W(),
    factor_warn(true)
{
  this->make_private();
}

libsnp::rfunc::rfunc(const rfunc& irf, INTEGER ily, INTEGER ilx,
    INTEGER sqsh, REAL infl,
    INTEGER ip, char iPtype, INTEGER iq, char iQtype,
    INTEGER iv, char iVtype, INTEGER iw, char iWtype)
  : ly(ily), lx(ilx), squasher(sqsh), inflection(infl),
    lR((ly*ly+ly)/2), R(lR,1,0.0), R0(lR,1,0.0), 
    p(ip), Ptype(iPtype), rowsP(0), colsP(0), P(), 
    q(iq), Qtype(iQtype), rowsQ(0), colsQ(0), Q(),
    v(iv), Vtype(iVtype), rowsV(0), colsV(0), V(),
    w(iw), Wtype(iWtype), rowsW(0), colsW(0), W(),
    factor_warn(true)
{

  if ( ly != irf.get_ly() ) {
    error("Error, rfunc, rfunc, y lengths differ");
  }
  
  if ( lx != irf.get_lx() ) {
    error("Error, rfunc, rfunc, x lengths differ");
  }
  
  this->make_private();

  bool nested = true;

  realmat old_R0 = irf.get_R0();
  for (INTEGER i=1; i<=lR; ++i) R0[i] = old_R0[i];
  
  INTEGER old_p = irf.get_p();
  
  if ( (p > 0) && (old_p > 0) ) {
    
    INTEGER min_p;
    if (p >= old_p) {
      min_p = old_p;
    }
    else {
      nested = false;
      min_p = p;
    }
      
    realmat old_P = irf.get_P();

    switch (Ptype) {

      case 's':
        switch (irf.get_Ptype()) {
          case 's':
             for (INTEGER k=1; k<=min_p; ++k) P[k]=old_P[k];
             break;
          case 'd':
            nested = false;
            for (INTEGER k=1; k<=min_p; ++k) P[k]=old_P(1,k);
            break;
          case 'f':
            nested = false;
            for (INTEGER k=1; k<=min_p; ++k) P[k]=old_P(1,ly*(k-1)+1);
            break;
        }
        break;

      case 'd':
        switch (irf.get_Ptype()) {
          case 's':
            for (INTEGER k=1; k<=min_p; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                P(i,k)=old_P[k];
              }
            }
            break;
          case 'd':
            for (INTEGER k=1; k<=min_p; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                P(i,k)=old_P(i,k);
              }
            }
            break;
          case 'f':
            nested = false;
            for (INTEGER k=1; k<=min_p; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                P(i,k)=old_P(i,ly*(k-1)+i);
              }
            }
            break;
        }
        break;

      case 'f':
        switch (irf.get_Ptype()) {
          case 's':
            for (INTEGER k=1; k<=min_p; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                P(i,ly*(k-1)+i)=old_P[k];
              }
            }
            break;
          case 'd':
            for (INTEGER k=1; k<=min_p; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                P(i,ly*(k-1)+i)=old_P(i,k);
              }
            }
            break;
          case 'f':
            for (INTEGER k=1; k<=min_p; ++k) {
              for (INTEGER j=1; j<=ly; ++j) {
                for (INTEGER i=1; i<=ly; ++i) {
                  P(i,ly*(k-1)+j)=old_P(i,ly*(k-1)+j);
                }
              }
            }
            break;
        }
        break;
    }
  }
  else if ( (p == 0) && (old_p != 0) ) {
    nested = false;
  }

  INTEGER old_q = irf.get_q();
  
  if ( (q > 0) && (old_q > 0) ) {
    
    INTEGER min_q;
    if (q >= old_q) {
      min_q = old_q;
    }
    else {
      nested = false;
      min_q = q;
    }
      
    realmat old_Q = irf.get_Q();

    switch (Qtype) {

      case 's':
        switch (irf.get_Qtype()) {
          case 's':
             for (INTEGER k=1; k<=min_q; ++k) Q[k]=old_Q[k];
             break;
          case 'd':
            nested = false;
            for (INTEGER k=1; k<=min_q; ++k) Q[k]=old_Q(1,k);
            break;
          case 'f':
            nested = false;
            for (INTEGER k=1; k<=min_q; ++k) Q[k]=old_Q(1,ly*(k-1)+1);
            break;
        }
        break;

      case 'd':
        switch (irf.get_Qtype()) {
          case 's':
            for (INTEGER k=1; k<=min_q; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                Q(i,k)=old_Q[k];
              }
            }
            break;
          case 'd':
            for (INTEGER k=1; k<=min_q; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                Q(i,k)=old_Q(i,k);
              }
            }
            break;
          case 'f':
            nested = false;
            for (INTEGER k=1; k<=min_q; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                Q(i,k)=old_Q(i,ly*(k-1)+i);
              }
            }
            break;
        }
        break;

      case 'f':
        switch (irf.get_Qtype()) {
          case 's':
            for (INTEGER k=1; k<=min_q; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                Q(i,ly*(k-1)+i)=old_Q[k];
              }
            }
            break;
          case 'd':
            for (INTEGER k=1; k<=min_q; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                Q(i,ly*(k-1)+i)=old_Q(i,k);
              }
            }
            break;
          case 'f':
            for (INTEGER k=1; k<=min_q; ++k) {
              for (INTEGER j=1; j<=ly; ++j) {
                for (INTEGER i=1; i<=ly; ++i) {
                  Q(i,ly*(k-1)+j)=old_Q(i,ly*(k-1)+j);
                }
              }
            }
            break;
        }
        break;
    }
  }
  else if ( (q == 0) && (old_q != 0) ) {
    nested = false;
  }

  INTEGER old_v = irf.get_v();
  
  if ( (v > 0) && (old_v > 0) ) {
    
    INTEGER min_v;
    if (v >= old_v) {
      min_v = old_v;
    }
    else {
      nested = false;
      min_v = v;
    }
      
    realmat old_V = irf.get_V();

    switch (Vtype) {

      case 's':
        switch (irf.get_Vtype()) {
          case 's':
             for (INTEGER k=1; k<=min_v; ++k) V[k]=old_V[k];
             break;
          case 'd':
            nested = false;
            for (INTEGER k=1; k<=min_v; ++k) V[k]=old_V(1,k);
            break;
          case 'f':
            nested = false;
            for (INTEGER k=1; k<=min_v; ++k) V[k]=old_V(1,ly*(k-1)+1);
            break;
        }
        break;

      case 'd':
        switch (irf.get_Vtype()) {
          case 's':
            for (INTEGER k=1; k<=min_v; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                V(i,k)=old_V[k];
              }
            }
            break;
          case 'd':
            for (INTEGER k=1; k<=min_v; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                V(i,k)=old_V(i,k);
              }
            }
            break;
          case 'f':
            nested = false;
            for (INTEGER k=1; k<=min_v; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                V(i,k)=old_V(i,ly*(k-1)+i);
              }
            }
            break;
        }
        break;

      case 'f':
        switch (irf.get_Vtype()) {
          case 's':
            for (INTEGER k=1; k<=min_v; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                V(i,ly*(k-1)+i)=old_V[k];
              }
            }
            break;
          case 'd':
            for (INTEGER k=1; k<=min_v; ++k) {
              for (INTEGER i=1; i<=ly; ++i) {
                V(i,ly*(k-1)+i)=old_V(i,k);
              }
            }
            break;
          case 'f':
            for (INTEGER k=1; k<=min_v; ++k) {
              for (INTEGER j=1; j<=ly; ++j) {
                for (INTEGER i=1; i<=ly; ++i) {
                  V(i,ly*(k-1)+j)=old_V(i,ly*(k-1)+j);
                }
              }
            }
            break;
        }
        break;
    }
  }
  else if ( (v == 0) && (old_v != 0) ) {
    nested = false;
  }


  INTEGER old_w = irf.get_w();
  
  if ( (w > 0) && (old_w > 0) ) {
    
    INTEGER min_w;
    if (w >= old_w) {
      min_w = old_w;
    }
    else {
      nested = false;
      min_w = w;
    }
      
    realmat old_W = irf.get_W();

    switch (Wtype) {

      case 's':
        switch (irf.get_Wtype()) {
          case 's':
             for (INTEGER k=1; k<=min_w; ++k) W[k]=old_W[k];
             break;
          case 'd':
            nested = false;
            for (INTEGER k=1; k<=min_w; ++k) W[k]=old_W(1,k);
            break;
          case 'f':
            nested = false;
            for (INTEGER k=1; k<=min_w; ++k) W[k]=old_W(1,lx*(k-1)+1);
            break;
        }
        break;

      case 'd':
        switch (irf.get_Wtype()) {
          case 's':
            for (INTEGER k=1; k<=min_w; ++k) {
              for (INTEGER i=1; i<=lx; ++i) {
                W(i,k)=old_W[k];
              }
            }
            break;
          case 'd':
            for (INTEGER k=1; k<=min_w; ++k) {
              for (INTEGER i=1; i<=lx; ++i) {
                W(i,k)=old_W(i,k);
              }
            }
            break;
          case 'f':
            nested = false;
            for (INTEGER k=1; k<=min_w; ++k) {
              for (INTEGER i=1; i<=lx; ++i) {
                W(i,k)=old_W(i,lx*(k-1)+i);
              }
            }
            break;
        }
        break;

      case 'f':
        switch (irf.get_Wtype()) {
          case 's':
            for (INTEGER k=1; k<=min_w; ++k) {
              for (INTEGER i=1; i<=lx; ++i) {
                W(i,lx*(k-1)+i)=old_W[k];
              }
            }
            break;
          case 'd':
            for (INTEGER k=1; k<=min_w; ++k) {
              for (INTEGER i=1; i<=lx; ++i) {
                W(i,lx*(k-1)+i)=old_W(i,k);
              }
            }
            break;
          case 'f':
            for (INTEGER k=1; k<=min_w; ++k) {
              for (INTEGER j=1; j<=lx; ++j) {
                for (INTEGER i=1; i<=lx; ++i) {
                  W(i,lx*(k-1)+j)=old_W(i,lx*(k-1)+j);
                }
              }
            }
            break;
        }
        break;
    }
  }
  else if ( (w == 0) && (old_w != 0) ) {
    nested = false;
  }

  if (!nested) warn("Warning, rfunc, rfunc, models not nested");
}

libsnp::rfunc::rfunc()
: ly(0), x(realmat()), lx(0), 
  squasher(0), inflection(2.0), 
  lR(0), R(realmat()), R0(realmat()),
  p(0), Ptype(' '), rowsP(0), colsP(0), P(realmat()),
  q(0), Qtype(' '), rowsQ(0), colsQ(0), Q(realmat()),
  v(0), Vtype(' '), rowsV(0), colsV(0), V(realmat()),
  w(0), Wtype(' '), rowsW(0), colsW(0), W(realmat()),
  lRparms(0), Rparms(realmat()), 
  Pstate(vector<sE_k>()), Qstate(vector<sH_k>()), 
  Vstate(vector<sE_k>()), Wstate(vector<sX_k>())
{ }

libsnp::rfunc::rfunc(const rfunc& irf)
: ly(irf.ly), x(irf.x), lx(irf.lx), 
  squasher(irf.squasher), inflection(irf.inflection), 
  lR(irf.lR), R(irf.R), R0(irf.R0),
  p(irf.p), Ptype(irf.Ptype), rowsP(irf.rowsP), colsP(irf.colsP), P(irf.P),
  q(irf.q), Qtype(irf.Qtype), rowsQ(irf.rowsQ), colsQ(irf.colsQ), Q(irf.Q),
  v(irf.v), Vtype(irf.Vtype), rowsV(irf.rowsV), colsV(irf.colsV), V(irf.V),
  w(irf.w), Wtype(irf.Wtype), rowsW(irf.rowsW), colsW(irf.colsW), W(irf.W),
  lRparms(irf.lRparms), Rparms(irf.Rparms), 
  Pstate(irf.Pstate), Qstate(irf.Qstate), 
  Vstate(irf.Vstate), Wstate(irf.Wstate)
{ }

rfunc& libsnp::rfunc::operator=(const rfunc& irf)
{
  if (this != &irf) {
    ly=irf.ly; x=irf.x; lx=irf.lx; 
    squasher=irf.squasher; inflection=irf.inflection; 
    lR=irf.lR; R=irf.R; R0=irf.R0;
    p=irf.p; Ptype=irf.Ptype; rowsP=irf.rowsP; colsP=irf.colsP; P=irf.P;
    q=irf.q; Qtype=irf.Qtype; rowsQ=irf.rowsQ; colsQ=irf.colsQ; Q=irf.Q;
    v=irf.v; Vtype=irf.Vtype; rowsV=irf.rowsV; colsV=irf.colsV; V=irf.V;
    w=irf.w; Wtype=irf.Wtype; rowsW=irf.rowsW; colsW=irf.colsW; W=irf.W;
    lRparms=irf.lRparms; Rparms=irf.Rparms;
    Pstate=irf.Pstate; Qstate=irf.Qstate; 
    Vstate=irf.Vstate; Wstate=irf.Wstate;
  }
  return *this;
}

const realmat& libsnp::rfunc::operator() 
  (const realmat& iy, const realmat& iu, const realmat& ix,
    const realmat& iduwb0, const kronprd& iduwB,
    realmat& dRwb0, realmat& dRwB,
    realmat& dRwR0, realmat& dRwP, realmat& dRwQ,
    realmat& dRwV,  realmat& dRwW)
{
  /* 
  The notation and ideas for the differentiation come from Gallant,
  A. Ronald, Nonlinear Statistical Models, Wiley, New York, p. 363, 377.
  */

  squash sqsh;
  switch (squasher) {
    case 0: 
      sqsh.set_identity();
      break;
    case 1:
      sqsh.set_spline(inflection);
      break;
    case 2:
      sqsh.set_logistic(inflection);
      break;
    default:
      sqsh.set_identity();
  }

  INTEGER lH = ly*ly;              //Length of vec(H)

  //Accumulators: H accumulates variance terms, the others accumulate
  //derivatives.  Initialization occurs at first use. 

  realmat H;
  realmat dHwR0;
  realmat dHwb0;
  realmat dHwB;
  realmat dHwP;
  realmat dHwQ;
  realmat dHwV;
  realmat dHwW;

  realmat H0(ly,ly);              //H0 = R0*T(R0) (after R0 expanded from vech)
  realmat dH0wR0;
  
  intvec vechH(lR);               //Selects vech(H) from vec(H)
  intvec vecTH(lH);               //Selects vec(T(H)) from vec(H)

  realmat I(ly,ly,0.0);           //Identity of order ly
  for (INTEGER i=1; i<=ly; ++i) { 
    I(i,i) = 1.0;
  }

  {
    realmat eR0(ly,ly,0.0) ;      //Expanded R0
    realmat deR0wR0(lH,lR,0.0);
  
    for (INTEGER j=1; j<=ly; ++j) {
      for (INTEGER i=1; i<=j; ++i) {
        INTEGER R_ij = (j*j-j)/2+i;
        INTEGER H_ij = ly*(j-1)+i;
        INTEGER H_ji = ly*(i-1)+j;
        eR0(i,j) = R0[R_ij];
        eR0(j,i) = R0[R_ij];
        deR0wR0(H_ij,R_ij) = 1.0;
        deR0wR0(H_ji,R_ij) = 1.0;
        vechH[R_ij] = H_ij;
        vecTH[H_ij] = H_ji;
        vecTH[H_ji] = H_ij;
      }
    }
    
    H0 = eR0*T(eR0);
    dH0wR0 = (kronprd(eR0,I) + kronprd(I,eR0))*deR0wR0;
  }                                

  H = H0;
  dHwR0 = dH0wR0;
    
  bool intercept = false;
  bool regression = false;

  if ( p > 0 ) { //There is an MA part and might be intercept or regression.

    if ( iy.get_rows() != ly || iu.get_rows() != ly ) {
      error("Error, rfunc, operator (), y or u has wrong number of rows.");
    }
  
    intercept = true;
    if (iduwb0.size() == 0) {
      intercept = false;
    }
    else if ( iduwb0.get_rows() != ly ) {
      error("Error, rfunc, operator (), iuwb0 has wrong row dimension");
    }
    INTEGER colsb0 = 0;
    if (intercept) colsb0 = iduwb0.get_cols();

    regression = true;
    if (iduwB.size() == 0) {
      regression = false;
    }
    else if (iduwB.get_rows() != ly) {
      error("Error, rfunc, operator (), iduwB has wrong row dimension");
    } 
    INTEGER colsB = 0;
    if (regression) colsB = iduwB.get_cols();

    if (intercept) dHwb0.resize(lH,colsb0,0.0);
    if (regression) dHwB.resize(lH,colsB,0.0);

    INTEGER ldk;
    INTEGER ld;
  
    for (INTEGER k=p-1; k>0; --k) {
      Pstate[k]=Pstate[k-1];
    }
    Pstate[0] = sE_k(iy, iu, iduwb0, iduwB, false);

    switch (Ptype) {

      case 's':
        ldk = 1;                       //Number of cols in dHwP_k
        ld = p*ldk;                    //Number of cols in dHwP

        dHwP.resize(lH,ld,0.0);

        for (INTEGER k=0; k<p; ++k) {
          
          if (Pstate[k].initial) break;

          intvec idxd = seq(ldk*k+1,ldk*k+ldk); //Selects dHwP_k from dHwP

          //Compute P_k*E_k*T(P_k), where E_k = e[k]*T(e[k]).
          //Accumulate in H.

          REAL P_k = P[k+1];            //Pstate indexes from zero, P from one
          REAL P_k_P_k = pow(P_k,2);
          realmat e_k = Pstate[k].y - Pstate[k].u;

          realmat dswe(ly,ly,0.0);            //Jacobian of squasher
          for (INTEGER i=1; i<=ly; ++i) {
            e_k[i] = sqsh(e_k[i],dswe(i,i));  //e_k replaced by squashed
          }                                   //residuals. Better notation
                                              //would have been s_k hereafter,
                                              //but only want to make minimal 
                                              //changes to legacy code.
          realmat E = e_k*T(e_k);

          H += P_k_P_k*E;

          //Compute derivatives of P_k*E_k*T(P_k) wrt b0, B, and P_k.
          //and accumulate in dHwb0, dHwB, and dHwP.

          realmat dPEPwe = P_k_P_k*(kronprd(e_k,I) + kronprd(I,e_k));
          dPEPwe = dPEPwe*dswe;

          realmat dHwP_k(lH,1);
          for (INTEGER i=1; i<=lH; ++i) {
            dHwP_k[i] = 2.0*P_k*E[i];
          }

          if (intercept) dHwb0 -= dPEPwe*Pstate[k].duwb0;
          if (regression) dHwB -= dPEPwe*Pstate[k].duwB;

          for (INTEGER j=1; j<=ldk; ++j) {
            for (INTEGER i=1; i<=lH; ++i) {
              dHwP(i,idxd[j]) += dHwP_k(i,j);
            }
          }
        } 
        break;

      case 'd':
        ldk = ly;                      //Number of cols in dHwP_k
        ld = p*ldk;                    //Number of cols in dHwP

        dHwP.resize(lH,ld,0.0);

        for (INTEGER k=0; k<p; ++k) {
          
          if (Pstate[k].initial) break;

          intvec idxd = seq(ldk*k+1,ldk*k+ldk); //Selects dHwP_k from dHwP

          //Compute P_k*E_k*T(P_k), where E_k = e[k]*T(e[k]).
          //Accumulate in H.

          realmat P_k = P("",k+1);   //Pstate indexes from zero,  P from one
          realmat e_k = Pstate[k].y - Pstate[k].u;

          realmat dswe(ly,ly,0.0);            //Jacobian of squasher
          for (INTEGER i=1; i<=ly; ++i) {
            e_k[i] = sqsh(e_k[i],dswe(i,i));  //e_k replaced by squashed
          }                                   //residuals. Better notation
                                              //would have been s_k hereafter,
                                              //but only want to make minimal 
                                              //changes to legacy code.

          for (INTEGER j=1; j<=ly; ++j) {
            for (INTEGER i=1; i<=ly; ++i) {
              H(i,j) += P_k[i]*e_k[i]*e_k[j]*P_k[j];
            }
          }

          //Compute derivatives of P_k*E_k*T(P_k) wrt b0, B, and P_k.
          //and accumulate in dHwb0, dHwB, and dHwP.

          realmat dPEPwe(lH,ly);
          realmat dHwP_k(lH,ly);
          for (INTEGER c=1; c<=ly; ++c) {
            for (INTEGER j=1; j<=ly; ++j) {
              for (INTEGER i=1; i<=ly; ++i) {
                INTEGER ij = ly*(j-1)+i;
                if ( i != c && j != c ) {
                  dPEPwe(ij,c) = 0.0;
                  dHwP_k(ij,c) = 0.0;
                }
                else if ( i == c && j != c ) {
                  dPEPwe(ij,c) = P_k[i]*e_k[j]*P_k[j];
                  dHwP_k(ij,c) = e_k[i]*e_k[j]*P_k[j];
                }
                else if ( i != c && j == c ) {
                  dPEPwe(ij,c) = P_k[i]*e_k[i]*P_k[j];
                  dHwP_k(ij,c) = P_k[i]*e_k[i]*e_k[j];
                }
                else {
                  dPEPwe(ij,c) = 2.0*P_k[i]*e_k[i]*P_k[j];
                  dHwP_k(ij,c) = 2.0*P_k[i]*e_k[i]*e_k[j];
                }
              }
            }
          }

          dPEPwe = dPEPwe*dswe;

          if (intercept) dHwb0 -= dPEPwe*Pstate[k].duwb0;
          if (regression) dHwB -= dPEPwe*Pstate[k].duwB;

          for (INTEGER j=1; j<=ldk; ++j) {
            for (INTEGER i=1; i<=lH; ++i) {
              dHwP(i,idxd[j]) += dHwP_k(i,j);
            }
          }
        } 
        break;

      case 'f':
        ldk = ly*ly;                   //Number of cols in dHwP_k
        ld = p*ldk;                    //Number of cols in dHwP

        dHwP.resize(lH,ld,0.0);

        for (INTEGER k=0; k<p; ++k) {
  
          if (Pstate[k].initial) break;

          intvec idxd = seq(ldk*k+1,ldk*k+ldk); //Selects dHwP_k from dHwP

          intvec idxPk = seq(ly*k+1,ly*k+ly);   //Selects P_k from P
          
          //Compute P_k*E_k*T(P_k), where E_k = e[k]*T(e[k]).
          //Accumulate in H.

          realmat P_k = P("",idxPk);
          realmat e_k = Pstate[k].y - Pstate[k].u;

          realmat dswe(ly,ly,0.0);            //Jacobian of squasher
          for (INTEGER i=1; i<=ly; ++i) {
            e_k[i] = sqsh(e_k[i],dswe(i,i));  //e_k replaced by squashed
          }                                   //residuals. Better notation
                                              //would have been s_k hereafter,
                                              //but only want to make minimal 
                                              //changes to legacy code.

          realmat Pe_k = P_k*e_k;

          H += Pe_k*T(Pe_k);
  
          //Compute derivatives of P_k*E_k*T(P_k) wrt b0, B, and P_k,
          //and accumulate in dHwb0, dHwB, and dHwP.

          realmat dPEPwe = kronprd(Pe_k,P_k) + kronprd(P_k,Pe_k);
          dPEPwe = dPEPwe*dswe;
          
          if (intercept) dHwb0 -= dPEPwe*Pstate[k].duwb0;
          if (regression) dHwB -= dPEPwe*Pstate[k].duwB;

          realmat PE = Pe_k*T(e_k);
          realmat dHwP_k = kronprd(PE,I) + (kronprd(I,PE))("",vecTH);

          for (INTEGER j=1; j<=ldk; ++j) {
            for (INTEGER i=1; i<=lH; ++i) {
              dHwP(i,idxd[j]) += dHwP_k(i,j);
            }
          }
  
        }
        break;
    }
  }
  
  if ( v > 0 ) { //There is a leverage effect; might be intercept or regression.

    INTEGER colsb0 = 0;
    if (intercept) colsb0 = iduwb0.get_cols();

    INTEGER colsB = 0;
    if (regression) colsB = iduwB.get_cols();

    INTEGER ldk;
    INTEGER ld;
  
    for (INTEGER k=v-1; k>0; --k) {
      Vstate[k]=Vstate[k-1];
    }
    Vstate[0] = sE_k(iy, iu, iduwb0, iduwB, false);

    switch (Vtype) {

      case 's':
        ldk = 1;                       //Number of cols in dHwV_k
        ld = v*ldk;                    //Number of cols in dHwV

        dHwV.resize(lH,ld,0.0);

        for (INTEGER k=0; k<v; ++k) {
          
          if (Vstate[k].initial) break;

          intvec idxd = seq(ldk*k+1,ldk*k+ldk); //Selects dHwV_k from dHwV

          //Compute max(0,V_k*s(e_k))*max'(0,V_k*s(e_k))
          //Accumulate in H.

          //Notation:
          //e_k = residual
          //s_k = s(e_k)         dswe is Jacobian   s_k is squashed residual
          //x_k = V_k*s_k        V_k*I and s_k are Jacobians
          //f_k = max(0,x_k)     dfwx is Jacobian
          //ff  = f_k*T(f_k)     dffwf is Jacobian

          REAL V_k = V[k+1];        //Vstate indexes from zero,  V from one
          realmat e_k = Vstate[k].y - Vstate[k].u;

          realmat s_k(ly,1);
          realmat dswe(ly,ly,0.0);            
          for (INTEGER i=1; i<=ly; ++i) {
            s_k[i] = sqsh(e_k[i],dswe(i,i));  //Squash residuals
          }
           
          realmat x_k = V_k*s_k;

          smomax f(0.1);                    //f.max0(x) returns max(0,x)
                                            //f.max1(x) returns derivative
          realmat f_k(ly,1);
          realmat dfwx(ly,ly,0.0);
          for (INTEGER i=1; i<=ly; ++i) {
            f_k[i] = f.max0(x_k[i]);
            dfwx(i,i) = f.max1(x_k[i]);
          }  

          H += f_k*T(f_k);
  
          //Compute derivatives of f_k*T(f_k) wrt b0, B, and V_k,
          //and accumulate in dHwb0, dHwB, and dHwV.

          //Note:
          //(kronprd(f_k,I) + kronprd(I,f_k))*dfwx 
          //  = kronprd(f_k,dfwx) + kronprd(dfwx,f_k)
          
          realmat dffwx = kronprd(f_k,dfwx) + kronprd(dfwx,f_k);
          realmat dffwe = V_k*dffwx*dswe;
          
          if (intercept) dHwb0 -= dffwe*Vstate[k].duwb0;
          if (regression) dHwB -= dffwe*Vstate[k].duwB;

          realmat dffwV_k = dffwx*s_k;
          
          for (INTEGER j=1; j<=ldk; ++j) {
            for (INTEGER i=1; i<=lH; ++i) {
              dHwV(i,idxd[j]) += dffwV_k(i,j);
            }
          }
        } 
        break;

      case 'd':
        ldk = ly;                      //Number of cols in dHwV_k
        ld = v*ldk;                    //Number of cols in dHwV

        dHwV.resize(lH,ld,0.0);

        for (INTEGER k=0; k<v; ++k) {
          
          if (Vstate[k].initial) break;

          intvec idxd = seq(ldk*k+1,ldk*k+ldk); //Selects dHwV_k from dHwV

          //Compute max(0,V_k*s(e_k))*max'(0,V_k*s(e_k))
          //Accumulate in H.

          //Notation:
          //e_k = residual
          //s_k = s(e_k)         dswe is Jacobian   s_k is squashed residual
          //x_k = diag(V_k)*s_k  diag(V_k) and diag(s_k) are Jacobians
          //f_k = max(0,x_k)     dfwx is Jacobian
          //ff  = f_k*T(f_k)     dffwf is Jacobian

          realmat V_k = V("",k+1);   //Vstate indexes from zero,  V from one
          realmat e_k = Vstate[k].y - Vstate[k].u;

          realmat s_k(ly,1);
          realmat dswe(ly,ly,0.0);            
          for (INTEGER i=1; i<=ly; ++i) {
            s_k[i] = sqsh(e_k[i],dswe(i,i));  //Squash residuals
          }
           
          realmat x_k = diag(V_k)*s_k;

          smomax f(0.1);                    //f.max0(x) returns max(0,x)
                                            //f.max1(x) returns derivative
          realmat f_k(ly,1);
          realmat dfwx(ly,ly,0.0);
          for (INTEGER i=1; i<=ly; ++i) {
            f_k[i] = f.max0(x_k[i]);
            dfwx(i,i) = f.max1(x_k[i]);
          }  

          H += f_k*T(f_k);
  
          //Compute derivatives of f_k*T(f_k) wrt b0, B, and V_k,
          //and accumulate in dHwb0, dHwB, and dHwV.

          //Note:
          //(kronprd(f_k,I) + kronprd(I,f_k))*dfwx 
          //  = kronprd(f_k,dfwx) + kronprd(dfwx,f_k)
          
          realmat dffwx = kronprd(f_k,dfwx) + kronprd(dfwx,f_k);
          realmat dffwe = dffwx*diag(V_k)*dswe;
          
          if (intercept) dHwb0 -= dffwe*Vstate[k].duwb0;
          if (regression) dHwB -= dffwe*Vstate[k].duwB;

          realmat dffwV_k = dffwx*diag(s_k);
          
          for (INTEGER j=1; j<=ldk; ++j) {
            for (INTEGER i=1; i<=lH; ++i) {
              dHwV(i,idxd[j]) += dffwV_k(i,j);
            }
          }
        } 
        break;

      case 'f':
        ldk = ly*ly;                   //Number of cols in dHwV_k
        ld = v*ldk;                    //Number of cols in dHwV

        dHwV.resize(lH,ld,0.0);

        for (INTEGER k=0; k<v; ++k) {
  
          if (Vstate[k].initial) break;

          intvec idxd = seq(ldk*k+1,ldk*k+ldk); //Selects dHwV_k from dHwV

          intvec idxVk = seq(ly*k+1,ly*k+ly);   //Selects V_k from V

          //Compute max(0,V_k*s(e_k))*max'(0,V_k*s(e_k))
          //Accumulate in H.

          //Notation:
          //e_k = residual
          //s_k = s(e_k)       dswe is Jacobian   s_k is squashed residual 
          //x_k = V_k*s_k      V_k and krnprd(T(s),I) are Jacobians
          //f_k = max(0,x_k)   dfwx is Jacobian
          //ff  = f_k*T(f_k)   dffwf is Jacobian

          realmat V_k = V("",idxVk);
          realmat e_k = Vstate[k].y - Vstate[k].u;

          realmat s_k(ly,1);
          realmat dswe(ly,ly,0.0);            
          for (INTEGER i=1; i<=ly; ++i) {
            s_k[i] = sqsh(e_k[i],dswe(i,i));  //Squash residuals
          }

          realmat x_k = V_k*s_k;

          smomax f(0.1);                    //f.max0(x) returns max(0,x)
                                            //f.max1(x) returns derivative
          realmat f_k(ly,1);
          realmat dfwx(ly,ly,0.0);
          for (INTEGER i=1; i<=ly; ++i) {
            f_k[i] = f.max0(x_k[i]);
            dfwx(i,i) = f.max1(x_k[i]);
          }  

          H += f_k*T(f_k);
  
          //Compute derivatives of f_k*T(f_k) wrt b0, B, and V_k,
          //and accumulate in dHwb0, dHwB, and dHwV.

          //Note:
          //(kronprd(f_k,I) + kronprd(I,f_k))*dfwx 
          //  = kronprd(f_k,dfwx) + kronprd(dfwx,f_k)
          
          realmat dffwx = kronprd(f_k,dfwx) + kronprd(dfwx,f_k);
          realmat dffwe = dffwx*V_k*dswe;

          if (intercept) dHwb0 -= dffwe*Vstate[k].duwb0;
          if (regression) dHwB -= dffwe*Vstate[k].duwB;

          //Note:
          //dffwx*kronprd(T(s_k),I) = kronprd(T(s_k),dffwx)

          realmat dffwV_k = kronprd(T(s_k),dffwx);

          for (INTEGER j=1; j<=ldk; ++j) {
            for (INTEGER i=1; i<=lH; ++i) {
              dHwV(i,idxd[j]) += dffwV_k(i,j);
            }
          }
  
        }
        break;
    }
  }
  
  if ( q > 0 ) { //There is an AR part. 
    
    INTEGER ldk;
    INTEGER ld;
        
    switch (Qtype) {

      case 's':
        ldk = 1;                       //Number of cols in dHwQ_k
        ld = q*ldk;                    //Number of cols in dHwQ
        
        dHwQ.resize(lH,ld,0.0);

        for (INTEGER k=0; k<q; ++k) {

          intvec idxd = seq(ldk*k+1,ldk*k+ldk); //Selects dHwQ_k from dHwQ
        
          REAL Q_k = Q[k+1];                    //Select Q_k from Q 
          REAL Q_k_Q_k = pow(Q_k,2); 

          if (Qstate[k].initial) {
          
            //Compute Q_k*H_k*T(Q_k), where H_k = H0.
            //Accumulate in H.

            H += Q_k_Q_k*H0;

            //Compute derivatives of Q_k*H_k*T(Q_k) wrt R0 and Q_k 
            //and accumulate in dHwR0 and dHwQ.
            
            dHwR0 += Q_k_Q_k*dH0wR0;

            realmat dHwQ_k(lH,1);
            for (INTEGER i=1; i<=lH; ++i) {
              dHwQ_k[i] = 2.0*Q_k*H0[i];
            }

            for (INTEGER j=1; j<=ldk; ++j) {
              for (INTEGER i=1; i<=lH; ++i) {
                dHwQ(i,idxd[j]) += dHwQ_k(i,j);
              }
            }
          }
          else {
          
            //Compute Q_k*H_k*T(Q_k), where H_k = Qstate[k].H.
            //Accumulate in H.

            H += Q_k_Q_k*Qstate[k].H;
  
            //Compute derivatives of Q_k*H_k*T(Q_k) 
            //wrt R0, b0, B, P, Q_k, Q, V, W
            //and accumulate in dHwb0, dHwB, dHwR0, dHP, dHwQ, dHwV, dHwW.
  
            if (intercept) dHwb0 += Q_k_Q_k*Qstate[k].dHwb0;
            if (regression) dHwB += Q_k_Q_k*Qstate[k].dHwB;
            dHwR0 += Q_k_Q_k*Qstate[k].dHwR0;
            dHwP += Q_k_Q_k*Qstate[k].dHwP;
            dHwQ += Q_k_Q_k*Qstate[k].dHwQ;
            if (v > 0) dHwV += Q_k_Q_k*Qstate[k].dHwV;
            
            realmat dHwQ_k(lH,1);
            for (INTEGER i=1; i<=lH; ++i) {
              dHwQ_k[i] = 2.0*Q_k*Qstate[k].H[i];
            }

            for (INTEGER j=1; j<=ldk; ++j) {
              for (INTEGER i=1; i<=lH; ++i) {
                dHwQ(i,idxd[j]) += dHwQ_k(i,j);
              }
            }
          }
        } // end for loop

        for (INTEGER k=q-1; k>0; --k) {
          Qstate[k]=Qstate[k-1];
        }
        Qstate[0] = sH_k(H, dHwb0, dHwB, dHwR0, dHwP, dHwQ, dHwV, dHwW, false);

        break;

      case 'd':
        ldk = ly;                      //Number of cols in dHwQ_k
        ld = q*ldk;                    //Number of cols in dHwQ
        
        dHwQ.resize(lH,ld,0.0);

        for (INTEGER k=0; k<q; ++k) {

          intvec idxd = seq(ldk*k+1,ldk*k+ldk); //Selects dHwQ_k from dHwQ
        
          realmat Q_k = Q("",k+1);              //Select Q_k from Q 

          if (Qstate[k].initial) {
          
            //Compute Q_k*H_k*T(Q_k), where H_k = H0.
            //Accumulate in H.

            for (INTEGER j=1; j<=ly; ++j) {
              for (INTEGER i=1; i<=ly; ++i) {
                H(i,j) += Q_k[i]*H0(i,j)*Q_k[j];
              }
            }

            //Compute derivatives of Q_k*H_k*T(Q_k) wrt R0 and Q_k 
            //and accumulate in dHwR0 and dHwQ.
            
            realmat dQHQwH(lH,lH,0.0);
            for (INTEGER j=1; j<=ly; ++j) {
              for (INTEGER i=1; i<=ly; ++i) {
                INTEGER ij = ly*(j-1)+i;
                dQHQwH(ij,ij) = Q_k[i]*Q_k[j];
              }
            }
  
            dHwR0 += dQHQwH*dH0wR0;

            realmat dHwQ_k(lH,ly);
            for (INTEGER c=1; c<=ly; ++c) {
              for (INTEGER j=1; j<=ly; ++j) {
                for (INTEGER i=1; i<=ly; ++i) {
                  INTEGER ij = ly*(j-1)+i;
                  if ( i != c && j != c ) {
                    dHwQ_k(ij,c) = 0.0;
                  }
                  else if ( i == c && j != c ) {
                    dHwQ_k(ij,c) = H0(i,j)*Q_k[j];
                  }
                  else if ( i != c && j == c ) {
                    dHwQ_k(ij,c) = Q_k[i]*H0(i,j);
                  }
                  else {
                    dHwQ_k(ij,c) = 2.0*H0(c,c)*Q_k[c];
                  }
                }
              }
            }

            for (INTEGER j=1; j<=ldk; ++j) {
              for (INTEGER i=1; i<=lH; ++i) {
                dHwQ(i,idxd[j]) += dHwQ_k(i,j);
              }
            }
          }
          else {
          
            //Compute Q_k*H_k*T(Q_k), where H_k = Qstate[k].H.
            //Accumulate in H.

            for (INTEGER j=1; j<=ly; ++j) {
              for (INTEGER i=1; i<=ly; ++i) {
                H(i,j) += Q_k[i]*Qstate[k].H(i,j)*Q_k[j];
              }
            }

            //Compute derivatives of Q_k*H_k*T(Q_k) wrt R0, b0, B, P, Q_k, Q
            //and accumulate in dHwb0, dHwB, dHwR0, dHP, dHwQ.
  
            realmat dQHQwH(lH,lH,0.0);
            for (INTEGER j=1; j<=ly; ++j) {
              for (INTEGER i=1; i<=ly; ++i) {
                INTEGER ij = ly*(j-1)+i;
                dQHQwH(ij,ij) = Q_k[i]*Q_k[j];
              }
            }
            
            if (intercept) dHwb0 += dQHQwH*Qstate[k].dHwb0;
            if (regression) dHwB += dQHQwH*Qstate[k].dHwB;
            dHwR0 += dQHQwH*Qstate[k].dHwR0;
            dHwP += dQHQwH*Qstate[k].dHwP;
            dHwQ += dQHQwH*Qstate[k].dHwQ;
            if (v > 0) dHwV += dQHQwH*Qstate[k].dHwV;

            realmat dHwQ_k(lH,ly);
            for (INTEGER c=1; c<=ly; ++c) {
              for (INTEGER j=1; j<=ly; ++j) {
                for (INTEGER i=1; i<=ly; ++i) {
                  INTEGER ij = ly*(j-1)+i;
                  if ( i != c && j != c ) {
                    dHwQ_k(ij,c) = 0.0;
                  }
                  else if ( i == c && j != c ) {
                    dHwQ_k(ij,c) = Qstate[k].H(i,j)*Q_k[j];
                  }
                  else if ( i != c && j == c ) {
                    dHwQ_k(ij,c) = Q_k[i]*Qstate[k].H(i,j);
                  }
                  else {
                    dHwQ_k(ij,c) = 2.0*Qstate[k].H(c,c)*Q_k[c];
                  }
                }
              }
            }

            for (INTEGER j=1; j<=ldk; ++j) {
              for (INTEGER i=1; i<=lH; ++i) {
                dHwQ(i,idxd[j]) += dHwQ_k(i,j);
              }
            }
          }

        }  // end for loop

        for (INTEGER k=q-1; k>0; --k) {
          Qstate[k]=Qstate[k-1];
        }
        Qstate[0] = sH_k(H, dHwb0, dHwB, dHwR0, dHwP, dHwQ, dHwV, dHwW, false);
          
        break;
        
      case 'f':

        ldk = ly*ly;                   //Number of cols in dHwQ_k
        ld = q*ldk;                    //Number of cols in dHwQ
        
        dHwQ.resize(lH,ld,0.0);

        for (INTEGER k=0; k<q; ++k) {

          intvec idxQk = seq(ly*k+1,ly*k+ly);   //Selects Q_k from Q
          intvec idxd = seq(ldk*k+1,ldk*k+ldk); //Selects dHwQ_k from dHwQ
        
          realmat Q_k = Q("",idxQk);            //Select Q_k from Q 

          if (Qstate[k].initial) {
          
            //Compute Q_k*H_k*T(Q_k), where H_k = H0.
            //Accumulate in H.

            H += Q_k*H0*T(Q_k);

            //Compute derivatives of Q_k*H_k*T(Q_k) wrt R0 and Q_k 
            //and accumulate in dHwR0 and dHwQ.
            
            kronprd dQHQwH = kronprd(Q_k,Q_k);

            dHwR0 += dQHQwH*dH0wR0;

            realmat QH = Q_k*H0;
            realmat dHwQ_k = kronprd(QH,I) + (kronprd(I,QH))("",vecTH);

            for (INTEGER j=1; j<=ldk; ++j) {
              for (INTEGER i=1; i<=lH; ++i) {
                dHwQ(i,idxd[j]) += dHwQ_k(i,j);
              }
            }
          }
          else {
          
            //Compute Q_k*H_k*T(Q_k), where H_k = Qstate[k].H.
            //Accumulate in H.

            H += Q_k*Qstate[k].H*T(Q_k);
  
            //Compute derivatives of Q_k*H_k*T(Q_k) wrt R0, b0, B, P, Q_k, Q
            //and accumulate in dHwb0, dHwB, dHwR0, dHP, dHwQ.
  
            kronprd dQHQwH = kronprd(Q_k,Q_k);

            if (intercept) dHwb0 += dQHQwH*Qstate[k].dHwb0;
            if (regression) dHwB += dQHQwH*Qstate[k].dHwB;
            dHwR0 += dQHQwH*Qstate[k].dHwR0;
            dHwP += dQHQwH*Qstate[k].dHwP;
            dHwQ += dQHQwH*Qstate[k].dHwQ;
            if (v > 0) dHwV += dQHQwH*Qstate[k].dHwV;
            
            realmat QH = Q_k*Qstate[k].H;
            realmat dHwQ_k = kronprd(QH,I) + (kronprd(I,QH))("",vecTH);
            
            for (INTEGER j=1; j<=ldk; ++j) {
              for (INTEGER i=1; i<=lH; ++i) {
                dHwQ(i,idxd[j]) += dHwQ_k(i,j);
              }
            }
          }
        } // end for loop

        for (INTEGER k=q-1; k>0; --k) {
          Qstate[k]=Qstate[k-1];
        }
        Qstate[0] = sH_k(H, dHwb0, dHwB, dHwR0, dHwP, dHwQ, dHwV, dHwW, false);
          
        break;
    }
  }
 
  if ( w > 0 ) {      //There is a level effect; the level effect is
                      //not run through the GARCH recursion although
                      //dHwW is in Qstate in case that changes in which
                      //case this code must be moved after code for V 
                      //and before code for Q
    INTEGER ldk;
    INTEGER ld;
  
    for (INTEGER k=w-1; k>0; --k) {
      Wstate[k]=Wstate[k-1];
    }
    Wstate[0] = sX_k(ix, false);

    switch (Wtype) {

      case 's':
        ldk = 1;                       //Number of cols in dHwW_k
        ld = w*ldk;                    //Number of cols in dHwW

        dHwW.resize(lH,ld,0.0);

        for (INTEGER k=0; k<w; ++k) {
          
          if (Wstate[k].initial) break;

          intvec idxd = seq(ldk*k+1,ldk*k+ldk); //Selects dHwW_k from dHwW

          //Compute W_k*X_k*T(W_k), where X_k = x[k]*T(x[k]).
          //Accumulate in H.

          REAL W_k = W[k+1];            //Wstate indexes from zero,  W from one
          REAL W_k_W_k = pow(W_k,2);
          realmat x_k = Wstate[k].x;

          realmat X = x_k*T(x_k);

          H += W_k_W_k*X;

          //Compute derivatives of W_k*X_k*T(W_k) wrt W_k.
          //and accumulate in dHwW.

          realmat dWXWwx = W_k_W_k*(kronprd(x_k,I) + kronprd(I,x_k));

          realmat dHwW_k(lH,1);
          for (INTEGER i=1; i<=lH; ++i) {
            dHwW_k[i] = 2.0*W_k*X[i];
          }

          for (INTEGER j=1; j<=ldk; ++j) {
            for (INTEGER i=1; i<=lH; ++i) {
              dHwW(i,idxd[j]) += dHwW_k(i,j);
            }
          }
        } 
        break;

      case 'd':
        ldk = lx;                      //Number of cols in dHwW_k
        ld = w*ldk;                    //Number of cols in dHwW

        dHwW.resize(lH,ld,0.0);

        for (INTEGER k=0; k<w; ++k) {
          
          if (Wstate[k].initial) break;

          intvec idxd = seq(ldk*k+1,ldk*k+ldk); //Selects dHwW_k from dHwW

          //Compute W_k*X_k*T(W_k), where X_k = x[k]*T(x[k]).
          //Accumulate in H.

          realmat W_k = W("",k+1);   //Wstate indexes from zero,  W from one
          realmat x_k = Wstate[k].x;

          for (INTEGER j=1; j<=lx; ++j) {
            for (INTEGER i=1; i<=lx; ++i) {
              H(i,j) += W_k[i]*x_k[i]*x_k[j]*W_k[j];
            }
          }

          //Compute derivatives of W_k*X_k*T(W_k) wrt W_k.
          //and accumulate in dHwW.

          realmat dWXWwx(lH,lx);
          realmat dHwW_k(lH,lx);
          for (INTEGER c=1; c<=lx; ++c) {
            for (INTEGER j=1; j<=lx; ++j) {
              for (INTEGER i=1; i<=lx; ++i) {
                INTEGER ij = lx*(j-1)+i;
                if ( i != c && j != c ) {
                  dWXWwx(ij,c) = 0.0;
                  dHwW_k(ij,c) = 0.0;
                }
                else if ( i == c && j != c ) {
                  dWXWwx(ij,c) = W_k[i]*x_k[j]*W_k[j];
                  dHwW_k(ij,c) = x_k[i]*x_k[j]*W_k[j];
                }
                else if ( i != c && j == c ) {
                  dWXWwx(ij,c) = W_k[i]*x_k[i]*W_k[j];
                  dHwW_k(ij,c) = W_k[i]*x_k[i]*x_k[j];
                }
                else {
                  dWXWwx(ij,c) = 2.0*W_k[i]*x_k[i]*W_k[j];
                  dHwW_k(ij,c) = 2.0*W_k[i]*x_k[i]*x_k[j];
                }
              }
            }
          }

          for (INTEGER j=1; j<=ldk; ++j) {
            for (INTEGER i=1; i<=lH; ++i) {
              dHwW(i,idxd[j]) += dHwW_k(i,j);
            }
          }
        } 
        break;

      case 'f':
        ldk = lx*lx;                   //Number of cols in dHwW_k
        ld = w*ldk;                    //Number of cols in dHwW

        dHwW.resize(lH,ld,0.0);

        for (INTEGER k=0; k<w; ++k) {
  
          if (Wstate[k].initial) break;

          intvec idxd = seq(ldk*k+1,ldk*k+ldk); //Selects dHwW_k from dHwW

          intvec idxWk = seq(lx*k+1,lx*k+lx);   //Selects W_k from W
          
          //Compute W_k*X_k*T(W_k), where X_k = x[k]*T(x[k]).
          //Accumulate in H.

          realmat W_k = W("",idxWk);
          realmat x_k = Wstate[k].x;

          realmat Wx_k = W_k*x_k;

          H += Wx_k*T(Wx_k);
  
          //Compute derivatives of W_k*X_k*T(W_k) wrt W_k,
          //and accumulate in dHwW.

          realmat dWXWwx = kronprd(Wx_k,W_k) + kronprd(W_k,Wx_k);

          realmat WX = Wx_k*T(x_k);
          realmat dHwW_k = kronprd(WX,I) + (kronprd(I,WX))("",vecTH);

          for (INTEGER j=1; j<=ldk; ++j) {
            for (INTEGER i=1; i<=lH; ++i) {
              dHwW(i,idxd[j]) += dHwW_k(i,j);
            }
          }
  
        }
        break;
    }
  }
 
  //Factor H
  
  realmat dRwS(lR,lR);               //S is the upper triangle of H

  for (INTEGER j=1; j<=ly; ++j) {
    for (INTEGER i=1; i<=j; ++i) {
      INTEGER R_ij = (j*j-j)/2+i;
      INTEGER H_ij = ly*(j-1)+i;
      R[R_ij] = H[H_ij];             //R now contains the upper triangle of H
    }
  }
 
  INTEGER ier = factor(R.get_x(),dRwS.get_x(),ly);

  if (factor_warn && (ier == -1)) {
    warn("Warning, rfunc, operator (), factorization failure");
    factor_warn = false;
  }
  
  intvec null;

  //dRwR0 = dRwS*dHwR0(vechH,"");          //Note application of selection
  dgmprd(dRwS,dHwR0,vechH,null,dRwR0);     //matrix here and below

  if (p > 0) {
    //dRwP = dRwS*dHwP(vechH,"");
    dgmprd(dRwS,dHwP,vechH,null,dRwP);
  }

  if (q > 0) {
    //dRwQ = dRwS*dHwQ(vechH,"");
    dgmprd(dRwS,dHwQ,vechH,null,dRwQ);
  }
  
  if (intercept) {
    //dRwb0 = dRwS*dHwb0(vechH,"");
    dgmprd(dRwS,dHwb0,vechH,null,dRwb0);
  }

  if (regression) {
    //dRwB = dRwS*dHwB(vechH,"");
    dgmprd(dRwS,dHwB,vechH,null,dRwB);
  }

  if (v > 0) {
    //dRwV = dRwS*dHwV(vechH,"");
    dgmprd(dRwS,dHwV,vechH,null,dRwV);
  }

  if (w > 0) {
    //dRwW = dRwS*dHwW(vechH,"");
    dgmprd(dRwS,dHwW,vechH,null,dRwW);
  }

  return R;
}

libsnp::rfunc::sE_k::sE_k() 
      : y(), u(), duwb0(), duwB(), initial(true) {}
    
libsnp::rfunc::sE_k::sE_k
       (const realmat& iy, const realmat& iu, 
        const realmat& iduwb0, const kronprd& iduwB,
        bool init)
      : y(iy), u(iu), duwb0(iduwb0), duwB(iduwB), 
        initial(init) {}

libsnp::rfunc::sE_k::sE_k(const sE_k& isE_k)
{
  /*
  y = isE_k.y;
  u = isE_k.u;
  duwb0 = isE_k.duwb0;
  duwB = isE_k.duwB;
  */
  dgmcpy(isE_k.y,y);
  dgmcpy(isE_k.u,u);
  dgmcpy(isE_k.duwb0,duwb0);
  duwB = isE_k.duwB;            //Don't have a dgmcpy for kronprd
  initial = isE_k.initial;
}

libsnp::rfunc::sE_k& rfunc::sE_k::operator=(const sE_k& isE_k)
{
  if (this != &isE_k) {
    /*
    y = isE_k.y;
    u = isE_k.u;
    duwb0 = isE_k.duwb0;
    duwB = isE_k.duwB;
    */
    dgmcpy(isE_k.y,y);
    dgmcpy(isE_k.u,u);
    dgmcpy(isE_k.duwb0,duwb0);
    duwB = isE_k.duwB;            //Don't have a dgmcpy for kronprd
    initial = isE_k.initial;
  }
  return *this;
}

libsnp::rfunc::sH_k::sH_k() 
      : H(), dHwb0(), dHwB(), dHwR0(), dHwP(), dHwQ(), dHwV(), dHwW(),
        initial(true) {}
    
libsnp::rfunc::sH_k::sH_k
       (const realmat& iH, const realmat& idHwb0, const realmat& idHwB,
        const realmat& idHwR0, const realmat& idHwP, const realmat& idHwQ, 
        const realmat& idHwV, const realmat& idHwW, bool init)
      : H(iH), dHwb0(idHwb0), dHwB(idHwB), 
        dHwR0(idHwR0), dHwP(idHwP), dHwQ(idHwQ), 
        dHwV(idHwV), dHwW(idHwW), initial(init) {}

libsnp::rfunc::sH_k::sH_k(const sH_k& isH_k)
{
  /*
  H = isH_k.H;
  dHwb0 = isH_k.dHwb0;
  dHwB = isH_k.dHwB;
  dHwR0 = isH_k.dHwR0;
  dHwP = isH_k.dHwP;
  dHwQ = isH_k.dHwQ;
  dHwV = isH_k.dHwV;
  dHwW = isH_k.dHwW;
  */
  dgmcpy(isH_k.H,H);
  dgmcpy(isH_k.dHwb0,dHwb0);
  dgmcpy(isH_k.dHwB,dHwB);
  dgmcpy(isH_k.dHwR0,dHwR0);
  dgmcpy(isH_k.dHwP,dHwP);
  dgmcpy(isH_k.dHwQ,dHwQ);
  dgmcpy(isH_k.dHwV,dHwV);
  dgmcpy(isH_k.dHwW,dHwW);
  initial = isH_k.initial;
}

libsnp::rfunc::sH_k& rfunc::sH_k::operator=(const sH_k& isH_k)
{
  if (this != &isH_k) {
    /*
    H = isH_k.H;
    dHwb0 = isH_k.dHwb0;
    dHwB = isH_k.dHwB;
    dHwR0 = isH_k.dHwR0;;
    dHwP = isH_k.dHwP;
    dHwQ = isH_k.dHwQ;
    dHwV = isH_k.dHwV;
    dHwW = isH_k.dHwW;
    */
    dgmcpy(isH_k.H,H);
    dgmcpy(isH_k.dHwb0,dHwb0);
    dgmcpy(isH_k.dHwB,dHwB);
    dgmcpy(isH_k.dHwR0,dHwR0);
    dgmcpy(isH_k.dHwP,dHwP);
    dgmcpy(isH_k.dHwQ,dHwQ);
    dgmcpy(isH_k.dHwV,dHwV);
    dgmcpy(isH_k.dHwW,dHwW);
    initial = isH_k.initial;
  }
  return *this;
}

libsnp::rfunc::sX_k::sX_k() 
      : x(), initial(true) {}
    
libsnp::rfunc::sX_k::sX_k (const realmat& ix, bool init)
      : x(ix), initial(init) {}

libsnp::rfunc::sX_k::sX_k(const sX_k& isX_k)
{
  //x = isX_k.x;
  dgmcpy(isX_k.x,x);
  initial = isX_k.initial;
}

libsnp::rfunc::sX_k& rfunc::sX_k::operator=(const sX_k& isX_k)
{
  if (this != &isX_k) {
    //x = isX_k.x;
    dgmcpy(isX_k.x,x);
    initial = isX_k.initial;
  }
  return *this;
}

const realmat& libsnp::rfunc::get_Rparms() const
{
  INTEGER k=0;

  for (INTEGER i=1; i<=lR; ++i) {
    ++k;
    Rparms[k] = R0[i];
  }

  if (p > 0 ) {
    for (INTEGER i=1; i<=rowsP*colsP; ++i) {
      ++k;
      Rparms[k] = P[i];
    }
  }

  if (q > 0 ) {
    for (INTEGER i=1; i<=rowsQ*colsQ; ++i) {
      ++k;
      Rparms[k] = Q[i];
    }
  }

  if (v > 0 ) {
    for (INTEGER i=1; i<=rowsV*colsV; ++i) {
      ++k;
      Rparms[k] = V[i];
    }
  }

  if (w > 0 ) {
    for (INTEGER i=1; i<=rowsW*colsW; ++i) {
      ++k;
      Rparms[k] = W[i];
    }
  }

  return Rparms;
}

void libsnp::rfunc::set_Rparms(const realmat& iRparms)
{
  INTEGER k=0;

  for (INTEGER i=1; i<=lR; ++i) {
    ++k;
    R0[i] = iRparms[k];
  }

  if (p > 0 ) {
    for (INTEGER i=1; i<=rowsP*colsP; ++i) {
      ++k;
      P[i] = iRparms[k];
    }
  }

  if (q > 0 ) {
    for (INTEGER i=1; i<=rowsQ*colsQ; ++i) {
      ++k;
      Q[i] = iRparms[k];
    }
  }

  if (v > 0 ) {
    for (INTEGER i=1; i<=rowsV*colsV; ++i) {
      ++k;
      V[i] = iRparms[k];
    }
  }

  if (w > 0 ) {
    for (INTEGER i=1; i<=rowsW*colsW; ++i) {
      ++k;
      W[i] = iRparms[k];
    }
  }
}
    
const realmat& libsnp::rfunc::operator() (
    const realmat& iy, const realmat& iu, const realmat& ix,
    const realmat& iduwb0, const kronprd& iduwB,
    realmat& dRwb0, realmat& dRwB,
    realmat& dRwRparms)
{
  realmat dRwR0;
  realmat dRwP;
  realmat dRwQ;
  realmat dRwV;  
  realmat dRwW;

  (*this)(iy,iu,ix,iduwb0,iduwB,dRwb0,dRwB,dRwR0,dRwP,dRwQ,dRwV,dRwW);
  
  if ( dRwRparms.get_rows()!=lR || dRwRparms.get_cols()!=lRparms) {
    dRwRparms.resize(lR,lRparms);
  }

  INTEGER k=0;

  for (INTEGER j=1; j<=lR; ++j) {
    ++k;
    for (INTEGER i=1; i<=lR; ++i) dRwRparms(i,k) = dRwR0(i,j);
  }

  if (p > 0 ) {
    for (INTEGER j=1; j<=rowsP*colsP; ++j) {
      ++k;
      for (INTEGER i=1; i<=lR; ++i) dRwRparms(i,k) = dRwP(i,j);
    }
  }

  if (q > 0 ) {
    for (INTEGER j=1; j<=rowsQ*colsQ; ++j) {
      ++k;
      for (INTEGER i=1; i<=lR; ++i) dRwRparms(i,k) = dRwQ(i,j);
    }
  }

  if (v > 0 ) {
    for (INTEGER j=1; j<=rowsV*colsV; ++j) {
      ++k;
      for (INTEGER i=1; i<=lR; ++i) dRwRparms(i,k) = dRwV(i,j);
    }
  }

  if (w > 0 ) {
    for (INTEGER j=1; j<=rowsW*colsW; ++j) {
      ++k;
      for (INTEGER i=1; i<=lR; ++i) dRwRparms(i,k) = dRwW(i,j);
    }
  }

  return R;
}

const realmat& libsnp::rfunc::operator() 
  (const realmat& iy, const realmat& iu, const realmat& ix) 
{

  /*
  This code was obtained by stripping out derivative computations from the
  full version.  The state vectors have derivative matrices in them.  The
  null matrices, e.g. those immediately below, become place keepers in those
  state vectors.  They should serve no other purpose if the stripping was
  done correctly.
  */
  
  realmat iduwb0;
  kronprd iduwB;
  realmat dRwb0;
  realmat dRwB;
  realmat dRwR0; 
  realmat dRwP; 
  realmat dRwQ;
  realmat dRwV;  
  realmat dRwW;

  squash sqsh;
  switch (squasher) {
    case 0: 
      sqsh.set_identity();
      break;
    case 1:
      sqsh.set_spline(inflection);
      break;
    case 2:
      sqsh.set_logistic(inflection);
      break;
    default:
      sqsh.set_identity();
  }

  INTEGER lH = ly*ly;              //Length of vec(H)

  //Accumulators: H accumulates variance terms, the others accumulate
  //derivatives.  Initialization occurs at first use. 

  realmat H;
  realmat dHwR0;
  realmat dHwb0;
  realmat dHwB;
  realmat dHwP;
  realmat dHwQ;
  realmat dHwV;
  realmat dHwW;

  realmat H0(ly,ly);              //H0 = R0*T(R0) (after R0 expanded from vech)
  realmat dH0wR0;
  
  intvec vechH(lR);               //Selects vech(H) from vec(H)
  intvec vecTH(lH);               //Selects vec(T(H)) from vec(H)

  realmat I(ly,ly,0.0);           //Identity of order ly
  for (INTEGER i=1; i<=ly; ++i) { 
    I(i,i) = 1.0;
  }

  {
    realmat eR0(ly,ly,0.0) ;      //Expanded R0
  
    for (INTEGER j=1; j<=ly; ++j) {
      for (INTEGER i=1; i<=j; ++i) {
        INTEGER R_ij = (j*j-j)/2+i;
        INTEGER H_ij = ly*(j-1)+i;
        INTEGER H_ji = ly*(i-1)+j;
        eR0(i,j) = R0[R_ij];
        eR0(j,i) = R0[R_ij];
        vechH[R_ij] = H_ij;
        vecTH[H_ij] = H_ji;
        vecTH[H_ji] = H_ij;
      }
    }
    
    H0 = eR0*T(eR0);
  }                                

  H = H0;
    
  if ( p > 0 ) { //There is an MA part and might be intercept or regression.

    if ( iy.get_rows() != ly || iu.get_rows() != ly ) {
      error("Error, rfunc, operator (), y or u has wrong number of rows.");
    }
  
    for (INTEGER k=p-1; k>0; --k) {
      Pstate[k]=Pstate[k-1];
    }
    Pstate[0] = sE_k(iy, iu, iduwb0, iduwB, false);

    switch (Ptype) {

      case 's':
        for (INTEGER k=0; k<p; ++k) {
          
          if (Pstate[k].initial) break;

          //Compute P_k*E_k*T(P_k), where E_k = e[k]*T(e[k]).
          //Accumulate in H.

          REAL P_k = P[k+1];            //Pstate indexes from zero,  P from one
          REAL P_k_P_k = pow(P_k,2);
          realmat e_k = Pstate[k].y - Pstate[k].u;

          realmat dswe(ly,ly,0.0);            //Jacobian of squasher
          for (INTEGER i=1; i<=ly; ++i) {
            e_k[i] = sqsh(e_k[i],dswe(i,i));  //e_k replaced by squashed
          }                                   //residuals. Better notation
                                              //would have been s_k hereafter,
                                              //but only want to make minimal 
                                              //changes to legacy code.
          realmat E = e_k*T(e_k);

          H += P_k_P_k*E;
        } 
        break;

      case 'd':
        for (INTEGER k=0; k<p; ++k) {
          
          if (Pstate[k].initial) break;

          //Compute P_k*E_k*T(P_k), where E_k = e[k]*T(e[k]).
          //Accumulate in H.

          realmat P_k = P("",k+1);   //Pstate indexes from zero,  P from one
          realmat e_k = Pstate[k].y - Pstate[k].u;

          realmat dswe(ly,ly,0.0);            //Jacobian of squasher
          for (INTEGER i=1; i<=ly; ++i) {
            e_k[i] = sqsh(e_k[i],dswe(i,i));  //e_k replaced by squashed
          }                                   //residuals. Better notation
                                              //would have been s_k hereafter,
                                              //but only want to make minimal 
                                              //changes to legacy code.

          for (INTEGER j=1; j<=ly; ++j) {
            for (INTEGER i=1; i<=ly; ++i) {
              H(i,j) += P_k[i]*e_k[i]*e_k[j]*P_k[j];
            }
          }
        } 
        break;

      case 'f':
        for (INTEGER k=0; k<p; ++k) {
  
          if (Pstate[k].initial) break;

          intvec idxPk = seq(ly*k+1,ly*k+ly);   //Selects P_k from P
          
          //Compute P_k*E_k*T(P_k), where E_k = e[k]*T(e[k]).
          //Accumulate in H.

          realmat P_k = P("",idxPk);
          realmat e_k = Pstate[k].y - Pstate[k].u;

          realmat dswe(ly,ly,0.0);            //Jacobian of squasher
          for (INTEGER i=1; i<=ly; ++i) {
            e_k[i] = sqsh(e_k[i],dswe(i,i));  //e_k replaced by squashed
          }                                   //residuals. Better notation
                                              //would have been s_k hereafter,
                                              //but only want to make minimal 
                                              //changes to legacy code.

          realmat Pe_k = P_k*e_k;

          H += Pe_k*T(Pe_k);
        }
        break;
    }
  }
  
  if ( v > 0 ) { //There is a leverage effect; might be intercept or regression.

    for (INTEGER k=v-1; k>0; --k) {
      Vstate[k]=Vstate[k-1];
    }
    Vstate[0] = sE_k(iy, iu, iduwb0, iduwB, false);

    switch (Vtype) {

      case 's':

        for (INTEGER k=0; k<v; ++k) {
          
          if (Vstate[k].initial) break;

          //Compute max(0,V_k*s(e_k))*max'(0,V_k*s(e_k))
          //Accumulate in H.

          //Notation:
          //e_k = residual
          //s_k = s(e_k)         dswe is Jacobian   s_k is squashed residual
          //x_k = V_k*s_k        V_k*I and s_k are Jacobians
          //f_k = max(0,x_k)     dfwx is Jacobian
          //ff  = f_k*T(f_k)     dffwf is Jacobian

          REAL V_k = V[k+1];        //Vstate indexes from zero,  V from one
          realmat e_k = Vstate[k].y - Vstate[k].u;

          realmat s_k(ly,1);
          realmat dswe(ly,ly,0.0);            
          for (INTEGER i=1; i<=ly; ++i) {
            s_k[i] = sqsh(e_k[i],dswe(i,i));  //Squash residuals
          }
           
          realmat x_k = V_k*s_k;

          smomax f(0.1);                    //f.max0(x) returns max(0,x)
                                            //f.max1(x) returns derivative
          realmat f_k(ly,1);
          for (INTEGER i=1; i<=ly; ++i) {
            f_k[i] = f.max0(x_k[i]);
          }  

          H += f_k*T(f_k);
        } 
        break;

      case 'd':
        for (INTEGER k=0; k<v; ++k) {
          
          if (Vstate[k].initial) break;

          //Compute max(0,V_k*s(e_k))*max'(0,V_k*s(e_k))
          //Accumulate in H.

          //Notation:
          //e_k = residual
          //s_k = s(e_k)         dswe is Jacobian   s_k is squashed residual
          //x_k = diag(V_k)*s_k  diag(V_k) and diag(s_k) are Jacobians
          //f_k = max(0,x_k)     dfwx is Jacobian
          //ff  = f_k*T(f_k)     dffwf is Jacobian

          realmat V_k = V("",k+1);   //Vstate indexes from zero,  V from one
          realmat e_k = Vstate[k].y - Vstate[k].u;

          realmat s_k(ly,1);
          realmat dswe(ly,ly,0.0);            
          for (INTEGER i=1; i<=ly; ++i) {
            s_k[i] = sqsh(e_k[i],dswe(i,i));  //Squash residuals
          }
           
          realmat x_k = diag(V_k)*s_k;

          smomax f(0.1);                    //f.max0(x) returns max(0,x)
                                            //f.max1(x) returns derivative
          realmat f_k(ly,1);
          for (INTEGER i=1; i<=ly; ++i) {
            f_k[i] = f.max0(x_k[i]);
          }  

          H += f_k*T(f_k);
        } 
        break;

      case 'f':
        for (INTEGER k=0; k<v; ++k) {
  
          if (Vstate[k].initial) break;

          intvec idxVk = seq(ly*k+1,ly*k+ly);   //Selects V_k from V

          //Compute max(0,V_k*s(e_k))*max'(0,V_k*s(e_k))
          //Accumulate in H.

          //Notation:
          //e_k = residual
          //s_k = s(e_k)       dswe is Jacobian   s_k is squashed residual 
          //x_k = V_k*s_k      V_k and krnprd(T(s),I) are Jacobians
          //f_k = max(0,x_k)   dfwx is Jacobian
          //ff  = f_k*T(f_k)   dffwf is Jacobian

          realmat V_k = V("",idxVk);
          realmat e_k = Vstate[k].y - Vstate[k].u;

          realmat s_k(ly,1);
          realmat dswe(ly,ly,0.0);            
          for (INTEGER i=1; i<=ly; ++i) {
            s_k[i] = sqsh(e_k[i],dswe(i,i));  //Squash residuals
          }

          realmat x_k = V_k*s_k;

          smomax f(0.1);                    //f.max0(x) returns max(0,x)
                                            //f.max1(x) returns derivative
          realmat f_k(ly,1);
          for (INTEGER i=1; i<=ly; ++i) {
            f_k[i] = f.max0(x_k[i]);
          }  

          H += f_k*T(f_k);
        }
        break;
    }
  }
  
  if ( q > 0 ) { //There is an AR part. 
    
    switch (Qtype) {

      case 's':
        for (INTEGER k=0; k<q; ++k) {

          REAL Q_k = Q[k+1];                    //Select Q_k from Q 
          REAL Q_k_Q_k = pow(Q_k,2); 

          if (Qstate[k].initial) {
          
            //Compute Q_k*H_k*T(Q_k), where H_k = H0.
            //Accumulate in H.

            H += Q_k_Q_k*H0;
          }
          else {
          
            //Compute Q_k*H_k*T(Q_k), where H_k = Qstate[k].H.
            //Accumulate in H.

            H += Q_k_Q_k*Qstate[k].H;
          }
        } // end for loop

        for (INTEGER k=q-1; k>0; --k) {
          Qstate[k]=Qstate[k-1];
        }
        Qstate[0] = sH_k(H, dHwb0, dHwB, dHwR0, dHwP, dHwQ, dHwV, dHwW, false);

        break;

      case 'd':
        for (INTEGER k=0; k<q; ++k) {

          realmat Q_k = Q("",k+1);              //Select Q_k from Q 

          if (Qstate[k].initial) {
          
            //Compute Q_k*H_k*T(Q_k), where H_k = H0.
            //Accumulate in H.

            for (INTEGER j=1; j<=ly; ++j) {
              for (INTEGER i=1; i<=ly; ++i) {
                H(i,j) += Q_k[i]*H0(i,j)*Q_k[j];
              }
            }
          }
          else {
          
            //Compute Q_k*H_k*T(Q_k), where H_k = Qstate[k].H.
            //Accumulate in H.

            for (INTEGER j=1; j<=ly; ++j) {
              for (INTEGER i=1; i<=ly; ++i) {
                H(i,j) += Q_k[i]*Qstate[k].H(i,j)*Q_k[j];
              }
            }
          }

        }  // end for loop

        for (INTEGER k=q-1; k>0; --k) {
          Qstate[k]=Qstate[k-1];
        }
        Qstate[0] = sH_k(H, dHwb0, dHwB, dHwR0, dHwP, dHwQ, dHwV, dHwW, false);
          
        break;
        
      case 'f':
        for (INTEGER k=0; k<q; ++k) {

          intvec idxQk = seq(ly*k+1,ly*k+ly);   //Selects Q_k from Q
        
          realmat Q_k = Q("",idxQk);            //Select Q_k from Q 

          if (Qstate[k].initial) {
          
            //Compute Q_k*H_k*T(Q_k), where H_k = H0.
            //Accumulate in H.

            H += Q_k*H0*T(Q_k);
          }
          else {
          
            //Compute Q_k*H_k*T(Q_k), where H_k = Qstate[k].H.
            //Accumulate in H.

            H += Q_k*Qstate[k].H*T(Q_k);
          }
        } // end for loop

        for (INTEGER k=q-1; k>0; --k) {
          Qstate[k]=Qstate[k-1];
        }
        Qstate[0] = sH_k(H, dHwb0, dHwB, dHwR0, dHwP, dHwQ, dHwV, dHwW, false);
          
        break;
    }
  }
 
  if ( w > 0 ) {      //There is a level effect; the level effect is
                      //not run through the GARCH recursion although
                      //dHwW is in Qstate in case that changes in which
                      //case this code must be moved after code for V 
                      //and before code for Q
  
    for (INTEGER k=w-1; k>0; --k) {
      Wstate[k]=Wstate[k-1];
    }
    Wstate[0] = sX_k(ix, false);

    switch (Wtype) {

      case 's':
        for (INTEGER k=0; k<w; ++k) {
          
          if (Wstate[k].initial) break;

          //Compute W_k*X_k*T(W_k), where X_k = x[k]*T(x[k]).
          //Accumulate in H.

          REAL W_k = W[k+1];            //Wstate indexes from zero,  W from one
          REAL W_k_W_k = pow(W_k,2);
          realmat x_k = Wstate[k].x;

          realmat X = x_k*T(x_k);

          H += W_k_W_k*X;
        } 
        break;

      case 'd':
        for (INTEGER k=0; k<w; ++k) {
          
          if (Wstate[k].initial) break;

          //Compute W_k*X_k*T(W_k), where X_k = x[k]*T(x[k]).
          //Accumulate in H.

          realmat W_k = W("",k+1);   //Wstate indexes from zero,  W from one
          realmat x_k = Wstate[k].x;

          for (INTEGER j=1; j<=lx; ++j) {
            for (INTEGER i=1; i<=lx; ++i) {
              H(i,j) += W_k[i]*x_k[i]*x_k[j]*W_k[j];
            }
          }
        } 
        break;

      case 'f':
        for (INTEGER k=0; k<w; ++k) {
  
          if (Wstate[k].initial) break;

          intvec idxWk = seq(lx*k+1,lx*k+lx);   //Selects W_k from W
          
          //Compute W_k*X_k*T(W_k), where X_k = x[k]*T(x[k]).
          //Accumulate in H.

          realmat W_k = W("",idxWk);
          realmat x_k = Wstate[k].x;

          realmat Wx_k = W_k*x_k;

          H += Wx_k*T(Wx_k);
        }
        break;
    }
  }
 
  //Factor H
  
  realmat dRwS(lR,lR);               //S is the upper triangle of H

  for (INTEGER j=1; j<=ly; ++j) {
    for (INTEGER i=1; i<=j; ++i) {
      INTEGER R_ij = (j*j-j)/2+i;
      INTEGER H_ij = ly*(j-1)+i;
      R[R_ij] = H[H_ij];             //R now contains the upper triangle of H
    }
  }
 
  INTEGER ier = factor(R.get_x(),dRwS.get_x(),ly);

  if (factor_warn && (ier == -1)) {
    warn("Warning, rfunc, operator (), factorization failure");
    factor_warn = false;
  }

  return R;
}

REAL libsnp::rfunc::stability() const
{
  REAL rv = 0.0;

  INTEGER max_lag = p;
  max_lag = q > max_lag ? q : max_lag;

  if (max_lag == 0) return rv;

  INTEGER  ly2 = ly*ly; 

  realmat companion(ly2*max_lag,ly2*max_lag,0.0);
  for (INTEGER i=1; i<=ly2*(max_lag-1); ++i) companion(ly2+i,i) = 1.0;

  for (INTEGER k=1; k<=p; ++k) {
    realmat Pk(ly,ly,0.0);
    if (Ptype == 's') {
      for (INTEGER i=1; i<=ly; ++i) Pk(i,i) = P[k];
    }
    else if (Ptype == 'd') {
      for (INTEGER i=1; i<=ly; ++i) Pk(i,i) = P[ly*(k-1) + i];
    }
    else if (Ptype == 'f') {
      for (INTEGER j=1; j<=ly; ++j) {
        for (INTEGER i=1; i<=ly; ++i) {
          Pk(i,j) = P[ly2*(k-1) + ly*(j-1) + i];
        }
      }
    }
    else {
      error("Error, rfunc, stability, this should never happen");
    }

    kronprd PPk(Pk,Pk);

    for (INTEGER j=1; j<=ly2; ++j) {
      for (INTEGER i=1; i<=ly2; ++i) {
        companion(i,ly2*(k-1)+j) += PPk(i,j);
      }
    }
  }

  for (INTEGER k=1; k<=q; ++k) {
    realmat Qk(ly,ly,0.0);
    if (Qtype == 's') {
      for (INTEGER i=1; i<=ly; ++i) Qk(i,i) = Q[k];
    }
    else if (Qtype == 'd') {
      for (INTEGER i=1; i<=ly; ++i) Qk(i,i) = Q[ly*(k-1) + i];
    }
    else if (Qtype == 'f') {
      for (INTEGER j=1; j<=ly; ++j) {
        for (INTEGER i=1; i<=ly; ++i) {
          Qk(i,j) = Q[ly2*(k-1) + ly*(j-1) + i];
        }
      }
    }
    else {
      error("Error, rfunc, stability, this should never happen");
    }

    kronprd QQk(Qk,Qk);

    for (INTEGER j=1; j<=ly2; ++j) {
      for (INTEGER i=1; i<=ly2; ++i) {
        companion(i,ly2*(k-1)+j) += QQk(i,j);
      }
    }
  }

  realmat A(ly2,lR,0.0);
  for (INTEGER j=1; j<=ly; ++j) {
    for (INTEGER i=1; i<=j; ++i) {
      A(ly*(j-1) + i, (j*j-j)/2 + i) = 1.0;
      A(ly*(i-1) + j, (j*j-j)/2 + i) = 1.0;
    }
  }

  realmat B(lR,ly2,0.0);
  for (INTEGER j=1; j<=ly; ++j) {
    for (INTEGER i=1; i<=j; ++i) {
      B((j*j-j)/2 + i, ly*(i-1) + j) = 1.0;
    }
  }

  realmat I(max_lag,max_lag,0.0);
  for (INTEGER i=1; i<=max_lag; ++i) {
    I(i,i) = 1.0;
  }

  kronprd IA(I,A);
  kronprd IB(I,B);

  companion = IB*companion*IA;

  realmat workspace;
  INTEGER ier;
  rv = eigen(companion,workspace,ier);
  if (ier != 0) warn("Warning, rfunc, eigen value computation failed");

  return rv;
}

