models.h

Go to the documentation of this file.

/*!
  \file   models.h
  \author Klaus Holst
  \date   Wed Oct 13 14:21:51 2004
  
  \brief  Stochastic model classes and corresponding Prediction-Based estimating function.
    
*/

#ifndef _MODELS_H_
#define _MODELS_H_

#include <cmath>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <ctime>
#include "linalg.h"       // various (matrix)functions
#include "linrdist.h"     // Various random number generators
#include "qscr.h"         // Screen output
#include "svd.h"          // Singular Value Decomposition (and pseudoinverse)
#include "util.h"         // misc. utilities



typedef safevector<double(*)(double)> functionvector;


//////////////////////////
//      CLASSES         //
/////////////////////////


//! Abstract model class
/*! 
  
*/
class MODEL {
 public:
  //! Constructor
  /*!
  
  \return 
  */
  MODEL() {}
  //! Destructor
  /*!
  
  \return 
  */
  virtual ~MODEL() {}
  //! Virtual simulation function
  /*! 
    
  \param n : number of variables simulated
 
  \return (VIRTUAL)
  */
  virtual matrix simulate(long n) =0;                    
  //! Virtual simulation function
  /*! 
    
  \param n : number of variables simulated
  
  \return vector of variables opposed to function simulate which returns matrix.
  */
  virtual unsigned getdim() =0;
  //! Returns dimension of parameterspace.
  safevector<double> simulatevector(long n);                    
  //! sets the seed.
  void setseed(long s){ seed = &seedvar; seedvar = s;  }
  //! sets the parameters.
  virtual void setpar(matrix par) {parameters=par;}
  //! returns the parameters.
  matrix getpar() {return(parameters);}
  //! returns the number of parameters.
  long noofpar() {return(parameters.row());}
  //! returns the data.
  matrix getdata() {return(data);}
  //! returns the datasize.
  long getdatasize() {return(data.row());}
  //! sets the data.
  void setdata(matrix Y) {data=Y;}
  //! the estimation function described in MS.
  /*! 
    
  \param par : the parameters of the model.
  
  \return 
  */
  virtual matrix estfunc(matrix par) =0;
 protected:
  //! matrix holding the data.
  matrix data;

  //! matrix holding the parameters.
  matrix parameters;
  //! variable helping to set the seed.
  long seedvar;
  //! the seed used in the simulations
  long *seed;

};


//! Pre-est. model
/*! 
  inherits from MODEL
*/
class PREESTMODEL: public MODEL {
 public:
  //! sets the data and calculates f_j(data).
  void setdata(matrix Y);
  //! the estimating function MS(2.1)
  matrix estfunc(matrix par);
  //! Returns two-norm of estimating function with parameters given by par.
  double twonorm_estfunc( matrix par ){ return(estfunc(par).twonorm()); }
  //! MS(3.11) disregarding A(theta).
  matrix sumH(matrix par);
  //! Returns two-norm of _approximated_ estimating function with parameters given by par
  double twonorm_sumH( matrix par ){ return(sumH(par).twonorm()); }
  //! calculates and sets Mbar_terms.
  void set_Mbar_terms(double lambda);
  //! Determines wether Mbar (time-consuming calculation) should be calculated.
  void set_iter_Mbar(bool val) { iter_Mbar = val; }
  //! Sets pre-estimate of sumH to x.
  void set_sumH_est(double x) { sumH_est = x; }

  //! Calculates matrix \f$\mathbf{X}\f$ with \f$j\f$th row containing observations \f$(x_{s},\ldots,x_{s+1-q}).\f$
  void calculate_X();

  matrix get_X() { return X; }
  //! Returns the number of inner-functions in pbef.
  unsigned get_N() { return N; }
  //! Z_{jk}^{(i)}, j=0,...,N-1, k=1,...,q, i
  virtual double Zval(unsigned j, long i, long k, matrix const& Y);
  //! Returns \f$\breve{\pi}_{j+1}^{(i-1)}(\theta) = \breve{\alpha}_{(j+1)0}(\theta) + \breve{\alpha}_{j+1}(\theta)^TZ_{j+1}^{(i-1)}\f$. Here \f$Z_j^{(i-1)}\f$ is our actual data.
  double prediction(unsigned j, long i, matrix const& Y);
  //! boolean to determine wether to print debug-info or not
  bool printdebug;
  
 protected:
  //! matrix holding f_j(data), j=1,...,N.
  safevector<matrix> fdata;
  //! The number of inner-functions, \f$f_1,...,f_N\f$, in the estimating function.
  unsigned N;
  //! inner-functions in estimating-function G.
  functionvector ifuncs;
  //! The lag-length in the prediction.
  long q;
  //! Summation start-index, \f$s\f$, of estimating function.
  long s;
  //! The number of simulations used in calculating Mbar.
  long MbarNOS;
  //! summation-start in the expresion of Mbar.
  long l;
  //! Number of simulated variables used in determination of \f$\breve{a}_j(\theta), C_j(\theta)\f$, etc.
  unsigned numsimul;
  
