models.cpp

#include "models.h"     

/////////////////////////////////////////////////////////
//  functions belonging to the class MODEL:            //
/////////////////////////////////////////////////////////


safevector<double> MODEL::simulatevector(long n) {
  safevector<double> res(n);
  matrix simdata = simulate(n);
  for (long i=0; i<n; i++) 
    res[i] = simdata.get(i);
  return res;
}

/////////////////////////////////////////////////////////
//  functions belonging to the class MARKOVMODEL:      //
/////////////////////////////////////////////////////////

matrix MARKOVMODEL::simulate(long n)  
{
  matrix X(n);
  long i;
  X.set(startvalue,0);
  for(i=1;i<=n-1;i++)
    X.set(next(X.get(i-1)),i);
  return(X);
}


///////////////////////////////////////////////////////
//  functions belonging to the class SVMODEL:  //
///////////////////////////////////////////////////////

matrix SVMODEL::simulateS(long n) 
{ 
  matrix S(n);
  matrix v;
  double tmp;
  long i,j;

  setinitvalue();
  
  for(i=0;i<=n-1;i++)
    {
      v = MARKOVMODEL::simulate(PointsPerDelta);
      tmp=0.5*(v.get(0)+v.get(PointsPerDelta-1));
      for(j=1;j<=PointsPerDelta-1;j++)
        {
          tmp += v.get(j); 
        }
      S.set(delta*tmp,i);
      MARKOVMODEL::setstartvalue(v.get(PointsPerDelta-1));
    }
  return(S);
}

matrix SVMODEL::simulate(long n)
{
  matrix S(n);
  matrix Y(n);
  long i;

  S = simulateS(n);

  for(i=0;i<=n-1;i++)
    Y.set(sqrt(S.get(i))*rnorm2(seed),i,0);
  return(Y);
}

////////////////////////////////////////////////////////////
//  functions belonging to the class PREESTMODEL:         //
////////////////////////////////////////////////////////////


//#############
//############# WE ONLY NEED TO CHANCE THEESE TWO FUNCTIONS IN ORDER TO EXTEND THE PROGRAM TO MORE GENEREAL
//############# CHOICES OF Z_j^i. PREESTMODEL::Zval & PREESTMODEL::calculate_Z
//############# 
//############# The smart way to do this is to create a subclass with theese two functions. 
//############# If for instance you create a new class MYMODEL inheriting from PREESTMODEL, and you have written the appropiate MYMODEL::next function (simulation), then it automatically inherits the canonic choice of Z-functions. If you however to choose to look at other Z-functions, then you just need to create a subclass MODEL1 that inherits from MYMODEL, and then simply reimplent MYMODEL1::Zval and MODEL1::calculate_Z. 

double PREESTMODEL::Zval(unsigned j, long i, long k, matrix const& Y) { // Z_{jk}^{(i)}, j=0,...,N-1, k=1,...,q, i
  double (*fj)(double) = ifuncs[j];
  double x = Y.get(i-k);
  return fj(x);
  //  return fdata[j].get(i-k); // f_j(X_{i+1-k})
}

void PREESTMODEL::calculate_Z() {
 
  Z.clear();
  fXnext.clear();
  fX1.clear();
  for (unsigned j=0; j<N; j++) {
    double (*fj)(double) = ifuncs[j];
    Z.push_back(X.map(fj)); // Should probably be calculated within Zval. Then we only need to chance Zval in the future.
    fXnext.push_back(Xnext.map(fj));
    fX1.push_back(X1.map(fj));
  }
}

//#############
//############# ^^^^^^^^^^^^^^ READ ABOVE COMMENT
//############# 
//#############


void PREESTMODEL::setdata(matrix Y) {
  MODEL::setdata(Y); // Setting data
  fdata.clear();
  for (unsigned j=0; j<N; j++) {
    double (*fj)(double) = ifuncs[j];
    fdata.push_back(Y.map(fj));
  }
}

void PREESTMODEL::calculate_X() {
  //  setseed(-12345);
  // Initializing matrices...
  X.newsize(numsimul,q);
  Xnext.newsize(numsimul);
  X1.newsize(numsimul);
  
  assert(s>=q+1);
  long start = s-1;
  for (unsigned k=0; k<numsimul; k++) {
    matrix Xs = simulate(s+1); // We simulate X_1,...,X_{s+1} each numsimul times.
    Xnext.set(Xs.get(s),k); // X_{s+1}
    X1.set(Xs.get(0),k);  // X_1
    for (long i=0; i<q; i++) { 
      X.set(Xs.get(start-i),k,i); // X_s,...,X_{s+1-q}
    }
  }
}