  //! For fixed \f$s\f$ the \f$numsimul\times q\f$ matrix with \f$j\f$th row given by observations \f$x_{s},\ldots, x_{s+1-q})\f$.
  matrix X;
  // A numsimul-vector with observations of \f$X_{s+1}\f$.
  matrix Xnext;
  // A numsimul-vector with observations of \f$X_{1}\f$.
  matrix X1;
  //! Boolean that determines wether to calculate the time-consuming Mbar in current iteration
  bool iter_Mbar;
  //! simulated Y for calculations.
  matrix simY;
  //! \f$A_n^*(\theta)\f$
  matrix Astar;

  //! On index j: holds simulations of \f$Z_{j1},\ldots,Z_{jq}\f$, i.e. \f$f_j\f$ mapped onto matrix X. 
  safevector <matrix> Z;
  // On index j thibxs vector holds q-vector \f$(1,EZ_{j1},...,E_{jq})\f$
  safevector <matrix> EtildeZ;
  //! On index j: Holds simulations of \f$f_j(X_{s+1})\f$.
  safevector <matrix> fXnext;
  //! On index j: Holds simulations of \f$f_j(X_{1})\f$.
  safevector <matrix> fX1;
  //!  Terms used in projection. \f$\breve{\alpha}_j\f$. MS(2.4)
  safevector <matrix> a;
  //! \f$\breve{\alpha}_{j0}(\theta)=E_\theta(f_j(X_1)) - \breve{\alpha}_j(\theta)^TE_\theta(Z_j^{(s)})\f$. MS(2.5)
  safevector <double> a0;
  //! defined (as b_j) in MS on the top of page 9.
  safevector <matrix> b;
  //! Covariance matrix of Z_j (as C_j). See MS on the top of page 9.
  safevector <matrix> C;
  //! EHH in Mbar (MS(3.14)).
  matrix EHH;
  //! columnvector containing the simulated value of var(H_i) in the i'th row.
  matrix varH;
  //! estimate of theta, based on sumH.
  double sumH_est;
  //! The number of terms in the expression of Mbar.
  int Mbar_terms;
  //! MS(3.14)
  matrix Mbar;

  virtual void calculate_Z();
  //! Calculates the terms \f$a_{j,0},\ldots,a_{j,q}\f$ used in calculation of projections \f$\pi_j\f$.
  void calculate_a();
  //! Calculates \f$b_j(\theta) = \left(COV_\theta(Z_{j1}^{(s)},f_j(X_{s+1})),\ldots, COV_\theta(Z_{jq}^{(s)},f_j(X_{s+1})))\right)^T \f$
  void calculate_b();
  //! Calculates covariance matrix of Z
  void calculate_C();
  //! The first term in the expression of \f$\overline{M}\f$
  void calculate_EHH();
  //! Calculation of \f$\overline{M}\f$
  void calculate_Mbar();
  //! Partial derivatives of \f$\breve{\alpha}\f$
  matrix calculate_Dalpha();
  //! \f$\overline{C} = diag(\tilde{C}_1,\ldots,\tilde{C}_N)\f$ where 
  // \\f$tilde{C}_j(\theta)\f$ = C_j' + E_\theta(\tilde{Z}_j^{(s)})^\otimes 2\f$
  // where Z_j' is the covariance matrix of Z_j^{(s)} with an extra column+row of zeros inserted at position 0.
  matrix calculate_Cbar();
  //! Calculates \f$A_n^*(\theta)\f$ derived from Mbar,Dalpha and Cbar.
  void calculate_Astar();

  //! 

};


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


//! Markov model class
/*! 
  inherits from class PREESTMODEL
*/
class MARKOVMODEL: public PREESTMODEL {
 public: 
  virtual void setinitvalue() =0;              
  //! updating in the simulation.        
  virtual double next(double oldvalue) const=0;
  //! simulates a process given by "next".
  /*! 
    
  \param n : number of variables to simulate.
  
  \return 
  */
  matrix simulate(long n);
  //! sets the startvalue for the simulation.
  void setstartvalue(double x) {startvalue = x;}
  //! sets delta.
  void setdelta(double d) {delta = d;}

 protected:
  //! variable used in the simulation.
  double delta;
  //! startvalue for the simulation.
  double startvalue;
};

//! Stochastic volatility model
/*! 
  inherits from MARKOVMODEL
*/
class SVMODEL: public MARKOVMODEL {
 protected:
  //! the timedistance between the observations. 
  double Delta;
  //! the number of subintervals for each Delta.
  long PointsPerDelta;
 public:
  //! simulates the SV-process.
  virtual matrix simulate(long n);
  //! simulates the volatility of the SV-process.
  virtual matrix simulateS(long n);
};






//////////////////////////////////////////////////////
//         GLOBAL DEFINITIONS                        //
//////////////////////////////////////////////////////
 
typedef  double (PREESTMODEL::*PREpointerFn)(matrix par);
#define callMemberFunction(object,ptrToMember)((object).*(ptrToMember))







//********************* REAL MODELS: **********************


//! Stochastic volatility model with CIR process
/*! 
  inherits from SVMODEL
*/
class SVCIR: public SVMODEL {
 public:
  //! Constructor of SV-CIR-process model
  /*! 
    
  \param FF : Inner functions
  \param d : The time-distance between the observations (delta)
  \param PPD : Points pr. delta (time-distance)
  \param N1 : Number of variables to simulate in Mbar-calculation
  \param N2 : Number of variables to simulate to estimate moments of process 
  \param qq : The lag-length in the prediction
  \param ss : Summation-start in estimating function
  \param ll :Summation-start in the expresion of Mbar
  \param sd : seed
  
  \return 
  */
  SVCIR(functionvector FF, double d = 1.0, long PPD = 200, long N1=150, long N2=150, long qq=5, long ss = 6, long ll=10, long sd = -123456) {
    ifuncs = FF;
    Delta=d; PointsPerDelta=PPD; setdelta(d/PPD); MbarNOS = N1; q = qq; 
    if (q>ss) 
      s = q;
    else
      s = ss; 
    l = ll; setseed(sd); 
    iter_Mbar = true; // Astar (and hence Mbar) is calculated
    numsimul = N2;
    N = FF.size();
  }
  //! Sets parameters of the model. Only used in simulations
  /*! 

  \param par : \f$ 3\times 1\f$-matrix containing parameters (alpha,theta,sigma).
  */
  virtual void setpar(matrix par);
  virtual unsigned getdim();

 protected:

  double alpha;                            // par[0]
  double theta;                            // par[1]
  double sigma;                            // par[2]

  virtual void setinitvalue() {MARKOVMODEL::setstartvalue(alpha);}
  //! updating in the simulation.
  virtual double next(double oldvalue) const;

};



//! Stochastic volatility model with Exponential Ornstein-Uhlenbeck-type process
/*! 
  inherits from SVMODEL
*/
class SVEOU: public SVMODEL 
{
 public:
  //! Constructor of SVEOU-process model
  /*! 

  \param FF : Inner functions
  \param d : The time-distance between the observations (delta)
  \param PPD : Points pr. delta (time-distance)
  \param N1 : Number of variables to simulate in Mbar-calculation
  \param N2 : Number of variables to simulate to estimate moments of process 
  \param qq : The lag-length in the prediction
  \param ss : Summation-start in estimating function
  \param ll :Summation-start in the expresion of Mbar
  \param sd : seed  
  
  \return 
  */
  SVEOU(functionvector FF, double d = 1.0, long PPD = 200, long N1=100, long N2=100, long qq=5, long ss = 6, long ll=10, long sd = -123456) {
    ifuncs = FF;
    Delta=d; PointsPerDelta=PPD; setdelta(d/PPD); MbarNOS = N1; q = qq; 
    if (q>ss) 
      s = q;
    else
      s = ss; 
    l = ll; setseed(sd); 
    iter_Mbar = true; 
    numsimul = N2;
    N = FF.size();
  }
  virtual inline ~SVEOU() {}

  //! Sets parameters of the model. Only used in simulations
  /*! 
   
  \param par : \f$ 3\times 1\f$-matrix containing parameters (alpha,theta,sigma).
  */
  virtual void setpar(matrix par);
  virtual unsigned getdim();

 protected:
  double alpha;  // Parameters[0]
  double theta;  // Parameters[1]
  double sigma;  // Parameters[2]
  virtual double next(double oldvalue) const; 
  virtual void setinitvalue() {MARKOVMODEL::setstartvalue(alpha);} ///exp(theta));} 
};



//! Stochastic volatility model with Exponential Ornstein-Uhlenbeck-type process with NEW Z_i
/*! 
  inherits from SVEOU. Simple example of implementation of new Z_i-fuunctions. In this case 
  the squareroot of each f_i is taken.
*/
class SVEOU1: public SVEOU
{
 public:
  //! Constructor of SVEOU1-process model
  /*! 

  \param FF : Inner functions
  \param d : The time-distance between the observations (delta)
  \param PPD : Points pr. delta (time-distance)
  \param N1 : Number of variables to simulate in Mbar-calculation
  \param N2 : Number of variables to simulate to estimate moments of process 
  \param qq : The lag-length in the prediction
  \param ss : Summation-start in estimating function
  \param ll :Summation-start in the expresion of Mbar
  \param sd : seed  
  
  \return 
  */
  SVEOU1(functionvector FF, double d = 1.0, long PPD = 200, long N1=100, long N2=100, long qq=5, long ss = 6, long ll=10, long sd = -123456): SVEOU(FF,d,PPD,N1,N2,qq,ss,ll,sd) {}
  virtual inline ~SVEOU1() {}
 
 protected:
  //! Z_{jk}^{(i)}, j=0,...,N-1, k=1,...,q, i
  virtual double Zval(unsigned j, long i, long k, matrix const& Y);
  virtual void calculate_Z();  
};





#endif // _MODELS_H_

---