void PREESTMODEL::calculate_C() {
  C.clear();
  for (unsigned j=0; j<N; j++) {
    C.push_back(Z[j].covar());
  }
}

void PREESTMODEL::calculate_b()
{
  b.clear();
  for (unsigned j=0; j<N; j++) {
    b.push_back(covar2(Z[j],fXnext[j]));
  }
}

void PREESTMODEL::calculate_a()
{
  setseed(-123456); // Fixed seed. Since we need to calculate the derivative of alpha.

  calculate_X();
  calculate_Z();
  calculate_b();
  calculate_C();
  
  a.clear();
  a0.clear();
  EtildeZ.clear();
  for (unsigned j=0; j<N; j++) {
    matrix aj = pseudoinverse(C[j])*b[j];
  
    double EfX1 = (fX1[j].mean())[0];
    safevector<double> EZj = Z[j].mean();
    
    matrix EtildeZj(q+1); // MS. bottom p.9 
    EtildeZj.set(1,0);
    for (long kk=1; kk<=q; kk++) {
      EtildeZj.set(EZj[kk-1], kk);
    }
    EtildeZ.push_back(EtildeZj);

    double prod = 0; // Calculating inner-product <a,EZj>
    for (int i=0; i<q; i++) {
      prod += aj.get(i)*EZj[i];
    }
    double aj0 = EfX1 - prod;  // a_j0
    a.push_back(aj);
    a0.push_back(aj0);
  }
  b.clear();  
}


double PREESTMODEL::prediction(unsigned j, long i, matrix const& Y)
{

  assert (j<N);
  assert (i>=1 && i<=Y.row());

  double res = a0[j];
  for(int k=1; k<=q; k++) {
    res -= a[j].get(k-1)*Zval(j,i,k, Y); // <a_j|Z_j^{i-1}>
  }
  return(res);
}


matrix PREESTMODEL::sumH(matrix par)

{
  double res;
  matrix sumH(N*(q+1));

  setpar(par);
  long n=getdatasize();
  
  calculate_a();

  for (long i=s+1; i<=n; i++) // sum_{i=s+1}^n H^i
    { 
      for (unsigned j=0; j<N; j++) {// calculate Z_{j0}, j=1,...,N
        double residual = fdata[j].get(i-1)-prediction(j,i-1,getdata()); // f_j(X_i)-pi_j(i-1)
        res = Zval(j,i-1,0,getdata())*residual; // Z_{j0}^{(i-1)}
        sumH.add(res,j*(q+1)); // Z_{j0} on row j(q+1)
        for(long k=1; k<=q; k++) {
          res = Zval(j,i-1,k,getdata())*residual;
          sumH.add(res,j*(q+1)+k); // Z_{jk}
        }
      }
    }    

  if (printdebug) {
    std::cerr << "\n\nsumH=\n";
    sumH.stdprint();
    std::cerr << "\n";    
  }

    
  return(sumH);
}


//! First term in \f$\overline{M}\f$ : \f$E_\theta(H^{(r)}(\theta)H^{(r+i)}(\theta))\f$. Under the assumption that our incoming process is ergodic, we can approximate this term by calculating the mean.
void PREESTMODEL::calculate_EHH() {

  long m=l;
  if(q>l) m=q;
  
  matrix H(N*(q+1)); // Simulation of H
  matrix M(N*(q+1), N*(q+1)); 

  simY=simulate(MbarNOS+1);

  for(long i=m; i<MbarNOS; i++) {
    //!   Iteration no. i

    for (unsigned j=0; j<N; j++) {
      double (*fj)(double) = ifuncs[j];
      for (long k=0; k<=q; k++) {
        //! OBS: \breve{\alpha} has already been calculated, so we can use prediction directly without calling calculate_a();
        double res = fj(simY.get(i))-prediction(j,i-1,simY);
        if (k>0) // Pr. defintion: Z_{j0}^{(i-1)} = 1
          res *= Zval(j,i-1,k,simY); 
        // index (j,k) is placed in H at location (j-1)*(q+1)+k (We index from 0).
        unsigned Hindex = j*(q+1)+k;
        H.set(res,Hindex);
      }
    } // Now we have calculated H^(i)(theta)


    // Next we calculate the mean*normalization.
    for (unsigned r=0; r<N*(q+1); r++) {
      for (unsigned s=0; s<N*(q+1); s++) {     
        M.add(H.get(r)*H.get(s),r,s); // Here we add (H*H^t)_(rs) = (H)_r*(H)_s to the matrix M.
      }
    }
  }

  double w = (double)1/(double)(MbarNOS - m +1);  
  EHH = M*w; 

  qmove(1,25);
  varH.newsize(N*(q+1));
  
  for (unsigned k=0; k<N*(q+1); k++)                                
    varH.set(EHH.get(k,k),k);                        // the diagonal of E((H^r)(H^r)^T) is the
                                                     // variances of the coordinates of H^r.
                                                     // This is the reason why we break the
                                                     // calculation of Mbar up into two pieces.
}

void PREESTMODEL::set_Mbar_terms(double lambda)
{
  double K;

  K=varH.maxvalue();

  for(long k=0; k<MbarNOS-q;k++){
    if(K/exp(lambda*(k-q))<0.0001){
      Mbar_terms=k;
      break;}
  }

}

void PREESTMODEL::calculate_Mbar()
{
  int m=l;

  matrix H1(N*(q+1)),H2(N*(q+1)),M(N*(q+1),N*(q+1));
  calculate_EHH();
  set_Mbar_terms(sumH_est);                         // lambda=theta is set to sumH_est. The estimate
                                                    // is based on sumH, and should be found prior to calling estfunc.
  Mbar=EHH;

  if(q>l) m=q;

  std::string ss = "Mbar_terms = " + numStr(Mbar_terms);
  if (printdebug)
    qmvaddstr(30,50, ss.c_str());  qclrtoeol();

  qmove(29,12);
  for(long k=1;k<=Mbar_terms;k++){
    M=M*0;

    for(long i=m; i<MbarNOS-k; i++) {
      
      
      for (unsigned j=0; j<N; j++) {
        double (*fj)(double) = ifuncs[j];
        for (long r=0; r<=q; r++) {
          // OBS: \breve{\alpha} has already been calculated, so we can use prediction directly without calling calculate_a();
          double res1 = fj(simY.get(i))-prediction(j,i-1,simY);
          double res2 = fj(simY.get(i+k))-prediction(j,i+k-1,simY);
          if (k>0) { // Pr. defintion: Z_{j0}^{(i-1)} = 1 
            res1 *= Zval(j,i-1,r,simY); 
            res1 *= Zval(j,i+k-1,r,simY); 
          }
          // index (j,k) is placed in H at location (j-1)*(q+1)+k (We index from 0).
          unsigned Hindex = j*(q+1)+r;
          H1.set(res1,Hindex);
          H2.set(res2,Hindex);
        }
      } // Now we have calculated H^(i)(theta) and H^(i+r)(theta)
          
      M = M+(H1*H2.trans());    
    }

    M=M+M.trans(); // See MS(3.14).
    double w1=(double) 1/(MbarNOS-k-m);                          
    double w2=(MbarNOS-q-k);
    double w3=(double) 1/(MbarNOS-q);
    M=M*w1;
    Mbar=Mbar+M*w2*w3; 
  }

  qmove(26,1); qclrtoeol();
  qmove(27,1); qclrtoeol();
  qmove(28,1); qclrtoeol();
  qmove(29,1); qclrtoeol();
}



matrix PREESTMODEL::calculate_Cbar() {
  matrix Cbar(N*(q+1), N*(q+1));

  for (unsigned j=0; j<N; j++) {
    matrix EZ2 = EtildeZ[j]*EtildeZ[j].trans();
    double val;
    for (long r=0; r<=q; r++) {
      for (long s=0; s<=q; s++) {
        long super_r = j*(q+1)+r;
        long super_s = j*(q+1)+s;
        if (r==0 || s==0) {
          val = EZ2.get(r,s);
        }
        else {
          val = C[j].get(r-1,s-1) + EZ2.get(r,s);
        }
        Cbar.set(val,super_r,super_s);
         
      }
    }
  }

  return Cbar;
}


matrix PREESTMODEL::calculate_Dalpha() { // Estimating derivatives of function alpha
  unsigned p = noofpar();
  matrix res(N*(q+1),p);
  double epsilon = 0.01;
  matrix origpar = getpar();

  safevector<matrix> a_orig = a;
  safevector<double> a0_orig = a0;

  
  for (unsigned i=0; i<p; i++) {
    matrix newpar = origpar;
    newpar.set(newpar.get(i)+epsilon,i); // used to estimate i'th partial derivative of alpha
    setpar(newpar);

    calculate_a();

    for (unsigned j=0; j<N; j++) {  
      matrix cur_aj = a[j];
      matrix orig_aj = a_orig[j];
      double Da0 = (a0[j]-a0_orig[j])/epsilon; // (f(x+e)-f(x))/e ~= f'(x)
      res.set(Da0,j*(q+1),i);
      for (int k=0; k<q; k++) {
        double Daj = (cur_aj.get(k)-orig_aj.get(k))/epsilon;
        unsigned rr= (unsigned)j*(q+1)+k+1;
        res.set(Daj,rr,i);
      }
    }
  }

  setpar(origpar);
  return res;
}



void PREESTMODEL::calculate_Astar() {
  if (printdebug)
    qmvaddstr(30,50, "Calculating Astar:");    
  matrix Cbar = calculate_Cbar();
  if (printdebug)
    qmvaddstr(30,50, "Calculating Mbar:");  qclrtoeol();  
  calculate_Mbar();
  if (printdebug)
    qmvaddstr(30,50, "Calculating MbarInv (using pseudoinverse):");
  matrix MbarInv = pseudoinverse(Mbar);
  if (printdebug)
    qmvaddstr(30,50, "Calculating Dalpha:"); qclrtoeol();
  matrix Dalpha = calculate_Dalpha();

  matrix U = Cbar*Dalpha;
  Astar = U.trans()*MbarInv;
  if (printdebug)
    qmvaddstr(30,50, ""); qclrtoeol();

  if (printdebug) {
    std::cerr << "\nMbarInv=\n";
    MbarInv.stdprint();
    std::cerr << "\nU=\n";
    U.stdprint();
    std::cerr << "\nDalpha=\n";
    Dalpha.stdprint();

  //   std::cerr << "Dalpha.trans = " << Dalpha.row() << "x" << Dalpha.column() << std::endl;
  //   std::cerr << "dim(Cbar) = " << Cbar.row() << "x" << Cbar.column()<< std::endl;
  //   std::cerr << "dim(MbarInv) = " << MbarInv.row() << "x" << MbarInv.column()<< std::endl;
  //   std::cerr << "Dalpha=\n";
  //   Dalpha.stdprint();
    std::cerr << std::endl << std::endl;
  }

}


matrix PREESTMODEL::estfunc(matrix par)
{
  matrix sumHH=sumH(par);
  setpar(par);  // Pre-estimate

  // We only calculate Astar (and hence Mbar) if iter_Mbar is set.
  if (iter_Mbar) {
    calculate_Astar();
  }
  return Astar*sumHH;
}



/////////////////////////////////////////////////////
//  functions belonging to the class SVCIR:        //
/////////////////////////////////////////////////////


void SVCIR::setpar(matrix par) {
  alpha=par.get(0); theta=par.get(1); sigma=par.get(2); MODEL::setpar(par);
}

unsigned SVCIR::getdim() { return(3); }

double SVCIR::next(double oldvalue) const
{ 
  double newvalue;
  double deltaW;

  deltaW=sqrt(delta)*rnorm2(seed);
  newvalue = oldvalue - theta*(oldvalue-alpha)*delta+ 
    sigma*sqrt(oldvalue)*deltaW + 
    ((double) sigma*sigma/4)*(deltaW*deltaW-delta);
  return(newvalue);
}

/////////////////////////////////////////////////////
//  functions belonging to the class SVEOU:        //
/////////////////////////////////////////////////////

void SVEOU::setpar(matrix par) {
  alpha=par.get(0); theta=par.get(1); sigma=par.get(2); MODEL::setpar(par);
} 

unsigned SVEOU::getdim() { return(3); }

double SVEOU::next(double oldvalue) const
{
  double newvalue;
  double DeltaW;
  
  DeltaW =  sqrt(sigma*sigma/(2*theta)*(1-exp(-2*theta*delta)))*rnorm2(seed);
  newvalue = exp(alpha*(1-exp(-theta*delta)) + exp(-theta*delta)*log(oldvalue) + DeltaW);

  return(newvalue);
}


/////////////////////////////////////////////////////
//  functions belonging to the class SVEOU1:        //
/////////////////////////////////////////////////////


double SVEOU1::Zval(unsigned j, long i, long k, matrix const& Y) { // Z_{jk}^{(i)}, j=0,...,N-1, k=1,...,q, i
  double (*fj)(double) = ifuncs[j];
  double x = Y.get(i-k);
  return sqrt(fj(x));
}

void SVEOU1::calculate_Z() {
 
  Z.clear();
  fXnext.clear();
  fX1.clear();
  for (unsigned j=0; j<N; j++) {
    double (*fj)(double) = ifuncs[j];
    Z.push_back(X.map(fj)); // Should probably be calculated within Zval. Then we only need to chance Zval in the future.
    fXnext.push_back(Xnext.map(fj));
    fX1.push_back(X1.map(fj));
  }
}

////////////////////////////////////////////////////////////////////////////////


#ifdef TEST

int main() {
  try {
    std::cerr << "Nothing to test in models.cpp\n";
  } catch(std::string const& err) {
    std::cerr << err << std::endl;
    return 1;
  }
  return 0;
}


#endif//TEST

---